复合形法大作业

优化理论与最优控制大作业(2013--2014年度第1学期)题目：复合形法大作业院系：控制与计算机工程学院小组成员：研控计1320班：范冠男1132227028习春苗1132227001程丕建1132227027王凯1132227013郭萍1132227020研控计1322班：赵亮1132227143成绩：日期：2013年12月8日一、作业题目利用复合形法求解Schaffer’sfunction2222122212))(001.01(5.0)(sin5.0)(maxxxxxXf（44ixi=1,2）注：本组各函数中的n值均取为2二、复合形法的基本原理及本文思路1、复合形法原理：复合形法的基本思路是在n维空间的可行域中选取K个设计点（通常取n+1K2n）作为初始复合形（多面体）的顶点。然后比较复合形各顶点目标函数的大小，其中目标函数值最大的点作为坏点，以坏点之外其余各点的中心为映射中心，寻找坏点的映射点，一般说来此映射点的目标函数值总是小于坏点的，也就是说映射点优于坏点。这时，以映射点替换坏点与原复合形除坏点之外其余各点构成K个顶点的新的复合形。如此反复迭代计算，在可行域中不断以目标函数值低的新点代替目标函数值最大的坏点从而构成新复合形，使复合形不断向最优点移动和收缩，直至收缩到复合形的各顶点与其形心非常接近、满足迭代精度要求时为止。最后输出复合形各顶点中的目标函数值最小的顶点作为近似最优点。2、本文思路：本文在理解复合形法的基础上，提出将定义域区域进行等分，分成m×m个小块。然后对每小块区域选取一个初始点进行寻优，最后比较这些初值点找到的最优值，并把最好的一个最优值作为最终的输出最优值。三、基本程序流程图图一、程序流程图四、求解寻优过程1、函数三维图形：-4-3-2-101234-4-3-2-10123400.10.20.30.40.50.60.70.80.91ZX:0Y:0Z:1图二、目标函数三维图2、理论结果：由函数的三维图形不难看出，该函数的理论最大值为1。即当X=0，Y=0，时，Z=f(x)取最大值为1。故理论解为X=[0，0]T，f(X)=1。3、寻优过程：本文在利用复合形法求解过程中，在平面区域内将区域等分成为64小块，并在每个小块中选取一个初始值作为复合形法的初始值进行寻优计算，并将最终的最优值作为寻优结果。在区域内初始点散点图如下，从中可以看到在每个网状线格子中都有一个初值。-4-3-2-101234-4-3-2-101234XY初始点在区域散点图图三、64个初始点在定义域散点图同时，本文选取部分初始值寻优结果列入下表：初始点序号X0(1)X0(2)迭代次数过程最优值1-3.464932895-3.464932895270.7829676572-2.03612987-3.464932895240.9902787243-1.884374121-3.464932895210.9902777784-0.948551707-3.464932895370.99028408950.304348946-3.464932895460.984610841…………………………221.580191833-1.884374121230.989391712232.530964452-1.884374121230.990026091243.901208093-1.884374121500.99909109625-3.464932895-0.948551707210.80354054326-2.03612987-0.948551707400.995184066…………………………60-0.9485517073.901208093390.990284089610.3043489463.901208093260.984379257621.5801918333.901208093220.99028409632.5309644523.901208093490.918694359643.9012080933.901208093470.990146661表一、64个初始点寻优记录表为了清楚地展示复合形的寻优过程，本文绘制了复合形法在迭代过程的寻优轨迹，也即最大值的寻找过程。下图为复合形法中找到最优值时的寻优轨迹图。0510152025303540455000.10.20.30.40.50.60.70.80.91迭代次数函数值寻优轨迹寻优轨迹(f(xh))形心点的轨迹图四、全局最优点的寻优轨迹图4、寻优结果：由该方法找到的最优解为X=[0.02955566；-0.00589362]，此时最大值f(X)=0.99909101。由此可知，该结果与理论值很接近的，证明了算法的有效性。五、寻优分析与探讨1、复合形法通过上面的求解过程，我们得出并不是任意给定的初始点都能找到全局最优点，也即函数的最大值。本文通过在定义域内选取大量的初始点来进行优化求解，并且在将区域分成64块时找到了最优点。但是实际上在这64个初始点中，能找到最优点的概率还是很低的。当然，通过仿真实验我们还发现随着在定义域内选取的初始点越多，也即分的区域块数越多，找到最优点的概率越大。2、Matlab工具箱求解为了验证分区选取初始点的有效性，本文还通过Matlab中自带的优化工具箱，即求解非线性规划的fmincon命令来求取该函数的最大值。实际上，对于给定函数解析式的非线性函数，该方法比复合形法要更有效。下表为结合初始点分区选择和非线性规划方法求解该函数最值的过程。初始点序号X0(1)X0(2)过程最优值1-3.447309367-3.4473093670.6468487762-0.75535277-3.4473093670.99028384431.176723369-3.4473093670.9902840943.926936932-3.4473093670.6468487765-3.447309367-0.755352770.9902838446-0.75535277-0.755352770.99999999271.176723369-0.755352770.99999998983.926936932-0.755352770.9902840859-3.4473093671.1767233690.9902840910-0.755352771.1767233690.999999989111.1767233691.1767233690.646848776123.9269369321.1767233690.9902840913-3.4473093673.9269369320.64684877614-0.755352773.9269369320.990284085151.1767233693.9269369320.99028409163.9269369323.9269369320.646848776表二、16个初始点寻优记录表由此表可知，本文仅将定义域分成16块即找到3次全局最优点（即表中绿色部分初始点）。同时，本文也做过仿真实验，当分的区域块数越大，找到全局最优的机率越大。而对于寻优而言，我们只需找到一次全局最优点即得到该函数的最大值，进一步验证了分区的有效性。六、总结本文在理解复合形法的基础上，针对复合形法寻优过程对初值的依赖性很大这一问题，提出将定义域区域进行等分，然后对每小块区域选取一个初始点进行寻优，然后比较这些初值点找到的最优值，把最好的一个最为最终的最优值。实验证明该方法很有效。同时，我们也认识到复合形法也存在一定的问题，运算比较慢。本文通过Matlab求解非线性规划的方法进一步对定义域分区的思想进行验证。从仿真结果中结果中可以看出，这种方法比复合形法更有效。七、附录1、目标函数的三维图形绘制程序x=-4:0.1:4;y=-4:0.1:4;[XY]=meshgrid(x,y);Z=0.5-(sin((sqrt(X.^2+Y.^2))).^2-0.5)./(1+0.001*(X.^2+Y.^2)).^2;mesh(X,Y,Z)xlabel('X');xlabel('Y');xlabel('Z');2、复合形法求解程序如下：symsx1x2f=-(0.5-(sin((sqrt(x1.^2+x2.^2))).^2-0.5)./(1+0.001*(x1.^2+x2.^2)).^2);%目标函数a=[-4;-4];b=[4;4];alpha=1;var=[x1;x2];e=1.0e-8;e1=1.0e-6;sita=0.5;M=[];%记录每个初值迭代后最优值的向量W=[];%记录各个初始值取最优值时的解向量t=[];%记录各个初始值取最优值时的解向量D=[];%记录各个初始值取最优值时的迭代次数m=8;%定义域分成8*8个块数a1=-4;b1=4;z=zeros(1,m);X0=zeros(2,m*m);fori=1:mz(i)=a1+(b1-a1)/m*(i-1)+rand(1,1)*(b1-a1)/m;endfori=1:mforj=1:mX0(1,i+m*(j-1))=z(i);endendfori=1:mW=[Wones(1,m)*z(i)];endX0(2,:)=W;fori=1:m*m[x,d,minf]=childfun(f,X0(:,i),a,b,alpha,sita,var,e,e1);M=[M,minf];t=[t,x];D=[D,d];end[maxfindex]=min(M)x=t(:,index)d1=D(index)function[x,d,minf]=childfun(f,x0,a,b,alpha,sita,var,e,e1)%f为目标函数%%g为约束条件%a为xi的下限a=[a1;a2;…;an]%b为xi的上限b=[b1;b2;…;bn]%alpha为反射系数o%var为自变量向量var=[x1;x2;…;xn]%e为运算中止精度%e1为反射系数收缩下限%sita为紧缩系数aa=a;bb=b;n=2;k=3;while1fx=zeros([1,k]);X=zeros([n,k]);g=[var-aa;-var+bb];%约束函数g(X)%产生初值X(:,1)=x0;fori=2:kr=abs(rand([2,1]));X(:,i)=aa+r.*(bb-aa);end%%寻优traceFXk=[0];%用来记录每次迭代所产生的最坏点Xhtracefxc0=[0];%用来记录每次迭代所产生的形心点Xc0FXk=[];while1fori=1:kfx(i)=subs(f,var,X(:,i));%计算复合形所有顶点的函数值end[FX,IX]=sort(fx);%对复合形所有顶点的函数值从小到大排序Xsorted=X(:,IX);%得到排序后的函数值所对应的x值traceFXk=[traceFXk,FX(k)];xc0=sum(Xsorted,2)/k;%复合形所有顶点的形心点fxc0=subs(f,var,xc0);tracefxc0=[tracefxc0,fxc0];Sum=0;fori=1:kSum=Sum+(FX(i)-fxc0)^2;endE=sqrt(Sum/k);%终止迭代条件ifE=ex=Xsorted(:,1);%令x=xLbreak;elsexc=sum(Xsorted(:,1:(k-1)),2)/(k-1);%除最坏点外其余K-1个顶点的形心点gxc=subs(g,var,xc);ifmin(gxc)=0%若形心点满足约束xr=xc+alpha*(xc-Xsorted(:,k));fxr=subs(f,var,xr);gxr=subs(g,var,xr);ifmin(gxr)=0iffxrFX(k)%若f(xr)f(xh),则令xh=xr，产生新的复合形Xsorted(:,k)=xr;elseifalpha=e1%若f(xr)f(xh),但此时反射系数alpha已经小于e1x0=Xsorted(