单纯形法C语言程序文档格式.docx
《单纯形法C语言程序文档格式.docx》由会员分享，可在线阅读，更多相关《单纯形法C语言程序文档格式.docx（16页珍藏版）》请在冰点文库上搜索。
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（n<
m）%判断是否为标准型
fprintf（'
不符合要求需引入松弛变量'
）
flag=0;
else
flag=1;
B=A（:
n-m+1:
n）;
%找基矩阵
cB=c（n-m+1:
%基矩阵对应目标值的c
whileflag
w=cB/B;
%计算单纯形乘子，cB/B=cB*inv（B）,用cB/B的目的是，为
了提高运行速度。
。
panbieshu=w*A-c%计算判别数,后面没有加分号，就是为了计算后能够显示出
来。
[z,k]=max（panbieshu）;
%k作为进基变量下标。
b'
'
./（B\\A（:
%d））为'
k）;
./（B\A（:
k））
if（z<
0.000000001）
%所有判别数都小于0时达到最优解。
已找到最优解!
\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
ifmark==0
X（i）=0;
%如果D中没有X（i）,则X（i）为非基变量，所以
4
X（i）,0。
基向量为:
X
目标函数值为:
）;
f
if（B\A（:
k）<
=0）%如果B\A（;
k）中的每一个分量都小于
零。
\n此问题不存在最优解!
%若B\A（:
k）
的第k列均不大于0，则该问题不存在最优解。
b1=B\b'
temp=inf;
if（（A（i,k）>
0）&
&
（b1（i）/（A（i,k）+eps））<
temp）
temp=b1（i）/A（i,k）;
%找退基变量
r=i;
x（%d）进基，x（%d）退基\n'
k,D（r,2））;
%显示进基
变量和退基变量
B（:
r）=A（:
cB（r）=c（k）;
%确定进基退基变量后，相应的基矩阵及新基对应的目
标值的c也相应改变
D（r,2）=k;
%改变D中的映射关系
程序保存为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
2））为;
ans=
-4
6
x
（1）进基，x（5）退基
0500-4
12/5
-3
x
（2）进基，x（4）退基
000-1-2
1/0
已找到最优解!
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]运行后的结果为:
1300
7
400-3
6/5
-1
x
（1）进基，x（3）退基
00-4/5-3/5
3/58/500
-27/5
3x,xmin3）12
3x，3x，x,30S.t123
4x,4x，x,16124
2x,x,1212
窗口输入
>
rundanchunxinA=[33100;
4-4010;
2-1001]
8
b=[301612]c=[-3-1000]运行后的结果为:
31000
10
x
（1）进基，x（4）退基
040-3/40
-16
12
x
（2）进基，x（3）退基
00-2/3-1/40
-1/0
16
9
73001
-24
xxx,，2min4）123
xxxx，,，,210st..1234
248xxx,，,123
，,,xxx244123
xj,,0,1,...,4j
rundanchunxinA=[1-1110;
-21-201]
b=[510]
c=[-31000]
运行后的结果为:
3-1000
-5
02-3-30
-10
此问题不存在最优解!
五:
心得体会:
通过本次实验对单纯形了解更深刻，此次实验中inf表示为一个无穷大的数。
本次做的只是最简单的线性规划问题，面对以后更大的、更复杂的问题，虽然起不了什么非常大的作用，但这是基础，所以我非常认真对待这次实验，做完本次实验，使我对单纯形方法，更加熟练，对matlab程序设计也更加熟悉。
单纯形法完全c语言程序,能运行
#include"
math.h"
stdio.h"
#defineN2
voidpaixu（p,n）
intn;
doublep[];
{intm,k,j,i;
doubled;
k=0;
m=n-1;
while（k<
m）
{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>
i--）
if（p[i-1]>
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）;
\n"
//给xx[N][N+1]赋值,每一行构成单纯形的一个定点
//***********************
for（i=0;
i<
N;
i++）
xx[0]=x0;
for（i=1;
N+1;
for（j=0;
j<
j++）
{if（j==i-1）
xx[j]=x0[j]+a;
xx[j]=x0[j]+b;
for（i=0;
{for（j=0;
printf（"
xx[%d][%d]%13.7e"
i,j,xx[j]）;
loop1:
//求单纯形的每个定点的函数值f0,f和x1是过渡数组
{for（j=0;
x1[j]=xx[j];
f0=mubiao（x1）;
f=mubiao（x1）;
f0[%d]=%13.7ef[%d]=%13.7e\n"
i,f0,i,f）;
}12printf（"
//比较f的大小，f[0]是最小值，f[N]是最大值
paixu（f,N+1）;
for（i=N;
i>
=0;
i--）
f[%d]=%13.7e\n"
i,f）;
//找最好点和最坏点分别是哪一个点，即xx[][]的行数
{if（f0==f[0]）
k=i;
if（f0==f[N]）
l=i;
最好点k=%d\n"
最坏点l=%d\n"
l）;
//终止判断条件
f[N]-f[0]=%13.7e\n"
f[N]-f[0]）;
if（（f[N]-f[0]）<
eps1+eps2*fabs（f[N]））
{printf（"
迭代次数m=%d\n"
m）;
optx[%d]=%13.7e\n"
j,xx[k][j]）;
fmin=%13.7e\n"
f[0]）;
{
m=m+1;
//把xx[][]中最好点移到第一行，最坏点移到最后一行
{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];
xx[%d][%d]=%13.7e"
//求除最坏点f[N]外其余点的中点xc[]
xa=0;
j++）13{{for（i=0;
xa[j]=xa[j]+xx[j];
xa[j]=xa[j]/N;
xa[%d]=%13.7exb[%d]=%13.7exw[%d]=%13.7e\n"
i,xa,i,xb,i,xw）;
//求xw[N]的反射点xr[N];
{xr=xa+aef*（xa-xw）;
xr[%d]=%13.7e"
i,xr）;
}
//求xr[N]的函数值fr
fr=mubiao（xr）;
fr=%13.7e\n"
fr）;
//判断xr与xb的好坏
if（fr<
=f[0]）
{for（i=0;
{xe=xr+r*（xr-xa）;
xe[%d]=%13.7e"
i,xe）;
fe=mubiao（xe）;
if（fe<
xx[N][j]=xe[j];
xx[N][j]=xr[j];
gotoloop1;
fw=f[N];
if（fr>
=fw）
xc=xa-bt*（xa-xw）;
fc=mubiao（xc）;
if（fc>
{for（i=1;
xx[j]=xx[0][j]-rou*（xx[j]-xx[0][j]）;
}14else
xx[N][j]=xc[j];
{if（fr>
=fe）
xo=xa+bt*（xa-xw）;
fo=mubiao（xo）;
if（fo>
=fr）
xx[N][j]=xo[j];
15
升级会员 