单纯形法C语言程序..doc

______________________________________________________________________________________________________________
实验：编制《线性规划》计算程序
一、实验目的：
(1)使学生在程序设计方面得到进一步的训练，掌握Matlab (C或VB)语言进行程序设计中一些常用方法。
(2)使学生对线性规划的单纯形法有更深的理解.
二、实验用仪器设备、器材或软件环境
计算机, Matlab R2009a
三、算法步骤、计算框图、计算程序等
本实验主要编写如下线性规划问题的计算程序:

其中初始可行基为松弛变量对应的列组成.
对于一般标准线性规划问题:

１．求解上述一般标准线性规划的单纯形算法（修正）步骤如下：
对于一般的标准形式线性规划问题(求极小问题),首先给定一个初始基本可行解。设初始基为B,然后执行如下步骤：
(1).解,求得，

(2).计算单纯形乘子w, ,得到,对于非基变量，计算判别数,可直接计算令
,R为非基变量集合
若判别数 ,则得到一个最优基本可行解，运算结束；否则，转到下一步
(3).解,得到；若,即的每个分量均非正数，
则停止计算，问题不存在有限最优解，否则，进行步骤(4).确定下标r,使

;
２、计算框图为：

开始

初始可行基B

是

否
得到最优解

是

否
不存在有限最优解
确定下标r,使得

3．计算程序(Matlab)：
A=input('A=');
b=input('b=');
c=input('c=');
format rat %可以让结果用分数输出
[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:n); %基矩阵对应目标值的c
while flag
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';
for i=1:n
mark=0;
for j=1:m
if (D(j,2)==i)
mark=1;
X(i)=xB(D(j,1)); %利用D找出xB与X之间的关系。。
end
end
if mark==0
X(i)=0; %如果D中没有X(i),则X(i)为非基变量，所以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;
for i=1:m
if ((A(i,k)>0) && (b1(i)/(A(i,k)+eps))<temp )
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)

窗口输入：
run danchunxin
A=[-1 2 1 0 0;2 3 0 1 0;1 -1 0 0 1]
b=[4 12 3]
c=[-4 -1 0 0 0]
运行后的结果为：

panbieshu =

4 1 0 0 0

b'./(B\A(:,2))为;
ans =

-4
6
3

x(1)进基，x(5)退基

panbieshu =

0 5 0 0 -4

b'./(B\A(:,2))为;
ans =

4
12/5
-3

x(2)进基，x(4)退基

panbieshu =

0 0 0 -1 -2

b'./(B\A(:,2))为;
ans =

1/0
1/0
3

已找到最优解!
基向量为:
X =

21/5 6/5 29/5 0 0

目标函数值为:
f =

-18

2)
S.t

窗口输入：
run danchunxin
A=[2 3 1 0;-1 1 0 1]
b=[6 1]
c=[-1 -3 0 0]
运行后的结果为：

panbieshu =

1 3 0 0

b'./(B\A(:,2))为;
ans =

2
1

x(2)进基，x(4)退基

panbieshu =

4 0 0 -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/5 8/5 0 0

目标函数值为:
f =

-27/5

3)
S.t

窗口输入
>> run danchunxin
A=[3 3 1 0 0;4 -4 0 1 0;2 -1 0 0 1]
b=[30 16 12]
c=[-3 -1 0 0 0]
运行后的结果为：
panbieshu =

3 1 0 0 0

b'./(B\A(:,2))为;
ans =

10
4
6

x(1)进基，x(4)退基

panbieshu =

0 4 0 -3/4 0

b'./(B\A(:,2))为;
ans =

5
-16
12

x(2)进基，x(3)退基

panbieshu =

0 0 -2/3 -1/4 0

b'./(B\A(:,2))为;
ans =

-1/0
16
1/0

已找到最优解!
基向量为:
X =

7 3 0 0 1

目标函数值为:
f =

-24

4)

窗口输入：
>> run danchunxin
A=[1 -1 1 1 0;-2 1 -2 0 1]
b=[5 10]
c=[-3 1 0 0 0]
运行后的结果为：
panbieshu =

3 -1 0 0 0

b'./(B\A(:,2))为;
ans =

5
-5

x(1)进基，x(4)退基

panbieshu =

0 2 -3 -3 0

