单纯形法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(nm) %判断是否为标准型 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(z0.) 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) & (b1(i)/(A(i,k)+eps) 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 (km) j=m-1; m=0;for (i=k; ipi+1) d=p; p=pi+1; pi+1=d; m=i;j=k+1; k=0;for (i=m; i=j; i-)if (pi-1p) 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;iN;i+)xx0=x0;for(i=1;iN+1;i+)for(j=0;jN;j+)if(j=i-1)xxj=x0j+a;else xxj=x0j+b;for (i=0;iN+1;i+)for (j=0;jN;j+)printf(xx%d%d%13.7e ,i,j,xxj);printf(n);loop1:/求单纯形的每个定点的函数值f0,f和x1是过渡数组printf(n);printf(n);for(i=0;iN+1;i+)for(j=0;j=0;i-)printf(f%d=%13.7e n,i,f);/找最好点和最坏点分别是哪一个点，即xx的行数for(i=0;iN+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;jN;j+)printf(optx%d=%13.7en,j,xxkj);printf(fmin=%13.7en,f0);elsem=m+1;/把xx中最好点移到第一行，最坏点移到最后一行for(j=0;jN;j+)xbj=xxkj;xxkj=xx0j;xx0j=xbj;/xwj=xxlj;xxlj=xxNj;xxNj=xwj;for (i=0;iN+1;i+)for (j=0;jN;j+)printf(xx%d%d=%13.7e ,i,j,xxj);printf(n);/求除最坏点fN外其余点的中点xcfor(i=0;iN;i+)xa=0;for(j=0;jN;j+)for(i=0;iN;i+)xaj=xaj+xxj;xaj=xaj/N;for(i=0;iN;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;iN;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;iN;i+)xe=xr+r*(xr-xa);/printf(xe%d=%13.7e ,i,xe);printf(n);fe=mubiao(xe);if(fe=f0)for(j=0;jN;j+)xxNj=xej;elsefor(j=0;j=fw)for(i=0;i=fw)for(i=1;iN+1;i+)for(j=0;jN;j+)xxj=xx0j-rou*