运筹学与最优化MATLAB 编程实验报告割平面法求解整数规划问题一、 引言：通过对MATLAB 实践设计的学习，学会使用MATLAB 解决现实生活中的问题。

该设计是在MATLAB 程序设计语言的基础上，对实际问题建立数学模型并设计程序，使用割平面法解决一个整数规划问题。

经实验，该算法可成功运行并求解出最优整数解。

二、 算法说明：割平面法有许多种类型，本次设计的原理是依据Gomory 的割平面法。

Gomory 割平面法首先求解非整数约束的线性规划，再选择一个不是整数的基变量，定义新的约束，增加到原来的约束中，新的约束缩小了可行域，但是保留了原问题的全部整数可行解。

算法具体设计步骤如下：1、首先，求解原整数规划对应的线性规划,*)(min x c x f =⎩⎨⎧≥≤0..x bAx t s ，设最优解为x*。

2、如果最优解的分量均为整数，则x*为原整数规划的最优解；否则任选一个x*中不为整数的分量，设其对应的基变量为x p ，定义包含这个基变量的切割约束方程con jj ij p b x r x =+∑，其中x p 为非基变量。

3、令][ij ij ij r r r -=，][con con con b b b -=，其中[]为高斯函数符号，表示不大于某数的最大整数。

将切割约束方程变换为∑∑-=-+jjij con con jj ij p x r b b x r x ][][，由于0<ij r <1，0<con b <1，所以有1<-∑jj ij con x r b ，因为自变量为整数，则∑-jj ij con x r b 也为整数，所以进一步有0≤-∑jj 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 ,090008.003.010005.002.070005.002.085003.003.003.001.001.001.0..635241654321 且为整数，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 ,090008.003.010005.002.070005.002.085003.003.003.001.001.001.0..10639528417654321 且为整数，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.050 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 01];>> 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 0intf =-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) ~= ndisp('基变量的个数要与约束矩阵的行数相等！');mx = NaN;mf = NaN;return;endM = 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:vrele = find(baseVector == l,1);if(isempty(ele))mx(l) = 0;elsemx(l)=xb(ele);endendif 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));endfor i=(index_x+1):scAtmp(1,i) = A(num,i) - floor(A(num,i));endAtmp(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:vrele = find(baseVector == l,1);if(isempty(ele))mx_1(l) = 0;elsemx_1(l)=xb(ele);endendintx = 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));endfor i=(index_x+1):scAtmp(1,i) = A(num,i) - floor(A(num,i));endAtmp(1,index_x) = 0; %下一次对偶单纯形迭代的初始表格A = [A zeros(sr,1);-Atmp(1,:) 1];xb = [xb;-fi];baseVector = [baseVector sc+1];sigma = [sigma 0];continue;endelse%如果右端向量不全大于0，则进行对偶单纯形法的换基变量过程 minb_1 = inf;chagB_1 = inf;sA = size(A);[br,idb] = min(xb);for j=1:sA(2)if A(idb,j)<0bm = sigma(j)/A(idb,j);if bm<minb_1minb_1 = bm;chagB_1 = j;endendendsigma = 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 ~= idbxb(i) = xb(i)-A(i,chagB_1)*xb(idb);A(i,:) = A(i,:) - A(i,chagB_1)*A(idb,:);endendbaseVector(idb) = chagB_1;end%------------------------------------------------------------end%--------------------对偶单纯形法的迭代过程--------------------- endelse%如果检验数有不小于0的，则进行单纯形算法的迭代过程minb = inf;chagB = inf;for j=1:nif A(j,ind)>0bz = xb(j)/A(j,ind);if bz<minbminb = bz;chagB = j;endendendsigma = sigma -A(chagB,:)*maxs/A(chagB,ind);xb(chagB) = xb(chagB)/A(chagB,ind);A(chagB,:) = A(chagB,:)/A(chagB,ind);for i =1:nif i ~= chagBxb(i) = xb(i)-A(i,ind)*xb(chagB);A(i,:) = A(i,:) - A(i,ind)*A(chagB,:);endendbaseVector(chagB) = ind;endM = M + 1;if (M == 1000000)disp('找不到最优解！');mx = NaN;minf = NaN;return;endend。