运筹学与最优化MATLAB编程
实验报告
院系：
专业：
姓名：
学号：
指导老师：
完成日期：
割平面法求解整数规划问题
一、引言：
通过对MATLAB实践设计的学习，学会使用MATLAB解决现实生活中的问题。该设计是在MATLAB程序设计语言的基础上，对实际问题建立数学模型并设计程序，使用割平面法解决一个整数规
划问题。经实验，该算法可成功运行并求解出最优整数解。 二、 算法说明：
割平面法有许多种类型，本次设计的原理是依据Gomory 的割平面法。Gomory 割平面法首先求解非整数约束的线性规划，再选择一个不是整数的基变量，定义新的约束，增加到原来的约束中，新的约束缩小了可行域，但是保留了原问题的全部整数可行解。
算法具体设计步骤如下：
1、首先，求解原整数规划对应的线性规划
,*)(m i n x c x f =???≥≤0
..x b
Ax t s ，设最优解为x*。
2、如果最优解的分量均为整数，则x*为原整数规划的最优解；否则任选一个x*中不为整数的分量，设其对应的基变量为x p ，定义包含这个基变量的切割约束方程con j
j ij p b x r x =+∑，其中x p 为非基变量。
3、令][ij ij ij r r r -=，][con con con b b b -=，其中[]为高斯函数符号，表示不大于某数的最大整数。将切割约束方程变换为
∑∑-=-+j
j
ij con con j
j ij p x r b b x r x ][][，由于0 1<-∑j j ij con x r b ，因为自变量为整数，则∑-j j ij con x r b 也为整数，所以进 一步有0≤-∑j j ij con x r b 。 4、将切割方程加入约束方程中，用对偶单纯形法求解线性规划 ?? ? ? ???≥≤-≤=∑00..,*)(min x x r b b Ax t s x c x f j j ij con ，然后在转入步骤2进行求解，直 到求出最优整数解停止迭代。 三、程序实现： 程序设计流程图如图1，具体设计代码(见附录)。 图1 四、算例分析： 已知AM工厂是一个拥有四个车间的玩具生产厂商，该厂商今年新设计出A、B、C、D、E、F六种玩具模型，根据以前的生产情况及市场调查预测，得知生产每单位产品所需的工时、每个车间在一季度的工时上限以及产品的预测价格，如下表所示。问：每种设计产品在这个季度各应生产多少，才能使AM工厂这个季度的生产总值达到最大？ 1、问题分析并建立模型： 由题意可知这是一个求解产量使产值最大的整数规划问题。根据上述问题和已知数据,可以假设每种产品在这个季度各应生产产量分别为：x 1、x 2、x 3、x 4、x 5、x 6，则有以下线性方程组 maxZ=20x 1+14x 2+16x 3+36x 4+32x 5+30x 6 ????? ????=≥≤+≤+≤+≤+++++6 ,,1i ,0900 08.003.010005.002.0700 05.002.085003.003.003.001.001.001.0..6 35241654321 且为整数，i x x x x x x x x x x x x x t s 2、实验步骤： 首先引入松弛变量x 7、x 8、x 9 、x 10，使其化为标准型 minZ=-20x 1-14x 2-16x 3-36x 4-32x 5-30x 6 ????? ????=≥=++=++=++=++++++10 ,,1i ,0900 08.003.010005.002.0700 05.002.085003.003.003.001.001.001.0..10 639528417654321 且为整数，i x x x x x x x x x x x x x x x x x t s 其次从标准型可表示出约束系数矩阵、右端项常数矩阵、目标系数矩阵分别为A、b、c。然后调用DividePlane函数，使用割平面法进行求解。 在MATLAB的命令窗口输入一下命令： >> A=[0.01 0.01 0.01 0.03 0.03 0.03 1 0 0 0;0.02 0 0 0.05 0 0 0 1 0 0;0 0.02 0 0 0.05 0 0 0 1 0;0 0 0.03 0 0 0.08 0 0 0 1]; >> c=[-20;-14;-16;-36;-32;-30]; >> b=[850;700;100;900]; >> [intx,intf]=DividePlane(A,c,b,[7 8 9 10]) 3、实验结果及分析： intx = 35000 5000 30000 0 0 0 intf = -1250000 实验结果求出的目标函数值是化为标准型的最小值，则转化为原问题的目标函数值应取相反数，所以从实验结果可知：生产各种产品的产量分别应为为，生产A 35000、生产B 5000、生产C 30000、生产D 、E、F均为0，此时的季度产值为最大即1250000元。该结果是可信的，故通过该实例说明该程序能够运用于实际，用来解决实际生活中求解整数规划的问题。 五、结束语： Matlab是个很强大的软件，提供了大量的函数来处理各种数学、工程、运筹等的问题，并且含有处理二维、三维图形的功能，使用matlab能够解决许多实际生活中的问题。通过这个学期的学习，仅是了解了matlab的部分函数功能和简单的GUI界面设计，掌握了一些基本的程序编写技能，同时，在老师的指导下简单了解了使用LinGo 和Excel解决线性和非线性规划问题的求解方法，收获相当丰富，同时认识到要学好matlab仍然需要一个长期的过程。 六、参考文献： [1] 龚纯，王正林.精通MATLAB最优化计算.北京：电子工业出版社，2009 [2]吴祈宗，郑志勇，邓伟等.运筹学与最优化MATLAB编程.北京：机械工业出版社，2009 [3]邓成梁.运筹学的原理和方法(第二版).武汉：华中科技大学出版社，2002 七、附录： function [intx,intf] = DividePlane(A,c,b,baseVector) %功能：用割平面法求解整数规划 %调用格式：[intx,intf]=DividePlane(A,c,b,baseVector) %其中，A：约束矩阵； % c：目标函数系数向量； % b：约束右端向量； % baseVector：初始基向量； % intx：目标函数取最小值时的自变量值； % intf：目标函数的最小值； sz = size(A); nVia = sz(2); n = sz(1); xx = 1:nVia; if length(baseVector) ~= n disp('基变量的个数要与约束矩阵的行数相等！'); mx = NaN; mf = NaN; return; end M = 0; sigma = -[transpose(c) zeros(1,(nVia-length(c)))]; xb = b; %首先用单纯形法求出最优解 while 1 [maxs,ind] = max(sigma); %--------------------用单纯形法求最优解--------------------------------------if maxs <= 0 %当检验数均小于0时，求得最优解。 vr = find(c~=0 ,1,'last'); for l=1:vr ele = find(baseVector == l,1); if(isempty(ele)) mx(l) = 0; else mx(l)=xb(ele); end end if max(abs(round(mx) - mx))<1.0e-7 %判断最优解是否为整数解，如果是整数解。 intx = mx; intf = mx*c; return; else%如果最优解不是整数解时，构建切割方程 sz = size(A); sr = sz(1); sc = sz(2); [max_x, index_x] = max(abs(round(mx) - mx)); [isB, num] = find(index_x == baseVector); fi = xb(num) - floor(xb(num)); for i=1:(index_x-1) Atmp(1,i) = A(num,i) - floor(A(num,i)); end for i=(index_x+1):sc Atmp(1,i) = A(num,i) - floor(A(num,i)); end Atmp(1,index_x) = 0; %构建对偶单纯形法的初始表格 A = [A zeros(sr,1);-Atmp(1,:) 1]; xb = [xb;-fi]; baseVector = [baseVector sc+1]; sigma = [sigma 0]; %-------------------对偶单纯形法的迭代过程---------------------- while 1 %---------------------------------------------------------- if xb >= 0 %判断如果右端向量均大于0，求得最优解 if max(abs(round(xb) - xb))<1.0e-7 %如果用对偶单纯形法求得了整数解，则返回最优整数解 vr = find(c~=0 ,1,'last'); for l=1:vr ele = find(baseVector == l,1); if(isempty(ele)) mx_1(l) = 0; else mx_1(l)=xb(ele); end end intx = mx_1; intf = mx_1*c; return; else%如果对偶单纯形法求得的最优解不是整数解，继续添加切割方程 sz = size(A); sr = sz(1); sc = sz(2); [max_x, index_x] = max(abs(round(mx_1) - mx_1)); [isB, num] = find(index_x == baseVector); fi = xb(num) - floor(xb(num)); for i=1:(index_x-1) Atmp(1,i) = A(num,i) - floor(A(num,i)); end for i=(index_x+1):sc Atmp(1,i) = A(num,i) - floor(A(num,i)); end Atmp(1,index_x) = 0; %下一次对偶单纯形迭代的初始表格 A = [A zeros(sr,1);-Atmp(1,:) 1]; xb = [xb;-fi]; baseVector = [baseVector sc+1]; sigma = [sigma 0]; continue; end else%如果右端向量不全大于0，则进行对偶单纯形法的换基变量过程 minb_1 = inf; chagB_1 = inf; sA = size(A); [br,idb] = min(xb); for j=1:sA(2) if A(idb,j)<0 bm = sigma(j)/A(idb,j); if bm minb_1 = bm; chagB_1 = j; end end end sigma = sigma -A(idb,:)*minb_1; xb(idb) = xb(idb)/A(idb,chagB_1); A(idb,:) = A(idb,:)/A(idb,chagB_1); for i =1:sA(1) if i ~= idb xb(i) = xb(i)-A(i,chagB_1)*xb(idb); A(i,:) = A(i,:) - A(i,chagB_1)*A(idb,:); end end baseVector(idb) = chagB_1; end %------------------------------------------------------------end %--------------------对偶单纯形法的迭代过程--------------------- end else%如果检验数有不小于0的，则进行单纯形算法的迭代过程 minb = inf; chagB = inf; for j=1:n if A(j,ind)>0 bz = xb(j)/A(j,ind); if bz minb = bz; chagB = j; end end end sigma = sigma -A(chagB,:)*maxs/A(chagB,ind); xb(chagB) = xb(chagB)/A(chagB,ind); A(chagB,:) = A(chagB,:)/A(chagB,ind); for i =1:n if i ~= chagB xb(i) = xb(i)-A(i,ind)*xb(chagB); A(i,:) = A(i,:) - A(i,ind)*A(chagB,:); end end baseVector(chagB) = ind; end M = M + 1; if (M == 1000000) disp('找不到最优解！'); mx = NaN; minf = NaN; return; end end 常用无约束最优化方法评价准则 方法算法特点适用条件 最速下降法属于间接法之一。方法简便，但要计算一阶偏导 数，可靠性较好，能稳定地使函数下降，但收敛 速度较慢，尤其在极点值附近更为严重 适用于精度要求不高或用于对 复杂函数寻找一个好的初始 点。 Newton法属于间接法之一。需计算一、二阶偏导数和Hesse 矩阵的逆矩阵，准备工作量大，算法复杂，占用 内存量大。此法具有二次收敛性，在一定条件下 其收敛速度快，要求迭代点的Hesse矩阵必须非 奇异且定型（正定或负定）。对初始点要求较高， 可靠性较差。 目标函数存在一阶\二阶偏导 数，且维数不宜太高。 共轭方向法属于间接法之一。具有可靠性好，占用内存少， 收敛速度快的特点。 适用于维数较高的目标函数。 变尺度法属于间接法之一。具有二次收敛性，收敛速度快。 可靠性较好，只需计算一阶偏导数。对初始点要 求不高，优于Newton法。因此，目前认为此法是 最有效的方法之一，但需内存量大。对维数太高 的问题不太适宜。 适用维数较高的目标函数 （n=10～50）且具有一阶偏导 数。 坐标轮换法最简单的直接法之一。只需计算函数值，无需求 导，使用时准备工作量少。占用内存少。但计算 效率低，可靠性差。 用于维数较低（n<5）或目标函 数不易求导的情况。 单纯形法此法简单，直观，属直接法之一。上机计算过程 中占用内存少，规则单纯形法终止条件简单，而 不规则单纯形法终止条件复杂，应注意选择，才 可能保证计算的可靠性。 可用于维数较高的目标函数。 常用约束最优化方法评价标准 方法算法特点适用条件 外点法将约束优化问题转化为一系列无约束优化问题。 初始点可以任选，罚因子应取为单调递增数列。 初始罚因子及递增系数应取适当较大值。 可用于求解含有等式约束或不等 式约束的中等维数的约束最优化 问题。 内点法将约束优化问题转化为一系列无约束优化问题。 初始点应取为严格满足各个不等式约束的内点， 障碍因子应取为单调递减的正数序列。初始障碍 因子选择恰当与否对收敛速度和求解成败有较大 影响。 可用于求解只含有不等式约束的 中等维数约束优化问题。 混合罚函数法将约束优化问题转化为一系列无约束优化问题， 用内点形式的混合罚函数时，初始点及障碍因子 的取法同上；用外点形式的混合罚函数时，初始 点可任选，罚因子取法同外点法相同。 可用于求解既有等式约束又有不 等式约束的中等维数的约束化问 题。 约束坐标轮换法由可行点出发，分别沿各坐标轴方向以加步探索 法进行搜索，使每个搜索点在可行域内，且使目 标函数值下降。 可用于求解只含有不等式约束， 且维数较低（n<5），目标函数的 二次性较强的优化问题。 复合形法在可行域内构造一个具有n个顶点的复合形，然 后对复合形进行映射变化，逐次去掉目标函数值 最大的顶点。 可用于求解含不等式约束和边界 约束的低维优化问题。 项目三 常用无约束最优化方法（一） [实验目的] 编写最速下降法、Newton 法（修正Newton 法）的程序。 [实验学时] 2学时 [实验准备] 1．掌握最速下降法的思想及迭代步骤。 2．掌握Newton 法的思想及迭代步骤； 3．掌握修正Newton 法的思想及迭代步骤。 [实验内容及步骤] 编程解决以下问题：【选作一个】 1．用最速下降法求 22120min ()25[22]0.01T f X x x X ε=+==，，，． 2．用Newton 法求 22121212min ()60104f X x x x x x x =--++-， 初始点 0[00]0.01T X ε==，，． 最速下降法 Matlab 程序： clc;clear; syms x1 x2; X=[x1,x2]; fx=X(1)^2+X(2)^2-4*X(1)-6*X(2)+17; fxd1=[diff(fx,x1) diff(fx,x2)]; x=[2 3]; g=0; e=0.0005; a=1; fan=subs(fxd1,[x1 x2],[x(1) x(2)]); g=0; for i=1:length(fan) g=g+fan(i)^2; end g=sqrt(g); step=0; while g>e step=step+1; dk=-fan; %点x(k)处的搜索步长 ak=((2*x(1)-4)*dk(1)+(2*x(2)-6)*dk(2))/(dk(1)*dk(2)-2*dk(1)^2-2*dk(2)^2); xu=x+ak*dk; x=xu; %输出结果 optim_fx=subs(fx,[x1 x2],[x(1) x(2)]); fprintf(' x=[ %d %d ] optim_fx=%d\n',x(1),x(2),optim_fx); %计算目标函数点x(k+1)处一阶导数值 fan=subs(fxd1,[x1 x2],[x(1) x(2)]); g=0; for i=1:length(fan) g=g+fan(i)^2; end g=sqrt(g); end %输出结果 optim_fx=subs(fx,[x1 x2],[x(1) x(2)]); fprintf('\n最速下降法\n结果：\n x=[ %d %d ] optim_fx=%d\n',x(1),x(2),optim_fx); c++程序 #include 第三章无约束最优化方法 本章内容及教学安排 第一节概述 第二节迭代终止原则 第三节常用的一维搜索方法 第四节梯度法 第五节牛顿法 第六节共轭方向法 第七节变尺度法 第八节坐标轮换法 第九节鲍威尔方法 第一节概述 优化问题可分为 无约束优化问题 有约束优化问题 无约束最优化问题求解基于古典极值理论的一种数值迭代方法，主要用来求解非线性规划问题 迭代法的基本思想： 所以迭代法要解决三个问题 1、如何选择搜索方向 2、如何确定步长 3、如何确定最优点（终止迭代） 第二节 迭代终止准则 1）1K K X X ε+-≤ 111/2 21K K K K n i i i X X X X ε++=??-=-≤???? ∑（） 2） 11()()()() () K K K K K f X f X f X f X or f X ε ε ++-≤-≤ 3）(1)()K f X ε+?≤ 第三节 常用的一维搜索方法 本节主要解决的是如何确定最优步长的问题。 从初始点(0)X 出发，以一定的步长沿某一个方向，可以找到一个新的迭代点，其公式如下： (1)(0)00(2)(1)11(1)() K K k k X X S X X S X X S ααα+=+=+= + 现在假设K S 已经确定，需要确定的是步长k α，就把求多维目标函数的极小值这个多维算过程中，当起步点和方向问题，变成求一个变量即步长的最优值的一维问题了。即 (1)()min ()min ()min ()K K K k k f X f X S f αα+=+= 由此可见，最佳步长*K α由一维搜索方法来确定 求*k α，使得()()()()()()min K K K K f f X S αα=+→ 一、一维搜索区间的确定 区间[,]a b 应满足 ()(*)()f a f f b α><常用最优化方法评价准则
常用无约束最优化方法(一)
第三章 无约束最优化方法
相关主题
文本预览