牛顿法的matlab实现,牛顿法、DFP以及Lagrange乘子法的MATLAB实现

一、设计方法选择适当的初始点，用牛顿法和DFP法求解下列模型，结合计算结果，对两种算法的优缺点进行比较，精度要求

.

(1) 牛顿法matlab实现

function newton()

syms x1 x2 x3 x4

f=(x1+2*x2)^2+3*(x3-x4)^2+(x2-2*x3)^2+2*(x1-x4)^2+9;

v=[x1,x2,x3,x4];

df=jacobian(f,v);

df=df.';

df

G=jacobian(df,v);

G

epson=1e-4;

xm=[1,1,1,1]';

gl=subs(df,{x1,x2,x3,x4},{xm(1,1),xm(2,1),xm(3,1),xm(4,1)});

Gl=subs(G,{x1,x2,x3,x4},{xm(1,1),xm(2,1),xm(3,1),xm(4,1)});

k=0;

while(norm(gl)>epson)

xm=xm-inv(Gl)*gl;

gl=subs(df,{x1,x2,x3,x4},{xm(1,1),xm(2,1),xm(3,1),xm(4,1)});

Gl=subs(G,{x1,x2,x3,x4},{xm(1,1),xm(2,1),xm(3,1),xm(4,1)});

k=k+1;

end

k

xm'

end

(2) DFP算法matlab实现

function dfp()

syms x1 x2 x3 x4 t;

f=(x1+2*x2)^2+3*(x3-x4)^2+(x2-2*x3)^2+2*(x1-x4)^2+9;

fx1=diff(f,x1);

fx2=diff(f,x2);

fx3=diff(f,x3);

fx4=diff(f,x4);

fi=[fx1 fx2 fx3 fx4];

x0=[1,1,1,1];

f0=subs(f,[x1 x2 x3 x4],x0);

g0=subs(fi,[x1 x2 x3 x4],x0);

H0=eye(4);

xk=x0;

fk=f0;

gk=g0;

Hk=H0;

k=1;

while(norm(gk)>1e-4)

pk=-Hk*gk';

xk=xk+t*pk';

f_t=subs(f,[x1 x2 x3 x4],xk);

df_t=diff(f_t,t);

tk=solve(df_t);

xk=subs(xk,t,tk);

fk=subs(f,[x1 x2 x3 x4],xk);

gk0=gk;

gk=subs(fi,[x1 x2 x3 x4],xk);

yk=gk-gk0;

sk=tk*pk';

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

k=k+1;

k

xk

end

% xk

end

二、用乘子法求解下列模型，初始点： (2, 2, 2) ；精度要求

.

function [X,FVAL]=AL_main(x_al,r_al,N_equ,N_inequ)

global r_al pena N_equ N_inequ;

pena=10;

c_scale=2;

cta=0.5;

e_al=1e-4;

max_itera=100;

out_itera=1;

while out_itera

x_al0=x_al;

r_al0=r_al;

compareFlag=compare(x_al0);

[X,FVAL]=fminunc(@AL_obj,x_al0);

x_al=X;

if compare(x_al)

break;

end

if compare(x_al)/compareFlag>=cta

pena=c_scale*pena;

end

[h,g]=constrains(x_al);

for i=1:N_equ

r_al(i)=r_al(i)-pena*h(i);

end

for i=1:N_inequ

r_al(i+N_equ)=max(0,(r_al(i+N_equ)-pena*g(i)));

end

out_itera=out_itera+1;

end

X=x_al;

FVAL=obj(X);

end

function f=AL_obj(x)

global r_al pena N_equ N_inequ;

h_equ=0;

h_inequ=0;

[h,g]=constrains(x);

for i=1:N_equ

h_equ=h_equ-h(i)*r_al(i)+(pena/2)*h(i).^2;

end

for i=1:N_inequ

h_inequ=h_inequ+(0.5/pena)*(max(0,(r_al(i)-pena*g(i))).^2-r_al(i).^2);

end

f=obj(x)+h_equ+h_inequ;

end

function f=compare(x)

global r_al pena N_equ N_inequ;

h_equ=0;

h_inequ=0;

[h,g]=constrains(x);

for i=1:N_equ

h_equ=h_equ+h(i).^2;

end

for i=1:N_inequ

h_inequ=h_inequ+(min(g(i),r_al(i+N_equ)/pena)).^2;

end

f=sqrt(h_equ+h_inequ);

end

function [h,g]=constrains(x)

h=x(1)^2+x(2)^2+x(3)^2-12;

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

g(2)=x(1);

g(3)=x(2);

g(4)=x(3);

end

function f=obj(x)

f=20-x(1)^2-2*x(2)^2-x(3)^2-x(1)*x(2)-x(1)*x(3);

end

function main()

clear all;

clc;

x_al=[2,2,2];

r_al=[1,1,1,1,1];

N_equ=1;

N_inequ=4;

[X,FVAL]=AL_main(x_al,r_al,N_equ,N_inequ)

end

Reference