为什么80%的码农都做不了架构师?>>>
一、设计方法选择适当的初始点,用牛顿法和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<max_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)<e_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
牛顿法求解最优化计算Matlab
DFP算法及Matlab程序
约束优化算法:拉格朗日乘子法