欧拉–拉格朗日系统在权重不平衡有向图下的分布式优化算法

Distributed Optimization Algorithm for Euler-Lagrange Systems over Weight-Unbalanced Digraph

摘要: 本文研究了欧拉–拉格朗日(EL)多智能体系统在权重不平衡有向图下的分布式优化问题，优化目标为通过智能体间的局部通讯最小化全局目标函数，该目标函数为智能体自身局部目标函数的和。为解决该问题，本文设计平衡补偿变量调节拓扑权重，并提出EL系统在权重不平衡有向图下的分布式优化算法，该算法使智能体状态达成一致的同时，协同最小化全局目标函数。最后，给出一个基于Simulink的数值仿真验证所提出算法的有效性。

Abstract: This paper investigates the distributed optimization problem of Euler-Lagrange (EL) multi-agent systems in weighted unbalanced digraphs, the optimization objective of this paper is to minimize the global cost function through local communications among agents, where the global cost function is summed up by local ones assigned to corresponding agent. To address this problem, the balanced compensation variables are designed to adjust the topology weights, then, a distributed optimization algorithm for EL systems over weight-unbalanced digraphs is proposed, which enables agents to achieve consensus while cooperatively minimizing the global cost function. Finally, numerical simulations based on Simulink are provided to verify the effectiveness of the proposed algorithm.

1. 引言

近年来，随着移动通信技术和人工智能等相关领域的不断发展，多智能体系统协同控制的实用性日益凸显。得益于其在智能电网[1]、大规模机器学习[2]、医学研究[3]等各领域的广泛应用，分布式优化问题逐步成为多智能体系统分布式协同控制领域研究的热点之一。在该问题背景下，系统中的每个智能体拥有一个仅自身可知的局部目标函数，系统的全局目标函数为所有局部目标函数的和，分布式优化的目标是每个智能体通过与邻居智能体间的信息交互协同最小化全局目标函数[4]。不同于集中式，分布式的方法意味着每个智能体通过局部信息交互与迭代达成一致性，克服了集中式算法在复杂问题上的局限。文献[4]-[11]提出了多种算法以分析求解该分布式优化问题。文献[4]提出了一种针对离散时间多智能体系统的分布式子梯度方法，通过提出投影共识算法将文献[5]中的工作扩展到解决代理被限制在给定集合内的问题。文献[6]提出了考虑离散时间条件、局部约束集和通信延迟的分布式次梯度投影算法。文献[7]针对非均匀梯度增益、有限时间收敛性和凸约束集，提出了一类分布式连续时间算法来处理连续时间分布优化问题。为了扩大适用性，文献[8]-[11]提出了一些基于拉格朗日的分布式算法。此外，结合物理动力学，分布式优化问题得到了越来越广泛的重视，当控制方案最优时，每个智能体的行为也是最优的。文献[12]为一般线性多智能体系统建立了两种自适应分布式最优算法，以最小化性能函数实现一致性。

当解决复杂的非线性质点系系统的动力学问题时，EL系统经常被使用，故其具有较强的通用性和显著的研究价值[13] [14]。同时，多EL系统的分布式优化问题也引起了各位研究者的关注。文献[15] [16]针对多EL系统建立了几种分布式优化算法，但其控制参数的选择与全局信息有关，导致了这些设计的局限性。为了满足实际场景的需要，文献[3]提出了一种针对无向图的不确定EL系统的自适应分布式优化算法，实现了隐私保护和完全分布式的需要。文献[17]研究了不确定多EL系统的约束分布式优化问题。文献[18]设计了一种具有局部和相对测量值的分布式非线性控制器，以实现多EL系统的分布式优化结果。

目前为止的研究结果大多是在通信拓扑被建模为无向图或加权平衡有向图[15]-[17]的条件下提出的。但双向通信是一种十分理想的信息交互方式，每个智能体都需要同时获取入邻信息和出邻信息，在实际应用中往往难以实现。因此，在权重不平衡的有向网络上多EL系统的分布式优化问题仍有待解决。为解决该问题，本文提出了一种EL系统在权重不平衡有向图下的分布式优化算法。与文献[15]-[18]不同，本文考虑的智能体间的通讯结构为更一般的非平衡有向拓扑结构，系统的总入度和总出度无需相等，拓扑的不对称性增加了理论分析的复杂度；此外，文献[15]-[17]以解耦的方式解决分布式优化问题，实际应用中难以实现，为解决该问题，本文的算法采用更为复杂的耦合结构，增加了所提出算法的实用性。

基于上述论述，本文的主要结构如下：第一章给出符号、图论和凸函数等相关预备知识；第二章给出问题的表述；第三章给出EL系统在权重不平衡有向图下的分布式优化算法；第四章利用数值仿真验证所提出算法的有效性；第五章给出最后的结论。

2. 预备知识

2.1. 符号

