最近偶尔翻阅一本写的不错的最优化理论教材，该书讲得很详细，很透彻。我对非线性规划理论又有了全新的认识，发现牛顿法可以说是是无约束优化中最重要的方法，其他方法：LM方法，高斯牛顿，拟牛顿法，共轭梯度法可以说是对牛顿法的扩展。准备闲来无事时将牛顿法的原理以及求解过程用数例好好再过一遍。

 

适用对象：二阶可微函数

 

1. 牛顿法的几何意义本质

       在原函数的某一点处用一个二次函数近似原函数，然后用这个二次函数的极小值点作为原函数的下一个迭代点。

上面这句话也说明，若原函数本身是一个二次函数，则牛顿法一步就能到达极小点或鞍点。若原函数本身是一个二次正定函数，则牛顿法一步到达最小值点。

 

2. 牛顿法的代数意义

 

梯度与黑塞矩阵分别由下列符号表示：

  

  

设任意点为 Xk，下个一迭代点位 Xk+Sk， 该迭代点在 Xk 处的泰勒展开式为：



用下个迭代点的值代替改点的值，即：



因此：



所以，迭代方向为：



该方向又称作牛顿方向。

 

3. 牛顿法的二阶收敛性

若初始点 x0 充分靠近极值点 x*，并且极值点 x* 的黑塞矩阵非奇异，并且黑塞矩阵在极值点附近 Lipschitz 连续，则牛顿法具有二阶收敛性。

注：Lipschitz 连续是一种比普通连续性更强的连续，它限制了函数的改变速度。对于函数可行域的任意两点，存在一个常数 K，使得：



证明：

由黑塞矩阵的非奇异性与连续性知道，在 x*附近，存在一个常数 M，对于任意的 k，使得



而：



上式右端成立，是因为 g(x*)=0。继续：



上式是因为：



继续，再利用到 Lipstchitz 连续，得到：



因此，牛顿迭代法二阶收敛。

一、牛顿法
    在博文“优化算法——牛顿法(Newton Method)”中介绍了牛顿法的思路，牛顿法具有二阶收敛性，相比较最速下降法，收敛的速度更快。在牛顿法中使用到了函数的二阶导数的信息，对于函数，其中表示向量。在牛顿法的求解过程中，首先是将函数在处展开，展开式为：



其中，，表示的是目标函数在的梯度，是一个向量。，表示的是目标函数在处的Hesse矩阵。省略掉最后面的高阶无穷小项，即为：



上式两边对求导，即为：



    在基本牛顿法中，取得最值的点处的导数值为，即上式左侧为。则：



求出其中的：

从上式中发现，在牛顿法中要求Hesse矩阵是可逆的。  
    当时，上式为：



此时，是否可以通过，，和模拟出Hesse矩阵的构造过程？此方法便称为拟牛顿法(QuasiNewton)，上式称为拟牛顿方程。在拟牛顿法中，主要包括DFP拟牛顿法，BFGS拟牛顿法。
二、DFP拟牛顿法
1、DFP拟牛顿法简介
        DFP拟牛顿法也称为DFP校正方法，DFP校正方法是第一个拟牛顿法，是有Davidon最早提出，后经Fletcher和Powell解释和改进，在命名时以三个人名字的首字母命名。
对于拟牛顿方程：



化简可得：



令，可以得到：



在DFP校正方法中，假设：



2、DFP校正方法的推导
    令：，其中均为的向量。，。
    则对于拟牛顿方程可以简化为：



将代入上式：



将代入上式：





已知：为实数，为的向量。上式中，参数和解的可能性有很多，我们取特殊的情况，假设，。则：



代入上式：





令，，则：






则最终的DFP校正公式为：



3、DFP拟牛顿法的算法流程
    设对称正定，由上述的DFP校正公式确定，那么对称正定的充要条件是。
    在博文“优化算法——牛顿法(Newton Method)”中介绍了非精确的线搜索准则：Armijo搜索准则，搜索准则的目的是为了帮助我们确定学习率，还有其他的一些准则，如Wolfe准则以及精确线搜索等。在利用Armijo搜索准则时并不是都满足上述的充要条件，此时可以对DFP校正公式做些许改变：



       DFP拟牛顿法的算法流程如下：

4、求解具体的优化问题
    求解无约束优化问题



其中，。

python程序实现：
function.py
#coding:UTF-8
'''
Created on 2015年5月19日

@author: zhaozhiyong
'''

from numpy import *

#fun
def fun(x):
    return 100 * (x[0,0] ** 2 - x[1,0]) ** 2 + (x[0,0] - 1) ** 2

