约束优化设计实验报告力学系型号:联想y470CPU:i5-2450M内存:2GB系统:win7-64位求解问题:如上是以下三个约束方法共同需要求解的问题,预估结果:在(x1,x2,x3)≈(23,13,12)点附近存在极值。
其中,每个方法对应的初始条件分别为:(1)随机试验法设计变量范围:随机试验点数:N=1000精度:eps=0.001(2)随机方向法初始点:x0=(25,15,5)初始步长:a0=0.5精度:eps=0.001(3)线性规划单纯形法初始复合形:X=[20 23 25 30;10 13 15 20;10 9 5 0]顶点个数:n=4精度:eps=0.01计算结果:程序说明:主程序为main,运行main后按提示即可得到相应约束方法的求解结果。
程序如下:1、主程序clear;global kk;kk=0;disp('1.随机试验法');disp('2.随机方向法');disp('3.线性规划单纯形法');while 1n0=input('请输入上面所想选择约束优化方法的编号(1、2、3):');if n0==1||n0==2||n0==3break;enddisp('此次输入无效.');enddisp(' ');disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');[xx,yy]=fmins(n0);fprintf('迭代次数为:%8.0f\n', kk);disp('所求极值点的坐标向量为:');fprintf(' %16.5f\n', xx);fprintf('所求函数的极值为:%16.5f\n', yy);2、调用函数function [xx,yy]=fmins(n0)if n0==1tic;[xx,yy]=suijishiyan();toc;elseif n0==2tic;[xx,yy]=suijifangxiang();toc;elseif n0==3tic;[xx,yy]=danchunxing();toc;endend3、单纯形法function [xx,yy]=danchunxing()clear;global kk;syms a b c;f=-a*b*c;g=[-a+2*b+2*c>=0;a+2*b+2*c<=72;abs(a-b-10)<=1e-3;b>=10;b<=20];X=[20 23 25 30;10 13 15 20;10 9 5 0];alpha=1.3;sita=0.5;gama=1;beta=0.7;var=[a;b;c];eps=0.001;N=size(X);n=N(2);FX=zeros(1,n);while 1for i=1:nFX(i)=double(subs(f,var,X(:,i)));end[XS,IX]=sort(FX);Xsorted=X(:,IX);px=sum(Xsorted(:,1:(n-1)),2)/(n-1);Fpx=double(subs(f,var,px));SumF=0;for i=1:nSumF=SumF+(FX(IX(i))-Fpx)^2;endSumF=sqrt(SumF/(n-1));if SumF<=epsxx=Xsorted(:,1);break;elsebcon_1=1;cof_alpha=alpha;while bcon_1x2=px+cof_alpha*(px-Xsorted(:,n));gx2=double(subs(g,var,x2));if min(gx2)>0bcon_1=0;elsecof_alpha=0.7*(cof_alpha);endendfx2=double(subs(f,var,x2));if fx2<XS(1)cof_gama=gama;bcon_2=1;while bcon_2x3=x2+cof_gama*(x2-px);gx3=double(subs(g,var,x3));fx3=double(subs(f,var,x3));if min(gx3)>0bcon_2=0;if fx3<XS(1)count=1;elsecount=2;endelsebcon_2=0;count=3;endendif count==1Xsorted(:,n)=x3;X=Xsorted;continueelseXsorted(:,n)=x2;X=Xsorted;continueendelseif fx2<XS(n-1)Xsorted(:,n)=x2;X=Xsorted;continueelseif fx2<XS(n)Xsorted(:,n)=x2;cof_beta=beta;bcon_3=1;while bcon_3<4x4=Xsorted(:,n)+cof_beta*(px-Xsorted(:,n));gx4=double(subs(g,var,x4));if min(gx4)>0bcon_3=5;elsecof_beta=cof_beta/2;bcon_3=bcon_3+1;endendif min(gx4)>0fx4=double(subs(f,var,x4));FNnew=double(subs(f,var,Xsorted(:,n)));if fx4<FNnewXsorted(:,n)=x4;X=Xsorted;continueelsex0=Xsorted(:,1);for i=1:nXsorted(:,i)=x0+sita*(Xsorted(:,i)-x0);endendelsex0=Xsorted(:,1);for i=1:nXsorted(:,i)=x0+sita*(Xsorted(:,i)-x0);X=Xsorted;continueendendelsex0=Xsorted(:,1);for i=1:nXsorted(:,i)=x0+sita*(Xsorted(:,i)-x0);X=Xsorted;continueendendendendendX=Xsorted;kk=kk+1;endyy=-double(subs(f,var,xx));end4、随机方向法function [xx,yy]=suijifangxiang()clear;global kk;x0 = [25;15;5];n = 3;eps = 0.001;a0 = 0.5;k = 10000;f0 = -x0(1)*x0(2)*x0(3);fl = -x0(1)*x0(2)*x0(3);i = 1;while 1a = a0;r = -1 + 2*rand(n,k);for j = 1:kss=(sum(abs(r(:,j))))^0.5;e = r(:,j)./ss;x = x0 + a*e;if (-x(1)+2*x(2)+2*x(3)>=0)&&(x(1)+2*x(2)+2*x(3)<=72)&&(abs(x(1)-x(2)-10)<=1e-3)&&(x(2)>=10)&&(x(2)<=20)f =-x(1)*x(2)*x(3);if f < flfl = f;d = e;xl = x;endendendwhile 1x = xl;a = 1.3*a;x = x+a*d;if (-x(1)+2*x(2)+2*x(3)>=0)&&(x(1)+2*x(2)+2*x(3)<=72)&&(abs(x(1)-x(2)-10)<=1e-3)&&(x(2)>=10)&&(x(2)<=20)f =-x(1)*x(2)*x(3);if f < flfl = f;elsebreakendelsebreakkk=kk+1;endfor b = 1:50x = x - a*d;a = 0.7*a;x = x + a*d;if (-x(1)+2*x(2)+2*x(3)>=0)&&(x(1)+2*x(2)+2*x(3)<=72)&&(abs(x(1)-x(2)-10)<=1e-3)&&(x(2)>=10)&&(x(2)<=20)f =-x(1)*x(2)*x(3);if f < flfl = f;breakendendendepsl = abs((f0-fl)/f0);if epsl < eps&&x(1)<23breakelsei = i+1;x0 = x;fl = f;endkk=kk+1;endxx = x;yy =-f;end5、随机试验法function [xx,yy]=suijishiyan()format longN=1000;x1=[];x2=[];x3=[];zmax=-inf;x10=unifrnd(20,30,N,1);x20=unifrnd(10,20,N,1);x30=unifrnd(-10,16,N,1);for i=1:Nfor j=1:Nfor k=1:Nif (-x10(i)+2*x20(j)+2*x30(k)>=0)&& (x10(i)+2*x20(j)+2*x30(k)<=72)&&(abs(x10(i)-x20(j)-10)<=1e-3)z=x10(i)*x20(j)*x30(k);if zmax<zzmax=z;x1=x10(i);x2=x20(j);x3=x30(k);endendendendendxx=[x1;x2;x3];yy=zmax;end。