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

为什么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程序

约束优化算法：拉格朗日乘子法

转载于:https://my.oschina.net/keyven/blog/804688