#gfun
def gfun(x):
    result = zeros((2, 1))
    result[0, 0] = 400 * x[0,0] * (x[0,0] ** 2 - x[1,0]) + 2 * (x[0,0] - 1)
    result[1, 0] = -200 * (x[0,0] ** 2 - x[1,0])
    return result

dfp.py
#coding:UTF-8
'''
Created on 2015年5月19日

@author: zhaozhiyong
'''

from numpy import *
from function import *

def dfp(fun, gfun, x0):
    result = []
    maxk = 500
    rho = 0.55
    sigma = 0.4
    m = shape(x0)[0]
    Hk = eye(m)
    k = 0
    while (k < maxk):
        gk = mat(gfun(x0))#计算梯度
        dk = -mat(Hk)*gk
        m = 0
        mk = 0
        while (m < 20):
            newf = fun(x0 + rho ** m * dk)
            oldf = fun(x0)
            if (newf < oldf + sigma * (rho ** m) * (gk.T * dk)[0,0]):
                mk = m
                break
            m = m + 1
        
        #DFP校正
        x = x0 + rho ** mk * dk
        sk = x - x0
        yk = gfun(x) - gk
        if (sk.T * yk > 0):
            Hk = Hk - (Hk * yk * yk.T * Hk) / (yk.T * Hk * yk) + (sk * sk.T) / (sk.T * yk)
        
        k = k + 1
        x0 = x
        result.append(fun(x0))
    
    return result

testDFP.py
#coding:UTF-8
'''
Created on 2015年5月19日

@author: zhaozhiyong
'''

from bfgs import *
from dfp import dfp

import matplotlib.pyplot as plt  

x0 = mat([[-1.2], [1]])
result = dfp(fun, gfun, x0)

n = len(result)
ax = plt.figure().add_subplot(111)
x = arange(0, n, 1)
y = result
ax.plot(x,y)

plt.show()

5、实验结果

第6章
回到目录
第8章



第7章-二阶时滞多智能体系统分组一致性
7.1 引言
7.2 预备知识
7.3 问题描述与分析
定理7.1
7.4 例子与数值仿真
7.5 本章小结

7.1 引言
分组一致性主要考虑两个关键因素：

分组一致性的具体收敛状态
分组收敛的速度

一般采用牛顿第二定律进行描述

牛顿第二定律又称为“加速度定律”。

牛顿第二定律被誉为经典力学的灵魂。在经典力学里，它能够主导千变万化的物体运动与精彩有序的物理现象。牛顿第二定律的用途极为广泛，它可以用来设计平稳地耸立于云端的台北101摩天大厦，也可以用来计算从地球发射火箭登陆月球的运动轨道。

牛顿第二定律是一个涉及到物体运动的理论，根据这定律，任意物体的运动所出现的改变，都是源自于外力的施加于这物体。这理论导致了经典力学的诞生，是科学史的一个里程碑，先前只是描述自然现象的理论不再被采纳，取而代之的是这个创立了一种理性的因果关系架构的新理论。

牛顿第二定律

相比于一阶时滞网络，二阶系统一致性需要考虑两个物理量，即速度和位置。
牵制算法
LaSalle 不变集原理
Lyapunov 稳定性理论
牵制控制方法

牵制控制的基本思想是: 通过有选择地对网络中的少部分节点施加控制而使得整个网络具有期望的行为。比起对所有节点进行控制的方式，由于较小的资源花费，对网络的小部分节点应用局部反馈控制的牵制控制是更有吸引力和令人期待的。

Lyapunov 第一定理
Hopf 分叉理论

当考虑系统动态环节后, 系统将可能遇到 Hopf 分岔 (Hopf bifurcation, HB) 现象, 相应的失稳为振荡型失稳。

微分方程理论在自动控制、航天技术、生态生物等方面一直有着广泛的应用，在这些实际应用中，系统通常都是一些含有参数的微分方程组。考虑如下形式的系统：

$dX/dt=f(X, λ,μ) \tag{1}$

系统（1）的解显然随参数μ的变化而变化。如果 $λ$ 在 $λ$ 等于一个确定值的一个小邻域内变化时，系统（1）在相空间的相图拓扑结构发生了变化，那么就称系统发生了分岔，称 $λ$ 为分岔参数，$λ$ 的确定值为分岔值。

而今在应用数学中，Hopf 分岔理论已经成为研究微分方程小振幅周期解产生和消亡的经典工具。因此，对 Hopf 分岔的研究是十分有意义的。

7.2 预备知识
包含 $n+m$ 个智能体的二阶连续系统如下：

