复合形法大作业

优化理论与最优控制大作业

( 2013-- 2014 年度第 1 学期)

题 目： 复合形法大作业 院 系： 控制与计算机工程学院 小组成员：

研控计1320班： 范冠男 1132227028

习春苗 1132227001 程丕建 1132227027 王凯 1132227013

郭萍 1132227020

研控计1322班： 赵亮 1132227143

成 绩：

日期：2013年12月8日

一、 作业题目

利用复合形法求解Schaffer’s function

2(sin2x12x2)0.5 maxf(X)0.5 （4xi4 i=1,2） 222(10.001(x1x2)) 注：本组各函数中的n值均取为2

二、 复合形法的基本原理及本文思路

1、复合形法原理：复合形法的基本思路是在n维空间的可行域中选取K个设计点（通常取n+1K2n）作为初始复合形（多面体）的顶点。然后比较复合形各顶点目标函数的大小，其中目标函数值最大的点作为坏点，以坏点之外其余各点的中心为映射中心，寻找坏点的映射点，一般说来此映射点的目标函数值总是小于坏点的，也就是说映射点优于坏点。这时，以映射点替换坏点与原复合形除坏点之外其余各点构成K个顶点的新的复合形。如此反复迭代计算，在可行域中不断以目标函数值低的新点代替目标函数值最大的坏点从而构成新复合形，使复合形不断向最优点移动和收缩，直至收缩到复合形的各顶点与其形心非常接近、满足迭代精度要求时为止。最后输出复合形各顶点中的目标函数值最小的顶点作为近似最优点。

2、本文思路：本文在理解复合形法的基础上，提出将定义域区域进行等分，分成m×m个小块。然后对每小块区域选取一个初始点进行寻优，最后比较这些初值点找到的最优值，并把最好的一个最优值作为最终的输出最优值。

三、 基本程序流程图

图一、程序流程图

四、求解寻优过程

1、函数三维图形：

10.90.80.70.60.50.40.30.20.1043210-1-2-3-4-4X: 0Y: 0Z: 14210-1-3-23Z

图二、目标函数三维图

2、理论结果：

由函数的三维图形不难看出，该函数的理论最大值为1。即当X=0，Y=0，时，Z=f(x)取最大值为1。故理论解为X=[0 ，0]T，f(X)=1。 3、寻优过程：

本文在利用复合形法求解过程中，在平面区域内将区域等分成为64小块，并在每个小块中选取一个初始值作为复合形法的初始值进行寻优计算，并将最终的最优值作为寻优结果。在区域内初始点散点图如下，从中可以看到在每个网状线格子中都有一个初值。

初始点在区域散点图43210Y-1-2-3-4-4-3-2-10X1234

图三、64个初始点在定义域散点图

同时，本文选取部分初始值寻优结果列入下表： 初始点序号 1 2 3 4 5 X0(1) -3.464932895 -2.03612987 -1.884374121 -0.948551707 0.304348946 X0(2) -3.464932895 -3.464932895 -3.464932895 -3.464932895 -3.464932895 迭代次数 27 24 21 37 46 过程最优值 0.782967657 0.990278724 0.990277778 0.990284089 0.984610841 …… 22 23 24 25 26 …… 1.580191833 2.530964452 3.901208093 -3.464932895 -2.03612987 …… -1.884374121 -1.884374121 -1.884374121 -0.948551707 -0.948551707 …… 23 23 50 21 40 …… 0.989391712 0.990026091 0.999091096 0.803540543 0.995184066 …… 60 61

…… -0.948551707 0.304348946 …… 3.901208093 3.901208093

…… 39 26 …… 0.990284089 0.984379257

62 63 64 1.580191833 2.530964452 3.901208093 3.901208093 3.901208093 3.901208093 22 49 47 0.99028409 0.918694359 0.990146661 表一、64个初始点寻优记录表

为了清楚地展示复合形的寻优过程，本文绘制了复合形法在迭代过程的寻优轨迹，也即最大值的寻找过程。下图为复合形法中找到最优值时的寻优轨迹图。

寻优轨迹10.90.80.70.6函数值 0.50.40.30.20.10 051015202530迭代次数35404550寻优轨迹(f(xh))形心点的轨迹

图四、全局最优点的寻优轨迹图

4、寻优结果：

由该方法找到的最优解为X=[0.02955566；-0.00589362]，此时最大值f(X)=

0.99909101。由此可知，该结果与理论值很接近的，证明了算法的有效性。

五、寻优分析与探讨

1、复合形法

通过上面的求解过程，我们得出并不是任意给定的初始点都能找到全局最优

点，也即函数的最大值。本文通过在定义域内选取大量的初始点来进行优化求解，并且在将区域分成64块时找到了最优点。但是实际上在这64个初始点中，能找到最优点的概率还是很低的。当然，通过仿真实验我们还发现随着在定义域内选取的初始点越多，也即分的区域块数越多，找到最优点的概率越大。 2、Matlab工具箱求解

为了验证分区选取初始点的有效性，本文还通过Matlab中自带的优化工具箱，即求解非线性规划的fmincon命令来求取该函数的最大值。实际上，对于给定函数解析式的非线性函数，该方法比复合形法要更有效。下表为结合初始点分区选择和非线性规划方法求解该函数最值的过程。 初始点序号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 X0(1) -3.447309367 -0.75535277 1.176723369 3.926936932 -3.447309367 -0.75535277 1.176723369 3.926936932 -3.447309367 -0.75535277 1.176723369 3.926936932 -3.447309367 -0.75535277 1.176723369 3.926936932 X0(2) -3.447309367 -3.447309367 -3.447309367 -3.447309367 -0.75535277 -0.75535277 -0.75535277 -0.75535277 1.176723369 1.176723369 1.176723369 1.176723369 3.926936932 3.926936932 3.926936932 3.926936932 过程最优值 0.646848776 0.990283844 0.99028409 0.646848776 0.990283844 0.999999992 0.999999989 0.990284085 0.99028409 0.999999989 0.646848776 0.99028409 0.646848776 0.990284085 0.99028409 0.646848776 表二、16个初始点寻优记录表

