优化方法MATLAB编程——大连理工大学

优化方法上机大作业

学 院： 姓 名： 学 号：

指导老师： 肖 现 涛

第一题

源程序如下：

function zy_x = di1ti(x)

%di1ti是用来求解优化作业第一题的函数。 x0=x; yimuxulong=0.000001; g0=g(x0);s0=-g0; A=2*ones(100,100); k=0; while k<100

lanmed=-(g0)'*s0/(s0'*A*s0); x=x0+lanmed*s0;

1

g=g(x); k=k+1;

if norm(g)fprintf('After %d iterations,obtain the optimal solution.\\n \\n The optimal solution is \\n %f.\\n\\nThe optimal \"x\" is \"ans\".',k,f(x) )

break; end

miu=norm(g)^2/norm(g0)^2; s=-g+miu*s0; g0=g; s0=s;x0=x; end

function f=f(x)

f=(x'*ones(100,1))^2-x'*ones(100,1); function g=g(x)

g=(2*x'*ones(100,1))*ones(100,1)-ones(100,1); 代入x0,运行结果如下： >> x=zeros(100,1); >> di1ti(x)

After 1 iterations,obtain the optimal solution.

2

The optimal solution is -0.250000.

The optimal \"x\" is \"ans\". ans =0.005*ones(100,1). 第二题

1. 最速下降法。 源程序如下：

function zy_x=di2titidu(x) %该函数用来解大作业第二题。 x0=x; yimuxulong=1e-5; k=0; g0=g(x0); s0=-g0; while k>=0

if norm(g0)3

break; else

lanmed=10;c=0.1;i=0; while i>=0&i<100 x=x0+lanmed*s0;

if f(x)>(f(x0)+c*lanmed*g0'*s0) lanmed=lanmed/2; i=i+1; else break; end end

x=x0+lanmed*s0; x0=x; g0=g(x); s0=-g0; k=k+1; end end zy_x=x;

4

zyj=f(x);

fprintf('after %d iterations,obtain the optimal solution.\\n\\nThe optimal solution is %f.\\n\\n The optimal \"x\" is \"ans\".\\n',k,zyj);

function f=f(x) x1=[1 0 0 0]*x; x2=[0 1 0 0]*x; x3=[0 0 1 0]*x; x4=[0 0 0 1]*x;

f=(x1-1)^2+(x3-1)^2+100*(x2-x1^2)^2+100*(x4-x3^2)^2; function g=g(x) x1=[1 0 0 0]*x; x2=[0 1 0 0]*x; x3=[0 0 1 0]*x; x4=[0 0 0 1]*x;

g=[2*(x1-1)-400*x1*(x2-x1^2);200*(x2-x1^2);2*(x3-1)-400*x3*(x4-x3^2);200*(x4-x3^2)];

>> x=[-1.2 1 -1.2 1]';

5

>> di2titidu(x)

after 5945 iterations,obtain the optimal solution. The optimal solution is 0.000000. The optimal \"x\" is \"ans\". ans = 1.0000 1.0000 1.0000 1.0000 2. 牛顿法 源程序如下：

function zy_x=di2tinewton(x) %该函数用来解大作业第二题。 x0=x; yimuxulong=1e-5; k=0; g0=g(x0); h0=h(x0);s0=-inv(h0)*g0; while k>=0&k<1000 if norm(g0)6

x0=x; g0=g(x); h0=h(x); s0=-inv(h0)*g0; k=k+1; end end zy_x=x; zyj=f(x);

fprintf('after %d iterations,obtain the optimal solution.\\n\\nThe optimal solution is %f.\\n\\n The optimal \"x\" is \"ans\".\\n',k,zyj);

function f=f(x) x1=[1 0 0 0]*x; x2=[0 1 0 0]*x; x3=[0 0 1 0]*x; x4=[0 0 0 1]*x;

f=(x1-1)^2+(x3-1)^2+100*(x2-x1^2)^2+100*(x4-x3^2)^2;

7

function g=g(x) x1=[1 0 0 0]*x; x2=[0 1 0 0]*x; x3=[0 0 1 0]*x; x4=[0 0 0 1]*x;

g=[2*(x1-1)-400*x1*(x2-x1^2);200*(x2-x1^2);2*(x3-1)-400*x3*(x4-x3^2);200*(x4-x3^2)];

function h=h(x) x1=[1 0 0 0]*x; x2=[0 1 0 0]*x; x3=[0 0 1 0]*x; x4=[0 0 0 1]*x;

h=[2+1200*x1^2-400*x2 -400*x1 0 0;-400*x1 200 0 0;0 0 2+1200*x3^2-400*x4 -400*x3;0 0 -400*x3 200];

代入初始值，运行结果如下： >> x=[-1.2 1 -1.2 1]'; >> di2tinewton(x)

after 6 iterations,obtain the optimal solution. The optimal solution is 0.000000. The optimal \"x\" is \"ans\".

8

ans = 1.0000 1.0000 1.0000 1.0000

可以看出，用Newton法经过6次迭代就能求出最优解。 3. BFGS法 源程序如下：

function zy_x=di2tiBFGS(x) %该函数用来解大作业第二题。 x0=x; yimuxulong=1e-5; k=0; g0=g(x0); H0=eye(4);s0=-H0*g0; while k>=0&k<100 if norm(g0)lanmed=10;c=0.1;i=0; while i>=0&i<100 x=x0+lanmed*s0;

if f(x)>(f(x0)+c*lanmed*g0'*s0)

9

lanmed=lanmed/2; i=i+1; else break; end end

x=x0+lanmed*s0; dete_x=x-x0; dete_g=g(x)-g0;

miu=1+dete_g'*H0*dete_g/(dete_x'*dete_g);

H=H0+(miu*dete_x*dete_x'-H0*dete_g*dete_x'-dete_x*dete_g'*H0)/(dete_x'*dete_g);

s=-H*g(x); x0=x; s0=s; H0=H; g0=g(x); k=k+1; end end

10

zy_x=x; zyj=f(x);

fprintf('after %d iterations,obtain the optimal solution.\\n\\nThe optimal solution is %f.\\n\\n The optimal \"x\" is \"ans\".\\n',k,zyj);

function f=f(x) x1=[1 0 0 0]*x; x2=[0 1 0 0]*x; x3=[0 0 1 0]*x; x4=[0 0 0 1]*x;

f=(x1-1)^2+(x3-1)^2+100*(x2-x1^2)^2+100*(x4-x3^2)^2; function g=g(x) x1=[1 0 0 0]*x; x2=[0 1 0 0]*x; x3=[0 0 1 0]*x; x4=[0 0 0 1]*x;

g=[2*(x1-1)-400*x1*(x2-x1^2);200*(x2-x1^2);2*(x3-1)-400*x3*(x4-x3^2);200*(x4-x3^2)];

11

代入初始值，计算结果如下： >> x=[-1.2 1 -1.2 1]'; >> di2tiBFGS(x)

after 53 iterations,obtain the optimal solution. The optimal solution is 0.000000. The optimal \"x\" is \"ans\". ans = 1.0000 1.0000 1.0000 1.0000 第三题

1. 惩罚函数法 源程序如下：

12

function zy_x=di3ti(x) %该函数用来解大作业第三题。 x0=x; M=100; c=4; m=1; while m>0

g0=g(x0,M); yimuxulong=1e-5;k=0;s0=-inv(H(x0,M))*g0; while k>=0

if norm(g0)x=x0+s0; %牛顿法； x0=x; g0=g(x,M);

s0=-inv(H(x0,M))*g0; k=k+1; end end

if max([abs(h(x)),g1(x),g2(x),g3(x)])<0.5 break; else M=M*c;

13

m=m+1; end end zy_x=x; zyj=f(x);

fprintf('after %d iterations,obtain the optimal solution.\\n\\nThe optimal solution is %f.\\n\\n The optimal \"x\" is \"ans\".\\n',m,zyj);

function F=F(x,M) x1=[1 0]*x; x2=[0 1]*x;

F=4*x1-x2^2-12+M*(h^2+g1^2+g2^2+g3^2); function g=g(x,M) x1=[1 0]*x; x2=[0 1]*x;

g=[4+M*(-4*(25-x1^2-x2^2)*x1+2*(10*x1-x1^2+10*x2-x2^2-34)*(10-2*x1)+2*x1);-2*x2+M*(-4*(25-x1^2-x2^2)*x2+2*(10*x1-x1^2+10*x2-x2^2-34)*(10-2*x2)+2*x2)];

function H=H(x,M) x1=[1 0]*x;

14

x2=[0 1]*x;

H=[M*(24*x1^2-120*x1+8*x2^2-40*x2+238),M*(16*x1*x2-40*x1-40*x2+200);M*(16*x1*x2-40*x1-40*x2+200),-2+M*(24*x2^2-120*x2+8*x1^2-40*x1+238)];

function f=f(x) x1=[1 0]*x; x2=[0 1]*x; f=4*x1-x2^2-12; function h=h(x) x1=[1 0]*x; x2=[0 1]*x; h=25-x1^2-x2^2; function g1=g1(x) x1=[1 0]*x; x2=[0 1]*x;

g=10*x1-x1^2+10*x2-x2^2-34; if g<0 g1=g; else g1=0; end

5

1

function g2=g2(x) x1=[1 0]*x; x2=[0 1]*x; if x1>=0 g2=0; else g2=x1; end

function g3=g3(x) x1=[1 0]*x; x2=[0 1]*x; if x2>=0 g3=0; else g3=x2; end

代入任意初始值，运算结果如下。 >> x=rand(2,1); >> di3ticf(x)

after 1 iterations,obtain the optimal solution.

16

The optimal solution is -31.490552. The optimal \"x\" is \"ans\". ans = 1.0024 4.8477 2. 乘子法 源程序如下：

function [x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0)

%功能: 用乘子法解一般约束问题: min f(x), s.t. h(x)=0, g(x).=0

%输入: x0是初始点, fun, dfun分别是目标函数及其梯度；

% hf, dhf分别是等式约束（向量）函数及其Jacobi矩阵的转置；

% gf, dgf分别是不等式约束（向量）函数及其Jacobi矩阵的转置；

%输出: x是近似最优点，mu, lambda分别是相应于等式约束和不等式约束的乘子向量;

% output是结构变量, 输出近似极小值f, 迭代次数, 内迭代次数等

maxk=500;

17

c=2.0;

eta=2.0;theta=0.8;

k=0;ink=0;

epsilon=0.00001;

x=x0;he=feval(hf,x);gi=feval(gf,x);

n=length(x);l=length(he);m=length(gi);

mu=zeros(l,1);lambda=zeros(m,1);

btak=10;btaold=10;

while(btak>epsilon&&k%调用BFGS算法程序求解无约束

ink=ink+ik;

he=feval(hf,x);gi=feval(gf,x);

18

子问题

[x,ival,ik]=bfgs('mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c);

btak=0;

for i=1:l

btak=btak+he(i)^2;

end

%更新乘子向量

for i=1:m

temp=min(gi(i),lambda(i)/c);

btak=btak+temp^2;

end

btak=sqrt(btak);

if btak>epsilon

if k>=2&&btak>theta*btaold c=eta*c;

9

1

end

for i=1:l

mu(i)=mu(i)-c*he(i);

end

for i=1:m

lambda(i)=max(0,lambda(i)-c*gi(i));

end

k=k+1;

btaold=btak;

x0=x;

end

end

f=feval(fun,x);

20

output.fval=f;

output.iter=k;

%增广拉格朗日函数

function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c)

f=feval(fun,x);he=feval(hf,x);gi=feval(gf,x);

l=length(he);m=length(gi);

psi=f;s1=0;

for i=1:l

psi=psi-he(i)*mu(i);

s1=s1+he(i)^2;

end

psi=psi+0.5*c*s1;

s2=0;

21

for i=1:m

s3=max(0,lambda(i)-c*gi(i));

s2=s2+s3^2-lambda(i)^2;

end

psi=psi+s2/(2*c);

%不等式约束函数文件g1.m

function gi=g1(x)

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

%目标函数的梯度文件df1.m

function g=df1(x)

g=[4, -2*x(2)]';

%等式约束（向量）函数的Jacobi矩阵（转置）文件dh1.m

function dhe=dh1(x)

22

dhe=[-2*x(1), -2*x(2)]'

%不等式约束（向量）函数的Jacobi矩阵（转置）文件dg1.m

function dgi=dg1(x)

dgi=[10-2*x(1), 10-2*x(2)]';

function [x,val,k]=bfgs(fun,gfun,x0,varargin)

maxk=500;

rho=0.55;sigma=0.4;epsilon=0.00001;

k=0;n=length(x0);

Bk=eye(n);

while(kgk=feval(gfun,x0,varargin{:});

if(norm(gk)break;

23

end

dk=-Bk\\gk;

m=0;mk=0;

while(m<20)

newf=feval(fun,x0+rho^m*dk,varargin{:});

oldf=feval(fun,x0,varargin{:});

if(newfmk=m;

break;

end

m=m+1;

end

x=x0+rho^mk*dk;

24

sk=x-x0;

yk=feval(gfun,x,varargin{:})-gk;

if(yk'*sk>0)

Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);

end

k=k+1;x0=x;

end

val=feval(fun,x0,varargin{:});

结果

x=[2 2]';

[x,mu,lambda,output]=multphr('fun','hf','gf1','df','dh','dg',x0)

x =

1.0013

25

4.87

mu =

0.7701

lambda =

0

0

0.9434

output =

fval: -31.9923iter: 4

26