实用最优化方法Matlab程序设计

实用最优化方法Matlab程序设计

程序如下：

function lambda = nonexact(x0,s0)

g0 = grad(x0);

f0=f(x0);

a=0;

c1=0.1;

c2=0.5;

lambdak=1;

sk=s0;

d=-c1*lambdak*g0'*sk;

xk=x0 + lambdak*sk;

f1=f(xk);

e=f0-f1;

while d>e

lambdak=(lambdak+a)/2;

xk=x0 + lambdak*sk;

fk=f(xk);

e=f0-fk;

d=-c1*lambdak*g0'*sk;

end

lambdak

function g=grad(x)

g=zeros(2,1);

g(1)=-4*x(1)*(x(2)-x(1)^2)-2*(1-x(1));

g(2)=2*(x(2)-x(1)^2);

function fa=f(x)

fa=(x(2)-x(1)^2)^2+(1-x(1))^2;

在命令窗口中输入x0=[0;1]; s0=[-1;1];nonexact(x0,s0) 输出结果为：

程序如下：

function x_star = cong(x0,eps)

g0 = grad(x0);

res0 = norm(g0);

resk=res0;

k = 0;

xk =x0;

while resk>eps

k=k+1;

if k==1

sk=-grad(xk);

else

sk=-grad(xk)+resk^2/res0^2*s0;

end

lambdak=step(xk,sk);

x0=xk;

xk=x0 + lambdak*sk;

res0=resk;

resk=norm(grad(xk));

s0=sk;

fprintf('------the%d-th iteration, the residual is %f, the lambdak is %f\n\n', k ,resk,lambdak);

end

x_star=xk;

disp('the optimal solution is');

x_star

%sub functions

function g=grad(x)

g=zeros(4,1);

g(1)=2*x(1)-2*x(2)+2;

g(2)=4*x(2)-2*x(1)-x(3)+3;

g(3)=2*x(3)-x(2)-1;

g(4)=2*x(4);

function lambda=step(x,s)

a=2*x(1)*s(1)-2*x(2)*s(1)-2*x(1)*s(2)+4*x(2)*s(2)+2*s(3)*x(3)+2*x(4)*s(4)-s(2)*x (3)-s(3)*x(2)+2*s(1)+3*s(2)-s(3);

b=2*s(1)^2+4*s(2)^2-4*s(1)*s(2)+2*s(3)^2+2*s(4)^2-2*s(2)*s(3);

lambda=-a/b;

在命令窗口中输入x0=[0;0;0;0];eps=1e-6;cong(x0,eps)

输出结果为

-

当

H= 2I时

k

程序如下：

function x_star = dfph(x0,eps)

g0 = grad(x0);

res0 = norm(g0);

res=res0;

k = 0;

xk =x0;

while res>eps

k=k+1;

H=eye(length(x0));

sk=-H*grad(xk);

lambdak = nonexact(xk,sk);

x0=xk;

xk=x0 + lambdak*sk;

res=norm(grad(xk));

fprintf('------the%d-th iteration, the residual is %f, the lambdak is %f\n\n', k ,res,lambdak);

end

x_star=xk;

disp('the optimal solution is');

x_star

function lambda = nonexact(x0,s0)

g0 = grad(x0);

f0=f(x0);

a=0;

c1=0.1;

c2=0.5;

lambdak=1;

d=-c1*lambdak*g0'*s0;

xk=x0 + lambdak*s0;

f1=f(xk);

e=f0-f1;

while d>e

lambdak=(lambdak+a)/2;

xk=x0 + lambdak*s0;

fk=f(xk);

e=f0-fk;

d=-c1*lambdak*g0'*s0;

end

lambda=lambdak;

function g=grad(x)

g=zeros(2,1);

g(1)=2*x(1)*exp(x(1)^2+x(2)^2)+1;

g(2)=4*x(2)+2*x(2)*exp(x(1)^2+x(2)^2);

function fa=f(x)

fa=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);

在命令窗口中输入x0=[0;1];eps=1e-3;dfph(x0,eps) 输出结果为：

当1

2()k

k H f x -??=???

时

程序如下：

function x_star = newton1(x0,eps)

g0 = grad(x0);

res0 = norm(g0);

res=res0;

k = 0;

xk =x0;

while res>eps

k=k+1;

H=inv(GRAND(xk));

sk=-H*grad(xk);

lambdak = nonexact(xk,sk);

x0=xk;

xk=x0 + lambdak*sk;

res=norm(grad(xk));

fprintf('------the%d-th iteration, the residual is %f, the lambdak is %f\n\n', k ,res,lambdak);

end

x_star=xk;

disp('the optimal solution is');

x_star

function lambda = nonexact(x0,s0)