由此表可知，本文仅将定义域分成16块即找到3次全局最优点（即表中绿色部分初始点）。同时，本文也做过仿真实验，当分的区域块数越大，找到全局最优的机率越大。而对于寻优而言，我们只需找到一次全局最优点即得到该函数的最大值，进一步验证了分区的有效性。

六、总结

本文在理解复合形法的基础上，针对复合形法寻优过程对初值的依赖性很大这一问题，提出将定义域区域进行等分，然后对每小块区域选取一个初始点进行寻优，然后比较这些初值点找到的最优值，把最好的一个最为最终的最优值。实验证明该方法很有效。同时，我们也认识到复合形法也存在一定的问题，运算比较慢。本文通过Matlab求解非线性规划的方法进一步对定义域分区的思想进行验证。从仿真结果中结果中可以看出，这种方法比复合形法更有效。

七、附录

1、目标函数的三维图形绘制程序

x=-4:0.1:4; y=-4:0.1:4; [X Y]=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、复合形法求解程序如下：

syms x1 x2

f=-(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); for i=1:m

z(i)=a1+(b1-a1)/m*(i-1)+rand(1,1)*(b1-a1)/m; end for i=1:m for j=1:m

X0(1,i+m*(j-1))=z(i); end end for i=1:m

W=[W ones(1,m)*z(i)]; end X0(2,:)=W; for i=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

[maxf index]=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; while 1 fx=zeros([1,k]); X=zeros([n,k]);

g=[var-aa;-var+bb]; %约束函数g(X) %产生初值 X(:,1)=x0; for i=2:k

r=abs(rand([2,1])); X(:,i)=aa+r.*(bb-aa); end %%寻优

traceFXk=[0]; %用来记录每次迭代所产生的最坏点Xh tracefxc0=[0]; %用来记录每次迭代所产生的形心点Xc0 FXk=[]; while 1 for i=1:k

fx(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; for i=1:k

Sum=Sum+(FX(i)-fxc0)^2; end

E=sqrt(Sum/k); %终止迭代条件 if E<=e

x=Xsorted(:,1); % 令x=xL break; else

xc=sum(Xsorted(:,1:(k-1)),2)/(k-1); %除最坏点外其余K-1个顶点的形心点 gxc=subs(g,var,xc); if min(gxc)>=0 %若形心点满足约束 xr=xc+alpha*(xc-Xsorted(:,k)); fxr=subs(f,var,xr); gxr=subs(g,var,xr); if min(gxr)>=0

if fxrelse

if alpha<=e1 %若f(xr)>f(xh),但此时反射系数alpha已经小于e1

x0=Xsorted(:,1);%则需将复合形中的所有顶点向最好点收缩以产生新的复合形 for i=1:k

Xsorted(:,i)=x0+sita*(Xsorted(:,i)-x0); end

else%若f(xr)>f(xh),但此时反射系数alpha仍大于e1

alpha=alpha/2;%则将反射系数减小为原来的一半，重新计算反射点 end end else

alpha=alpha/2; end else

for i=1:n %若形心点不可行，则分别用形心点xc(i)和最好点xL(i)代替变量xi的上下限 if xc(i)aa(i)=Xsorted(i,1); bb(i)=xc(i); end end break; end end

X=Xsorted; %得到新的复合形 end

if aa==a&bb==b break; else

continue; end end format long minf=subs(f,var,x)

minf1=minf*ones(1,length(traceFXk)); d=length(traceFXk); xx=1:length(traceFXk); figure

plot(xx,-traceFXk,'b',xx,-tracefxc0,'r--') legend('寻优轨迹(f(xh))','形心点的轨迹',0); xlabel('迭代次数'); ylabel('函数值');

title('寻优轨迹'); end

3、Matlab工具箱求解命令如下： %各参数初始化

M=[];%记录每个初值迭代后最优值的向量 W=[];

t=[];%记录各个初始值取最优值时的解向量 m=4;%定义域分成4*4个块数 a=-4;%区域下限 b=4;%区域上限 z=zeros(1,m);

X0=zeros(2,m*m);%记录所有初值矩阵

%初值的产生过程 for i=1:m

z(i)=a+(b-a)/m*(i-1)+rand(1,1)*(b-a)/m; end for i=1:m for j=1:m

X0(1,i+m*(j-1))=z(i); end end for i=1:m

W=[W ones(1,m)*z(i)]; end X0(2,:)=W;

%全局最优值求解过程 for i=1:m*m VLB=[-4;-4]; VUB=[4;4];

[x,fval,exitflag]=fmincon(@fun,X0(:,i),[],[],[],[],VLB,VUB); %优化工具箱函数调用 if exitflag>0 %判断迭代是否收敛 M=[M,fval]; t=[t,x]; end end

[maxf index]=min(M) x=t(:,index)

function f=fun(x) %目标函数

f=-(0.5-(sin((sqrt(x(1).^2+x(2).^2))).^2-0.5)./(1+0.001*(x(1).^2+x(2).^2)).^2); end