单纯形法C语言程序.docx
《单纯形法C语言程序.docx》由会员分享，可在线阅读，更多相关《单纯形法C语言程序.docx（18页珍藏版）》请在冰点文库上搜索。
单纯形法C语言程序
实验：
编制《线性规划》计算程序
一、实验目的：
⑴使学生在程序设计方面得到进一步的训练，掌握Matlab（C或VB）语言进行
程序设计中一些常用方法。
（2）使学生对线性规划的单纯形法有更深的理解.
二、实验用仪器设备、器材或软件环境
计算机,MatlabR2009a
三、算法步骤、计算框图、计算程序等
本实验主要编写如下线性规划问题的计算程序：
minex
”Ax兰b
st.«xZ0,b30
其中初始可行基为松弛变量对应的列组成•
对于一般标准线性规划问题：
minex
-Ax=b
st.*
xK0,bK0
1•求解上述一般标准线性规划的单纯形算法（修正）步骤如下：
对于一般的标准形式线性规划问题（求极小问题）,首先给定一个初始基本可行解。
设初始基为B,然后执行如下步骤：
（1）解Bxb=b求得XB=B°b令Xn=0,计算目标函数值f=qxB
以b（i=1,2,...,m）记B」b的第i个分量
（2）.计算单纯形乘子w,wB=CB,得到w二CbB」，对于非基变量，计算判别
数5=乙-G=cBB」Pi-c,可直接计算石=CbB丄A—c令
=max^},R为非基变量集合
若判别数孤空0,则得到一个最优基本可行解，运算结束；否则，转到下一步
1
（3）.解Byk二pk,得到yk=Bpk；若0,即yk的每个分量均非正数，
则停止计算，问题不存在有限最优解，否则，进行步骤（4）.确定下标r,使
=鲁,且ytjo}xBr为离基变量,
xk为进基变量，用Pk替换PB,得到新的基矩阵B,还回步骤
（1）;
2、计算框图为:
最优解
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）;
if（n%创建一个一一映射，为了结果能够标准输出
%初始化X
%判断是否为标准型
fprintf（'不符合要求需引入松弛变量'）
flag=0;
else
flag=1;
%k作为进基变量下标为',k）;
panbieshu=w*A-c
来。
。
%计算判别数,后面没有加分号，就是为了计算后能够显示出
[z,k]=max（panbieshu）;
fprintf（'b''./（B\\A（:
%d））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）为非基变量，所以
X（i）=0。
。
end
end
flag=O;
fprintf（'\n此问题不存在最优解！
\n'）;%若B\A（:
k）
的第k列均不大于0,则该问题不存在最优解。
。
else
b仁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文件
四、数值实验及其结果:
1）min_4x^X2
s・t_X^2X^X^4
2X13X2x^12
XXj-0,j二1,…，5.
窗口输入:
rundanchunxin
A=[-12100;23010;1-1001]b=[4123]
c=[-4-1000]运行后的结果为：
panbieshu=
4100
b'./（B\A（:
2））为;
ans=
-4
6
3
x
（1）进基，x（5）退基
panbieshu=
-4
-2
0500
b'./（B\A（:
2））为;
ans=
4
12/5
-3
x
（2）进基，x（4）退基
panbieshu=
000-1
b'./（B\A（:
2））为;
ans=
1/0
1/0
3
已找到最优解
基向量为:
21/56/529/500
目标函数值为：
f=
-18
2）minx13x2
s.t2x13x2x36
-x1x2x4=1
Xj-0,j二1,...,4.
窗口输入：
rundanchunxin
A=[2310;-1101]
b=[61]
c=[-1-300]
运行后的结果为：
panbieshu=
1300
b'./（B\A（:
2））为；
ans=
2
1
x
（2）进基，x（4）退基
panbieshu=
400-3
b'./（B\A（:
2））为；
ans=
6/5
-1
x
（1）进基，
x（3）
退基
panbieshu=
0
0
-4/5
-3/5
b'./（B\A（:
2））
为；
ans=
6
1/0
已找到最优解
基向量为：
X=
3/58/50
目标函数值为
-27/5
3）min-3xi-X2
s.t3x13x2x330
4x〔-4x2x4=16
2xXj-0,j二h…，4-
窗口输入
>>rundanchunxin
A=[33100;4-4010;2-1001]
b=[301612]
c=[-3-1000]
运行后的结果为：
panbieshu=
310
b'./（B\A（:
2））为;
ans=
10
4
6
x
（1）进基，x（4）退基panbieshu=
-3/40
-1/40
040
b'./（B\A（:
2））为;
ans=
5
-16
12
x
（2）进基，x（3）退基panbieshu=
00-2/3
b'./（B\A（:
2））为;
ans=
-1/0
16
1/0
已找到最优解!
基向量为:
X=
730目标函数值为:
-24
4）
minx1-2x2X3
stX[+X2-2x3+X4=10
2x^-x24x3乞8
-X2x2-4x3一4
Xj-0,j二1,…,4
窗口输入：
>>rundanchunxin
A=[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
此问题不存在最优解
五：
心得体会:
通过本次实验对单纯形了解更深刻，此次实验中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（kp[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;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;ifor（i=1;ixx[j]=x0[j]+a;
elsexx[j]=x0[j]+b;
}
for（i=0;i{for（j=0;jprintf（"\n"）;
}
loop1:
//求单纯形的每个定点的函数值f0,f和x1是过渡数组printf（"\n"）;
printf（"\n"）;
for（i=0;i{for（j=0;jf0=mubiao（x1）;
f=mubiao（x1）;
printf（"f0[%d]=%13.7ef[%d]=%13.7e\n",i,f0,i,f）;
}
printf（"\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;j{{for（i=0;ixa[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]的函数值frfr=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;ifc=mubiao（xc）;
if（fc>=fw）
{for（i=1;ifor（j=0;jgotoloop1;
}
else
{for（j=0;jxx[N][j]=xc[j];
gotoloop1;
}
}
else
{if（fr>=fe）
{
for（i=0;iif（fo>=fr）
{for（i=1;ifor（j=0;j}
else
{for（j=0;jgotoloop1;
}
}
else
{for（j=0;jgotoloop1;
}
}
}
}
}
升级会员 