在本文中，$ℝ$ 表示一个全体实数集，${ℝ}^n$ 表示一个n维的实向量空间，${ℝ}^{+}$ 表示正实数空间。${ℝ}^{m\times n}$ 表示$m\times n$ 维实矩阵组成的集合。$\Vert ·\Vert$ 表示向量的欧几里得范数和矩阵的诱导2-范数；$1_n$ 表示元素为1的n维列向量；$I_m$ 表示m维的单位阵；乘子$\otimes$ 表示Kronecker积。

2.2. 图论

$\mathcal{G}=\left\{\mathcal{N},ℰ,\mathcal{A}\right\}$ 表示智能体之间互相传递信息的加权有向网络，其中，$\mathcal{N}=\left\{1,\cdot \cdot \cdot ,n\right\}$ 表示智能体的集合；$ℰ\subseteq \mathcal{N}\times \mathcal{N}$ 表示有向边的集合；$\mathcal{A}={\left[a_{ij}\right]}_{{}_{n\times n}}$ 表示加权的邻接矩阵。$\left(i,j\right)$ 表示从i到j的有向边，其表示智能体i向智能体j发送信息，但反之不然。$a_{ij}$ 表示从j到i的有向边的权重。若$\left(j,i\right)\in ℰ$ ，则$a_{ij}=1$ ，否则$a_{ij}=0$ 。网络中的有向路径是由有向边连接的智能体序列。若任意智能体之间均有一条有向路径，则称$\mathcal{G}$ 是强连通的有向网络。

定义${\mathcal{N}}_{in}^i=\left\{j\in \mathcal{N}|a_{ij}>0\right\}$ 表示智能体i的入邻集。相应的，$d_i^{in}={\displaystyle {\sum }_{j\in {\mathcal{N}}_i}{a}_{ij}}$ 表示智能体i的入度。有向网络$\mathcal{G}=\left\{\mathcal{N},ℰ\right\}$ 的拉普拉斯矩阵定义为$L={\left[l_{ij}\right]}_{n\times n}$ ，其中，$l_{ii}={\displaystyle {\sum }_{j\ne i}{a}_{ij}}$ ；$l_{ij}=-a_{ij}$ ，$j\ne i$ 。若网络中有${\displaystyle {\sum }_{j=1}^n{a}_{ij}}\ne {\displaystyle {\sum }_{j=1}^n{a}_{ji}}$ 成立，则称$\mathcal{G}$ 是权重不平衡的有向网络。

引理1 [10]：权重不平衡的强连通有向网络$\mathcal{G}=\left\{\mathcal{N},ℰ,\mathcal{A}\right\}$ 满足如下性质：

对于拉普拉斯矩阵L的0特征值，$\xi ={\left({\xi }_1,\cdots ,{\xi }_n\right)}^{{\rm T}}$ ，${\xi }_i\ge 0$ ，$i=1,\cdots ,n$ 是其对应的左特征向量，并且$\xi$ 满足${\xi }^{{\rm T}}L=0_n^{{\rm T}}$ 及${\displaystyle {\sum }_{i=1}^n{\xi }_i}=1$ ；

拉普拉斯矩阵L的0特征值对应的右特征向量为$1_n$ ；

定义$\Xi =\text{diag}\left({\xi }_1,\cdots ,{\xi }_n\right)$ ，则$\Xi L$ 是对应于权重平衡的强连通有向网络的拉普拉斯矩阵。令$\widehat{L}=\frac{\Xi L+\Xi L^{{\rm T}}}{2}$ ，则有$x^{{\rm T}}\widehat{L}x=2x^{{\rm T}}Lx$ 。

2.3. 凸函数

对一个可微函数$f\left(x\right)$ ，$\nabla f\left(x\right)$ 表示其梯度。若对于$\forall 0<a<1$ ，$\forall x,y\in R^n$ ，有$f\left(ax+\left(1-a\right)y\right)\le af\left(x\right)+\left(1-a\right)f\left(y\right)$ 成立，则称函数$f\left(x\right)$ 是凸函数。若存在常数$\omega >0$ ，对于$\forall x,y\in R^n$ ，有${\left(x-y\right)}^{{\rm T}}\left(\nabla f\left(x\right)-\nabla f\left(y\right)\right)\ge \omega {\Vert x-y\Vert }^2$ 成立，则称函数$f\left(x\right)$ 是$\omega$ -强凸的。若存在常数$\theta >0$ ，对于$\forall x,y\in R^n$ ，$\Vert f\left(x\right)-f\left(y\right)\Vert \le \theta \Vert x-y\Vert$ ，则称函数$f\left(x\right)$ 是$\theta$ -Lipschitz的。若梯度$\nabla f\left(x\right)$ 在一个连通的开集$\Omega$ 上是局部Lipschitz的，则其在$x\in \Omega$ 处的Hessian广义矩阵为${\partial }_x^{2}f\left(x\right)={\nabla }^{2}f\left(x\right)$ 。