g0 = grad(x0);

f0=f(x0);

a=0;

c1=0.1;

c2=0.5;

lambdak=1;

d=-c1*lambdak*g0'*s0;

xk=x0 + lambdak*s0;

f1=f(xk);

e=f0-f1;

while d>e

lambdak=(lambdak+a)/2;

xk=x0 + lambdak*s0;

fk=f(xk);

e=f0-fk;

d=-c1*lambdak*g0'*s0;

end

lambda=lambdak;

function g=grad(x)

g=zeros(2,1);

g(1)=2*x(1)*exp(x(1)^2+x(2)^2)+1;

g(2)=4*x(2)+2*x(2)*exp(x(1)^2+x(2)^2);

function fa=f(x)

fa=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);

function G=GRAND(x)

G=zeros(2,2);

G(1)=2*exp(x(1)^2+x(2)^2)+4*x(1)^2*exp(x(1)^2+x(2)^2);

G(2)=4*x(1)*x(2)*exp(x(1)^2+x(2)^2);

G(3)=4*x(1)*x(2)*exp(x(1)^2+x(2)^2);

G(4)=4+(2+4*x(2)^2)*exp(x(1)^2+x(2)^2);

在命令窗口中输入x0=[0;1];eps=1e-3;newton1(x0,eps) 输出结果为：

取

H为BFGS公式时：

k

程序如下：

function x_star = bfgs1(x0,eps)

g0 = grad(x0);

res0 = norm(g0);

res=res0;

k = 0;

xk =x0;

while res>eps

k=k+1;

if k==1

H=eye(length(x0));

else

H0=H;

mu=1+dg'*H0*dg/(dx'*dg);

H=H0+(mu*dx*dx'-H0*dg*dx'-dx*dg'*H0)/(dx'*dg);

end

sk=-H*grad(xk);

lambdak = nonexact(xk,sk);

x0=xk;

xk=x0 + lambdak*sk;

dx=xk-x0;

dg=grad(xk)-grad(x0);

res=norm(grad(xk));

fprintf('------the%d-th iteration, the residual is %f, the lambdak is %f\n\n', k ,res,lambdak);

end

x_star=xk;

fm=f(xk);

disp('the optimal solution is');

x_star

fm

function lambda = nonexact(x0,s0)

g0 = grad(x0);

f0=f(x0);

a=0;

c1=0.1;

c2=0.5;

lambdak=1;

d=-c1*lambdak*g0'*s0;

xk=x0 + lambdak*s0;

f1=f(xk);

e=f0-f1;

while d>e

lambdak=(lambdak+a)/2;

xk=x0 + lambdak*s0;

fk=f(xk);

e=f0-fk;

d=-c1*lambdak*g0'*s0;

end

lambda=lambdak;

function g=grad(x)

g=zeros(2,1);

g(1)=2*x(1)*exp(x(1)^2+x(2)^2)+1;

g(2)=4*x(2)+2*x(2)*exp(x(1)^2+x(2)^2);

function fa=f(x)

fa=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);

在命令窗口中输入x0=[0;1];eps=1e-3;bfgs1(x0,eps) 输出结果为：

程序如下：

function x=qpact(H,c,Ae,be,Ai,bi,x0)

epsilon=1.0e-9;

err=1.0e-6;

k=0;

x=x0;

n=length(x);

kmax=1.0e3;

ne=length(be);

ni=length(bi);

lamk=zeros(ne+ni,1);

index=ones(ni,1);

for (i=1:ni)

if(Ai(i,:)*x>bi(i)+epsilon)

index(i)=0;

end

end

while (k

Aee=[];

if(ne>0)

Aee=Ae;

end

for(j=1:ni)

if(index(j)>0)

Aee=[Aee; Ai(j,:)];

end

end

gk=H*x+c;

[m1,n1] = size(Aee);

[dk,lamk]=qsubp(H,gk,Aee,zeros(m1,1));

if(norm(dk)<=err)

y=0.0;

if(length(lamk)>ne)

[y,jk]=min(lamk(ne+1:length(lamk)));

end

if(y>=0)

exitflag=0;

else

exitflag=1;

for(i=1:ni)

if(index(i)&(ne+sum(index(1:i)))==jk)

index(i)=0;

break;

end

end

end

k=k+1;

else

exitflag=1;

alpha=1.0; %求步长%

tm=1.0;

for(i=1:ni)

if((index(i)==0)&(Ai(i,:)*dk<0))

tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);

