单纯形法C语言程序.docx

1、单纯形法C语言程序单纯形法C语言程序实验:编制线性规划计算程序 一、实验目的: (1)使学生在程序设计方面得到进一步的训练;，掌握Matlab (C或VB)语言进行程序设计中一些常用方法。 (2)使学生对线性规划的单纯形法有更深的理解. 二、实验用仪器设备、器材或软件环境 计算机, Matlab R2009a 三、算法步骤、计算框图、计算程序等 本实验主要编写如下线性规划问题的计算程序: mincxAx,b ,s.t.,x,0,b,0,其中初始可行基为松弛变量对应的列组成. 对于一般标准线性规划问题: mincxAx,b ,s.t.,x,0,b,0,(求解上述一般标准线性规划的单纯形算法(修正

2、)步骤如下: 对于一般的标准形式线性规划问题(求极小问题),首先给定一个初始基本可行解。设初始基为B,然后执行如下步骤: ,1Bxb,xBb,令计算目标函数值xfcx,0,BNBBB(1).解,求得， ,1以b(1,2,.,)imBbi,记的第个分量i ,1wBC,wCB,BB(2).计算单纯形乘子w, ,得到,对于非基变量，计算判别,1,1,Ac,数,可直接计算令 ,zccBpciiiBii,cBB,R为非基变量集合 ,maxk,iR,0k若判别数 ,则得到一个最优基本可行解，运算结束;否则，转到下一步 ,1yy,0Byp,yBp,kkkkkk(3).解,得到;若,即的每个分量均非正数， 则

3、停止计算，问题不存在有限最优解，否则，进行步骤(4).确定下标r,使 bbtr,min,0且y,x为离基变量,tkyyrktkB:0,tyrtk xp为进基变量,用p替换得到新的基矩阵,B,还回步骤(1)kBkr; 1 2 ,、计算框图为: 开始 初始可行基B,1 令x,Bb,x,0,f,cxBNBB1 ,计算单纯性乘子w,cB,计算判别数,wp,c,j,R(非基变量)Bjjj令,max,j,Rkj是 ,0?k否 得到最优解 ,1 解方程By,p,得到y,Bp,kkkk是 y,0? k否 不存在有限最优解 确定下标r,使得 ,bbir,min|,0 y,ikyyrkik,x为退基变量，x进基变

4、量，以Bkr p代替p,得到新的基矩阵BkBr3 : 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(nm) %判断是否为标准型 fprintf(不符合要求需引入松弛变量) flag=0; else flag=1; B=A(:,n-m+1:n); %找基矩阵 cB=c(n-m+1:n); %基矩阵对应目标值的c w

5、hile 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.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)=

6、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)为非基变量，所以4 X(i),0。 end end fprintf(基向量为:); X fprintf(目标函数值为:) ; f else if(BA(:,k)0) & (b1(i)/(A(i,k)+eps) run danchunxin A=3 3 1 0 0;4 -4 0 1 0;2 -1 0 0 1 8 b=30 16 12 c=-3 -1 0 0 0 运行后的结果为: panbieshu = 3 1 0 0 0 b

7、./(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 = 9 7 3 0 0 1 目标函数值为: f = -24 xxx,，2min4) 123xxxx，,，,210 st.1234248xxx,，, 123,，,xxx244 123xj,0,1,.,4 j窗口输入: run danc

8、hunxin 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./(BA(:,2)为; ans = 5 -5 x(1)进基，x(4)退基 panbieshu = 0 2 -3 -3 0 b./(BA(:,2)为; ans = -5 -10 10 此问题不存在最优解! 五:心得体会: 通过本次实验对单纯形了解更深刻，此次实验中inf表示为一个无穷大的数。 本次做的只是最简单的线性规划问题，面对以后更大的、更复杂的问题，虽然起不了什么非常大的作用，但这是基础，所以我非常认真对待这次实验，做完

9、本次实验，使我对单纯形方法，更加熟练，对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 (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;

10、return; double mubiao(double *x) double y; y=x1-x0*x0; y=100.0*y*y; 11 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.

11、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

12、(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) pri

13、ntf(迭代次数m=%dn,m); for(j=0;jN;j+) printf(optx%d=%13.7en,j,xxkj); printf(fmin=%13.7en,f0); else m=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外其余点的中点xc f

14、or(i=0;iN;i+) xa=0; for(j=0;jN;j+) 13 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的函数值fr fr=mubiao(xr); printf(fr=%13.7e n,fr); /判断xr与xb的

15、好坏 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; else for(j=0;j=fw) for(i=0;i=fw) for(i=1;iN+1;i+) for(j=0;jN;j+) xxj=xx0j-rou*(xxj-xx0j); goto loop1; 14 else for(j=0;j=fe) for(i=0;i=fr) for(i=1;iN+1;i+) for(j=0;jN;j+) xxj=xx0j-rou*(xxj-xx0j); goto loop1; else for(j=0;jN;j+) xxNj=xoj; goto loop1; else for(j=0;jN;j+) xxNj=xrj; goto loop1; 15