运筹学实验报告.docx
《运筹学实验报告.docx》由会员分享，可在线阅读，更多相关《运筹学实验报告.docx（16页珍藏版）》请在冰豆网上搜索。
运筹学实验报告
（此文档为word格式，下载后您可任意编辑修改！
）
运筹学实验报告
班级：
数电四班姓名：
刘文搏学号：
一、实验目的
运用MATLAB程序设计语言完成单纯性算法求解线性规划问题。
二、实验内容
编写一个MATLAB的函数文件：
linp.m用于求解标准形的线性规划问题：
minf=c*xsubjectto：
A*x=b；x>=0;
1、函数基本调用形式：
[x,minf,optmatrx,flag]=linp（A,b,c）
2、参数介绍：
A:
线性规划问题的约束A*x=b且x>=0中变量的系数组成的矩阵，是一个m*n的矩阵。
c:
线性规划问题的目标函数f=c*x中各变量的系数向量，是一个n维的行
向量。
b:
线性规划问题的约束A*x=b且x>=0中的常数向量，是一个m维的列
向量。
x:
输出线性规划问题的最优解，当线性规划问题没有可行解或有可行解无
最优解时x=[].
minf:
输出线性规划问题的最优值，当线性规划问题没有可行解时minf=[],
当线性规划问题有可行解无最优解时minf=-Inf。
flag:
线性规划问题的求解结果标志值，当线性规划问题有最优解时flag=1,
当线性规划问题有可行解无最优解时flag=0,当线性规划问题没有可行解时flag=-1.
cpt:
输出最优解对应的单纯性表，当线性规划问题没有可行解或有可
行解无最优解时cpt=[].
三、Linp函数
%此函数是使用两阶段算法求解线性规划问题
function[x,minf,flag,cpt]=linp（A,b,c）;
fori=1:
p%判断b是否<=0;将b转换成大于0；
ifb（i）<0
A（i,:
）=-1*A（i,:
）;
b（i）=-1*b（i）;
end
end
%返回值：
x，第一张单纯形表，基，标志参数A，c，b
%********第一张单纯形表的初始化
[m,n]=size（A）;%获得矩阵A的维数
[p,q]=size（b）;
dcxb=zeros（m+2,m+n+1）;%确定第一张单纯形表的大小
dcxb（1,:
）=[-c,zeros（1,m+1）];%¸给表的第一行赋值
dcxb（2,:
）=[zeros（1,n）,-1*ones（1,m）,0];%¸给表的第二行赋值
dcxb（[3:
m+2],:
）=[A,eye（m,m）,b];%添A和b到表中
jxl=[n+1:
n+m];
fori=3:
m+2
dcxb（2,:
）=dcxb（2,:
）+dcxb（i,:
）;
fori=3:
m+2
dxcb（2,:
）=dxcb1（2,:
）+dxcb1（i,:
）;
end
dxcb;
%************辅助问题换基迭代**********************
dyl=find（dcxb（2,[1:
m+n]）>0）;
while~isempty（dyl）
firstnum=dyl
（1）;
dll=dcxb（[3:
m+2],firstnum）;
youduanb=dcxb（[3:
m+2],m+n+1）;
look=find（dll>0）;
ifisempty（look）
dcxb（2,firstnum）=0;
else
min=Inf;
fori=3:
m+2
ifdll（i-2）>0&youduanb（i-2）dll（i-2）min=youduanb（i-2）dll（i-2）;
line1=i;
end
end
dcxb（line1,:
）=dcxb（line1,:
）dcxb（line1,firstnum）;
fori=1:
m+2
ifi~=line1
dcxb（i,:
）=dcxb（i,:
）+（-1*dcxb（i,firstnum）*dcxb（line1,:
））;
end
end
jxl（line1-2）=firstnum;
end
dyl=find（dcxb（2,[1:
m+n]）>0）;
dcxb
end
ifdcxb（2,m+n+1）>0
fprintf（'g>0,´此问题没有可行解'）;
x=[];
minf=inf;
cpt=[];
flag=-1;
return
end
look1=find（jxl>n）;
ifdcxb（2,m+n+1）==0%等于0，判断基变量中是否有人工变量；
if~isempty（look1）%´存在时进行处理
while~isempty（look1）
line2=look1
（1）+2;
chdy0=find（dcxb（line2,[1:
n]）~=0）;
ifisempty（chdy0）%´存在人工变量都为零的那一行，去掉该行
dcxb（line2,:
）=[];
look1
（1）=[];
jxl（line2-2）=[];
else%否则进行换基迭代
secondnum=chdy0
（1）;
dcxb（line2,:
）=dcxb（line2,:
）dcxb（line2,secondnum）;
jxl（line2-2）=secondnum;
fori=1:
m+2
ifi~=line2
dcxb（i,:
）=dcxb（i,:
）+（-1*dcxb（i,secondnum）*dcxb（line2,:
））;
end;
end
look1
（1）=[];
end;
end
end
end
%去掉人工变量，得到单纯性的第一张表
dcxb（2,:
）=[];
dcxb（:
[n+1:
n+m]）=[];
%有可行解，判断z
dcxb2=dcxb;
look2=find（dcxb2（1,[1:
n]）>0）;
while~isempty（look2）
thirdnum=look2
（1）;
duilie=dcxb2（[2:
m+1],thirdnum）;
youduanb1=dcxb2（[2:
m+1],n+1）;
look3=find（duilie>0）;
ifisempty（look3）
fprintf（'´此问题有可行解，但没有最优解'）;
x=zeros（n,1）;
[mi,n1]=size（jxl）;
fori=1:
n1
x（jxl（i））=dcxb2（i+1,n+1）;
end
fprintf（'可行解为'）;
x
minf=-Inf
cpt=[]
flag=0
return
end
min1=Inf;
fori=1:
m
ifduilie（i）>0&youduanb1（i）duilie（i）min1=youduanb1（i）duilie（i）;
line=i+1;%记录行数
end
end
dcxb2（line,:
）=dcxb2（line,:
）dcxb2（line,thirdnum）;
fori=1:
m+1
ifi~=line
dcxb2（i,:
）=dcxb2（i,:
）+（-1*dcxb2（i,thirdnum）*dcxb2（line,:
））;
end
end
jxl（line-1）=thirdnum;
dcxb2
look2=find（dcxb2（1,[1:
n]）>0）;
end
minf=dcxb2（1,n+1）;
x=zeros（n,1）;
[p,q]=size（jxl）;
fprintf（'\最优解已找到n'）;
fori=1:
q
x（jxl（i））=dcxb2（i+1,n+1）;
end
fprintf（'最优可行解为:
'）;
x
fprintf（'最优值为:
'）;
minf
cpt=dcxb2;
fprintf（'最优解对应的单纯形表为:
'）;
cpt
flag=1
return
例题1.
A=[12112-23;320340];
b=[2;3];
c=[4030];
运行结果：
>>[x,minf,flag,cpt]=linp（A,b,c）
请一次输入系数矩阵A;输入右端向量b;输入所求问题的向量c
A=[12112-23;320340];
b=[2;3];
c=[4030];
dcxb=
-4.00000-3.00000000
2.00001.00001.2500-0.6667005.0000
0.50001.00000.5000-0.66671.000002.0000
1.500000.7500001.00003.0000
dcxb=
00-1.0000002.66678.0000
01.00000.2500-0.66670-1.33331.0000
01.00000.2500-0.66671.0000-0.33331.0000
1.000000.5000000.66672.0000
dcxb=
00-1.0000002.66678.0000
0000-1.0000-1.00000
01.00000.2500-0.66671.0000-0.33331.0000
1.000000.5000000.66672.0000
最优解已找到!
最优可行解为:
x=
2
1
0
0
最优值为:
minf=
8
最优解对应的单纯形表为:
cpt=
00-1.000008.0000
01.00000.2500-0.66671.0000
1.000000.500002.0000
flag=
1
ans=
2
1
0
0
例题2
A=[12112-23;320340;3-604];
b=[2;3;0];
c=[4030];
运行结果
>>[x,minf,flag,cpt]=linp（A,b,c）
请一次输入系数矩阵A;输入右端向量b;输入所求问题的向量c
A=[12112-23;320340;3-604];
b=[2;3;0];
c=[4030];
dcxb=
Columns1through6
-4.00000-3.0000000
5.0000-5.00001.25003.333300
0.50001.00000.5000-0.66671.00000
1.500000.7500001.0000
3.0000-6.000004.000000
Columns7through8
00
05.0000
02.0000
03.0000
1.00000
dcxb=
Columns1through6
0-8.0000-3.00005.333300
05.00001.2500-3.333300
02.00000.5000-1.33331.00000
03.00000.7500-2.000001.0000
1.0000-2.000001.333300
Columns7through8
1.33330
-1.66675.0000
-0.16672.0000
-0.50003.0000
0.33330
dcxb=
Columns1through6
00-1.000004.00000
0000.0000-2.50000
01.00000.2500-0.66670.50000
0000-1.50001.0000
1.000000.500001.00000
Columns7through8
0.66678.0000
-1.25000
-0.08331.0000
-0.25000
0.16672.0000
dcxb=
Columns1through6
00-1.000004.00000
0000-2.50000
01.00000.2500-0.66670.50000
0000-1.50001.0000
1.000000.500001.00000
Columns7through8
0.66678.0000
-1.25000
-0.08331.0000
-0.25000
0.16672.0000
最优解已找到!
最优可行解为:
x=
2
1
0
0
最优值为:
minf=
8
最优解对应的单纯形表为:
cpt=
00-1.000008.0000
01.00000.2500-0.66671.0000
1.000000.500002.0000
flag=
1
ans=
2
1
0
0
例题3
A=[12112-23;320120;3-604];
b=[2;3;0];
c=[4030];
运行结果:
>>[x,minf,flag,cpt]=linp（A,b,c）
请一次输入系数矩阵A;输入右端向量b;输入所求问题的向量c
A=[12112-23;320120;3-604];
b=[2;3;0];
c=[4030];
dcxb=
Columns1through6
-4.00000-3.0000000
5.0000-5.00001.00003.333300
0.50001.00000.5000-0.66671.00000
1.500000.5000001.0000
3.0000-6.000004.000000
Columns7through8
00
05.0000
02.0000
03.0000
1.00000
dcxb=
Columns1through6
0-8.0000-3.00005.333300
05.00001.0000-3.333300
02.00000.5000-1.33331.00000
03.00000.5000-2.000001.0000
1.0000-2.000001.333300
Columns7through8
1.33330
-1.66675.0000
-0.16672.0000
-0.50003.0000
0.33330
dcxb=
Columns1through6
00-1.000004.00000
00-0.25000.0000-2.50000
01.00000.2500-0.66670.50000
00-0.25000-1.50001.0000
1.000000.500001.00000
Columns7through8
0.66678.0000
-1.25000
-0.08331.0000
-0.25000
0.16672.0000
dcxb=
Columns1through6
00-1.000004.00000
00-0.25000-2.50000
01.00000.2500-0.66670.50000
00-0.25000-1.50001.0000
1.000000.500001.00000
Columns7through8
0.66678.0000
-1.25000
-0.08331.0000
-0.25000
0.16672.0000
最优解已找到!
最优可行解为:
x=
2
1
0
0
最优值为:
minf=
8
最优解对应的单纯形表为:
cpt=
00008.0000
01.00000-0.66671.0000
001.000000
1.00000002.0000
flag=
1
ans=
2
1
0
0
例题4
A=[-11-10;-1-10-1];
b=[1;2];
c=[2200];
运行结果
>>[x,minf,flag,cpt]=linp（A,b,c）
请一次输入系数矩阵A;输入右端向量b;输入所求问题的向量c
A=[-11-10;-1-10-1];
b=[1;2];
c=[2200];
dcxb=
-2-200000
-20-1-1003
-11-10101
-1-10-1012
g>0,此问题没有可行解
ans=
[]
例题5
A=[-10100;01010;-12001];
b=[4;3;8];
c=[-2-5000];
运行结果:
>>[x,minf,flag,cpt]=linp（A,b,c）
请一次输入系数矩阵A;输入右端向量b;输入所求问题的向量c
A=[-10100;01010;-12001];
b=[4;3;8];
c=[-2-5000];
dcxb=
250000000
-2311100015
-101001004
010100103
-120010018
dcxb=
200-500-50-15
-201-210-306
-101001004
010100103
-100-210-212
dcxb=
200-500-50-15
-100-21-1-302
-101001004
010100103
-100-210-212
dcxb=
200-500-50-15
00000-1-1-10
-101001004
010100103
-100-210-212
此问题有可行解，但是没有最优解可行解为：
x=
0
3
4
0
2
minf=
-Inf
cpt=
[]
flag=
0
ans=
0
3
4
0
2
升级会员 