if(tm1

tm=tm1;

ti=i;

end

end

end

alpha=min(alpha,tm);

x=x+alpha*dk;

if(tm<1) %修正有效集%

index(ti)=1;

end

end

if(exitflag==0)

break;

end

k=k+1;

end

fn=f(x);

x_star=x;

x_star

fn

% 求解子问题%

function [x,lambda]=qsubp(H,c,Ae,be)

ginvH=pinv(H);

[m,n]=size(Ae);

if(m>0)

rb=Ae*ginvH*c + be;

lambda=pinv(Ae*ginvH*Ae')*rb;

x=ginvH*(Ae'*lambda-c);

else

x=-ginvH*c;

lambda=0;

end

function fn=f(x)

fn=(x(1)-1)^2+(x(2)-2.5)^2;

在命令窗口中输入

H=[2 0;0 2];

c=[-2;-5];

Ae=[ ]; be=[ ];

Ai=[1 -2;-1 -2;-1 2;1 0;0 1]; bi=[-2;-6;-2;0 ;0];

x0=[0;0];

qpact(H,c,Ae,be,Ai,bi,x0)

输出结果为：

程序如下：

function x_star = lagrange(x0,epsn)

v=1;

u=0;

m=h(u,x0);

while m>epsn

if (u+4*d(x0))>0

u=u+4*d(x0);

else

u=0;

end

v=v+4*e(x0);

x0=bfgs(v,u,x0);

m=h(u,x0);

end

x_tar=x0;

x_tar

fn=fn(x0);

fn

function star = bfgs(v0,u0,x0)

eps=1e-4;

g0 = grad(v0,u0,x0);

res0 = norm(g0);

res=res0;

k = 0;

xk =x0;

while res>eps

k=k+1;

sk=-grad(v0,u0,x0);

lambdak = nonexact(v0,u0,xk,sk);

x0=xk;

xk=x0 + lambdak*sk;

res=norm(grad(v0,u0,xk));

end

star=xk;

function lambda = nonexact(v0,u0,x0,s0) g0 = grad(v0,u0,x0);

f0=f(v0,u0,x0);

a=0;

c1=0.1;

c2=0.5;

lambdak=1;

p=-c1*lambdak*g0'*s0;

xk=x0 + lambdak*s0;

f1=f(v0,u0,xk);

q=f0-f1;

while p>q

lambdak=(lambdak+a)/2;

xk=x0 + lambdak*s0;

fk=f(v0,u0,xk);

q=f0-fk;

p=-c1*lambdak*g0'*s0;

end

lambda=lambdak;

function g=grad(v,u,x)

g=zeros(2,1);

g(1)=4-2*x(1)*v-8*x(1)*(25-x(1)^2-x(2)^2)+(u+4*(10*x(1)-x(1)^2+10*x(2)-x(2)^2-34))*(10-2*x(1));

g(2)=-2*x(2)-2*x(2)*v-8*x(2)*(25-x(1)^2-x(2)^2)+(u+4*(10*x(1)-x(1)^2+10*x(2)-x (2)^2-34))*(10-2*x(2));

function f=f(v,u,x)

f=4*x(1)-x(2)^2-12+v*(25-x(1)^2-x(2)^2)+2*(25-x(1)^2-x(2)^2)^2+0.125*((u+4*(10 *x(1)-x(1)^2+10*x(2)-x(2)^2-34))^2-u^2);

function h=h(u,x)

if (10*x(1)-x(1)^2+10*x(2)-x(2)^2-34)>(-u/4)

h=(25-x(1)^2-x(2)^2)^2+(10*x(1)-x(1)^2+10*x(2)-x(2)^2-34)^2;

else

h=(25-x(1)^2-x(2)^2)^2+(-u/4)^2;

end

function d=d(x)

d=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;

function e=e(x)

e=25-x(1)^2-x(2)^2;

function fn=fn(x)

fn=4*x(1)-x(2)^2-12;

在命令窗口中输入x0=[1;1];epsn=1e-4;lagrange(x0,epsn)

输出结果为：

先写出求解二次规划子问题的程序

程序如下：

function [d,mu,lam,val,k]=qpsubp(dfk,Bk,Ae,hk,Ai,gk)

n=length(dfk);

l=length(hk);

m=length(gk);

gamma=0.05;

epsilon=1.0e-6;

rho=0.5;

sigma=0.2;

ep0=0.05;

mu0=0.05*zeros(l,1);

lam0=0.05*zeros(m,1);

d0=ones(n,1);

u0=[ep0;zeros(n+l+m,1)];

z0=[ep0; d0; mu0;lam0,];

k=0;

z=z0;

ep=ep0;

d=d0;

mu=mu0;

lam=lam0;

while (k<=150)

dh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);

if(norm(dh)

break;

end

A=JacobiH(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);

b=beta(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk,gamma)*u0-dh;

dz=A\b;

if(l>0&m>0)

de=dz(1); dd=dz(2:n+1); du=dz(n+2:n+l+1); dl=dz(n+l+2:n+l+m+1);

end

if(l==0)

de=dz(1); dd=dz(2:n+1); dl=dz(n+2:n+m+1);

end

if(m==0)

de=dz(1); dd=dz(2:n+1); du=dz(n+2:n+l+1);

end

i=0; mk=0;

while (mk<=20)

if(l>0&m>0)

dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk);

end

if(l==0)

dh1=dah(ep+rho^i*de,d+rho^i*dd,mu,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk);

end

if(m==0)

dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam,dfk,Bk,Ae,hk,Ai,gk);

end

if(norm(dh1)<=(1-sigma*(1-gamma*ep0)*rho^i)*norm(dh))

mk=i;

break;

end

i=i+1;

if(i==20)

mk=10;

end

end

alpha=rho^mk;

if(l>0&m>0)

ep=ep+alpha*de; d=d+alpha*dd;

mu=mu+alpha*du; lam=lam+alpha*dl;

end

if(l==0)

ep=ep+alpha*de; d=d+alpha*dd;

lam=lam+alpha*dl;

end

if(m==0)

ep=ep+alpha*de; d=d+alpha*dd;

mu=mu+alpha*du;

end

k=k+1;

end

function p=phi(ep,a,b)

p=a+b-sqrt(a^2+b^2+2*ep^2);

function dh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk)