$\dot{x}_i(t) = v_i(t)\\
\dot{v}(t) = u_i(t),\quad i=1,2,\cdots,m+n \tag{7.1}$
7.3 问题描述与分析
不同时滞影响下二阶系统的一致性协议：

$\begin{aligned}
u_i(t) = &\red{\alpha}\sum_{v_j\in N_i} a_{ij} (\red{x}_j(t-T_{ij}) - \red{x}_i(t-T)) + \\ &\red{\beta}\sum_{v_j\in N_i} a_{ij} (\red{v}_j(t-T_{ij}) - \red{v}_i(t-T)),\quad i=1,2,\cdots,N \tag{7.2}
\end{aligned}$
$\alpha,\beta$ 分别为该系统的耦合强度。
耦合强度
定理7.1
考虑包含 $n+m(n,m>1)$ 个智能体的多智能体系统（7.4），基于入度平衡的假设条件（A1）和（A2），当网络拓扑结构中包含全局可达节点时，若满足 $\tilde{d}_i(\alpha \cos(\omega_{i0} T_i) + \beta \omega_{i0} \sin(\omega_{i0} T_i)) < \frac{1}{2} \omega_{i0}^2$ 成立，则称系统能渐进实现分组一致。其中 $\tilde{d}_i = \sum_{k=1,k\ne i}^{m+n} a_{ik}$，$\omega_{i0}$ 表示圆盘中心 $G_{i0}(j\omega)=\tilde{d}_i \frac{e^{-j\omega T_i}}{j\omega} \frac{\alpha+\beta j \omega}{j\omega}$ 的奈奎斯特曲线与复平面实轴的交点，同时 $\tan(\omega_{i0} T_i) = \frac{\beta}{\alpha} \omega_{i0}$。
7.4 例子与数值仿真


以上工作，几经尝试，未出现一次收敛的情况。
第二天，换一种新思路来进行


当未进行任何更改时（即权重全都保持1，未添加任何时滞，并且为删除掉不同组的节点状态），直接出现了收敛状态。
结果如下图所示



在此基础上更改权重后，

状态信息如下：



速度信息如下：


更改权重之后，再将不同组节点的状态值删除掉



更改 aij 为入度
对比之前的发现一个有趣的现象，虽然状态同样还是发散，但是分组情况发生了变化。
之前使用出度时，1 4 5 为一组，2 3 为一组；

现在使用入度时，1 2 3 为一组，4 5 为一组。



更改拓扑图后，结果正确



在此基础更改入度为出度




在此做个小总结。

首先，系统的一致性研究主要分为以下几个主要因素：

系统拓扑关系图
权重的设定
控制协议
输入时滞的影响
通信时滞的影响

在此基础上，搭建仿真模型，进行控制协议的验证。
通过以上几组实验对比，可基本得出时由于系统拓扑关系图不满足一致性要求所引起的。
书中所给的定理7.1，系统若想一致，需要满足以下三个条件：


基于入度平衡的假设条件（A1）和（A2）；


当网络拓扑结构中包含全局可达节点时；


若满足 $\tilde{d}_i(\alpha \cos(\omega_{i0} T_i) + \beta \omega_{i0} \sin(\omega_{i0} T_i)) < \frac{1}{2} \omega_{i0}^2$ 成立，则称系统能渐进实现分组一致。


通过观察书中所给的系统拓扑关系图，均满足以上要求，并且未设置任何时滞，也一定满足条件（3）。
既然已经满足上述三个条件，但是还是不能收敛。接下来的实验，先验证之前搭建的系统模型没有问题，之后再研究拓扑图为什么不能收敛的问题。
之前的模型，将拓扑关系图更改之后，结果也都正常，包括节点分组情况也正确。

位置信息如下：



速度状态如下：


保持时滞为零，其他条件与采用第6章实验2所用拓扑图条件一致的情况下，仅更改为第7章拓扑图，还是无法收敛。算是再次验证了拓扑图是错误的。



7.5 本章小结
18. Ren W, Beard R. Consensus seeking in multiagent systems under dynamically changing interaction topologies[J]. IEEE Transactions on Automatic Control, 2005, 50(5):655-661.

转载地址：http://blog.csdn.net/google19890102/article/details/45848439
一、牛顿法
    在博文“优化算法——牛顿法(Newton
 Method)”中介绍了牛顿法的思路，牛顿法具有二阶收敛性，相比较最速下降法，收敛的速度更快。在牛顿法中使用到了函数的二阶导数的信息，对于函数，其中表示向量。在牛顿法的求解过程中，首先是将函数在处展开，展开式为：