以下的引理给出判断分布式优化问题是否达到最优的依据。

引理2 [19]：若$f\left(x\right):{ℝ}^m\to ℝ$ 是一个连续可微的凸函数，则$f\left(x\right)$ 在$x^{\ast }$ 取最小值当且仅当$\nabla f\left(x^{\ast }\right)=0$ 。

3. 问题描述

考虑有向网络$\mathcal{G}=\left\{\mathcal{N},ℰ\right\}$ 包含n个智能体。任一智能体$i\in \mathcal{N}$ 的动力学建模成如下欧拉–拉格朗日系统：

$M_i\left(q_i\right){\ddot{q}}_i+C_i\left(q_i,{\dot{q}}_i\right){\dot{q}}_i={\tau }_i$ (1)

其中，$q_i\in R^m$ 表示广义坐标；${\dot{q}}_i\in R^m$ 表示广义坐标导数；$M_i\left(q_i\right)\in R^{m\times m}$ 、$C_i\left(q_i,{\dot{q}}_i\right){\dot{q}}_i$ 、${\tau }_i$ 分别表示对称的惯性矩阵、科里奥利力向量和广义控制力向量。EL系统被广泛地用于机械臂系统、无人车、航天器等实际对象的建模。研究EL系统的分布式最优协同控制算法具有普遍的实际意义。EL系统满足如下经典性质：

${\dot{M}}_i\left(q_i\right)-2C_i\left(q_i,{\dot{q}}_i\right)$ 是反对称矩阵；

存在常数$\underset{\_}{M}$ ，$\overline{M}$ ，$k_c>0$ 使得对于$\forall q_i,{\dot{q}}_i$ ，有$\underset{\_}{M}{I}_m<M_i\left(q_i\right)<\overline{M}{I}_m$ 和$C_i\left(q_i,{\dot{q}}_i\right)\le k_c\Vert {\dot{q}}_i\Vert$ 成立。

对于动力学建模为(1)的智能体，每个智能体拥有一个仅自己可得的局部目标函数$f_i\left(q_i\right):{ℝ}^m\to ℝ$ 。对于非平衡有向网络中的多EL系统(1)，分布式优化问题的目标为最小化如下所示的全局目标函数：

$f\left(q\right)={\displaystyle {\sum }_{i=1}^n{f}_i}\left(q_i\right),\text{ s}\text{.t}\text{. }{q}_i=q_j,\forall i,j\in \mathcal{N}$ (2)

为了解决上述问题我们作出如下假设：

假设1：非平衡有向图$\mathcal{G}=\left\{\mathcal{N},ℰ\right\}$ 是强连通图。

假设2：任一局部目标函数$f_i\left(q_i\right)$ 是可微且$\omega$ -强凸的，其梯度$\nabla f_i\left(q_i\right)$ 是满足全局Lipshitz性，Lipschitz系数为$\theta$ 。

注释1：上述两个假设为解决分布式优化问题(2)的经典假设，且由假设2可知，全局目标函数$f\left(q\right)={\displaystyle {\sum }_{i=1}^n{f}_i\left(q_i\right)}$ 是强凸的，这保证了分布式优化问题(2)解的唯一性。

4. 主要内容

为解决多EL系统(1)在权重非平衡有向图下的分布式优化问题(1)，本文首先明确需为智能体设计恰当的广义控制力${\tau }_i$ ，使得每个智能体收敛到分布式优化问题(1)的全局最优解，并达成一致性。为此，对于智能体$\forall i\in \mathcal{N}$ ，本文给出如下所示的广义控制力：

${\tau }_i=-\left(c_1+c_2{\left|{\dot{q}}_i\right|}^2\right)\left({\dot{q}}_i-q_i^d\right)$ (3)

其中$c_1$ 、$c_2$ 分别为取值范围待定的正参数，$q_i^d$ 为中间变量，用于估计全局最优解。进一步，对$\forall i\in \mathcal{N}$ ，分布式优化算法及平衡补偿变量如下所示：

$q_i^d=-c_3{\left(z_i^i\right)}^{-1}\nabla f_i\left(q_i\right)-c_4{\displaystyle {\sum }_{j\in {\mathcal{N}}_i}{a}_{ij}\left(q_i-q_j\right)}-\mu {\zeta }_i$ (4)

${\dot{\zeta }}_i=c_3{c}_4{\displaystyle {\sum }_{j\in {\mathcal{N}}_i}{a}_{ij}}\left(q_i-q_j\right)$ (5)

${\dot{z}}_i=-{\displaystyle {\sum }_{j\in {\mathcal{N}}_i}{a}_{ij}\left(z_i-z_j\right)}$ (6)

其中$c_3$ 、$c_4$ 、$\mu$ 分别为取值范围待定的正参数，${\zeta }_i$ 为智能体i的内部状态，$z_i$ 是第i个平衡补偿器，$z_i^j$ 是$z_i$ 的第j个分量，$j=1,\cdot \cdot \cdot ,n$ 。此外，$z_i$ 的初始状态需满足如下条件：

