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

Four short words sum up what has lifted most successful individuals above the crowd: a little bit more.-author-date单纯形法C语言程序.单纯形法完全c语言程序,能运行实验：编制线性规划计算程序一、实验目的：(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); %初始化Xif(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''./(BA(:,%d)为',k); b'./(BA(:,k) if(z<0.000000001) flag=0; %所有判别数都小于0时达到最优解。 fprintf(' 已找到最优解!n'); xB=(Bb')' 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(BA(:,k)<=0) % 如果BA(;,k)中的每一个分量都小于零。 flag=0; fprintf(' n 此问题不存在最优解!n'); %若BA(:,k)的第k列均不大于0，则该问题不存在最优解。 else b1=Bb' 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 endend程序保存为 danchunxin.m文件四、 数值实验及其结果： 1) 窗口输入：run danchunxinA=-1 2 1 0 0;2 3 0 1 0;1 -1 0 0 1b=4 12 3c=-4 -1 0 0 0运行后的结果为：panbieshu = 4 1 0 0 0 b'./(BA(:,2)为;ans = -4 6 3 x(1)进基，x(5)退基panbieshu = 0 5 0 0 -4 b'./(BA(:,2)为;ans = 4 12/5 -3 x(2)进基，x(4)退基panbieshu = 0 0 0 -1 -2 b'./(BA(:,2)为;ans = 1/0 1/0 3 已找到最优解!基向量为:X = 21/5 6/5 29/5 0 0 目标函数值为:f = -18 2)S.t 窗口输入：run danchunxinA=2 3 1 0;-1 1 0 1b=6 1c=-1 -3 0 0运行后的结果为：panbieshu = 1 3 0 0 b'./(BA(:,2)为;ans = 2 1 x(2)进基，x(4)退基panbieshu = 4 0 0 -3 b'./(BA(:,2)为;ans = 6/5 -1 x(1)进基，x(3)退基panbieshu = 0 0 -4/5 -3/5 b'./(BA(:,2)为;ans = 6 1/0 已找到最优解!基向量为:X = 3/5 8/5 0 0 目标函数值为:f = -27/5 3) S.t 窗口输入>> run danchunxinA=3 3 1 0 0;4 -4 0 1 0;2 -1 0 0 1b=30 16 12c=-3 -1 0 0 0运行后的结果为：panbieshu = 3 1 0 0 0 b'./(BA(:,2)为;ans = 10 4 6 x(1)进基，x(4)退基panbieshu = 0 4 0 -3/4 0 b'./(BA(:,2)为;ans = 5 -16 12 x(2)进基，x(3)退基panbieshu = 0 0 -2/3 -1/4 0 b'./(BA(:,2)为;ans = -1/0 16 1/0 已找到最优解!基向量为:X = 7 3 0 0 1 目标函数值为:f = -24 4) 窗口输入：>> run danchunxinA=1 -1 1 1 0;-2 1 -2 0 1b=5 10c=-3 1 0 0 0运行后的结果为：panbieshu = 3 -1 0 0 0 b'./(BA(:,2)为;ans = 5 -5 x(1)进基，x(4)退基panbieshu = 0 2 -3 -3 0 b'./(BA(:,2)为;ans = -5 -10 此问题不存在最优解!五：心得体会：通过本次实验对单纯形了解更深刻，此次实验中inf表示为一个无穷大的数。本次做的只是最简单的线性规划问题，面对以后更大的、更复杂的问题，虽然起不了什么非常大的作用，但这是基础，所以我非常认真对待这次实验，做完本次实验，使我对单纯形方法，更加熟练，对matlab程序设计也更加熟悉。单纯形法完全c语言程序,能运行#include "math.h"#include "stdio.h"#define N 2void 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>pi+1) d=p; p=pi+1; pi+1=d; m=i;j=k+1; k=0;for (i=m; i>=j; i-)if (pi-1>p) d=p; p=pi-1; pi-1=d; k=i;return;double mubiao(double *x) double y;y=x1-x0*x0; y=100.0*y*y;y=y+(1.0-x0)*(1.0-x0);return(y); main() int i,j,k,l,m=0;double c,xxN+1N,f0N+1,fN+1,x0N=1.2,1,x1N,s=0.0;double a,b;double xaN,xbN,xcN,xeN,xwN,xrN,xoN;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");/给xxNN+1赋值,每一行构成单纯形的一个定点/*for(i=0;i<N;i+)xx0=x0;for(i=1;i<N+1;i+)for(j=0;j<N;j+)if(j=i-1)xxj=x0j+a;else xxj=x0j+b;for (i=0;i<N+1;i+)for (j=0;j<N;j+)printf("xx%d%d%13.7e ",i,j,xxj);printf("n");loop1:/求单纯形的每个定点的函数值f0,f和x1是过渡数组printf("n");printf("n");for(i=0;i<N+1;i+)for(j=0;j<N;j+)x1j=xxj;f0=mubiao(x1);f=mubiao(x1);printf("f0%d=%13.7e f%d=%13.7en",i,f0,i,f);printf("n");/比较f的大小，f0是最小值，fN是最大值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=f0)k=i;if(f0=fN)l=i;printf("最好点k=%dn",k);printf("最坏点l=%dn",l);/终止判断条件printf("fN-f0=%13.7e n",fN-f0);if(fN-f0)<eps1+eps2*fabs(fN)printf("迭代次数m=%dn",m);for(j=0;j<N;j+)printf("optx%d=%13.7en",j,xxkj);printf("fmin=%13.7en",f0);elsem=m+1;/把xx中最好点移到第一行，最坏点移到最后一行for(j=0;j<N;j+)xbj=xxkj;xxkj=xx0j;xx0j=xbj;/xwj=xxlj;xxlj=xxNj;xxNj=xwj;for (i=0;i<N+1;i+)for (j=0;j<N;j+)printf("xx%d%d=%13.7e ",i,j,xxj);printf("n");/求除最坏点fN外其余点的中点xcfor(i=0;i<N;i+)xa=0;for(j=0;j<N;j+)for(i=0;i<N;i+)xaj=xaj+xxj;xaj=xaj/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);/求xwN的反射点xrN;for(i=0;i<N;i+)xr=xa+aef*(xa-xw);printf("xr%d=%13.7e ",i,xr);printf("n");/求xrN的函数值frfr=mubiao(xr);printf("fr=%13.7e n",fr);/判断xr与xb的好坏if(fr<=f0)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<=f0)for(j=0;j<N;j+)xxNj=xej;elsefor(j=0;j<N;j+)xxNj=xrj;goto loop1;elsefw=fN;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+)xxj=xx0j-rou*(xxj-xx0j);goto loop1;elsefor(j=0;j<N;j+)xxNj=xcj;goto loop1;elseif(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+)xxj=xx0j-rou*(xxj-xx0j);goto loop1;elsefor(j=0;j<N;j+)xxNj=xoj;goto loop1;elsefor(j=0;j<N;j+)xxNj=xrj;goto loop1;-