【最优化方法】实验四 约束最优化方法的MATLAB实现

实验的目的和要求：通过本次实验使学生较为熟练使用MATLAB软件，并能利用该软件进行约束最优化方法的计算。

实验内容：

１、罚函数法的MATLAB实现

２、可行方向法的MATLAB实现

学习建议：

本次实验就是要通过对一些具体问题的分析进一步熟悉软件的操作并加深对理论知识的理解。

重点和难点：

可行点和辅助函数选取。

一 罚函数法

用罚函数法解min f(x)=(x1-1)^2+x22

​ S.T. g(x)=x2-1>=0

编写下面m文件

fahanshu.m

syms x1 x2 

f=(x1-1)^2+x2^2;

g=x2-1;

x=[x1;x2];

x0=[0;0];

e=0.0001;M=1;

while abs(subs(g,x,x0))>e

  if subs(g,x,x0)<0

​    Q=f+M*g^2;

  else 

​    Q=f;

  end

  x0=xzNewton(Q,x,x0,0.0001);

  M=10*M;

end

x0

运行得：fahanshu

x0 =

1.0000

0.9999

与理论值x=[1;1]很接近。

二 投影梯度法

1.梯度投影法基本原理和步骤

思想：当迭代点是可行域的内点时，将目标函数负梯度作为搜索方向，当迭代点在可行域边界上时，将目标函数负梯度在可行域边界上的投影作为搜索方向。无论何种情况，所构造的方向都是可行下降方向。然后在可行域内沿该方向进行最优一维搜索得到新的迭代点。

MATLAB实现:

2.代码及数值算例:

(1) 程序源代码:

function [ X,FMIN,K ] = tidutouying( f,A,b,x1,x,e )

%  [ X,FMIN,K ] = tidutouying( f,A,b,x1,e ) 梯度投影法

%  f  目标函数

%  A  约束矩阵  b 右端项

%  x1 初始点 x 自由变量

%  e  精度要求

%  X  极小点

%  FMIN 极小值

%  K  迭代次数

%  张超编写与2014/5/3

count=1;

n=length(x1);

tf=jacobian(f,x)';

while 1

  [A1,A2,b1,b2,k]=fenjie(A,b,x1);

  while 1

  M=A1;

  if isempty(M)

​    P=eye(n);

  else

​    P=eye(n)-M'*(M*M')^(-1)*M;

  end

  Pk=-P*subs(tf,x,x1);

  if norm(Pk)<=e

​    if isempty(M)

​      x1;break;

​    else

​      W=(M*M')^(-1)*M*subs(tf,x,x1);

​      u=W;

​      if min(u)<0

​        for i=1:length(u)

​          if u(i)==min(u)

​            j=i;

​          end

​        end

​        A1(j,:)=[];

​      else 

​        x1;break;

​      end

​    end

  else

​    b_=b2-A2*x1;

​    P_=A2*Pk;

​    for i=1:length(P_)

​      if P_(i)<0

​        r(i)=b_(i)/P_(i);

​      else 

​        r(i)=10000;

​      end

​    end

​    rmax=min(r);

​    syms t

​    y=x1+t*Pk;

​    ft(t)=subs(f,x,y);    

​    [r1]=find0618(ft,0,double(rmax),0.00001);

​    x1=x1+r1*Pk;   

​    break;

  end

  end

  count=count+1;

  if isempty(M)    

   break;

  end

  if min(u)>=0

​    break;

  end

end

 X=x1;

 FMIN=subs(f,x,X);

 K=count;

end

 

function [A1,A2,b1,b2,k ] = fenjie( A,b,x )

% 分解起作用约束

A=A;

 b=b;

 x0=x;

 k=0;q=0;

s=size(A);

A1=zeros(s(1),s(2));

 A2=zeros(s(1),s(2));

 b1=zeros(s(1),1);

 b2=zeros(s(1),1);

 for i=1:s(1)

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

  if abs(gi)<0.0000001 

​    k=k+1;

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

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

  else 

​    q=q+1;

​    A2(q,:)=A(i,:);

​    b2(q,1)=b(i);

  end

 end

if k>0

  A1=A1(1:k,:);

  b1=b1(1:k,:);

else

  A1=[];

end

if q>0

  A2=A2(1:q,:);

  b2=b2(1:q,:);

end 

end

(2) 数值算例:

Min f(x)= 2*x1^2+2*x2^2 – 2*x1*x2 – 4*x1 – 6*x2

S.T. –x1 – x2>= –2

–x1 – 5*x2>= –5

x1>=0

x2>=0

初值x0=[0;0]

在matlab command window里输入

syms x1 x2 

f=2*x1^2+2*x2^2-2*x1*x2-4*x1-6*x2;

A=[-1 -1;-1 -5;1 0;0 1];

b=[-2;-5;0;0];

x1=[0;0];

x=[x1;x2];

e=0.01;

[X,FMIN,K]=tidutouying(f,A,b,x1,x,e)

 

X =

  1.1292

  0.7742

FMIN =

  -7.1613

N =

  16