$\begin{array}{l}{z}_i^i\left(0\right)=1,\forall i\in {\mathcal{N}}_i\\ z_i^j\left(0\right)=0,\forall i\ne j\end{array}$ (7)

注释2：分布式优化问题(2)的全局最优解由引入的中间变量$q_i^d,i=1,\cdots ,n$ 估计。分布式优化算法(4)~(6)中，每个智能体在梯度$\nabla f_i\left(q_i\right)$ 驱使下最小化其局部目标函数$f_i\left(q_i\right)$ ，一致性项${\displaystyle {\sum }_{j\in {\mathcal{N}}_i}{a}_{ij}\left(q_i-q_j\right)}$ 保证不同智能体间的状态$q_i,i=1,\cdots ,n$ 达成一致，${\zeta }_i$ 平衡梯度项和一致性项的关系，保证智能体最终达成的一致性状态收敛到(2)的全局最优解上。同时，为处理非平衡有向图对分布式优化算法的影响，本文引入了平衡补偿变量$z_i$ ，对任意$i=1,\cdots ,n$ ，$z_i$ 的第i个元素将收敛到拉普拉斯矩阵L的0特征值对应左特征向量$\xi$ 的第i个分量上。

注释3：在广义控制力(3)的作用下，每个智能体的广义速度${\dot{q}}_i$ 将追踪上相对应的中间变量$q_i^d$ 。换言之，多EL系统(1)的分布式优化问题(2)等价于一个跟踪问题。此外，由分布式优化算法(4)~(6)可知，智能体利用自身真实信息估计全局最优解，即对全局最优解的估计与追踪同步进行，这样的效果使算法更具实际意义。

基于上述论述，本文主要定理如下：

定理1：当假设1、2、3成立时，给定初值$\zeta \left(0\right)=0$ ，分布式优化算法(3)~(6)会驱使多EL系统(1)的广义坐标渐近收敛到分布式优化问题(2)的全局最优解上。

证明：证明共分为四步，第一步证明平衡点最优，第二步给出中间变量$q_i^d$ 与分布式优化问题(2)的全局最优解的收敛性分析，第三步给出EL系统(1)的广义坐标导数与$q_i^d$ 的收敛性分析，第二步与第三步耦合，因此在第四步中给出整个系统的稳定性条件。首先给出分布式优化算法(4)~(6)平衡点与非平衡有向图下的分布式优化问题(2)全局最优解之间的关系。

第一步：由文献[10]可知，拉普拉斯矩阵L的0特征值对应的左特征向量$\xi$ 可以被平衡补偿器(6)估计，换言之，${\lim }_{n\to \infty }{Z}_n=\Xi$ ，$\Xi =\text{diag}\left({\xi }_1,\cdots ,{\xi }_n\right)$ 。

令$e=\dot{q}-q^d$ ，可将算法改写为：

$\dot{q}=-c_3{\left(Z_n\otimes I_m\right)}^{-1}\nabla f\left(q\right)-c_4\left(L\otimes I_m\right)q-\mu \zeta +e$ (8)

$\dot{\zeta }=c_3{c}_4\left(L\otimes I_m\right)q$ (9)

$\dot{z}=-\left(L\otimes I_m\right)z$ (10)

对(9)两边左乘${\xi }^{{\rm T}}\otimes I_m$ ，可得$\left({\xi }^{{\rm T}}\otimes I_m\right)\dot{\zeta }=0$ ，即${\displaystyle {\sum }_{i=1}^n{\xi }_i{\zeta }_i\left(t\right)=}{\displaystyle {\sum }_{i=1}^n{\xi }_i{\zeta }_i\left(0\right)=0_n}$ 。设$\left(\overline{q},\dot{\overline{q}},\overline{\zeta }\right)$ 为稳定点，系统稳定时$Z_n=\Xi$ ，故：

$0=\dot{\overline{q}}$ (11)

$0=-c_3{\left(\Xi \otimes I_m\right)}^{-1}\nabla f\left(\overline{q}\right)-c_4\left(L\otimes I_m\right)\overline{q}-\mu \overline{\zeta }+e$ (12)

$0=c_3{c}_4\left(L\otimes I_m\right)\overline{q}$ (13)

对(12)两边左乘$\Xi \otimes I_m$ ，可得：

$0=-c_3\nabla f\left(\overline{q}\right)-c_4\left[\left(\Xi L\right)\otimes I_m\right]\overline{q}-\left(\Xi \otimes I_m\right)\mu \overline{\zeta }+\left(\Xi \otimes I_m\right)e$ (14)

