大连理工大学庞丽萍最优化方法MATLAB程序

班级：优化1班授课老师：庞丽萍姓名：学号：

第二章

12.（1）用修正单纯形法求解下列LP问题:

>> clear

>> A=[1 2 1 1 0 0;1 2 3 0 1 0;2 1 5 0 0 1];

[m,n]=size(A);

b=[10;15;20];

r=[-1 -2 -3 1];

c=[-1 -2 -3 1];

bs=[3:3];

nbs=[1:4];

a1=A(:,3);

T=A(:,bs);

a2=inv(T)*a1;

b=inv(T)*b;

A=[eye(m),a2];

B=eye(m);

xb=B\b;

cb=c(bs);

cn=c(nbs);

con=1;

M=zeros(1);

while con

M=M+1;

t=cb/B;

r=c-t*A;

if all(r>=0)

x(bs)=xb;

x(nbs)=0;

fx=cb*xb;

disp(['当前解是最优解，minz=',num2str(fx)])

disp('对应的最优解为,x=')

disp(x)

break

end

rnbs=r(nbs);

kk=find(rnbs==min(rnbs));

k=kk(1);

Anbs=A(:,nbs);

yik=B\Anbs(:,k);

xb=B\b;%yi0

if all(yik<=0)

disp('此LP问题无有限的最优解，计算结束',x)

disp(xb)

break

else

i=find(yik>0);

w=abs(xb(i,1)./yik(i,1));

l=find(w==min(w));

rr=min(l);

yrrk=yik(rr,1);

Abs=A(:,bs);

D=Anbs(:,k);

Anbs(:,k)=Abs(:,rr);

Abs(:,rr)=D;

F=bs(rr);

bs(rr)=nbs(k);

nbs(k)=F;

AA=[Anbs,Abs];

EE=eye(m);

EE(:,rr)=-yik./yrrk;

Errk=EE;

Errk(rr,rr)=1/yrrk;

BB=Errk/B;

B=inv(BB);

cb=c(:,bs);

xb=Errk*xb;

x(bs)=xb;

x(nbs)=0;

fx=cb*xb;

end

if M>=1000

disp('此问题无有限最优解')

break

end

end

%结果

当前解是最优解，minz=-15

对应的最优解为,x=

2.5000 2.5000 2.5000 0

第三章

30题DFP算法求函数极小点的计算程序

function[x,val,k]=dfp(fun,gfun,x0)

%功能:用DFP算法求解无约束问题:minf(x)

%输入:x0是初始点,fun,gfun分别是目标函数及其梯度

%输出:x,val分别是近似最优点和最优值,k是迭代次数.

maxk=1e5;%给出最大迭代次数

rho=0.55;sigma=0.4;epsilon=1e-5;

k=0;n=length(x0);

Hk=inv(feval('Hess',x0));%Hk=eye(n);

