然后比较复合形各顶点目标函数的大小，其中目标函数值最大的点作为坏点，以坏点之外其余各点的中心为映射中心，寻找坏点的映射点，一般说来此映射点的目标函数值总是小于坏点的，也就是说映射点优于坏点。
这时，以映射点替换坏点与原复合形除坏点之外其余各点构成K个顶点的新的复合形。
如此反复迭代计算，在可行域中不断以目标函数值低的新点代替目标函数值最大的坏点从而构成新复合形，使复合形不断向最优点移动和收缩，直至收缩到复合形的各顶点与其形心非常接近、满足迭代精度要求时为止。
最后输出复合形各顶点中的目标函数值最小的顶点作为近似最优点。
补充：
关于复合形法定点数目的选取
数目多少的选取，要视具体情况而定，一般说来，为了防止迭代过程中产生降维，顶点数目取多一些较好。
因为只要在k个顶点中有n+1个顶点所构成的n个矢量线性无关，搜索就不会在降维的空间里进行。
所以k值大些，降维的可能性就小些。
但是从另一方面看，顶点数目多，显然会降低计算速度。
为此，对于优化问题维数n<6时通常取k=2n；对于n>5的优化问题，一般应适当减少顶点数目，而取k=（1.25——1.5）n（取整）。
当然，顶点的最少数目不得低于n+1.
二．复合形法的优缺点
复合形法不需要计算目标函数的导数，也不进行一维搜索，因此对目标函数和约束函数都没有特殊的要求，适用范围较广。
复合形法的收敛速度较慢，特别当目标函数的维数较高和约束条件的数目增多时，这一缺点尤为突出。
另外，复合形法不能用于求解具有等式约束的优化问题。
三．问题求解
下面分别用复合形法和matlab工具箱分别进行求解并比较计算结果
问题一：
1-1函数的三维立体图
1-2.复合形法求解寻优趋势图
1-3.求解结果对照表
方法
复合形法
Matlab工具箱
x1
1.0004
1
x2
1.0002
1
f
1.8048e-07
0
问题二：
2-1函数的三维立体图
2-2.复合形法求解寻优趋势图
2-3.求解结果对照表
方法
复合形法
Matlab工具箱
x1
0.9950
-0.1431*e-05
x2
0
-0.1431*e-05
f
0.9950
8.1197e-10
问题三：
2-1函数的三维立体图
3-2.复合形法求解寻优趋势图
3-3.求解结果对照表
方法
复合形法
Matlab工具箱
x1
-2.9561
1.5345
x2
1.0558
1.5345
f
0.9903
1.9860
四．结果分析
对于求解此类问题，做出目标函数的大致图形（当然只限于三维以内）有利于我们判断函数的极值点位置以及估计函数值，同时也可以用来检验计算结果的正确性。
对于问题一，两种求解方法的计算结果基本相同，但对于问题二和问题三的求解结果为何不同，我们猜想应该是matlab工具箱的求解方法对求解函数有着特殊的要求所导致的，例如，要求函数可导或者连续等。
附录：
程序（只给出问题三的求解程序，其他问题的求解类似）
一．复合形法求解程序如下：
symsst
f=-（（sin（sqrt（s^2+t^2）））^2-0.5）/（1+0.001*（s^2+t^2））^2+0.5;
g=[s+44-st+44-t];
r=rand（1,6）;
m=-4+（4-（-4））*r;
X=reshape（m,2,3）
[x,maxf,trace_value,trace_meanvalue]=Fuhexing（f,g,X,1.3,0.7,1,0.7,[st]）;
x
maxf
[tvx,tvy]=size（trace_value）;
fxx=[1:
tvy-1];
fyy=ones（1,tvy-1）;
fyy=maxf*fyy;
plot（fxx,trace_value（1,2:
tvy）,'green',fxx,trace_meanvalue（1,2:
tvy）,'blue',fxx,fyy,'red'）
legend（'最大值','平均值','终值'）
grid
function[x,maxf,trace_value,trace_meanvalue]=Fuhexing（f,g,X,alpha,sita,gama,beta,var,eps）
%f目标函数
%g约束函数
%X初始复合形
%alpha反射系数
%sita压缩系数
%gama扩展系数
%beta收缩系数
%var自变量向量
%eps精度
%x目标函数取最小值时的自变量
%minf目标函数的最小值
trace_value=[0];
%用于记录最大值的轨迹
trace_meanvalue=[0];
%用于记录平均值的轨迹
N=size（X）;
n=N
（2）;
%n为复合形的顶点个数
FX=zeros（1,n）;
%用于存放复合形的函数值
ifnargin==8
eps=1.0e-6;
end
N=size（X）;
n=N
（2）;
Fx=zeros（1,n）;
while1
%一直执行，知道找到满足收敛条件后用break跳出循环
fori=1:
n
FX（i）=subs（f,var,X（:
i））;
%求出复合形的函数值
end
[XS,IX]=sort（FX）;
%对FX从大到小排序
Xsorted=X（:
IX）;
%Xsorted为排好序的X
px=（sum（Xsorted,2）-Xsorted（:
1））/（n-1）;
%求的除最小值以为的平均值
trace_meanvalue=[trace_meanvalue,sum（FX）/n];
Fpx=subs（f,var,px）;
%求的平均值的函数值
aaa=subs（f,var,Xsorted（:
n））;
trace_value=[trace_value,aaa];
SumF=sqrt（sum（（FX-Fpx）.^2）/（n-1））;
%计算收敛值
ifSumF<=eps
x=Xsorted（:
n）;
break;
%如果收敛，则退出循环
else
bcon_1=1;
cof_alpha=alpha;
%反射
whilebcon_1
x2=px+cof_alpha*（px-Xsorted（:
1））;
%x2为反射点坐标
gx2=subs（g,var,x2）;
%计算约束值
ifmin（gx2）>=0
%在约束条件内
bcon_1=0;
else
cof_alpha=0.5*cof_alpha;
%如不在约束条件内，则循环，直到到达约束条件
end
end
fx2=subs（f,var,x2）;
%反射点函数值反射部分结束
iffx2>XS（n）
%如果反射部分优于最优值则进行扩张操作
cof_gama=gama;
x3=x2+cof_gama*（x2-px）;
%扩张
gx3=subs（g,var,x3）;
fx3=subs（f,var,x3）;
ifmin（gx3）>=0
%如果符合边界，
iffx3>XS（n）
count=1;
%优于最优值
else
count=2;
%比最优值差
end
else
count=3;
%不符合边界
end
%扩张结束
ifcount==1
Xsorted（:
1）=x3;
%用扩张点代替最差点
X=Xsorted;
continue
else
Xsorted（:
1）=x2;
%用反色点代替最差点
X=Xsorted;
continue
end
else
%如果反射部分比最优值差
iffx2>XS
（2）
%如果反射部分大于次坏值，即反射部分优于次坏值则反射点代替
Xsorted（:
1）=x2;
X=Xsorted;
continue
else
iffx2>XS
（1）
%如果反射点比最差点好，比次差点差这用反射点代替最差点
Xsorted（:
1）=x2;
cof_beta=beta;
bcon_3=1;
whilebcon_3<4
x4=Xsorted（:
1）+cof_beta*（px-Xsorted（:
1））;
gx4=subs（g,var,x4）;
ifmin（gx4）>=0
%符合边界值
bcon_3=5;
%退出收缩
else
cof_beta=0.5*cof_beta;
bcon_3=bcon_3+1;
%最多收缩系数调整4次
end
end
ifmin（gx4）>=0
%符合边界值
fx4=subs（f,var,x4）;
FNnew=subs（f,var,Xsorted（:
1））;
%计算最差点函数值
iffx4>FNnew
%如果收缩值优于最差值
Xsorted（:
1）=x4;
%收缩值代替最差值
X=Xsorted;
continue
else
%收缩值比最差值差,则对所有点进行压缩
x0=Xsorted（:
n）;
fori=1:
n
Xsorted（:
i）=x0+sita*（Xsorted（:
i）-x0）;
end
end
else
%如果收缩后的值不符合边界
x0=Xsorted（:
n）;
fori=1:
n
Xsorted（:
i）=x0+sita*（Xsorted（:
i）-x0）;
X=Xsorted;
%用压缩值代替原复合形
continue
end
end
else
%如果反射值比最差点还要差
x0=Xsorted（:
n）;
fori=1:
n
Xsorted（:
i）=x0+sita*（Xsorted（:
i）-x0）;
X=Xsorted;
%用压缩值代替原复合形
continue
end
end
end
end
end
X=Xsorted;
end
maxf=subs（f,var,x）;
二．matlab工具箱求解命令如下：
x0=[1;1];
VLB=[-4;-4];
VUB=[4;4];
[x,fval]=fmincon（'fun',x0,[],[],[],[],VLB,VUB）
[XY]=meshgrid（-4:
0.01:
4,-4:
0.01:
4）;
Z=0.5-（（sin（（X.^2+Y.^2）.^0.5））.^2-0.5）./（1+0.001.*（X.^2+Y.^2））.^2;
mesh（X,Y,Z）;
shadinginterp;
functionf=fun（x）
f=-0.5+（（sin（x
（1）^2+x
（2）^2）.^（1/2））.^2-0.5）/（1+0.001*（x
（1）^2+x
（2）^2））^2;
升级会员 