对(14)两边同乘$1_n^{{\rm T}}\otimes I_m$ ，可得$0=-c_3{\displaystyle {\sum }_{i=1}^n\nabla f_i\left({\overline{q}}_i\right)}+\left(1_n^{{\rm T}}\Xi \right)\otimes e$ ，由于整个系统稳定时$e=0$ ，故${\displaystyle {\sum }_{i=1}^n\nabla f_i\left({\overline{q}}_i\right)}=0$ ，又$\mathcal{G}$ 为强连通图，故${\overline{q}}_i={\overline{q}}_j$ ，即分布式优化算法(4)~(6)的稳定点是问题(2)的全局最优解。

下一步，基于李雅普诺夫理论，给出分布式优化算法(4)~(6)对全局最优解的收敛性分析。

第二步：将算法(4)~(6)改写为如下的紧凑形式：

$\begin{array}{l}\dot{q}=-c_3{\left(\Xi \otimes I_m\right)}^{-1}\nabla f\left(q\right)-c_4\left(L\otimes I_m\right)q-\mu \zeta +e-c_3{\left(\left(Z_n-\Xi \right)\otimes I_m\right)}^{-1}\nabla f\left(q\right)\\ \dot{\zeta }=c_3{c}_4\left(L\otimes I_m\right)q\\ \dot{z}=-\left(L\otimes I_m\right)z\end{array}$ (15)

由于${\lim }_{n\to \infty }{Z}_n^{-1}={\Xi }^{-1}$ ，证明(15)的稳定性等价于分析如下系统的收敛性：

$\begin{array}{l}\dot{q}=-c_3{\left(\Xi \otimes I_m\right)}^{-1}\nabla f\left(q\right)-c_4\left(L\otimes I_m\right)q-\mu \zeta +e\\ \dot{\zeta }=c_3{c}_4\left(L\otimes I_m\right)q\\ \dot{z}=-\left(L\otimes I_m\right)z\end{array}$ (16)

令$X=q-\overline{q}$ ，$Y=\zeta -\overline{\zeta }$ ，可得：

$\begin{array}{l}\dot{X}=-c_3{\left(\Xi \otimes I_m\right)}^{-1}\varphi \left(X\right)-c_4\left(L\otimes I_m\right)X-\mu Y+e\\ \dot{Y}=c_3{c}_4\left(L\otimes I_m\right)X\end{array}$ (17)

其中，$\varphi \left(X\right)=\nabla f\left(q\right)-\nabla f\left(\overline{q}\right)$ 。定义如下李雅普诺夫函数：

$V_1=\frac{c_3}{2}{X}^{{\rm T}}\left(\Xi \otimes I_m\right)X+\frac{1}{2c_3}{\left(c_{3}X+Y\right)}^{{\rm T}}\left(\Xi \otimes I_m\right)\left(c_{3}X+Y\right)$ (18)

根据(18)可以得出$V_1$ 的导数为：

$\begin{array}{c}{\dot{V}}_1=c_3{X}^{{\rm T}}\left(\Xi \otimes I_m\right)\dot{X}+\frac{1}{c_3}{\left(c_{3}X+Y\right)}^{{\rm T}}\left(\Xi \otimes I_m\right)\left(c_3\dot{X}+\dot{Y}\right)\\ =-2c_3^2{X}^{{\rm T}}\varphi \left(X\right)-c_3{c}_4{X}^{{\rm T}}\left(\Xi L\otimes I_m\right)X-2c_3\mu X^{{\rm T}}\left(\Xi \otimes I_m\right)Y\\ -c_3{Y}^{{\rm T}}\varphi \left(X\right)-\mu Y^{{\rm T}}\left(\Xi \otimes I_m\right)Y+2c_3{X}^{{\rm T}}\left(\Xi \otimes I_m\right)e+Y^{{\rm T}}\left(\Xi \otimes I_m\right)e\end{array}$(19)

由假设2和杨氏不等式可得：

$\begin{array}{c}{\dot{V}}_1\le -2c_3^2\omega {\Vert X\Vert }^2-c_3{c}_4{\lambda }_2\left(\widehat{L}\right){\Vert X\Vert }^2+c_3^2{\mu }^2{\xi }_{\max }^2{\epsilon }_1{\Vert X\Vert }^2+c_3^2{\epsilon }_1{\Vert X\Vert }^2+\frac{c_3^2{\theta }^2{\epsilon }_1}{4}{\Vert X\Vert }^2\\ +\frac{1}{{\epsilon }_1}{\Vert Y\Vert }^2+\frac{1}{{\epsilon }_1}{\Vert Y\Vert }^2-\mu {\xi }_{\min }{\Vert Y\Vert }^2+2{\epsilon }_1{\Vert Y\Vert }^2+\frac{{\xi }_{\max }^2}{{\epsilon }_1}{\Vert e\Vert }^2+\frac{{\xi }_{\max }^2}{8{\epsilon }_1}{\Vert e\Vert }^2\\ =-\left[2c_3^2\omega -c_3^2{\epsilon }_1\left({\mu }^2{\xi }_{\max }^2+\frac{{\theta }^2}{4}+1\right)+c_3{c}_4{\lambda }_2\left(\widehat{L}\right)\right]{\Vert X\Vert }^2\\ -\left[\mu {\xi }_{\min }-2\left(\frac{{\epsilon }_1+1}{{\epsilon }_1}\right)\right]{\Vert Y\Vert }^2+\frac{9{\xi }_{\max }^2}{8{\epsilon }_1}{\Vert e\Vert }^2\end{array}$(20)

