单纯形法C语言程序.docx
《单纯形法C语言程序.docx》由会员分享，可在线阅读，更多相关《单纯形法C语言程序.docx（16页珍藏版）》请在冰豆网上搜索。
单纯形法C语言程序
单纯形法C语言程序
实验:
编制《线性规划》计算程序
一、实验目的:
（1）使学生在程序设计方面得到进一步的训练;，掌握Matlab（C或VB）语言进行程序设计中一些常用方法。
（2）使学生对线性规划的单纯形法有更深的理解.
二、实验用仪器设备、器材或软件环境
计算机,MatlabR2009a
三、算法步骤、计算框图、计算程序等
本实验主要编写如下线性规划问题的计算程序:
mincx
Ax,b,s.t.,x,0,b,0,
其中初始可行基为松弛变量对应的列组成.
对于一般标准线性规划问题:
mincx
Ax,b,s.t.,x,0,b,0,
（求解上述一般标准线性规划的单纯形算法（修正）步骤如下:
对于一般的标准形式线性规划问题（求极小问题）,首先给定一个初始基本可行解。
设初始基为B,然后执行如下步骤:
1Bxb,xBb,令计算目标函数值xfcx,,0,BNBBB
（1）.解,求得，
1以b（1,2,...,）imBbi,记的第个分量i
1wBC,wCB,BB
（2）.计算单纯形乘子w,,得到,对于非基变量，计算判别
1,1,Ac,数,可直接计算令,,,,,zccBpciiiBii,cBB
R为非基变量集合,,max{}k,,iR
,0k若判别数,则得到一个最优基本可行解，运算结束;否则，转到下一步
1yy,0Byp,yBp,kkkkkk（3）.解,得到;若,即的每个分量均非正数，
则停止计算，问题不存在有限最优解，否则，进行步骤（4）.确定下标r,使
bbtr,,min,0且y,,x为离基变量,tkyyrktkB:
0,tyrtk
xp为进基变量,用p替换得到新的基矩阵,B,还回步骤
（1）kBkr;
1
2
、计算框图为:
开始
初始可行基B
1令x,Bb,x,0,f,cxBNBB
1,,计算单纯性乘子w,cB,计算判别数,wp,c,j,R（非基变量）Bjjj
令,,max{,,j,R}kj
是
,0?
k
否
得到最优
解
1解方程By,p,得到y,Bp,kkkk
是
y,0?
k
否
不存在有限
最优解
确定下标r,使得
,bbir,min|,0y,,ikyyrkik,,
x为退基变量，x进基变量，以Bkr
p代替p,得到新的基矩阵BkBr
3
:
3（计算程序（Matlab）
A=input（'A='）;
b=input（'b='）;
c=input（'c='）;
formatrat%可以让结果用分数输出
[m,n]=size（A）;
E=1:
m;E=E';
F=n-m+1:
n;F=F';
D=[E,F];%创建一个一一映射，为了结果能够标准输出X=zeros（1,n）;%初始化X
if（nfprintf（'不符合要求需引入松弛变量'）
flag=0;
else
flag=1;
B=A（:
n-m+1:
n）;%找基矩阵
cB=c（n-m+1:
n）;%基矩阵对应目标值的c
whileflag
w=cB/B;%计算单纯形乘子，cB/B=cB*inv（B）,用cB/B的目的是，为
了提高运行速度。
。
panbieshu=w*A-c%计算判别数,后面没有加分号，就是为了计算后能够显示出
来。
。
[z,k]=max（panbieshu）;%k作为进基变量下标。
。
fprintf（'b''./（B\\A（:
%d））为',k）;
b'./（B\A（:
k））
if（z<0.000000001）
flag=0;%所有判别数都小于0时达到最优解。
。
fprintf（'已找到最优解!
\n'）;
xB=（B\b'）';
f=cB*xB';
fori=1:
n
mark=0;
forj=1:
m
if（D（j,2）==i）
mark=1;
X（i）=xB（D（j,1））;%利用D找出xB与X之间的关系。
。
end
end
ifmark==0
X（i）=0;%如果D中没有X（i）,则X（i）为非基变量，所以
4
X（i）,0。
。
end
end
fprintf（'基向量为:
'）;X
fprintf（'目标函数值为:
'）;f
else
if（B\A（:
k）<=0）%如果B\A（;,k）中的每一个分量都小于
零。
。
flag=0;
fprintf（'\n此问题不存在最优解!
\n'）;%若B\A（:
k）
的第k列均不大于0，则该问题不存在最优解。
。
else
b1=B\b';
temp=inf;
fori=1:
m
if（（A（i,k）>0）&&（b1（i）/（A（i,k）+eps））temp=b1（i）/A（i,k）;%找退基变量
r=i;
end
end
fprintf（'x（%d）进基，x（%d）退基\n',k,D（r,2））;%显示进基
变量和退基变量
B（:
r）=A（:
k）;
cB（r）=c（k）;%确定进基退基变量后，相应的基矩阵及新基对应的目
标值的c也相应改变
D（r,2）=k;%改变D中的映射关系
end
end
end
end
程序保存为danchunxin.m文件
四、数值实验及其结果:
,4min1）xx12
st..,，，,24xxx123
2312，，,xxx124
，,3xxx125
5
x,0,j,1,...,5.j
窗口输入:
rundanchunxinA=[-12100;23010;1-1001]
b=[4123]
c=[-4-1000]运行后的结果为:
panbieshu=
41000
b'./（B\A（:
2））为;
ans=
-4
6
3
x
（1）进基，x（5）退基
panbieshu=
0500-4
b'./（B\A（:
2））为;
ans=
4
12/5
-3
x
（2）进基，x（4）退基
panbieshu=
000-1-2
b'./（B\A（:
2））为;
ans=
1/0
6
1/0
3
已找到最优解!
基向量为:
X=
21/56/529/500
目标函数值为:
f=
-18
minx，3x2）12
2x，3x，x,6S.t123
x，x，x,1124
x,0,j,1,...,4.j
窗口输入:
rundanchunxinA=[2310;-1101]
b=[61]
c=[-1-300]运行后的结果为:
panbieshu=
1300
b'./（B\A（:
2））为;
ans=
2
1
x
（2）进基，x（4）退基
panbieshu=
7
400-3
b'./（B\A（:
2））为;
ans=
6/5
-1
x
（1）进基，x（3）退基
panbieshu=
00-4/5-3/5
b'./（B\A（:
2））为;
ans=
6
1/0
已找到最优解!
基向量为:
X=
3/58/500
目标函数值为:
f=
-27/5
3x,xmin3）12
3x，3x，x,30S.t123
4x,4x，x,16124
2x,x,1212
x,0,j,1,...,4.j
窗口输入
>>rundanchunxinA=[33100;4-4010;2-1001]
8
b=[301612]c=[-3-1000]运行后的结果为:
panbieshu=
31000
b'./（B\A（:
2））为;
ans=
10
4
6
x
（1）进基，x（4）退基
panbieshu=
040-3/40
b'./（B\A（:
2））为;
ans=
5
-16
12
x
（2）进基，x（3）退基
panbieshu=
00-2/3-1/40
b'./（B\A（:
2））为;
ans=
-1/0
16
1/0
已找到最优解!
基向量为:
X=
9
73001
目标函数值为:
f=
-24
xxx,，2min4）123
xxxx，,，,210st..1234
248xxx,，,123
，,,xxx244123
xj,,0,1,...,4j
窗口输入:
>>rundanchunxinA=[1-1110;-21-201]
b=[510]
c=[-31000]
运行后的结果为:
panbieshu=
3-1000
b'./（B\A（:
2））为;
ans=
5
-5
x
（1）进基，x（4）退基
panbieshu=
02-3-30
b'./（B\A（:
2））为;
ans=
-5
-10
10
此问题不存在最优解!
五:
心得体会:
通过本次实验对单纯形了解更深刻，此次实验中inf表示为一个无穷大的数。
本次做的只是最简单的线性规划问题，面对以后更大的、更复杂的问题，虽然起不了什么非常大的作用，但这是基础，所以我非常认真对待这次实验，做完本次实验，使我对单纯形方法，更加熟练，对matlab程序设计也更加熟悉。
单纯形法完全c语言程序,能运行
#include"math.h"
#include"stdio.h"
#defineN2
voidpaixu（p,n）
intn;
doublep[];
{intm,k,j,i;
doubled;
k=0;m=n-1;
while（k{j=m-1;m=0;
for（i=k;i<=j;i++）
if（p>p[i+1]）
{d=p;p=p[i+1];p[i+1]=d;m=i;}j=k+1;k=0;
for（i=m;i>=j;i--）
if（p[i-1]>p）
{d=p;p=p[i-1];p[i-1]=d;k=i;}
}
return;
}
doublemubiao（double*x）
{doubley;
y=x[1]-x[0]*x[0];y=100.0*y*y;
11
y=y+（1.0-x[0]）*（1.0-x[0]）;
return（y）;
}
main（）
{inti,j,k,l,m=0;
doublec,xx[N+1][N],f0[N+1],f[N+1],x0[N]={1.2,1},x1[N],s=0.0;
doublea,b;
doublexa[N],xb[N],xc[N],xe[N],xw[N],xr[N],xo[N];
doublefr,fe,fw,fc,fo;
doubleaef=1.0,r=1.0,eps1=1.0e-30,eps2=1.0e-30,bt=0.5,rou=0.5;
c=1.0;
b=（c/（N*sqrt
（2）））*（sqrt（N+1）-1）;
a=b+c/sqrt
（2）;
//printf（"a=%13.7eb=%13.7e",a,b）;
//printf（"\n"）;
//给xx[N][N+1]赋值,每一行构成单纯形的一个定点
//***********************
for（i=0;ixx[0]=x0;
for（i=1;ifor（j=0;j{if（j==i-1）
xx[j]=x0[j]+a;
else
xx[j]=x0[j]+b;
}
for（i=0;i{for（j=0;jprintf（"xx[%d][%d]%13.7e",i,j,xx[j]）;
printf（"\n"）;
}
loop1:
//求单纯形的每个定点的函数值f0,f和x1是过渡数组
printf（"\n"）;
printf（"\n"）;
for（i=0;i{for（j=0;jx1[j]=xx[j];
f0=mubiao（x1）;
f=mubiao（x1）;
printf（"f0[%d]=%13.7ef[%d]=%13.7e\n",i,f0,i,f）;
}12printf（"\n"）;
//比较f的大小，f[0]是最小值，f[N]是最大值
paixu（f,N+1）;
for（i=N;i>=0;i--）
printf（"f[%d]=%13.7e\n",i,f）;//找最好点和最坏点分别是哪一个点，即xx[][]的行数
for（i=0;i{if（f0==f[0]）
k=i;
if（f0==f[N]）
l=i;
}
printf（"最好点k=%d\n",k）;
printf（"最坏点l=%d\n",l）;
//终止判断条件
printf（"f[N]-f[0]=%13.7e\n",f[N]-f[0]）;
if（（f[N]-f[0]）{printf（"迭代次数m=%d\n",m）;
for（j=0;jprintf（"optx[%d]=%13.7e\n",j,xx[k][j]）;
printf（"fmin=%13.7e\n",f[0]）;
}
else
{
m=m+1;
//把xx[][]中最好点移到第一行，最坏点移到最后一行
for（j=0;j{xb[j]=xx[k][j];
xx[k][j]=xx[0][j];
xx[0][j]=xb[j];
//
xw[j]=xx[l][j];
xx[l][j]=xx[N][j];
xx[N][j]=xw[j];
}
for（i=0;i{for（j=0;jprintf（"xx[%d][%d]=%13.7e",i,j,xx[j]）;
printf（"\n"）;
}
//求除最坏点f[N]外其余点的中点xc[]
for（i=0;ixa=0;
for（j=0;jxa[j]=xa[j]+xx[j];}
xa[j]=xa[j]/N;
}
for（i=0;iprintf（"xa[%d]=%13.7exb[%d]=%13.7exw[%d]=%13.7e\n",i,xa,i,xb,i,xw）;
//求xw[N]的反射点xr[N];
for（i=0;i{xr=xa+aef*（xa-xw）;
printf（"xr[%d]=%13.7e",i,xr）;}
printf（"\n"）;
//求xr[N]的函数值fr
fr=mubiao（xr）;
printf（"fr=%13.7e\n",fr）;
//判断xr与xb的好坏
if（fr<=f[0]）
{for（i=0;i{xe=xr+r*（xr-xa）;
//printf（"xe[%d]=%13.7e",i,xe）;
}
printf（"\n"）;
fe=mubiao（xe）;
if（fe<=f[0]）
for（j=0;jxx[N][j]=xe[j];
else
for（j=0;jxx[N][j]=xr[j];
gotoloop1;
}
else
{
fw=f[N];
if（fr>=fw）
{for（i=0;ixc=xa-bt*（xa-xw）;
fc=mubiao（xc）;
if（fc>=fw）
{for（i=1;ifor（j=0;jxx[j]=xx[0][j]-rou*（xx[j]-xx[0][j]）;
gotoloop1;
}14else
{for（j=0;jxx[N][j]=xc[j];
gotoloop1;
}
}
else
{if（fr>=fe）
{
for（i=0;ixo=xa+bt*（xa-xw）;
fo=mubiao（xo）;
if（fo>=fr）
{for（i=1;ifor（j=0;jxx[j]=xx[0][j]-rou*（xx[j]-xx[0][j]）;
gotoloop1;
}
else
{for（j=0;jxx[N][j]=xo[j];
gotoloop1;
}
}
else
{for（j=0;jxx[N][j]=xr[j];
gotoloop1;
}
}
}
}
}
15
升级会员 