一、设计方法选择适当的初始点，用牛顿法和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 % xkend
二、用乘子法求解下列模型，初始点： （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=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