其中$\widehat{L}=\left(\Xi L+L^{{\rm T}}\Xi \right)/2$ ，将(20)简写为：

${\dot{V}}_1\le -b_1{\Vert X\Vert }^2-b_2{\Vert Y\Vert }^2+b_3{\Vert e\Vert }^2$ (21)

其中$b_1=2c_3^2\omega -c_3^2{\epsilon }_1\left({\mu }^2{\xi }_{\max }^2+{\theta }^2/4+1\right)+c_3{c}_4{\lambda }_2\left(\widehat{L}\right)$ ，$b_2=\mu {\xi }_{\min }-2\left({\epsilon }_1+1\right)/{\epsilon }_1$ ，$b_3=9{\xi }_{\max }^2/8{\epsilon }_1$ 。

下一步，基于李雅普诺夫理论，给出智能体广义坐标导数与$q_i^d$ 间误差的稳定性分析。

第三步：定义如下李雅普诺夫函数：

$V_2=\frac{1}{2}{e}^{{\rm T}}Me$ (22)

根据(22)可以得出$V_2$ 的导数为：

$\begin{array}{c}{\dot{V}}_2=\frac{1}{2}{e}^{{\rm T}}\dot{M}e+e^{{\rm T}}M\dot{e}\\ =e^{{\rm T}}\tau -e^{{\rm T}}\left(M{\dot{q}}^d-C{q}^d\right)\end{array}$ (23)

由杨氏不等式可得：

${\dot{V}}_2\le e^{{\rm T}}\tau +\frac{1}{2{\epsilon }_2}\left({\Vert {\dot{q}}^d\Vert }^2+{\Vert q^d\Vert }^2\right)+\frac{{\epsilon }_2}{2}\left({\Vert e^{{\rm T}}M\Vert }^2+{\Vert e^{{\rm T}}C\Vert }^2\right)$ (24)

由(3)式可以得出：

$e^{{\rm T}}\tau \le \left(-c_1-c_2\times \underset{i}{\max }{\left|{\dot{q}}_i\right|}^2\right){\Vert e\Vert }^2$ (25)

此外：

$q^d=-c_3{\left(\Xi \otimes I_m\right)}^{-1}\varphi \left(X\right)-c_4\left(L\otimes I_m\right)X-\mu Y$ (26)

故由(26)可得：

$\begin{array}{c}{\Vert q^d\Vert }^2={\Vert -c_3{\left(\Xi \otimes I_m\right)}^{-1}\varphi \left(X\right)-c_4\left(L\otimes I_m\right)X-\mu Y\Vert }^2\\ \le k{\Vert X\Vert }^2+{\mu }^2{\Vert Y\Vert }^2\end{array}$ (27)

其中$k=c_3^2{\xi }_{\min }^{-2}{\theta }^2+c_4{\Vert L\Vert }^2$ ，此外：

$\begin{array}{c}{\dot{q}}^d=-c_3{\left(\Xi \otimes I_m\right)}^{-1}{\partial }_X\varphi \left(X\right)\dot{X}-c_4\left(L\otimes I_m\right)\dot{X}-\mu \dot{Y}\\ \le -c_3{\left(\Xi \otimes I_m\right)}^{-1}\theta \dot{X}-c_4\left(L\otimes I_m\right)\dot{X}-\mu \dot{Y}\end{array}$ (28)

故由(28)可得：

$\begin{array}{c}{\Vert {\dot{q}}^d\Vert }^2\le {\Vert -c_3{\left(\Xi \otimes I_m\right)}^{-1}\theta \dot{X}-c_4\left(L\otimes I_m\right)\dot{X}-\mu \dot{Y}\Vert }^2\\ \le k{\Vert \dot{X}\Vert }^2+{\mu }^2{\Vert \dot{Y}\Vert }^2\end{array}$ (29)

由(17)可得：

$\begin{array}{l}{\Vert \dot{X}\Vert }^2\le k{\Vert X\Vert }^2+\mu {\Vert Y\Vert }^2+{\Vert e\Vert }^2\\ {\Vert \dot{Y}\Vert }^2\le c_3^2{c}_4^2{\Vert L\Vert }^2{\Vert X\Vert }^2\end{array}$ (30)

故由(29)~(30)可得：

${\Vert {\dot{q}}^d\Vert }^2\le \left(k^2+c_3^2{c}_4^2{\mu }^2{\Vert L\Vert }^2\right){\Vert X\Vert }^2+{\mu }^{2}k{\Vert Y\Vert }^2+k{\Vert e\Vert }^2$ (31)

