牛顿法
使用条件：目标函数具有二阶导数，且海塞矩阵正定。
优缺点： 收敛速度快、计算量大、很依赖初始点的选择。
算法的基本步骤：
算法流程图：
阻尼牛顿法
与牛顿法基本相同，
只是加入了一维精确搜索：
优缺点：改善了局部收敛性。
我们假设要求f=(x-1)*(x-1)+y*y的最小值，具体算法实现如下，只需要运行NTTest.m文件，其它函数文件放在同一目录下即可：
1、脚本文件NTTest.m
clear all
clc
syms x y
f=(x-1)*(x-1)+y*y;
var=[x y];
x0=[1 1];eps=0.000001;
disp('牛顿法:')
minNT(f,x0,var,eps)
disp('阻尼牛顿法:')
minMNT(f,x0,var,eps)
2、minNT.m
function [x,minf]=minNT(f,x0,var,eps)
%目标函数：f
%初始点：x0
%自变量向量：var
%精度：eps
%目标函数取最小值时的自变量值：x;
%目标函数的最小值：minf
format long;
if nargin==3
eps=1.0e-6;
end
tol=1;
syms L
% x0=transpose(x0);
while tol>eps %不满足精度要求
gradf=jacobian(f,var); %梯度方向
jacf=jacobian(gradf,var); %雅克比矩阵
v=Funval(gradf,var,x0);%梯度的数值解
tol=norm(v);%计算梯度（即一阶导）的大小
pv=Funval(jacf,var,x0);%二阶导的数值解
p=-inv(pv)*transpose(v); %搜索方向
x1=x0+p';%进行迭代
x0=x1;
end
x=x1;
minf=Funval(f,var,x);
format short;
3、minMNT.m
function [x,minf]=minMNT(f,x0,var,eps)
%目标函数：f
%初始点：x0
%自变量向量：var
%精度：eps
%目标函数取最小值时的自变量值：x;
%目标函数的最小值：minf
format long;
if nargin==3
eps=1.0e-6;
end
tol=1;
syms L
% x0=transpose(x0);
while tol>eps %不满足精度要求
gradf=jacobian(f,var); %梯度方向
jacf=jacobian(gradf,var); %雅克比矩阵
v=Funval(gradf,var,x0);%梯度的数值解
tol=norm(v);%计算梯度（即一阶导）的大小
pv=Funval(jacf,var,x0);%二阶导的数值解
p=-inv(pv)*transpose(v); %搜索方向
%%%%寻找最佳步长%%%
y=x0+L*p';
yf=Funval(f,var,y);
[a,b]=minJT(yf,0,0.1);
xm=minHJ(yf,a,b); %黄金分割法进行一维搜索最佳步长
x1=x0+xm*p';%进行迭代
x0=x1;
end
x=double(x1);
minf=double(Funval(f,var,x));
format short;
4、minHJ.m
function [x,minf]=minHJ(f,a,b,eps)
%目标函数：f
%极值区间左端点：a
%极值区间右端点：b
%精度：eps
%目标函数取最小值时自变量的值：x
%目标函数所取的最小值：minf
format long;
if nargin==3
eps=1.0e-6;
end
l=a+0.382*(b-a); %试探点
u=a+0.618*(b-a); %试探点
k=1;
tol=b-a;
while tol>eps&&k<100000
fl=subs(f,findsym(f),l); %试探点函数值
fu=subs(f,findsym(f),u); %试探点函数值
if fl>fu
a=1; %改变区间左端点
l=u;
u=a+0.618*(b-a); %缩短搜索区间
else
b=u; %改变区间右端点
u=l;
l=a+0.382*(b-a); %缩短搜索区间
end
k=k+1;
tol=abs(b-a);
end
if k==100000
disp('找不到最小值！');
x=NaN;
minf=NaN;
return;
end
x=(a+b)/2;
minf=subs(f,findsym(f),x);
format short;
5、minJT.m
function [minx,maxx]=minJT(f,x0,h0,eps)
%目标函数：f
%初始点：x0
%初始步长：h0
%精度：eps
%目标函数取包含极值的区间左端点：minx
%目标函数取包含极值的区间右端点：maxx
format long
if nargin==3
eps=1.0e-6;
end
x1=x0;
k=0;
h=h0;
while 1
x4=x1+h; %试探步
k=k+1;
f4=subs(f,findsym(f),x4);
f1=subs(f,findsym(f),x1);
if f4<f1
x2=x1;
x1=x4;
f2=f1;
f1=f4;
h=2*h; %加大步长
else
if k==1
h=-h; %方向搜索
x2=x4;
f2=f4;
else
x3=x2;
x2=x1;
x1=x4;
break;
end
end
end
minx=min(x1,x3);
maxx=x1+x3-minx;
format short;
6、Funval.m
function fv=Funval(f,varvec,varval)
var=findsym(f);
varc=findsym(varvec);
s1=length(var);
s2=length(varc);
m=floor((s1-1)/3+1);
varv=zeros(1,m);
if s1~=s2
for i=0:((s1-1)/3)
k=findstr(varc,var(3*i+1));
index=(k-1)/3;
varv(i+1)=varval(index+1);
end
fv=subs(f,var,varv);
else
fv=subs(f,varvec,varval);
end
运行结果如下图：
单纯形法
单纯形法的理论还有点复杂，而本文主要针对算法的基本实现，因此，理论部分就此略过，详情可以参考网上的相关资料。
下面给出具体的实现：
我们以具体实例来说明：
具体的Matlab实现如下：
1、脚本文件：
clear all clc % A=[2 2 1 0 0 0 % 1 2 0 1 0 0 % 4 0 0 0 1 0 % 0 4 0 0 0 1]; % c=[-2 -3 0 0 0 0]; % b=[12 8 16 12]'; % baseVector=[3 4 5 6]; A=[1 1 -2 1 0 0 2 -1 4 0 1 0 -1 2 -4 0 0 1]; c=[1 -2 1 0 0 0]; b=[12 8 4]'; baseVector=[4 5 6]; [x y]=ModifSimpleMthd(A,c,b,baseVector)
2、ModifSimpleMthd.m文件
function [x,minf]=ModifSimpleMthd(A,c,b,baseVector)
%约束矩阵：A
%目标函数系数向量：c
%约束右端向量：b
%初始基向量：baseVector
%目标函数取最小值时的自变量值：x
%目标函数的最小值：minf
sz=size(A);
nVia=sz(2);
n=sz(1);
xx=1:nVia;
nobase=zeros(1,1);
m=1;
if c>=0
vr=find(c~=0,1,'last');
rgv=inv(A(:,(nVia-n+1):nVia))*b;
if rgv>=0
x=zeros(1,vr);
minf=0;
else
disp('不存在最优解');
x=NaN;
minf=NaN;
return;
end
end
for i=1:nVia %获取非基变量下标
if(isempty(find(baseVector==xx(i),1)))
nobase(m)=i;
m=m+1;
else
;
end
end
bCon=1;
M=0;
B=A(:,baseVector);
invB=inv(B);
while bCon
nB=A(:,nobase); %非基变量矩阵
ncb=c(nobase); %非基变量系数
B=A(:,baseVector); %基变量矩阵
cb=c(baseVector); %基变量系数
xb=invB*b;
f=cb*xb;
w=cb*invB;
for i=1:length(nobase) %判别
sigma(i)=w*nB(:,i)-ncb(i);
end
[maxs,ind]=max(sigma); %ind为进基变量下标
if maxs<=0 %最大值小于零，输出解最优
minf=cb*xb;
vr=find(c~=0,1,'last');
for l=1:vr
ele=find(baseVector==l,1);
if(isempty(ele))
x(l)=0;
else
x(l)=xb(ele);
end
end
bCon=0;
else
y=inv(B)*A(:,nobase(ind));
if y<=0 %不存在最优解
disp('不存在最优解！');
x=NaN;
minf=NaN;
return;
else
minb=inf;
chagB=0;
for j=1:length(y)
if y(j)>0
bz=xb(j)/y(j);
if bz<minb
minb=bz;
chagB=j;
end
end
end %chagB为基变量下标
tmp=baseVector(chagB); %更新基矩阵和非基矩阵
baseVector(chagB)=nobase(ind);
nobase(ind)=tmp;
for j=1:chagB-1 %基变量矩阵的逆矩阵变换
if y(j)~=0
invB(j,:)=invB(j,:)-invB(chagB,:)*y(j)/y(chagB);
end
end
for j=chagB+1:length(y)
if y(j)~=0
invB(j,:)=invB(j,:)-invB(chagB,:)*y(j)/y(chagB);
end
end
invB(chagB,:)=invB(chagB,:)/y(chagB);
end
end
M=M+1;
if(M==1000000) %迭代步数限制
disp('找不到最优解！');
x=NaN;
minf=NaN;
return;
end
end
运行结果如下图：
关于最优化的更多算法实现，请访问
http://download.csdn.net/detail/tengweitw/8434549，里面有每个算法的索引说明，当然也包括上述算法。
作者：nineheadedbird
