• 最速下降法

最速下降法又称为梯度法,是1847年由著名数学家Cauchy给出的。他是解析法中最古老的一种,其他解析方法或是它的变形,或是受它的启发而得到的,因此它是最优化方法的基础。

设无约束问题中的目标函数

image.png

一阶连续可微

最速下降法的基本思想是:从当前点

 

image.png

出发,取函数

image.png

在点处下降最快的方向作为我们的搜索方向

由Taylor展式知

image.png

略去t的高阶无穷小项不计,可见取image.png

函数值下降得最多。

  • 最速下降法步骤

image.png

由以上计算步骤可知,最速下降法迭代终止时,求得的是目标函数驻点的一个近似点。

image.png

 

function [ X,FMIN,K ] = zuisuxiajiang( f,x,x0,e )
%    [ X,FMIN,N ] =zuisuxiajiang()法求解无约束问题
%    X         极小点
%    FMIN      极小值
%    K         迭代次数
%    f         问题函数
%    x         变量
%    x0        初始点
%    e         终止误差
count=0;
td=jacobian(f,x)';
while norm(subs(td,x,x0))>e
    P=-subs(td,x,x0);
    syms r
    y=x0+r*P;
    ft(r)=subs(f,x,y);        
    [r0]=fibonacci(ft,0,100,0.01);
    x0=x0+r0*P;
    count=count+1;
end
X=x0;
FMIN=subs(f,x,x0);
K=count;
end
  • 牛顿法

 

对于一维搜索函数,假定已给出极小点的一个较好的近似点a0,因为一个连续可微的函数在极小点附近与一个二次函数很接近,所以可以在a0点附近用一个二次函数来逼近函数,即在点a0将f(a)进行泰勒展开,并保留到二次项,有

image.png

然后以二次函数的极小点作为极小点的一个新近似点,根据极值必要条件  

image.png


image.png

  • 牛顿法计算步骤

⑴给定初始点a0,控制误差e,令k=0

⑵计算f(x)在ak点的一阶和二阶导数

⑶利用牛顿法迭代公式求ak+1

⑷若|ak+1-ak|≤e,则求得近似解a*=ak+1,停止计算,否则作第⑸步

⑸令k=k+1,然后转第⑵步

function [ X,FMIN,K ] = ysNewton( f,x,x0,e )
%    [ X,FMIN,N ] =ysNewton()原始牛顿法求解无约束问题
%    X         极小点
%    FMIN      极小值
%    K         迭代次数
%    f         问题函数
%    x         变量
%    x0        初始点
%    e         终止误差
count=0;
td=jacobian(f,x)';
H=jacobian(td',x);
while norm(subs(td,x,x0))>e
    P=-subs(H,x,x0)^(-1)*subs(td,x,x0);
    x0=x0+P;
    count=count+1;
end
X=x0;
FMIN=subs(f,x,x0);
K=count;
end

 

 

  • 共轭梯度法

基本原理

 

设A是n阶对称正定矩阵

image.png

是k个A共轭的非零向量,则这个向量组线性无关。

 

设有函数

image.png

中A是n阶对称正定矩阵

image.png

是一组A共轭向量。以任意的

image.png

为初始点,依次沿

image.png

进行搜索,得到点

image.png

image.png

是函数image.png

上的极小点,其中image.png

是由

image.png

生成的子空间。特别地,当k=n时

image.png

上唯一的极小点。

推论:在上述定理条件下,必有

image.png

定理3:对于正定二次型函数

image.png

FR算法在次一维搜索后即终止,并且对所有的

image.png

下列关系成立

image.png

  • 共轭梯度法算法步骤

image.png

 

function [ X,FMIN,K ] = gongetidu( f,x,x0,e )
%    [ X,FMIN,N ] =gongetidu()共轭梯度法求解无约束问题
%    X         极小点
%    FMIN      极小值
%    K         迭代次数
%    f         问题函数
%    x         变量
%    x0        初始点
%    e         终止误差
count=1;
td=jacobian(f,x)';
H=jacobian(td',x);
if norm(subs(td,x,x0))>e
    P=-subs(td,x,x0);
    r0=-subs(td,x,x0)'*P/(P'*H*P);
    x0=x0+r0*P;
else 
    x0;
end
while norm(double(subs(td,x,x0)))>e
    b0=subs(td,x,x0)'*subs(td,x,x0)/(P'*P);
    P=-subs(td,x,x0)+b0*P;
    r0=-subs(td,x,x0)'*P/(P'*H*P);
    x0=x0+r0*P;
    count=count+1;
end
X=x0;
FMIN=subs(f,x,x0);
K=count;
end

实例:

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

初始点x0=(0,3)T

允许误差e=0.1

键入:

 

syms x1 x2
>> f=(x1-2)^2+(x1-2*x2)^2;
>> x=[x1;x2];
>> x0=[0;3];
>> e=0.1;

 

1.最速下降法;

 

   [X,FMIN,N]=zuisuxiajiang(f,x,x0,e)
X =
    1.9763
    0.9818
FMIN =
   7.2076e-04
N =
    10

 

2.牛顿法

 

>> [X,FMIN,N]=ysNewton(f,x,x0,e)
X =
     2
     1
FMIN =
     0
N =
     1

 

3.共轭梯度法

 

[X,FMIN,N]=gongetidu(f,x,x0,e)
X =
 2
 1
FMIN =
0
N =
     2

可知对于二次正定函数newton法只需一次迭代就得到正确结果,共轭梯度法只需进行两次迭代就得出正确结果。但最速下降法却迭代了10次,虽然一维搜索存在误差,但实际上最速下降法也需迭代多次。

 

  • 外罚函数法

 

考虑约束非线性最优化问题

image.png

其中

image.png

都是定义在上的实值连续函数.记可行域为S

对于等式约束问题

image.png

定义罚函数

image.png

其中参数image.png

是很大的正常数.这样,可将问题转化为无约束问题

image.png

对于不等式约束问题

image.png

定义罚函数

image.png

 

其中参数image.png

是很大的正常数.这样,可将问题转化为无约束问题

image.png

对于一般的约束最优化问题,定义罚函数

image.png

其中参数image.png

是很大的正数,即有下列形式

image.png

image.png是满足下列条件的实值连续函数

image.png

典型取法为

image.png

image.png

是给定的常数,通常取作12

这样,可将约束问题(10.1.1)转化为无约束问题

image.png

  • 外罚函数算法步骤

image.png

  • 内罚函数原理

对于不等式约束问题(10.1.4),其可行域S的内部

image.png

为了保持迭代点始终含于IntS,我们引入另一种罚函数

image.png

显然,罚函数的作用是对企图脱离可行域的点给予惩罚,相当于在可行域的边界设置了障碍,不让迭代点穿越到可行域之外,因此也称为障碍函数.

两种常用的形式为

image.png

分别称为倒数障碍函数和对数障碍函数.

  • 内罚函数步骤

image.png

用罚函数法解

min f(x)=(x1-1)^2+x2^

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

 

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]很接近。

投影梯度法

投影梯度法原理

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

f

unction [ 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   迭代次数
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
        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
   A1=A1(1:k,:);
   b1=b1(1:k,:);
else
    A1=[];
end
if q>0
   A2=A2(1:q,:);
   b2=b2(1:q,:);
end 
end

 数值算例:

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