优化理论与最优控制作业--复合形法

maxf(x)100(xi1xi)2(1xi)2i1n122.048xi2.048

2.GeneralizedRastrigin'sFunction

maxf(x)(xi10cos(2xi)10)

i1n125.12xi5.12

3.Schaffer’sFunction

maxf(x1,x2)0.5

(sinx1x2)20.5(10.001(x1x2))

222224xi4

二、复合形的算法思想复合型算法主要思想是：在满足一系列约束条件的情况下，确定初始顶点，计算每个顶点的初始函数值，进行比较，依次得出最好点，次差点和最差点。计算除最差点以外的所有点的形心点，以该点为反射中心，计算最差点的反射点，若反射点处的函数值优于最差点，则用反射点代替最差点；反之迭代进行直到满足反射点函数值优于最差点，若迭代多次直到反射系数小于某一给定精度时还未满足条件，则说明反射方向不利，此时取次坏点的反射方向为寻优继续进行的方向。

按照以上方法计算迭代，直到达到设定的最优收敛条件，此时对应的最小函数值点为最优点。

三、复合形算法实现具体迭代步骤1)给定设计变量数目n，变量界限范围ai,bi(i1,2,n)，复合形顶点数目K，精

度要求，。

2)产生初始复合形，利用Xiaii(biai)(其中j为复合形顶点标号，i为

调整变量的标号,ai,bi为调整变量的上下界，为[01]区间内服从均匀分布的伪随机数)，得初始复合形K个顶点Xj(j1,2,,K)。

3)计算复合形各顶点的目标函数值f(X(j)),j1,2,K,，在各顶点中找出最坏