while(k

gk=feval(gfun,x0);%计算梯度

if(norm(gk)

dk=-Hk*gk;%解方程组,计算搜索方向

m=0;mk=0;

while(m<20)%用Armijo搜索求步长

if(feval(fun,x0+rho^m*dk)

mk=m;break;

end

m=m+1;

end

%DFP校正

x=x0+rho^mk*dk;

sk=x-x0;

yk=feval(gfun,x)-gk;

if(sk'*yk>0)

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

end

k=k+1;

x0=x;

end

val=feval(fun,x0);

%习题26的程序调用方式及结果：

function y=fun(x)

%UNTITLED Summary of this function goes here

% Detailed explanation goes here

y=(x(1)-1)^2+5*(x2-x(1)^2)^2

end

function y=gfun(x)

%UNTITLED Summary of this function goes here

% Detailed explanation goes here

y=[diff(y,x1) diff(y,x2)]

end

x0=[2 0]’;

[x,val,k]=dfp(fun,gfun,x0)

%结果

x =

1.00000 1.00000

val=

k =

6

%习题27的程序调用方式及结果：

function y=fun(x)

%UNTITLED Summary of this function goes here % Detailed explanation goes here

y=x1+2*x（2）^2+exp(x(1)^2+x(2)^2)

end

function y=gfun(x)

%UNTITLED Summary of this function goes here % Detailed explanation goes here

y=[diff(y,x1) diff(y,x2)]

end

x0=[1 0]’;

[x,val,k]=dfp(fun,gfun,x0)

%结果

x =

-0.41936 0

val=

0.77291

k =

5

36题编写Hooke-Jeeves 方法求函数极小点的计算程序。

clear;

syms x1d alpha yibxlou y1k j e1e2...e(n);

f(x)=4*x1^2+x2^2-40*x1-12*x2+136;

x(1)=[4 8]’;alpha=1;d=1;yibxlou=0.2;

y1=x1;k=j;j=1;[e1 e2 ... e(n)]=[1 0 ... 0;0 1 0 ... 0;...;0 ... 1] while (f(y(j)+de(j))

y(j+1)=y(j)+de(j);

if (j

j=j+1;

else j=n;

if (f(y(n+1)>=f(x(k))

if(d<=yibuxlou)

x*=x(k);

else

d=d/2;

y1=x(k);

x(k+1)=x(k);

k=k+1;

j=1;

else

if(f(y(j)-de(j))

y(j+1)=y(j)-de(j);

else

y(j+1)=y(j);

end;

%结果

y=

x*=

5.00000

6.00000

第四章

24题，编写梯度投影法的MA TLAB程序

function [x,minf]=minRosen(f,A,b,x0,var,eps) %目标函数：f;

%约束矩阵：A;

%约束右端力量：b;

%初始可行点：x0;

%自变量向量：var;

%精度：eps;

%目标函数取最小值时的自变量值：x;

%目标函数的最小值：minf;

format long;

if nargin == 5

eps=1.0e-6;

end

syms l;

x0=transpose(x0);

n=length(var);

sz=size(A);

m=sz(l);

gf=jacobian(f,var);

bConti=l;

while bConti

k=0;

s=0;

A1=A;

A2=A;

b1=b;

b2=b;

for i=1:m

dfun=A(i,:)*x0-b(i);

if abs(dfun)<0.000000001

k=k+1;

A1(k,:)=A(i,:);

b1(k,1)=b(i);

else

s=s+1;

A2(s,:)=A(i,:);

b2(s,1)=b(i);

end

if k>0

A1=A1(1:k,:);

b1=b1(1:k,:);

end

if s>0

A2=A2(1:s,:);

b2=b2(1:s,:);

end

while l

P=eye(n,n);

if k>0

tM=transpose(A1);

P=P-tM*inv(A1*tM)*A1;

end

gv=Funval(gf,var,x0);

gv=transpose(gv);

d=-P*gv;

if d==0

if k==0

x=x0;

bConti=0;

break;

else

w=inv(A1*tM)*A1*gv;

if w>=0

x=x0;

bConti=0;

break;

else

[u,index]=min(w);

sA1=size(A1);

if sA1(1)==1

k=0;

else

k=sA1(2);

A1=[A1(1:(index-1),:);A1((index+1):sA1(2),:)]; %去掉A1对应的行

end

end

end

else

break;

end

end

y1=x0+l*d

tmpf=Funval(f,val,y1);

bb=b2-A2*x0;

dd=A2*d;

if dd>=0

[tmpI,lm]=minJT(tmpf,0,0.1); % else

lm=inf;

for i=1:length(dd)

if dd(i)<0

if bb(i)/dd(i)

lm=bb(i)/dd(i);

end

end

end

end

[xm,minf]=minHJ(tmpf,0,lm,1.0e-14);

tol=norm(xm*d);

if tol

x=x0;

break;

end

x0=x0+xm*d;

end

minf=Funval(f,var,x)

function fv = Funval(f,varvec,varval)

var = findsym(f);

varc = findsym(varvec);

s1 = length(var);

s2 = length(varc);

m =floor((s1-1)/3+1);

varv = zeros(1,m);

if s1 ~= s2

for i=0: ((s1-1)/3)

k = findstr(varc,var(3*i+1));

index = (k-1)/3;

varv(i+1) = varval(index+1);

end

fv = subs(f,var,varv);

else

fv = subs(f,varvec,varval);

function [x,minf] = minHJ(f,a,b,eps) format long;

if nargin == 3

eps = 1.0e-6;

end

l = a + 0.382*(b-a);

u = a + 0.618*(b-a);

k=1;

tol = b-a;

while tol>eps && k<100000

fl = subs(f , findsym(f), l);

fu = subs(f , findsym(f), u);

if fl > fu

a = l;

l = u;

u = a + 0.618*(b - a);

else

b = u;

u = l;

l = a + 0.382*(b-a);

end

k = k+1;

tol = abs(b - a);

end

if k == 100000

disp('找不到最小值！');

x = NaN;

minf = NaN;

return;

end

x = (a+b)/2;

minf = subs(f, findsym(f),x);

format short;

function [minx,maxx] = minJT(f,x0,h0,eps) format long;

if nargin == 3

eps = 1.0e-6;

end

x1 = x0;

k = 0;

h = h0;

while 1

x4 = x1 + h;

k = k+1;

f4 = subs(f, findsym(f),x4);

f1 = subs(f, findsym(f),x1);

if f4 < f1

x2 = x1;

x1 = x4;

f2 = f1;

f1 = f4;

h = 2*h;

else

if k==1

h = -h;

x2 = x4;

f2 = f4;

else

x3 = x2;

x2 = x1;

x1 = x4;

break;

end

end

end

minx = min(x1,x3);

maxx = x1+x3 - minx;

format short;

syms x1 x2 x3;

f=x1^2+x1*x2+2*x2^2+2*x3^2+2*x2*x3+4*x1+6*x2+12*x3; [x,mf]=minRosen(f,[1 1 1 ;1 1 -2],[6;-2],[1 1 3],[t s])

x =

0 1.00000 1.00000

mf =

14.00000

35题序列二次规划法的MA TLAB程序

function[x,mu,lam,val,k]=sqpm(x0,mu0,lam0)

%功能:用基于拉格朗日函数Hesse阵的SQP方法求解约束优化问题:

%minf(x)s.t.h˙i(x)=0,i=1,...,l.

%输入:x0是初始点,mu0是乘子向量的初始值

%输出:x,mu分别是近似最优点及相应的乘子,

%val是最优值,mh是约束函数的模,k是迭代次数.

maxk=100;%最大迭代次数

n=length(x0);l=length(mu0);m=length(lam0);

rho=0.5;eta=0.1;B0=eye(n);

x=x0;mu=mu0;lam=lam0;

Bk=B0;sigma=0.8;

epsilon1=1e-6;

epsilon2=1e-5;

[hk,gk]=cons(x);dfk=df1(x);

[Ae,Ai]=dcons(x);Ak=[Ae;Ai];

k=0;

while(k

[dk,mu,lam]=qpsubp(dfk,Bk,Ae,hk,Ai,gk);%求解子问题

mp1=norm(hk,1)+norm(max(-gk,0),1);

if(norm(dk,1)

break;end%检验终止准则

deta=0.05;%罚参数更新

tau=max(norm(mu,inf),norm(lam,inf));

if(sigma*(tau+deta)<1)

sigma=sigma;

else

sigma=1.0/(tau+2*deta);

end

im=0;%Armijo搜索

while(im<=20)

if(phi1(x+rho^im*dk,sigma)-phi1(x,sigma)

break;

end

im=im+1;

if(im==20),mk=10;end

end

alpha=rho^mk;

x1=x+alpha*dk;

[hk,gk]=cons(x1);dfk=df1(x1);

[Ae,Ai]=dcons(x1);Ak=[Ae;Ai];

lamu=pinv(Ak)’*dfk;%计算最小二乘乘子

if(l>0&m>0)

mu=lamu(1:l);

lam=lamu(l+1:l+m);

end

if(l==0),mu=[];lam=lamu;end

if(m==0),mu=lamu;lam=[];end

sk=alpha*dk;%更新矩阵Bk

yk=dlax(x1,mu,lam)-dlax(x,mu,lam);

if(sk’*yk>0.2*sk’*Bk*sk)

theta=1;

else

theta=0.8*sk’*Bk*sk/(sk’*Bk*sk-sk’*yk);

end

zk=theta*yk+(1-theta)*Bk*sk;

Bk=Bk+zk*zk’/(sk’*zk)-(Bk*sk)*(Bk*sk)’/(sk’*Bk*sk);

x=x1;k=k+1;

end

val=f1(x);

%p=phi1(x,sigma)

%dd=norm(dk)

%%%%%%%%l1精确价值函数%%%%%%%

functionp=phi1(x,sigma)f=f1(x);

[h,g]=cons(x);gn=max(-g,0);

l0=length(h);m0=length(g);

if(l0==0),p=f+1.0/sigma*norm(gn,1);end

if(m0==0),p=f+1.0/sigma*norm(h,1);end

if(l0>0&m0>0)

p=f+1.0/sigma*(norm(h,1)+norm(gn,1));

end

%%%%%价值函数的方向导数%%%%%

Function dp=dphi1(x,sigma,d)

df=df1(x);[h,g]=cons(x);gn=max(-g,0);

l0=length(h);m0=length(g);

if(l0==0),dp=df’*d-1.0/sigma*norm(gn,1);end

if(m0==0),dp=df’*d-1.0/sigma*norm(h,1);end

if(l0>0&m0>0)

dp=df’*d-1.0/sigma*(norm(h,1)+norm(gn,1));

end

%%%%%%%%%拉格朗日函数L(x,mu)%%%%%%%%%%%%% functionl=la(x,mu,lam)

f=f1(x);%调用目标函数文件

[h,g]=cons(x);%调用约束函数文件

l0=lemgth(h);m0=length(g);

if(l0==0),l=f-lam*g;endif(m0==0),l=f-mu’*h;

end

if(l0>0&m0>0)

l=f-mu’*h-lam’*g;

end

%%%%%%%%%拉格朗日函数的梯度%%%%%%%%%%%%%

Function dl=dlax(x,mu,lam)

df=df1(x);%调用目标函数梯度文件

[Ae,Ai]=dcons(x);%调用约束函数Jacobi矩阵文件

[m1,m2]=size(Ai);[l1,l2]=size(Ae);

if(l1==0),dl=df-Ai’*lam;end

if(m1==0),dl=df-Ae’*mu;end

if(l1>0&m1>0),dl=df-Ae’*mu-Ai’*lam;

end

首先编写四个m函数:

%%%%%%%%目标函数f(x)%%%%%%%%%%%

Function f=f1(x)%f1.m

f=x1*x2*x3*x4*x5;

%%%%%%%目标函数f(x)的梯度%%%%%%%%

Function df=df1(x)%df1.m

df=[ x2*x3*x4*x5, x1 *x3*x4*x5, x1*x2 *x4*x5, x1*x2*x3 *x5, x1*x2*x3*x4]’; %%%%%%%%%约束函数%%%%%%%%%%%%%%

Function [h,g]=cons(x) %cons.m

h=[x(1)^2+x(2)^2+x(3)^2+x(4)^2+x(5)^2-10;x(2)*x(3)-5*x(4)*x(5);x(1)^3+x(2)^3+1]; %%%%%%%%约束函数Jacobi矩阵%%%%%%%%

Function [dh,dg]=dcons(x)%dcons.m

dh=[2*x1,2*x2,2*x3,2*x4,2*x5;x3,x2,-5*x5,-5*x4;3*x1^2,3*x2^2];

在Matlab命令窗口依次输入下列命令:

x0=[3 3]’;

mu0=[0]’;

lam0=[0 0]’;

[x,mu,lam,val,k]=sqpm(x0,mu0,lam0)

得到计算结果:

x=

-1.71719

1.59575

1.82721

-0.763621

-0.763622

mu=

-3.9894

lam=

1.0e-007*

0.4248

0.0039

val=

-2.91963 k=

8