其中，，表示的是目标函数在的梯度，是一个向量。，表示的是目标函数在处的Hesse矩阵。省略掉最后面的高阶无穷小项，即为：






上式两边对求导，即为：






    在基本牛顿法中，取得最值的点处的导数值为，即上式左侧为。则：






求出其中的：



从上式中发现，在牛顿法中要求Hesse矩阵是可逆的。
  
    当时，上式为：






此时，是否可以通过，，和模拟出Hesse矩阵的构造过程？此方法便称为拟牛顿法(QuasiNewton)，上式称为拟牛顿方程。在拟牛顿法中，主要包括DFP拟牛顿法，BFGS拟牛顿法。
二、DFP拟牛顿法
1、DFP拟牛顿法简介
        DFP拟牛顿法也称为DFP校正方法，DFP校正方法是第一个拟牛顿法，是有Davidon最早提出，后经Fletcher和Powell解释和改进，在命名时以三个人名字的首字母命名。
对于拟牛顿方程：






化简可得：






令，可以得到：






在DFP校正方法中，假设：






2、DFP校正方法的推导
    令：，其中均为的向量。，。
    则对于拟牛顿方程可以简化为：






将代入上式：






将代入上式：










已知：为实数，为的向量。上式中，参数和解的可能性有很多，我们取特殊的情况，假设，。则：






代入上式：










令，，则：











则最终的DFP校正公式为：






3、DFP拟牛顿法的算法流程
    设对称正定，由上述的DFP校正公式确定，那么对称正定的充要条件是。
    在博文“优化算法——牛顿法(Newton
 Method)”中介绍了非精确的线搜索准则：Armijo搜索准则，搜索准则的目的是为了帮助我们确定学习率，还有其他的一些准则，如Wolfe准则以及精确线搜索等。在利用Armijo搜索准则时并不是都满足上述的充要条件，此时可以对DFP校正公式做些许改变：






       DFP拟牛顿法的算法流程如下：


4、求解具体的优化问题
    求解无约束优化问题






其中，。


python程序实现：

function.py



[python] view
 plain copy


 


#coding:UTF-8  

''''' 

Created on 2015年5月19日 

 

@author: zhaozhiyong 

'''  

  

from numpy import *  

  

#fun  

def fun(x):  

    return 100 * (x[0,0] ** 2 - x[1,0]) ** 2 + (x[0,0] - 1) ** 2  

  

#gfun  

def gfun(x):  

    result = zeros((2, 1))  

    result[0, 0] = 400 * x[0,0] * (x[0,0] ** 2 - x[1,0]) + 2 * (x[0,0] - 1)  

    result[1, 0] = -200 * (x[0,0] ** 2 - x[1,0])  

    return result  

dfp.py



[python] view
 plain copy


 


#coding:UTF-8  

''''' 

Created on 2015年5月19日 

 

@author: zhaozhiyong 

'''  

  

from numpy import *  

from function import *  

  

def dfp(fun, gfun, x0):  

    result = []  

    maxk = 500  

    rho = 0.55  

    sigma = 0.4  

    m = shape(x0)[0]  

    Hk = eye(m)  

    k = 0  

    while (k < maxk):  

        gk = mat(gfun(x0))#计算梯度  

        dk = -mat(Hk)*gk  

        m = 0  

        mk = 0  

        while (m < 20):  

            newf = fun(x0 + rho ** m * dk)  

            oldf = fun(x0)  

            if (newf < oldf + sigma * (rho ** m) * (gk.T * dk)[0,0]):  

                mk = m  

                break  

            m = m + 1  

          

        #DFP校正  

        x = x0 + rho ** mk * dk  

        sk = x - x0  

        yk = gfun(x) - gk  

        if (sk.T * yk > 0):  

            Hk = Hk - (Hk * yk * yk.T * Hk) / (yk.T * Hk * yk) + (sk * sk.T) / (sk.T * yk)  

          

        k = k + 1  

        x0 = x  

        result.append(fun(x0))  

      

    return result  

testDFP.py



[python] view
 plain copy


 


#coding:UTF-8  

''''' 

Created on 2015年5月19日 

 

@author: zhaozhiyong 

'''  

  

from bfgs import *  

from dfp import dfp  

  

import matplotlib.pyplot as plt    

  

x0 = mat([[-1.2], [1]])  

result = dfp(fun, gfun, x0)  

  

n = len(result)  

ax = plt.figure().add_subplot(111)  

x = arange(0, n, 1)  

y = result  

ax.plot(x,y)  

  

plt.show()