故由(27)、(31)可得：

$\frac{1}{2{\epsilon }_2}\left({\Vert {\dot{q}}^d\Vert }^2+{\Vert q^d\Vert }^2\right)\le \frac{1}{2{\epsilon }_2}\left(k^2+k+c_3^2{c}_4^2{\mu }^2{\Vert L\Vert }^2\right){\Vert X\Vert }^2+\frac{{\mu }^2}{2{\epsilon }_2}\left(k+1\right){\Vert Y\Vert }^2+\frac{k}{2{\epsilon }_2}{\Vert e\Vert }^2$ (32)

由EL系统的性质可得：

$\frac{{\epsilon }_2}{2}\left({\Vert e^{{\rm T}}M\Vert }^2+{\Vert e^{{\rm T}}C\Vert }^2\right)\le \frac{{\epsilon }_2{\overline{M}}^2}{2}{\Vert e\Vert }^2+\frac{{\epsilon }_2{k}_c^2}{2}\times \max {\Vert {\dot{q}}_i\Vert }^2{\Vert e\Vert }^2$ (33)

故由(25)、(32) (33)可得：

$\begin{array}{c}{\dot{V}}_2\le \frac{1}{2{\epsilon }_2}\left(k^2+k+c_3^2{c}_4^2{\mu }^2{\Vert L\Vert }^2\right){\Vert X\Vert }^2+\frac{{\mu }^2}{2{\epsilon }_2}\left(k+1\right){\Vert Y\Vert }^2\\ +\left[\frac{k}{2{\epsilon }_2}+\frac{{\epsilon }_2{\overline{M}}^2}{2}-c_1-\left(c_2-\frac{{\epsilon }_2{k}_c^2}{2}\right)\times \max {\Vert {\dot{q}}_i\Vert }^2\right]{\Vert e\Vert }^2\end{array}$(34)

由于分布式优化算法(4)~(6)与智能体的动力学(1)耦合，因此最后一步给出整个系统的收敛性条件。

第四步：定义如下李雅普诺夫函数：

$V=\sigma V_1+V_2$ (35)

其中$\sigma$ 为取值待定的正实数。根据(21)、(34)可知V的导数为：

$\begin{array}{c}\dot{V}=\sigma {\dot{V}}_1+{\dot{V}}_2\\ \le -\left[b_1\sigma -\frac{1}{2{\epsilon }_2}\left(k^2+k+c_3^2{c}_4^2{\mu }^2{\Vert L\Vert }^2\right)\right]{\Vert X\Vert }^2\\ -\left[b_2\sigma -\frac{{\mu }^2}{2{\epsilon }_2}\left(k+1\right)\right]{\Vert Y\Vert }^2\\ +\left[-c_1+b_3\sigma +\frac{k}{2{\epsilon }_2}+\frac{{\epsilon }_2{\overline{M}}^2}{2}\right]{\Vert e\Vert }^2\end{array}$ (36)

对于式(36)给定$c_3>0$ ，选取合适的$\mu >0$ 使得$b_2>0$ ；通过选取合适的$c_4>0$ 和${\epsilon }_1>0$ 使得$b_1>0$ ；接

下来，根据(34)、(36)，对任意的${\epsilon }_2>0$ 、$\sigma >\max \left\{k^2+k+c_3^2{c}_4^2{\mu }^2{\Vert L\Vert }^2/2{\epsilon }_1{b}_1,{\mu }^2\left(k+1\right)/2{\epsilon }_2{b}_2\right\}$ 、

$c_1>b_3\sigma +k/2{\epsilon }_2+{\epsilon }_2{\overline{M}}^2/2$ 、$c_2>{\epsilon }_2{k}_c^2/2$ 可得$\dot{V}<0$ 。因此由Barbalat’s引理[20]可得$X,Y$ 和追踪误差e的有界性。进而根据LaSalle不变集定理[21]可得${\lim }_{t\to \infty }X=0$ ，${\lim }_{t\to \infty }Y=0$ ，${\lim }_{t\to \infty }e=0$ 。因此，智能体真实状态的导数$\dot{q}$ 追踪参考轨迹$q_i^d$ ，同时智能体的状态q渐近收敛到分布式优化问题(2)的全局最优解。

5. 数值仿真

本节利用MATLAB中的Simulink模块，搭建一个数值实验来验证EL系统在权重非平衡有向图下的分布式优化算法(3)~(6)的有效性。给出n个被建模成欧拉–拉格朗日系统(1)的智能体，这些智能体协同定位一个被m个参考信号$r^i\in {ℝ}^2$ ，$i=1,\cdots ,m$ 围绕的未知动态信号源[16]，其中每个智能体只能感知对应的参考信号。因此系统的全局目标函数选取为$f\left(q\right)={\displaystyle {\sum }_{i=1}^n{f}_i\left(q_i\right)}$ ，其中$f_i\left(q_i\right)={\Vert q_i-r_i\Vert }^2$ 这表明每个智能体都倾向于靠近各自探测能力所及的参考信号，并且该目标函数满足假设2。智能体最终达成的一致性点是$f\left(q\right)$ 的最优解，该最优解可以看作是对未知信号源的位置的估计。