点X(H)和最好点X(L)X(H）:f(X(H））maxf(X(j))(j1,2,K）X(L）:f(X(L）(j)jj）minf(X)(j1,2,K）

转入第(8)步。

4)计算除最坏点X(H)外其余各顶点的中心X(C)X

(C)1K(j)XjHK1j15)检查X(C)点的可行性。若X(C)不在可行域D内，则D域可能是一个非凸集。

这时可在以X(L)点为起点、X(C)点为端点的超立方体中(二维则为长方形)利用随机数产生新复合形的各个顶点，即

Xi(L)Xi(L）(C)(i1,2,K)

则取：

(C)biXiaiXi(i1,2,K)

否则相反，然后转回第(2)步;若X(C)在可行域，则进行下一步。6)寻求映射点X(R)X(R)X(C)(X(C)X(H））

式中映射系数通常取1.3。亦须检查点X(R)是否在可行域内。若在可行域内，则进行第(7)步；否则将映射系数减半，即0.5，再按上式计算新的映射点使其向可行域方向靠拢；若新的映射点仍在可行域外，则将再次减半继续重复迭代直到映射点X(R)进入可行域为止。而后进行下一步。7)比较映射点X(R)与最坏点X(H)的目标函数值，构成新复合形。

计算映射点X(R)的目标函数值fX(R)并与最坏点的目标函数值fX(H)相

比较。若fX(R)fX(H)，则用X(R)替换点X(H)，即X(R)X(H)，构成一个新复合形，返回第(3)步；否则将步长减半，即0.5，返回第(6)则以新的X(R)替换X(H)构成新复合形，返回第(3)步。

如果经过若干次的减半映射系数，直到值小于一个预先给定的很小正数步，重新计算新的回缩后的映射点X(R)，循环迭代，直至fX(R)fX(H)，

(通常取105)，仍不能使映射点优于最坏点，这说明该映射方向不利。为了改变映射方向，找出复合形各顶点中的次坏点X(SH)，即

X(SH）:f(X(SH））maxf(X(j))(j1,2,K,jH）

并以次坏点X(SH)替换最坏点X(H)，即X(SH)XH，返回第(4)步。8)检验是否满足迭代终止条件

反复执行上述过程，复合形随之收缩。在每一个新复合形构成之时，均应检验终止条件以判别是否结束达代。当复合形缩得很小时，各顶点的目标函数值必然近于相等。常用各顶点与好点的目标函数值之差的均方根值小于误差限作为终止迭代条件。即

1KfXf(X

(j)j1K(L)

)

212如满足迭代终止条件，输出最优解，结束迭代；否则返回步骤（3），继续进行下一次迭代。

四、复合形算法流程图五、优化求解结果问题一：

n1i1maxf(x)100(xi1xi)2(1xi)2目标函数最优值:Fmax=3905.926223最优顶点：X1=-2.048000,X2=-2.048000>>

22.048xi2.0484000300020001000X: -2.048Y: -2.048Z: 3906X: 1.952Y: -2.048Z: 3433020-2-4-3

-2-1012

图1-1函数的三维立体图

寻优轨迹

4000

最大值曲线平均值曲线3500

3000寻优值2500

2000

1500

1000 0

10203040

5060迭代次数

708090100

图1-2复合形法寻优趋势图

寻优演化轨迹图

1.510.50-0.5-1-1.5-2-2.5-2.2

-2-1.8-1.6-1.4-1.2-1-0.8-0.6-0.4

图1-3复合形法寻优演变轨迹图

结合函数三维立体图可知，此函数有两个峰点，其中一个是全局最优点，另一个是局部最优点，尝试多次运行程序，可知有时确实会陷入局部最优点。程序多次运行结果如下：

x12.0480002.048000-2.048000-2.048000-2.0480002.048000-2.048000-2.048000

x1-2.048000-2.048000-2.048000-2.048000-2.048000-2.048000-2.048000-2.048000

f(x1,x2)

37.73421437.7342113905.9262253905.9242613905.92622337.7342263905.9260423905.926223

问题二：

maxf(x)(xi10cos(2xi)10)

目标函数最优值:Fmax=80.706580最优顶点：X1=4.522967,X2=-4.523016>>

i1n125.12xi5.12

1008060402001050X: -4.52Y: 4.58Z: 80.09X: 4.58Y: 4.58Z: 79.48X: 4.48Y: -4.52X: -4.52Z: 80.34Y: -4.52Z: 80.710

5-5-10-10

0-5图2-1函数的三维立体图

寻优轨迹

85807570寻优值6560555040 0

最大值曲线平均值曲线 10203040

5060迭代次数

708090100

图2-2复合形法寻优趋势图

寻优演化轨迹

0

-1

-2

-3

-4

-5

-60.5

11.522.533.4.55

图2-3复合形法寻优演变轨迹图

此函数有很多个“山峰”，程序运行时很容易陷入局部最优点中，多次运行如下：

x14.523004-3.2223-4.522960-4.5229523.517838-3.5176-4.52308.522906

x14.5230324.607229-4.522960-2.512599-3.5177684.5230074.523002-4.522946

f(x1,x2)

80.70658068.57329780.70658066.635127.62488072.66573180.70657980.706578

问题三：

max

f(x1,x2)0.5

(sin

x1x2)20.5

2222(10.001(x1x2))24xi4

目标函数最优值:Fmax=0.990284最优顶点：X1=-3.097152,X2=-0.510336>>

10.80.60.40.20420-2-4-4

X: 0Y: 0Z: 1X: 3.1Y: 0.7Z: 0.98874

20-2图3-1函数的三维立体图

寻优轨迹1最大值曲线平均值曲线0.95 0.9寻优值0.850.80.75 0510152025迭代次数303045图3-2复合形法寻优趋势图

寻优演化轨迹

0.5

0

-0.5

-1

-1.5

-2-4

-3.5-3-2.5-2-1.5-1-0.50

图3-3复合形法寻优演变轨迹图

此函数有无限个局部极大点，只有一个全局最大值点f(0,0)=1，此函数最大值峰周围有一圈脊，它们的取值均为0.990283，多次运行结果如下：

x1-1.0433723.044291-1.8401202.1583771.0510790.630734-2.0878570.326701

x1-2.9605240.761077-2.2395-2.2785382.957379-3.074066-2.34293.121485

f(x1,x2)

0.9902840.9902840.9902840.9902840.9902840.9902840.9902840.990284

六、总结

本文用复合形算法对三个函数进行寻优，学习总结如下：

（1）对于存在局部最优点的函数，可以根据三维图形，首先确定其全局最

优点的范围，在此范围内进行寻优，但该方法对函数图形依赖过高，具有一定的局限性。

（2）反射系数要取适当，过大则造成不稳定性，使寻优过程很难收敛。（3）复合形法虽然具有不需要计算目标函数导数，也不进行一维搜索，对

目标函数和约束函数均无特殊要求，适用范围较广的优点。但复合形法的收敛速度较慢，特别当目标函数的维数较高和约束条件的数目增多时，这一缺点尤为突出。此外，当函数图形本身比较复杂时，复合形法容易陷入局部最优点，并且无法跳出，比如问题二中的函数，全局最优值为1,然后它有无数个局部最优值，均为0.990283，用复合形寻优时，尝试多次也没有找到全局最优，都陷入局部最优无法跳出。

附录：

clearallclc

%f=@(x1,x2)(-(100*(x2-x1^2)^2+(1-x1)^2));

%f=@(x1,x2)(-(x1^2-10*cos(2*pi*x1)+x2^2-10*cos(2*pi*x2)+20));;

f=@(x1,x2)(((sin(sqrt(x1^2+x2^2)))^2-0.5)/(1+0.001*(x1^2+x2^2))^2-0.5)

k=4;alpha=1.3;delta=1.0e-5;epsilon=1.0e-6;alpha0=alpha;%随机数的产生

random1=rand(1,4);random2=rand(1,4);%顶点和顶点函数值的产生%a=-2.048;b=2.048;%a=-5.12;b=5.12;a=-4.0;b=4.0;X1=zeros(1,k);X2=zeros(1,k);F=zeros(1,k);fori=1:k

%顶点

%暂存各顶点函数值

X1(i)=a+random1(i)*(b-a);X2(i)=a+random2(i)*(b-a);x1(i)=X1(i);x2(i)=X2(i);

%存储顶点，画寻优演化轨迹用

X(:,i)=[X1(i);X2(i)];end

F(i)=f(X1(i),X2(i));

%计算获得最优最差点和次最差点[FX_p]=sort(F);X=X(:,X_p);%迭代eps=0;

%顶点值排序

%顶点对应位置随排序变化,最优，最差，次差点由序号即可确定

fori=1:kend

eps=eps+(F(i)-F(1))^2;

%各顶点与好点的函数值之差的均方根%反射点函数值%存储最大值%迭代次数

eps=sqrt(eps/k);Fr=0;

iter=0;Max=0;Mean=0;

%存储平均值

while(eps>epsilon)%收敛准则，各顶点与好点的函数值之差的均方根小于给定的误差限ifiter>250end

break

alpha=alpha0;

Xc=(sum(X,2)-X(:,k))/(k-1);%除最差点外的其他各点的形心

%判断形心是否在可行域内,不在，则可以最优点Xl点和Xc为界（若Xlwhile(Xc(1)b||Xc(2)b)

ifXl(1)a1=Xl(1);b1=Xc(1);a2=Xl(2);b2=Xc(2);a1=Xl(1);b1=Xc(1);a2=Xc(2);b2=Xl(2);a1=Xc(1);b1=Xl(1);a2=Xl(2);b2=Xc(2);a1=Xc(1);b1=Xl(1);a2=Xc(2);b2=Xl(2);

elseifXl(1)Xc(2)elseifXl(1)>Xc(1)&&Xl(2)fori=1:k

X1(i)=a1+random1(i)*(b1-a1);X2(i)=a2+random2(i)*(b2-a2);F=f(X1(i),X2(i));

[FX_p]=sort(F);X=X(:,X_p);end

end

Xc=(sum(X,2)-X(:,k))/(k-1);%除最差点外的其他各点的形心

%顶点对应位置随排序变化

%顶点值排序

Xr=Xc+alpha*(Xc-X(:,k));

while(Xr(1)b||Xr(2)b)不在，则将alpha减半至反射点在可行域内

alpha=alpha/2;

end

Xr=Xc+alpha*(Xc-X(:,k));

%判断反射点是否在可行域内,

Fr=f(Xr(1),Xr(2));ifFrX(:,k)=Xr;else满足要求

%反射点函数值大于最差点的函数值，反射点替换最差点%反射点函数值小于或等于最差点函数值，重新计算反射点直到

%反射点函数值

while(Fr>F(k)||Xr(1)b||Xr(2)b)%若Fr且Xr为可行点，则结束循环；否则将alpha值减半，如此反复

ifalpha>delta

alpha=alpha/2;

Xr=Xc+alpha*(Xc-X(:,k));else

Fr=f(Xr(1),Xr(2));

%当alpha持续减半到小于预先给定的标准还没使反射点函数值大于

最差点，则改换映射方向，取对次坏点的映射

break

endend

end

ifFrX(:,k)=Xr;F(k)=Fr;

[FX_p]=sort(F);X=X(:,X_p);eps=0;fori=1:kend

%可能直接满足要求，也可能经过反射点重新计算后满足要求%将暂存顶点值的数组F中的最差点更新%顶点对应位置随排序变化%顶点值排序

%更新顶点,反射点替换最差点

%修改误差

eps=eps+(F(i)-F(1))^2;

else

eps=(eps/k)^(1/2);

%当alpha持续减半到小于预先给定的标准还没使反射点函数

值大于最差点，则改换映射方向，取对次坏点的映射次坏点替换最坏点

X(:,k)=X(:,k-1);%次差点替代最差点F(k)=f(X(1,k),X(2,k));

%每循环一次，将新点（反射点）加入顶点存储数组中

end

x1(iter+3)=Xr(1);x2(iter+3)=Xr(2);iter=iter+1;Max(iter)=-F(1);end

Mean(iter)=-mean(F);

fprintf('Fmax=%f\\n',-F(1));

fprintf('最优顶点：X1=%f,X2=%f\\n',X(:,1));i=1:1:iter;figure

plot(i,Max,i,Mean)xlabel('迭代次数');ylabel('寻优值');title('寻优轨迹');

legend('最大值曲线','平均值曲线')%添加图例%寻优演化轨迹图figureplot(x1,x2)

title('寻优演化轨迹');%函数图形figure

x1=a:0.1:b;x2=a:0.1:b;

[x1,x2]=meshgrid(x1,x2);

%z=100.*(x2-x1.^2).^2+(1-x1).^2;

%z=x1.^2-10.*cos(2.*pi.*x1)+x2.^2-10.*cos(2.*pi.*x2)+20;mesh(x1,x2,z)

z=-(sin(sqrt(x1.^2+x2.^2)).^2-0.5)./(1+0.001.*(x1.^2+x2.^2)).^2+0.5;