阻尼牛顿法 matlab程序,【工程优化】最优化算法--牛顿法、阻尼牛顿法及单纯形法...

牛顿法

使用条件：目标函数具有二阶导数，且海塞矩阵正定。

优缺点： 收敛速度快、计算量大、很依赖初始点的选择。

算法的基本步骤：

算法流程图：

阻尼牛顿法

与牛顿法基本相同，只是加入了一维精确搜索：

优缺点：改善了局部收敛性。

我们假设要求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

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=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

运行结果如下图：