在上述问题场景中，每一个智能体被建模为如下的EL系统：

$\left[\begin{array}{cc}{M}_{11,i}& M_{12,i}\\ M_{21,i}& M_{22,i}\end{array}\right]\left[\begin{array}{c}{\ddot{q}}_{ix}\\ {\ddot{q}}_{iy}\end{array}\right]+\left[\begin{array}{cc}{C}_{11,i}& C_{12,i}\\ C_{21,i}& C_{22,i}\end{array}\right]\left[\begin{array}{c}{\dot{q}}_{ix}\\ {\dot{q}}_{iy}\end{array}\right]=\left[\begin{array}{c}{\tau }_{ix}\\ {\tau }_{iy}\end{array}\right]$ (37)

其中$M_{11,i}={\theta }_{1i}+2{\theta }_{2i}\cos q_{iy}$ ，$M_{12,i}=M_{21,i}={\theta }_{3i}+{\theta }_{2i}\cos q_{ix}$ ，$M_{22,i}={\theta }_{3i}$ ，此外$C_{11,i}=-2{\theta }_{2i}\sin q_{iy}{\dot{q}}_{iy}$ ，$C_{12,i}=C_{21,i}=-{\theta }_{2i}\sin q_{ix}{\dot{q}}_{ix}$ ，$C_{22,i}=0$ ，上述元素中的参数分别选取为${\theta }_{1i}=1.301$ ，${\theta }_{2i}=0.056$ ，${\theta }_{3i}=0.296$ 。智能体之间的通讯拓扑结构如图1所示，该拓扑结构是非平衡有向图，并且是强连通的，由图1可知，该图拉普拉斯矩阵0特征值对应的左特征向量$\xi ={\left[0.4,0.2,0.2,0.1,0.1\right]}^{{\rm T}}$ 。

考虑$n=5$ 和$m=5$ 的情形。每个智能体的动力学建模为(37)式。假设参考信号的位置为$r_1={\left(-3,4\right)}^{{\rm T}}$ ，$r_2={\left(4,4\right)}^{{\rm T}}$ ，$r_3={\left(-4,-3\right)}^{{\rm T}}$ ，$r_4={\left(3,-5\right)}^{{\rm T}}$ ，$r_5={\left(2,-3\right)}^{{\rm T}}$ 。智能体的广义坐标初始化为$q_1\left(0\right)={\left(5,1\right)}^{{\rm T}}$ ，$q_2\left(0\right)={\left(2,3\right)}^{{\rm T}}$ ，$q_3\left(0\right)={\left(-5,3\right)}^{{\rm T}}$ ，$q_4\left(0\right)={\left(-4,-3\right)}^{{\rm T}}$ ，$q_5\left(0\right)={\left(4,-3\right)}^{{\rm T}}$ 。

选取$c_1=68$ ，$c_2=0.27$ ，$c_3=2$ ，$c_4=5$ 和$\mu =1$ 。实验结果如图2~4所示。从图2可以看出，在分布式优化算法(4)~(6)的驱动下，所有智能体的广义坐标收敛到分布式优化问题(2)的全局最优解上。图3表示全局目标函数梯度的轨迹，当$t\to \infty$ 时轨迹收敛到0，由引理2可知全局目标函数取到了最小值。图4表示由(6)式生成且初值满足(7)式的平衡补偿器$z_i\left(t\right)$ 收敛到拉普拉斯矩阵L的0特征值对应的左特征向量$\xi$ 上。

Figure 1. Communication topology among agents

图1. 智能体通讯拓扑

Figure 2. Trajectories of each agent’s state

图2. 智能体状态轨迹

Figure 3. Trajectories of the gradient of global cost function

图3. 全局目标函数梯度轨迹

Figure 4. Trajectories of balanced compensation variables

图4. 平衡补偿变量轨迹

6. 总结

本文研究了欧拉–拉格朗日系统在权重非平衡有向图下的分布式优化问题，通过代数图论中有向拉普拉斯矩阵的性质，设计平衡补偿器以调节拓扑权重，将问题转化为权重平衡图下的分布式优化问题，进而设计相应的分布式优化算法，在本文的理论研究基础上，未来对于多欧拉–拉格朗日系统的分布式优化问题可以继续进行深入的研究。但本文所研究的工作，均为理论研究阶段。虽然通过实验验证了算法的有效性，但是展示性不够强。能否通过构建任务模型进行实物验证，使其具有更强的展示性是之后研究中值得考虑的问题。

基金项目

本论文由国家自然科学基金(项目号12161063)支持。

NOTES

^*通讯作者。