运筹学与最优化MATLAB编程.docx
《运筹学与最优化MATLAB编程.docx》由会员分享，可在线阅读，更多相关《运筹学与最优化MATLAB编程.docx（11页珍藏版）》请在冰豆网上搜索。
运筹学与最优化MATLAB编程
运筹学与最优化MATLAB编程
实验报告
院系：
专业：
姓名：
学号：
指导老师：
完成日期：
割平面法求解整数规划问题
一、引言：
通过对MATLAB实践设计的学习，学会使用MATLAB解决现实生活中的问题。
该设计是在MATLAB程序设计语言的基础上，对实际问题建立数学模型并设计程序，使用割平面法解决一个整数规划问题。
经实验，该算法可成功运行并求解出最优整数解。
二、算法说明：
割平面法有许多种类型，本次设计的原理是依据Gomory的割平面法。
Gomory割平面法首先求解非整数约束的线性规划，再选择一个不是整数的基变量，定义新的约束，增加到原来的约束中，新的约束缩小了可行域，但是保留了原问题的全部整数可行解。
算法具体设计步骤如下：
1、首先，求解原整数规划对应的线性规划
，设最优解为x*。
2、如果最优解的分量均为整数，则x*为原整数规划的最优解；否则任选一个x*中不为整数的分量，设其对应的基变量为xp，定义包含这个基变量的切割约束方程
，其中xp为非基变量。
3、令
，
，其中[]为高斯函数符号，表示不大于某数的最大整数。
将切割约束方程变换为
，由于0<
<1，0<
<1，所以有
，因为自变量为整数，则
也为整数，所以进一步有
。
4、将切割方程加入约束方程中，用对偶单纯形法求解线性规划
，然后在转入步骤2进行求解，直到求出最优整数解停止迭代。
三、程序实现：
程序设计流程图如图1，具体设计代码（见附录）。
图1
四、算例分析：
已知AM工厂是一个拥有四个车间的玩具生产厂商，该厂商今年新设计出A、B、C、D、E、F六种玩具模型，根据以前的生产情况及市场调查预测，得知生产每单位产品所需的工时、每个车间在一季度的工时上限以及产品的预测价格，如下表所示。
问：
每种设计产品在这个季度各应生产多少，才能使AM工厂这个季度的生产总值达到最大？
产品
车间
A
B
C
D
E
F
每个车间每季度工时上限
甲
乙
丙
丁
0.01
0.02
0.01
0.02
0.01
0.03
0.03
0.05
0.03
0.05
0.03
0.08
850
700
100
900
单价（元）
20
14
16
36
32
30
1、问题分析并建立模型：
由题意可知这是一个求解产量使产值最大的整数规划问题。
根据上述问题和已知数据,可以假设每种产品在这个季度各应生产产量分别为：
x1、x2、x3、x4、x5、x6，则有以下线性方程组
maxZ=20x1+14x2+16x3+36x4+32x5+30x6
2、实验步骤：
首先引入松弛变量x7、x8、x9、x10，使其化为标准型
minZ=-20x1-14x2-16x3-36x4-32x5-30x6
其次从标准型可表示出约束系数矩阵、右端项常数矩阵、目标系数矩阵分别为A、b、c。
然后调用DividePlane函数，使用割平面法进行求解。
在MATLAB的命令窗口输入一下命令：
>>A=[0.010.010.010.030.030.031000;0.02000.05000100;00.02000.0500010;000.03000.080001];
>>c=[-20;-14;-16;-36;-32;-30];
>>b=[850;700;100;900];
>>[intx,intf]=DividePlane（A,c,b,[78910]）
3、实验结果及分析：
intx=
35000500030000000
intf=
-1250000
实验结果求出的目标函数值是化为标准型的最小值，则转化为原问题的目标函数值应取相反数，所以从实验结果可知：
生产各种产品的产量分别应为为，生产A35000、生产B5000、生产C30000、生产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;
iflength（baseVector）~=n
disp（'基变量的个数要与约束矩阵的行数相等！
'）;
mx=NaN;
mf=NaN;
return;
end
M=0;
sigma=-[transpose（c）zeros（1,（nVia-length（c）））];
xb=b;
%首先用单纯形法求出最优解
while1
[maxs,ind]=max（sigma）;
%--------------------用单纯形法求最优解--------------------------------------
ifmaxs<=0%当检验数均小于0时，求得最优解。
vr=find（c~=0,1,'last'）;
forl=1:
vr
ele=find（baseVector==l,1）;
if（isempty（ele））
mx（l）=0;
else
mx（l）=xb（ele）;
end
end
ifmax（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））;
fori=1:
（index_x-1）
Atmp（1,i）=A（num,i）-floor（A（num,i））;
end
fori=（index_x+1）:
sc
Atmp（1,i）=A（num,i）-floor（A（num,i））;
end
Atmp（1,index_x）=0;%构建对偶单纯形法的初始表格
A=[Azeros（sr,1）;-Atmp（1,:
）1];
xb=[xb;-fi];
baseVector=[baseVectorsc+1];
sigma=[sigma0];
%-------------------对偶单纯形法的迭代过程----------------------
while1
%----------------------------------------------------------
ifxb>=0%判断如果右端向量均大于0，求得最优解
ifmax（abs（round（xb）-xb））<1.0e-7%如果用对偶单纯形法求得了整数解，则返回最优整数解
vr=find（c~=0,1,'last'）;
forl=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））;
fori=1:
（index_x-1）
Atmp（1,i）=A（num,i）-floor（A（num,i））;
end
fori=（index_x+1）:
sc
Atmp（1,i）=A（num,i）-floor（A（num,i））;
end
Atmp（1,index_x）=0;%下一次对偶单纯形迭代的初始表格
A=[Azeros（sr,1）;-Atmp（1,:
）1];
xb=[xb;-fi];
baseVector=[baseVectorsc+1];
sigma=[sigma0];
continue;
end
else%如果右端向量不全大于0，则进行对偶单纯形法的换基变量过程
minb_1=inf;
chagB_1=inf;
sA=size（A）;
[br,idb]=min（xb）;
forj=1:
sA
（2）
ifA（idb,j）<0
bm=sigma（j）/A（idb,j）;
ifbmminb_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）;
fori=1:
sA
（1）
ifi~=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;
forj=1:
n
ifA（j,ind）>0
bz=xb（j）/A（j,ind）;
ifbzminb=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）;
fori=1:
n
ifi~=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
升级会员 