n=length(dfk); l=length(hk); m=length(gk);

dh=zeros(n+l+m+1,1);

dh(1)=ep;

if(l>0&m>0)

dh(2:n+1)=Bk*d-Ae'*mu-Ai'*lam+dfk;

dh(n+2:n+l+1)=hk+Ae*d;

for(i=1:m)

dh(n+l+1+i)=phi(ep,lam(i),gk(i)+Ai(i,:)*d);

end

end

if(l==0)

dh(2:n+1)=Bk*d-Ai'*lam+dfk;

for(i=1:m)

dh(n+1+i)=phi(ep,lam(i),gk(i)+Ai(i,:)*d);

end

end

if(m==0)

dh(2:n+1)=Bk*d-Ae'*mu+dfk;

dh(n+2:n+l+1)=hk+Ae*d;

end

dh=dh(:);

function bet=beta(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk,gamma)

dh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);

bet=gamma*norm(dh)*min(1,norm(dh));

function [dd1,dd2,v1]=ddv(ep,d,lam,Ai,gk)

m=length(gk);

dd1=zeros(m,m); dd2=zeros(m,m); v1=zeros(m,1);

for(i=1:m)

fm=sqrt(lam(i)^2+(gk(i)+Ai(i,:)*d)^2+2*ep^2);

dd1(i,i)=1-lam(i)/fm;

dd2(i,i)=1-(gk(i)+Ai(i,:)*d)/fm;

v1(i)=-2*ep/fm;

end

function A=JacobiH(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk)

n=length(dfk); l=length(hk); m=length(gk);

A=zeros(n+l+m+1,n+l+m+1);

[dd1,dd2,v1]=ddv(ep,d,lam,Ai,gk);

if(l>0&m>0)

A=[1, zeros(1,n), zeros(1,l), zeros(1,m);

zeros(n,1), Bk, -Ae', -Ai';

zeros(l,1), Ae, zeros(l,l), zeros(l,m) ;

v1, dd2*Ai, zeros(m,l), dd1];

end

if(l==0)

A=[1, zeros(1,n), zeros(1,m);

zeros(n,1), Bk, -Ai';

v1, dd2*Ai, dd1];

end

if(m==0)

A=[1, zeros(1,n), zeros(1,l);

zeros(n,1), Bk, -Ae';

zeros(l,1), Ae, zeros(l,l)];

end

序列二次规划程序为：

function [x,mu,lam,val,k]=sqpm(x0,mu0,lam0)

maxk=100; %最大迭代次数

n=length(x0);

l=length(mu0);

m=length(lam0);

rho=0.5;

eta=0.1;

B0=eye(n);

x=x0;

mu=mu0;

lam=lam0;

Bk=B0;

sigma=0.8;

epsilon1=1e-6;

epsilon2=1e-5;

[hk,gk]=cons(x);

dfk=df1(x);

[Ae,Ai]=dcons(x);

Ak=[Ae; Ai];

k=0;

while(k

[dk,mu,lam]=qpsubp(dfk,Bk,Ae,hk,Ai,gk); %求解子问题

mp1=norm(hk,1)+norm(max(-gk,0),1);

if(norm(dk,1)

break;

end %检验终止准则

deta=0.05; %罚参数更新

tau=max(norm(mu,inf),norm(lam,inf));