大連理工優(yōu)化方法增廣拉格朗日方法MATLAB程序.doc

《大連理工優(yōu)化方法增廣拉格朗日方法MATLAB程序.doc》由會員分享，可在線閱讀，更多相關(guān)《大連理工優(yōu)化方法增廣拉格朗日方法MATLAB程序.doc（5頁珍藏版）》請在裝配圖網(wǎng)上搜索。

上機大作業(yè)II 定義目標(biāo)函數(shù)fun function f=fun(x) x1=x(1); x2=x(2); f=4*x1-x2^2-12; 定義目標(biāo)函數(shù)梯度函數(shù)dfun function f=dfun(x) x2=x(2); f=[4;-2*x2]; 定義等式約束函數(shù)hf function qua=hf(x) qua=25-x(1)^2-x(2)^2; 定義等式約束函數(shù)梯度函數(shù)dhf function qua=dhf(x) qua=[-2*x(1);-2*x(2)]; 定義不等式約束函數(shù)gfun function inq=gfun(x) inq=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34; 定義不等式約束梯度數(shù)dgf function inq=dgf(x) inq=[10-2*x(1);10-2*x(2)]; 定義增廣拉格朗日函數(shù)mpsi function psi=mpsi(x,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma) f=feval(fun,x); he=feval(hf,x); gi=feval(gfun,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*sigma*s1; s2=0.0; for i=1:m s3=max(0.0, lambda(i) - sigma*gi(i)); s2=s2+s3^2-lambda(i)^2; end psi=psi+s2/(2.0*sigma); 定義增廣拉格朗日函數(shù)梯度函數(shù)dmpsi function dpsi=dmpsi(x,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma) dpsi=feval(dfun,x); he=feval(hf,x); gi=feval(gfun,x); dhe=feval(dhf,x); dgi=feval(dgf,x); l=length(he); m=length(gi); for i=1:l dpsi=dpsi+(sigma*he(i)-mu(i))*dhe(:,i); end for i=1:m dpsi=dpsi+(sigma*gi(i)-lambda(i))*dgi(:,i); end 定義BFGS法函數(shù)函數(shù)bfgs function [x,val,k]=bfgs(mpsi,dmpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma) maxk=1000; rho=0.5; sigma1=0.4; epsilon1=1e-4; k=0; n=length(x0); Bk=eye(n); while(k0) Bk=Bk-((Bk*sk)*sk*Bk)/(sk*Bk*sk)+(yk*yk)/(yk*sk); end k=k+1; x0=x; end val=feval(mpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma); 定義增廣拉格朗日乘子法函數(shù)multphr function answer=multphr(fun,hf,gfun,dfun,dhf,dgf,x0) maxk=5000; sigma=2.0; eta=2.0; theta=0.8; k=0; ink=0; epsilon=1e-4; x=x0; he=feval(hf,x); gi=feval(gfun,x); l=length(he); m=length(gi); mu=0.1*ones(l,1); lambda=0.1*ones(m,1); btak=10; btaold=10; while(btak>epsilon&&kepsilon if(k>=2&&btak > theta*btaold) sigma=eta*sigma; end for i=1:l mu(i)=mu(i)-sigma*he(i); end for i=1:m lambda(i)=max(0.0,lambda(i)-sigma*gi(i)); end end k=k+1; btaold=btak; x0=x; end f=feval(fun,x); x f mu lambda k 運行求解 >> x0=[0;0] x0 = 0 0 >> multphr(fun,hf,gfun,dfun,dhf,dgf,x0) x = 1.00128148956437 4.89871784708758 f = -31.9923105871169 mu = 1.01559644571312 lambda = 0.754451167977228 k = 4