数值优化——单纯形法
之前过冷水和分享了几期优化算法的方法后就没有再更新相关类推文了，最近有接触单纯形法的学习，本期就和大家分享一下用单纯形法的思想来来求函数的极值。
对问题min f(x) ,在n维空间中适当选取n+1个点x0,x1,....xn,xn+1构成一个单纯形。一般可以要求这n+1 个点使向量组x1-x0，x2-x0，.... xn-x0线性无关。
(1)计算函数值f(xi) ,i=0,1,....n f.0.1 决定坏点xh和好点xl,于是
（2）计算除去点xh外的n个点的中心，并且求出反射点
(3)若fr=f(xr)≥fh，则进行压缩，即令xs=xh+λ(xr-xh) ,并求出fs，然后转步骤(5),其中λ为给定的压缩系数，若 fr< fh ,则转步骤(4)。 (4)进行扩张,即令
计算fe,若fe≤fr,则令xs=xe否则,令xs=xr，fs=fr
(5)若fs<< span="">fh则用xs替换xh,fs替换fh，把这样得到的新点xs和其他n个点构成一个新的单纯形，重新确定xh和xl，然后返回步骤(2);若fs≥fh则转步骤(6)。
(6)若
则计算结束。取x=xl
过冷水在整理文稿的时候就觉得理解起来好抽象，还不如自己给我一段代码让我自己看，为了让读者理解起来比较容易一点，以
为案例演示单纯形法的思路
（1）任取一组解，初始值1(x=9,y=9)，解得f1=192;
（2）从(x=9,y=9)任意出发改变到其它两个点2(x=11,y=9),3(x=9,y=11);解得f2=228;f2=264;
（3）以上三个点(f1，f2，f3)从函数的最大点向对边中点做平行四边形，得到新的点4(11,7),这里过冷水插一句已知三角形三点(1 、2、3)，平行四边形的第四个点是怎么求出来的？总不能一眼就看出来是(11,7)吧1
（4）弃去步骤3中函数值最大的点2(11,9)，将(1,3,4)组建成新的新的三角形，然后重复步骤三寻找第四点的方法不断操作就可以得到函数的极小值点。
根据上述方法可以编程
clear
triangle=[9 9;9 11;11 9];
f=(triangle(:,1)-1).^2+2*(triangle(:,2)-1).^2;
a=0;
while sum(f)>=sum((triangle(:,1)-1).^2+2*(triangle(:,2)-1).^2);
f=(triangle(:,1)-1).^2+2*(triangle(:,2)-1).^2;
[m,~]=find(max(f)==f);[n,~]=find(max(f)~=f);
dot=[sum(triangle(n,:))-triangle(m,:)];
triangle(m,:)=dot;
a=a+1;
Dot(a,:)=dot;
end
ans= Dot(a-1,:)
ans =
1 1看上去该代码很简单，实际过冷水在编程化的过程中遇到了很多问题，具这个案例只是为了让读者理解单纯形法是怎么实现最优解的计算的，该方法不仅仅局限于二维，多维也是可以的，只不过不建议使用过冷水自编程序，实战水平怎么样，谁用谁知道。
再次过冷水和大家分享比较完整的单纯形法原程序求解
clc;clear
xx.x1=[8,9];
xx.x2=[10,11];
xx.x3=[8,11];
xx.alpha = 2;
xx.beta = 0.5;
xx.epsilon = 1e-2;
[ xx ] = NeldSearch( xx );
function [ param ] = NeldSearch( param )
initPts = [param.x3;param.x2;param.x1];
disp(initPts);
fvals = [func(param.x3) func(param.x2) func(param.x1)];
if(max(fvals) - min(fvals) < param.epsilon)
param.solve = mean(initPts);
disp('最优化结果');
param.solve
return;
end
%先排序；
[~,index]= sort(fvals); %默认升序；
%确定b,g,w；
param.x3 = initPts(index(1),:);%3
param.x2 = initPts(index(2),:);%2
param.x1 = initPts(index(3),:);%1
param.x4 = (param.x3 + param.x2)/2;%4
param.x5 = param.x4 + (param.x4 - param.x1);%5
if(func(param.x5) < func(param.x3))
param.x6 = param.x4 + param.alpha*(param.x4-param.x1);%6
if(param.x6 < param.x5)
param.x1 = param.x6;
else
param.x1 = param.x5;
end
elseif(func(param.x3) < func(param.x5) && func(param.x5) <= func(param.x2))
param.x1 = param.x5;
elseif(func(param.x2) < func(param.x5) && func(param.x5) <= func(param.x1))
param.x7 = param.x4 + param.beta*(param.x5-param.x4);%7
param.x1 = param.x7;
else
param.x8 = param.x4 - param.beta*(param.x4-param.x1);%8
if(func(param.x8) < func(param.x1))
param.x1 = param.x8;
else
param.x9 = (param.x1 + param.x3)/2;%9
param.x1 = param.x9;
param.x10 = (param.x2 + param.x3)/2;%10
param.x2 = param.x10;
end
end
plot(param.x3(1),param.x3(2),'*',param.x2(1),param.x2(2),'*',param.x1(1),param.x1(2),'*');
param = NeldSearch(param); %递归；
end
function f = func(point)
x1= point(1);
x2= point(2);
f = 4*(x1-5)^2+(x2-6)^2;
end腾讯云开发者