b'./(B\A(:,2))为;
ans =

-5
-10

此问题不存在最优解!

五：心得体会：
通过本次实验对单纯形了解更深刻，此次实验中inf表示为一个无穷大的数。
本次做的只是最简单的线性规划问题，面对以后更大的、更复杂的问题，虽然起不了什么非常大的作用，但这是基础，所以我非常认真对待这次实验，做完本次实验，使我对单纯形方法，更加熟练，对matlab程序设计也更加熟悉。

单纯形法完全c语言程序,能运行
#include "math.h"
#include "stdio.h"
#define N 2

void paixu(p,n)
int n;
double p[];
{ int m,k,j,i;
double d;
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>=j; i--)
if (p[i-1]>p)
{ d=p; p=p[i-1]; p[i-1]=d; k=i;}
}
return;
}
double mubiao(double *x)
{ double y;
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()
{ int i,j,k,l,m=0;
double c,xx[N+1][N],f0[N+1],f[N+1],x0[N]={1.2,1},x1[N],s=0.0;
double a,b;
double xa[N],xb[N],xc[N],xe[N],xw[N],xr[N],xo[N];
double fr,fe,fw,fc,fo;
double aef=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.7e b=%13.7e ",a,b);
//printf("\n");
//给xx[N][N+1]赋值,每一行构成单纯形的一个定点
//***********************
for(i=0;i<N;i++)
xx[0]=x0;
for(i=1;i<N+1;i++)
for(j=0;j<N;j++)
{if(j==i-1)
xx[j]=x0[j]+a;
else
xx[j]=x0[j]+b;
}
for (i=0;i<N+1;i++)
{for (j=0;j<N;j++)
printf("xx[%d][%d]%13.7e ",i,j,xx[j]);
printf("\n");
}
loop1:
//求单纯形的每个定点的函数值f0,f和x1是过渡数组
printf("\n");
printf("\n");
for(i=0;i<N+1;i++)
{for(j=0;j<N;j++)
x1[j]=xx[j];
f0=mubiao(x1);
f=mubiao(x1);
printf("f0[%d]=%13.7e f[%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<N+1;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])<eps1+eps2*fabs(f[N]))
{printf("迭代次数m=%d\n",m);
for(j=0;j<N;j++)
printf("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<N;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<N+1;i++)
{for (j=0;j<N;j++)
printf("xx[%d][%d]=%13.7e ",i,j,xx[j]);
printf("\n");
}
//求除最坏点f[N]外其余点的中点xc[]
for(i=0;i<N;i++)
xa=0;
for(j=0;j<N;j++)
{{for(i=0;i<N;i++)
xa[j]=xa[j]+xx[j];}
xa[j]=xa[j]/N;
}
for(i=0;i<N;i++)
printf("xa[%d]=%13.7e xb[%d]=%13.7e xw[%d]=%13.7e \n",i,xa,i,xb,i,xw);
//求xw[N]的反射点xr[N];
for(i=0;i<N;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<N;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;j<N;j++)
xx[N][j]=xe[j];
else
for(j=0;j<N;j++)
xx[N][j]=xr[j];
goto loop1;
}
else
{
fw=f[N];
if(fr>=fw)
{for(i=0;i<N;i++)
xc=xa-bt*(xa-xw);
fc=mubiao(xc);
if(fc>=fw)
{for(i=1;i<N+1;i++)
for(j=0;j<N;j++)
xx[j]=xx[0][j]-rou*(xx[j]-xx[0][j]);
goto loop1;
}
else
{for(j=0;j<N;j++)
xx[N][j]=xc[j];
goto loop1;
}
}
else
{if(fr>=fe)
{
for(i=0;i<N;i++)
xo=xa+bt*(xa-xw);
fo=mubiao(xo);
if(fo>=fr)
{for(i=1;i<N+1;i++)
for(j=0;j<N;j++)
xx[j]=xx[0][j]-rou*(xx[j]-xx[0][j]);
goto loop1;
}
else
{for(j=0;j<N;j++)
xx[N][j]=xo[j];
goto loop1;
}
}
else
{for(j=0;j<N;j++)
xx[N][j]=xr[j];
goto loop1;
}
}
}
}
}

Welcome To
Download !!!

欢迎您的下载，资料仅供参考！

精品资料