随着控制技术的发展, 将控制技术与优化算法相结合, 以分布式方式求解网络优化问题成为近年来的研究热点^[1-7]. 作为网络优化领域里一个重要分支, 资源分配在网络系统中有很广泛的应用, 例如能源系统中的经济调度和通讯网络中的网络利用率最优化等问题.

近些年来, 基于离散时间算法的分布式优化问题已经得到了一些很好的结果. 例如, 文献[8-10]等针对有约束或无约束最优协调问题设计了一种基于梯度的算法. 而对于强连通有向网络下的资源分配问题, 文献[11-13]等基于多智能体一致性理论, 分别设计一种基于余量和基于Push-sum (推和)的分布式优化算法. 更多离散时间算法可参考文献[14-19]及它们的引用文献, 其中收敛速度最快的为文献[17]的线性收敛, 但是依然无法预知其收敛时间. 另一方面, 随着信息物理系统和连续时间控制技术的发展^[20-25], 利用基于连续时间的分布式控制策略也已被广泛应用于解决资源分配问题, 特别是其中的有限/固定/预定时间收敛理论, 更符合实际系统需求. 为了解决无向网络下的无约束资源分配问题, 基于多智能体一致性理论, 文献[26]设计了一种指数收敛速度二阶连续时间算法. 不同于文献[26]的渐近式收敛, 文献[27-30]分别设计一种基于固定时间和预定时间收敛的连续时间算法. 在此基础上, 文献[31-35]提出了一种考虑不等式约束的分布式算法, 该算法利用映射算子将节点状态限制在不等式约束集合内, 从而得到了一种无初值约束分布式算法, 并在文献[35]中实现了固定时间收敛. 与映射法不同, 文献[36-38]和文献[39]利用精确罚函数法将局部等式约束转移到目标函数中, 分别设计了一阶和二阶连续时间算法. 而对于有向网络下的资源分配问题, 也有很多学者提出了解决方法, 例如奇异摄动算法^[40]、原始对偶算法^[41]、映射法^[42]、二阶连续时间算法^[43]以及微分映射法^[44].

从前面的讨论可以看出, 目前只有一阶分布式优化算法在连续时间和离散时间情况下都取得了显著的效果. 而有关二阶多智能体分布式优化算法的文献不是很多. 其中文献[26, 45-46]研究了无约束资源分配问题, 并没有将节点局部不等式约束考虑进来. 但在许多工程应用中, 节点的决策会受到局部约束, 这些约束可能来自安全考虑、生产能力、执行能力等, 如电力系统中发电机的输出功率的边界限制^[1]. 文献[39]和文献[43]分别利用精确罚函数法和微分映射法解决了资源分配问题中的局部不等式约束. 但是上述两篇文献所提出的算法都对初始状态进行限定. 而这种初始化协调过程可能只适应于静态网络. 然而对于一个动态网络来说, 一旦网络发生变化或者节点生产约束发生突变时, 就必须重新分配资源配置更改. 因此这些优化算法重新启动时, 必须重新进行初始化协调, 这大大降低了它们的适用性. 此外, 文献[26, 39, 45-46]中考虑的都是无向网络.

然而, 由于通信延迟等造成的通信中断会使得无向网络变成有向网络, 进而可能会使得一些基于无向网络的分布式算法失效. 另外, 无向网络是有向平衡网络的一个特例, 在工程应用中, 有向平衡网络比无向网络也更易于实现, 因此研究有向平衡网络上的资源分配问题是有意义的. 由于二阶多智能体系统的广泛存在, 如: 移动机器人、发电机等, 本文考虑二阶多智能体的分布式资源分配问题. 针对有向平衡网络下的资源分配问题, 本文首先利用局部不等式及其梯度设计第一种算法, 其次利用固定时间理论和映射算子设计第二种算法, 其中映射算子用于处理不等式约束, 且确保此不等式在固定时间内即可满足, 此方法可推广到同类型带局部不等式约束的分布式优化问题中.

本文所提算法创新点如下:

1) 相比于文献[26, 45-46], 本文将二阶多智能体节点的局部不等式约束考虑进来;

2) 与文献[39]的无向网络相比, 我们研究有向平衡网络下的资源分配问题;

3) 相比于文献[43], 本文所提两种算法均对初始值无要求, 且分别实现渐近和指数收敛;

4) 针对节点局部不等式约束, 本文结合了固定时间收敛理论和映射算子, 使得收敛速度上会比单纯使用映射算子方法的算法要快, 并通过案例仿真可以验证此性质.

本文后续内容结构如下: 在第1节中, 将完成问题描述以及预备知识, 其中预备知识包括代数图论、固定时间收敛理论等; 在第2节中, 首先进行算法设计, 分别是渐近收敛算法和指数收敛算法, 其次给出相应的最优性引理以及各自算法的收敛性分析; 在第3节中, 给出3个仿真案例;本文的总结与展望放在第4节中.

1.   问题描述及预备知识

在后续描述中, 令$ \mathbf{R} $、$ \mathbf{R}_{+} $和$ \mathbf{R}^m $分别表示实数、非负实数和$ m $维实空间. 对于闭集合$\Omega\in\mathbf{R}^m,$ 定义映射算子$P_{\Omega}(\boldsymbol x) = \arg\min_{\forall \boldsymbol y\in\Omega}\{\Vert\boldsymbol x-\boldsymbol y\Vert\}$, 且有性质$(\boldsymbol x- $$ P_{\Omega}(\boldsymbol x))^{\text{T}}(\boldsymbol z-P_{\Omega}(\boldsymbol x))\le0,\forall \boldsymbol z\in\Omega,\boldsymbol x\in\mathbf{R}^m$. $ \text{sign}(\cdot) $表示符号函数.$ (x)^+ = \max\{x,0\} $. $ \otimes $表示克罗内克积, 对于矩阵$ \boldsymbol A,\boldsymbol B\in\mathbf{R}^{n \times n} $, $\boldsymbol A\otimes\boldsymbol B := [a_{ij}\boldsymbol B]\in $$ \mathbf{R}^{n^2\times n^2}$.$ \partial f(\boldsymbol x) $表示函数$ f(\boldsymbol x) $的梯度.

1.1   问题描述

本文所研究的二阶多智能体系统包含$ n $个智能体节点, 每个智能体$ i $都采用如下动力学模型:

其中, $ \boldsymbol x_i\in\mathbf{R}^m $表示第 $ i $个智能体的资源分配量, $\boldsymbol u_i\in $$ \mathbf{R}^m$表示节点$ i $的控制输入. 并且每个智能体节点$ i $都拥有一个决策变量$ \boldsymbol x_i $和私有代价函数$ f_i(\boldsymbol x_i):\mathbf{R}^m\to\mathbf{R} $, 其中$ \boldsymbol x_i $满足私有约束$ \boldsymbol x_i\in\Omega_i $. 智能体系统的目标就是在满足自身节点约束的同时, 如何以最小的代价通过邻居节点信息交互把有限的资源$ \sum_{i = 1}^n\boldsymbol d_i $分配给所有智能体节点, 且各节点只能利用自己的私有信息$\Omega_i$和$ f_i(\boldsymbol x_i) $, 其中资源量$ \boldsymbol d_i $为各节点的虚拟分配量或期望资源量. 上述多智能体系统的资源分配问题数学模型表示如下:

其中, $ f(\boldsymbol x) $为全局目标函数: 集合$\Omega_i$表示每个智能体的局部约束, 集合$ V $表示智能体(编号)集合. 此外, 局部约束也可描述为$\Omega_i = \{\boldsymbol x_i|\varphi_i(\boldsymbol x_i)\le{\bf 0}_{s_i}\}$, 其中$\varphi_i(\boldsymbol x_i) = $$ [\varphi_{i,1}(\boldsymbol x_i),\varphi_{i,2}(\boldsymbol x_i),\cdots,\varphi_{i,s_i}(\boldsymbol x_i)]^{\text{T}},\;s_i\in\mathbf{N}$表示节点$ i $不等式约束的个数, 令$ s = \sum_{i = 1}^ns_i $. 以智能电网中的经济调度问题为例, $ f_i(\boldsymbol x_i) $表示发电机组的发电成本函数, 变量$ \boldsymbol d_i $表示负荷侧/用户侧的用电需求, 此物理量可以是各发电公司利用神经网络等技术, 结合用户侧的历史数据做出的短期负荷预测.

假设 1. 目标函数$ f(\boldsymbol x) $是$\omega-$强凸函数, 二次可微, 且其梯度满足$ \theta- $李普希茨性质. 即存在$\omega > $$ 0,\theta > 0$, 使得$(\boldsymbol x-\boldsymbol y)^{\text{T}}(\partial f(\boldsymbol x)-\partial f(\boldsymbol y))\ge \omega\Vert\boldsymbol x-\boldsymbol y\Vert^2, \Vert\partial f(\boldsymbol x)- $$ \partial f(\boldsymbol y)\Vert\le\theta\Vert\boldsymbol x-\boldsymbol y\Vert.$

假设 2. 局部不等式约束中$ \varphi(\boldsymbol x) $为凸函数且存在一点$ \boldsymbol x $使得 $ \varphi(\boldsymbol x)\le0 $和$ \sum_{i = 1}^n\boldsymbol x_i = \sum_{i = 1}^n\boldsymbol d_i $ 且存在某一智能体$ i $ 满足$ \varphi_i(\boldsymbol x_i)<0 $, 其中 $\varphi(\boldsymbol x) = [(\varphi_i (\boldsymbol x_i))^{\text{T}},\cdots, $$ (\varphi_n(\boldsymbol x_n))^{\text{T}}]^{\text{T}}\in\mathbf{R}^s,\boldsymbol x = [\boldsymbol x_1^{\text{T}},\cdots,\boldsymbol x_n^{\text{T}}]^{\text{T}}\in\mathbf{R}^{mn} .$

注 1. 假设1中的强凸性保证了问题 (2)最优解的唯一性. 此假设已广泛应用于此类问题的算法设计中, 如离散时间算法^{[12-13, 17, 19]}和连续时间算法^{[26, 31-32, 34-35, 39, 43]}, 有助于作者分析所设计算法的收敛性. 假设2是约束优化问题中常用的Slater条件, 此条件确保了后续引理2的有效性.

引理 1^[39]. 若假设1和2成立, 则$ \boldsymbol x^*\in\mathbf{R}^{nm} $是问题 (2)的最优解, 当且仅当存在$ \boldsymbol y^*\in\mathbf{R}^m $, 使得$ (\boldsymbol x^*,\boldsymbol y^*) $ 满足

引理 2^[41]. (KKT) 若假设1和2成立, 则$\boldsymbol x^*\in $$ \mathbf{R}^{nm}$是问题 (2)的最优解, 当且仅当存在$\boldsymbol y^*\in $$ \mathbf{R}^m,\;\boldsymbol{\mu}^*\in\mathbf{R}^s$, 使得$ (\boldsymbol x^*,\boldsymbol y^*,\boldsymbol{\mu}^*) $满足

1.2   代数图论

针对式 (1)所示的多智能体系统, 智能体节点间的通信链路可以用图$ G = (V,E,\boldsymbol A) $表示, 其中$V = \{1,2,\cdots, $$ n\}$ 表示节点集合, $ E $表示节点间的通信链路集合, $ \boldsymbol A $表示各节点间的通信权重. 若存在有向对$ (i,j)\in E $, 则表示节点$ i $能够接收到来自节点$ j $的信息, 此时就有$ a_{ij}>0 $. 此外, 本文中不考虑节点存在自环, 即$ a_{ii} = 0 $. 若网络$ G $中存在一组有向对$(i,i_1),(i_1,i_2),\cdots,(i_k,j)$, 则称节点$ j $到节点$ i $是可达的. 若网络$ G $任意一点到其他节点都是可达的, 则称网络$ G $是强连通的. 定义网络$ G $的拉普拉斯矩阵为$ \boldsymbol L = [l_{ij}]\in\mathbf{R}^{n\times n} $, 其中$l_{ii} = \sum_{j = 1}^na_{ij}, \; l_{ij} = -a_{ij}, $$ \forall j\not = i$. 所以必有结论$ \boldsymbol {L1}_n = {\bf 0}_n $. 对于有向网络$ G $, 若满足$ \boldsymbol 1_n^{\text{T}}\boldsymbol L = {\bf 0}_n^{\text{T}} $, 则称$ G $为平衡网络.

假设 3. 本文假设通信网络$ G $是有向平衡网络且强连通.

引理 3^[42]. 若假设3成立, 则有

其中$\hat{\boldsymbol L} = \dfrac{{\boldsymbol L}+{\boldsymbol L}^{\text{T}}}{2}$.

引理 4^[47]. 考虑如下微分方程

其中, $ \boldsymbol x\in\mathbf{R}^{n},f(\boldsymbol x):\mathbf{R}^n\to\mathbf{R}^n $在实数域 $ \mathbf{R}^n $上连续, 且满足$ f({\bf 0}_n) = {\bf 0}_n $. 若存在一连续正定函数$V(\boldsymbol x): \mathbf{R}^n\to $$ \mathbf{R}$和正实数$ \alpha,\beta>0,0<\mu<1,\nu>1 $使得不等式$\dot{V}\le $$ -\alpha V^{\mu}-\beta V^{\nu}$. 则存在一正实数 $T > 0,\;$使得对任意$ t\ge T $, 有$ \boldsymbol x(t) = {\bf 0}_n $, 即系统 (6)在固定时间$ T $收敛到源点, 且固定时间$ T $满足$T\le\dfrac{1}{\alpha(1-\mu)}+ $$ \dfrac{1}{\beta(\nu-1)}$.

引理 5^[47]. 对于实数$x_1,x_2,\cdots,x_n\in\mathbf{R}$, 有

引理 6^[8]. 令$\Vert P_{\Omega}(\boldsymbol x)-\boldsymbol x\Vert$表示点$ \boldsymbol x $到集合$\Omega$的距离, 则$\Vert P_{\Omega}(\boldsymbol x)-\boldsymbol x\Vert$关于变量$ \boldsymbol x $ 连续, 且其关于变量$ \boldsymbol x $的微分为$\nabla(\Vert P_{\Omega}(\boldsymbol x)\!-\!\boldsymbol x\Vert^2) \!=\! -\!2(P_{\Omega} (\boldsymbol x)-\boldsymbol x).$

2.   分布式连续时间控制算法

针对问题 (2), 基于多智能体固定时间理论, 本节将设计两种基于映射的分布式优化算法, 首先设计一类渐近收敛算法, 其次设计一类指数收敛算法, 并给出算法的收敛性分析.

2.1   算法设计

由引理1可知, 解决问题 (2)的关键是找到最优拉格朗日乘子. 为实现算法分布式演化, 本文为每个智能体都分配一个本地拉格朗日乘子$ \boldsymbol{\lambda}_i\in\mathbf{R}^m $, 基于状态反馈和梯度下降法, 利用局部不等式约束的梯度信息, 每个节点将通过以下算法自主收敛到问题 (2)的最优分配方案:

其中$ \boldsymbol{\tau}_i = (\boldsymbol{\mu}_i+\varphi_i(\boldsymbol x_i))^+ $, 控制参数 $ \alpha,\beta,\gamma>0 $, $\boldsymbol y_i\in $$ \mathbf{R}^m$表示为节点资源变量$ \boldsymbol x_i $的估计, 其最终目的是收敛到问题 (2)的最优值, 辅助变量$ \boldsymbol z_i\in\mathbf{R}^m $是为了平衡资源变量估计值$ \boldsymbol y_i $ 与期望资源量$ \boldsymbol d_i $之间的差值.

与上述映射机制不同, 利用映射算子和固定时间收敛理论, 第二种算法设计如下:

其中固定时间收敛理论是为了确保智能体可在固定时间内将状态收敛到约束范围内.

注 2. 由引理1可知, 求解问题 (2)的最优解关键在于求解最优朗格朗日乘子, 与之对应的等式约束是全局性的, 所以要给每个智能体分配一个局部拉格朗日乘子$ \lambda_i $, 利用经典控制理论的PI (Proportional integral, 比例积分)控制思想, 对其建立基于状态反馈的PI微分方程. 为了将拉格朗日乘子$ \lambda_i $ 与多智能体状态$ x_i $ 联系起来, 本文通过对辅助变量$ y_i $建立一阶迭代算法, 使其沿着局部目标函数的负梯度$ -\nabla f_i(x_i) $方向下降, 然后对$ x_i $ 建立二阶跟踪算法, 最终$ x_i $ 跟踪$ y_i $同步收敛到问题 (2)的最优解. 针对问题 (2), 本文设计了两种连续时间算法, 两者的区别在于采用不同的方法处理局部约束. 其中, 算法 (8)利用局部约束函数的梯度信息将节点状态引导到局部约束$\Omega_i$内, 基于固定时间收敛理论, 而算法 (9)采用映射法将节点状态限制在局部约束$\Omega_i$内部. 在案例仿真中, 算法 (9) 的优势在于高效处理简单的局部约束, 例如电力资源分配问题中的发电机组生产能力约束$ x_i\in[x_i^{\min},x_i^{\max}] $, 其中$ x_i^{\min}, $$ x_i^{\max} $分别表示发电机组$ i $的最小和最大输出功率^[13]. 而算法 (8)的优势在处理局部不等式凸约束, 例如在二次布局问题中, 节点$ i $ 的约束范围为$x_{i,1}^2+x_{i,2}\le10, \; -4\le x_{i,j}\le 4,\; $$ j = 1,2$.

注 3. 在现有二阶多智能体连续时间算法中, 文献[39]使用精确罚函数法将节点局部约束项转移到目标函数中, 即令 $f_i^{\varsigma} = f_i(x_i)+\frac{1}{\varsigma}[(x_{i}^{\min}-x_i)^++ $$ (x_i-x_{i}^{\max})^+]$, 其中$ \varsigma $为惩罚因子, $ x_{i}^{\min}\le x_i\le x_{i}^{\max} $. 进而将原问题等价于解决无约束资源分配问题. 且定义端点$ x_i = x_i^{\max} $ 处的梯度为$\nabla f_i^{\varsigma}\in [\nabla f_i(x_{i}^{\max}), $$ \nabla f_i(x_{i}^{\max})+\frac{1}{\varsigma}]$, 此方法会导致目标函数梯度在端点处不唯一, 且惩罚因子过小容易引起算法在稳定点时产生波动, 过大无法确保算法收敛到最优值. 此外, 该文献只提供了解决一维资源分配问题的算法, 此外相比于文献[39]算法需要交互3个变量信息, 本文设计算法减少了节点间的通信负担. 针对问题 (2), 文献[43]通过微分映射算子$\Pi_{\Omega_i}(\cdot)$解决其中的局部不等式约束, 为了确保算法收敛, 该文献通过奇异摄动法进行收敛性分析, 这就需要一个很小正实数$ \epsilon $来确保算法的收敛, 但并没有给出其收敛范围. 此外该文献的微分映射算子是在本文所设计映射算子上继续求微分得到, 即将变量$ -\nabla f_i(x_i)-\lambda_i $映射到正切锥 $T_{\Omega_i}(x_i)$. 所以其增加了计算约束$\Omega_i$正切锥的额外负担^[31]. 为了确保微分映射算子的有效性, 该文献对节点初始值进行约束.

2.2   算法收敛性分析

在给出算法收敛性分析之前, 先给出一个有关问题 (2)最优解的引理.

引理 7. 假定系统 (9)的平衡点为$(\boldsymbol x^*,\boldsymbol v^*,\boldsymbol y^*, \boldsymbol{\lambda}^*, $$ \boldsymbol{z}^*)$, 则$ \boldsymbol x^* $为问题 (2)的最优解. 反之若$ \boldsymbol x^* $为问题 (2)的最优解, 则存在$ (\boldsymbol x^*,\boldsymbol v^*,\boldsymbol y^*,\boldsymbol{\lambda}^*,\boldsymbol{z}^*) $为系统 (9)的平衡点.

证明 . 由于$ (\boldsymbol x^*,\boldsymbol v^*,\boldsymbol y^*,\boldsymbol{\lambda}^*,\boldsymbol{z}^*) $为系统 (9)的平衡点, 则有

其中最后一个等式可由后续定理1给出, 结合$ \boldsymbol x^* = \boldsymbol y^* $, 所以有$\boldsymbol x^*\in\Omega$, 即引理1中的不等式约束成立. 所以有$ (\boldsymbol{L}\otimes\boldsymbol{I}_m)\boldsymbol{\lambda}^* = 0_{nm} $, 结合矩阵$ \boldsymbol L $ 的定义和假设3, 可知$ \boldsymbol{\lambda}^* = \boldsymbol 1_{n}\otimes\boldsymbol c,\boldsymbol c\in\mathbf{R}^m $. 进而有 $ {\bf 0}_{nm}\in\partial f(\boldsymbol x)+\boldsymbol c\otimes\boldsymbol 1_n $. 紧接着由于有向平衡网络中矩阵$ \boldsymbol L $ 满足 $\boldsymbol 1_{mn}^{\text{T}}(\boldsymbol{L}\otimes $$ \boldsymbol{I}_m) = {\bf 0}_{nm}$, 所以有结论 $ \sum_{i = 1}^n\boldsymbol x_i^* = \sum_{i = 1}^n\boldsymbol d_i^* $. 综上分析, 由引理1可知$ \boldsymbol x^* $为问题 (2)的最优解. 反之亦成立. □

对于系统 (8)也有同样的结论. 可以发现系统 (9)与系统 (8)的区别在如何解决局部不等式约束, 假设$(\boldsymbol x^*, $$ \boldsymbol v^*,\boldsymbol y^*,\boldsymbol{\lambda}^*,\boldsymbol{z}^*,\boldsymbol{\mu}^*,\boldsymbol{\tau}^*)$为系统 (8)的平衡点, 则有$\boldsymbol{\mu}^* = $$ (\boldsymbol{\mu}^*+\varphi(\boldsymbol x))^+ ,\;$ 则必然有结论 $\varphi(\boldsymbol x^*)\le 0,\boldsymbol{\mu}^*\ge0, (\boldsymbol{\mu}^*)^{\text{T}}\times $$ \varphi(\boldsymbol x^*) = 0.$ □

现在开始收敛性分析, 首先给出系统 (8)的收敛定理.

定理 1. 若假设1、2和3成立, 且控制参数满足$\alpha > \dfrac{1}{\rho_2\omega}+\dfrac{1}{2\omega},\beta>\dfrac{\Vert\boldsymbol L\Vert^2\alpha^2}{\rho_2^2},\gamma>\dfrac{\alpha\theta\sqrt{2\varepsilon}}{2}$, 其中$ \varepsilon>0 $. 多智能体系统 (8)渐近收敛到问题 (2)的最优分配方案.

证明. 为简便分析, 不失一般性, 这里只考虑$ m = 1 $情形. 首先对系统 (8)进行坐标变换, 令对上述系统作坐标变换 ${\bar{\boldsymbol x}} = \boldsymbol x-\boldsymbol y,\;{\bar{\boldsymbol v}} = \boldsymbol v- {\dot {{y}}}, \;{\bar{\boldsymbol y}} = \boldsymbol y-\boldsymbol y^*, \; {\bar{\boldsymbol{\lambda}}} = \boldsymbol{\lambda}-\boldsymbol{\lambda}^*, $$ \;{\bar{z}} = z-z^*$. 可得

其中$ \boldsymbol e = \partial f(\boldsymbol x)+\partial\varphi(\boldsymbol y)^{\text{T}}\boldsymbol{\tau}-\partial f(\boldsymbol y^*)+\partial\varphi(\boldsymbol y^*) $. 第二个等式推导如下: $ {\dot{\bar{\boldsymbol{v}}}} = \dot{\boldsymbol v}-\ddot{y} = -\gamma\bar{\boldsymbol{v}}-\gamma\bar{\boldsymbol{x}} $. 在上述第四个等式中, 利用到公式 (10)中的第四个等式. 继续作如下正交变换:

其中, $\xi_1,\;\eta_1,\;\delta_1\in\mathbf{R},\;\boldsymbol{\xi}_2,\;\boldsymbol{\eta}_2,\;\boldsymbol{\delta}_2\in\mathbf{R}^{n-1},\; \boldsymbol r = \dfrac{1}{\sqrt{n}}\boldsymbol 1_n,\; $$ \boldsymbol r^{\text{T}}\boldsymbol R = {\bf 0}_{n-1}^{\text{T}},\;\boldsymbol R^{\text{T}}\boldsymbol R = \boldsymbol I_{n-1},\;\boldsymbol R\boldsymbol R^{\text{T}} = \boldsymbol I_n-\dfrac{1}{n}\boldsymbol 1_n\boldsymbol 1_n^{\text{T}}$. 正交变换之后的系统如下:

其中$,\boldsymbol{\varrho} = \left[\boldsymbol r\;\; \boldsymbol R\right]^{\text{T}}(\boldsymbol{\mu}-\boldsymbol{\mu}^*)$.

定义李雅普诺夫函数如下:

其中,$ 0<\varepsilon<\dfrac{2\gamma^2}{\alpha^2\theta^2} $. 上述李雅普诺夫函数的有效性通过凸函数$f(\boldsymbol x)+\dfrac{1}{2}\Vert\boldsymbol{\tau}\Vert^2$的性质得以保证. 进而函数$ V_1 $在系统 (13)关于时间$ t $的微分为

结合引理3和Young不等式, 我们有结论:

接下来分析$\phi: = \alpha\boldsymbol{\xi}^{\text{T}}[\boldsymbol r,\boldsymbol R]^{\text{T}}\boldsymbol e+\dfrac{\alpha}{2}(\boldsymbol{\varrho}+[\boldsymbol r,\boldsymbol R]^{\text{T}}\bar{\boldsymbol{\tau}})^{\text{T}}.\; $$ (-\boldsymbol{\varrho}+[\boldsymbol r,\boldsymbol R]^{\text{T}}\bar{\boldsymbol{\tau}})$, 我们有

其中式 (17)中的最后两项为

上式中第一个不等式利用了结论$(\boldsymbol{\mu}+\varphi(\boldsymbol y)- $$ \boldsymbol{\tau})^{\text{T}}(\boldsymbol{\tau}^*-\boldsymbol{\tau})\le0$, 其中我们考虑了$ \boldsymbol{\tau},\boldsymbol{\tau}^* $在集合$ \mathbf{R}_{+} $上的映射$ P_{\mathbf{R}_{+}}(\cdot) $. 最后一个不等式利用了函数 $ \varphi(\boldsymbol x) $ 的凸性、$ \boldsymbol{\tau},\boldsymbol{\tau}^* $的非负性、$ (\boldsymbol{\tau}^*)^{\text{T}}\varphi(\boldsymbol x^*) = 0 $和$ \varphi(\boldsymbol x^*)\le $$ 0 $. 将结论 (18)代回式 (17), 同时利用函数$ f(\boldsymbol x) $的李普希茨性质可得

将式(16)和(19)代入式(15), 我们有

由李雅普诺夫定理可知, 系统(8)渐近收敛到问题(2)的最优解. □

定理 2. 若假设1、2和3成立, 且控制参数满足$\alpha > \max\{\dfrac{(\omega+1)}{\rho_2\omega^2}+\dfrac{1}{2\omega},1\},\;\beta>\dfrac{\Vert\boldsymbol L\Vert^2\alpha^2}{\rho_2^2},\;\gamma>\dfrac{\sqrt{\varepsilon}\alpha\theta(1+\omega)}{\omega},$$ 0<\mu<1,\;1<\nu $, 其中$ \varepsilon>0 $, 多智能体系统 (9)指数收敛到问题 (2) 的最优解.

证明. 为简便分析, 不失一般性, 这里只考虑$ m = 1 $情形. 证明将分两步进行. 1) 首先证明算法 (9)在时间$ t\ge T_1 $后节点变量$ y_i $满足条件$y_i\in\Omega_i$, 即$P_{\Omega_i}( y_i)- y_i = 0$, 其中$T_1 = \dfrac{2^{\frac{1-\mu}{2}}}{1-\mu}+\dfrac{2^{\frac{1-\nu}{2}}}{n^{\frac{1-\nu}{2}}(\nu-1)}$; 2) 其次证明当时间$ t\ge T_1 $ 后, 算法(9)指数收敛到问题(2)的最优分配方案.

步骤 1. 令$V_2 = \dfrac{1}{2}\sum_{i = 1}^n\Vert P_{\Omega_i}( y_i)- y_i\Vert^2$, 很容易验证$ V_2 $是有效的李雅普诺夫函数. 且当且仅当$y_i\in\Omega_i,\; \forall i\in $$ V$ 时有$ V_2 = 0 $.$ V_2 $在系统(9)上关于时间$ t $的微分为

其中第一个不等式利用了结论$(P_{\Omega_i}( y_i)- y_i)^{\text{T}} (\nabla f_i $$ ( x_i)+ \lambda_i)\!-\!|P_{\Omega_i}( y_i)\!-\! y_i|\cdot|\nabla f_i( x_i)+ \lambda_i|\le0.$ 由于$\dfrac{1+\mu}{2} < 1$, 结合引理5, 可得

同理, 由于$\dfrac{1+\nu}{2} > 1$, 所以有

将式(22)和(23)代入式(21), 可得

由引理4可知, 在时间$ t\ge T_1 $后, 变量 $y_i\in\Omega_i$, 且固定时间为 $T_1 = \dfrac{2^{\frac{1-\mu}{2}}}{1-\mu}+\dfrac{2^{\frac{1-\nu}{2}}}{n^{\frac{1-\nu}{2}}(\nu-1)}$.

步骤 2. 当$ t\ge T_1 $后, 系统 (9)演变为如下系统

对上述系统做类似于定理1中的坐标变换, 此时有

其中$ \bar{\boldsymbol h} = \nabla f(\boldsymbol x)-\nabla f(\boldsymbol x^*) $. 求解系统 (25)的平衡点等价于求解系统 (26)的平衡点, 且原点就是系统 (26)的平衡点. 令李雅普诺夫函数为

其中, $0 < \varepsilon<\dfrac{\omega^2\gamma^2}{\theta^2(\omega+1)^2\alpha^2}$. $ V_3 $在系统 (26) 上关于时间$ t $的微分为

其中

利用目标函数的强凸性, 我们有

利用目标函数梯度的李普希茨性质, 可得

将上述不等式代回式(29), 可得

同理, 有

由于$ \boldsymbol 1_n^{\text{T}}\boldsymbol{R\eta}_2 = {\bf 0}_n $, 利用引理3, 我们有

进一步, 利用Young不等式, 我们有

将式(32) ~ (35)代入式(28), 可得

其中$ \boldsymbol{\chi} = \text{col}(\bar{\boldsymbol x},\bar{\boldsymbol v},\boldsymbol{\xi},\boldsymbol{\eta},\boldsymbol{\delta}_2) $, $ \kappa $定义如下

从式 (36)可知, 当$ t\ge T_1 $后, 系统(9)指数收敛到问题(2)的最优解. □

3.   案例仿真

现在我们提供仿真案例来验证之前所提算法的有效性. 第一个案例数据来自文献[43], 第二个仿真案例数据来自文献[44].

案例 1. 针对6发电机系统, 每个发动机的代价函数采用如下二次函数形式: $f_i(x_i) = a_{i1}x_i^2+ a_{i2}x_i+ $$ a_{i3}$. 各发电机代价函数系数、初始值$ x_i(0) $、期望功率$ d_i $和输出功率约束总结在表1中. 令算法(8)和算法(9)中控制参数为$\alpha = 2.75,\;\beta = 380,\;\gamma = 30, $$ \;\mu = 0.5,\;\nu = 2$. 从任意初始状态出发, 详细的仿真轨迹分别如如图1和图2所示, 从图中可知问题的最优值为$\boldsymbol x^* = [23.01;\;35;50; \;31.01; $$ \;45.78;\; 30.19]^{\rm{T}}$, 此结果与文献[43]一致. 相比于算法(8)中的梯度法, 由于算法(9)采用映射算子和固定时间收敛理论处理问题(2)中的局部不等式约束, 所以算法(9)可以很快地将超出约束范围的节点快速拉回到约束范围内. 所以图1相比图2, 在算法到达约束边界时, 有更大的波动性. 基于案例1中的数据, 接下来从同一初始状态出发, 将本文所设计算法与文献[43]和文献[39]进行收敛速度对比, 结果如图3所示(文献[43]中控制参数设置与本文相同). 其中算法误差定义为$e = \ln(\dfrac{1}{2}\sum_{i = 1}^n(x_i-x_i^*)^2)$, 其中最优值$ x_i^* $由集中式拉格朗日乘子法求得. 则 $ e $ 表示算法运行解与理论最优解的差值. 从图3可以看出, 本文所设计算法下误差下降较快, 即拥有更快的收敛速度, 且算法在稳定点时有更好的稳定性. 而文献[39]由于在端点处梯度不唯一, 会导致算法在稳定点时有轻微波动. 此外, 算法 (9)和文献[39]在边界点的波动性导致其误差$ e $要比算法 (8)和文献[43]中的误差稍大.

表 1  案例1各节点参数

Table 1  System parameters in case 1

Agent	$a_{i1}$	$a_{i2}$	$a_{i3}$	功率约束	$d_i$	$x_i(0)$
1	2	3	0.5	$[20,40]$	45	40
2	1	4	1.5	$[25,35]$	40	24
3	0.5	5	3	$[35,50]$	$25$	35
4	1.5	2	1	$[25,45]$	35	45
5	1	3.5	2.5	$[30,47]$	30	28
6	1.5	4.5	2	$[28,42]$	40	50

下载: 导出CSV 

| 显示表格

图 1  案例1中算法(8)的各发电机曲线图

Fig. 1  The trajectories of each generator by algorithm (8) in case 1

下载: 全尺寸图片 幻灯片

图 2  案例1中算法(9)的各发电机曲线图

Fig. 2  The trajectories of each generator by algorithm (9) in case 1

下载: 全尺寸图片 幻灯片

图 3  基于案例1数据的算法性能对比

Fig. 3  The comparison between the existence algorithms and ours based on case 1

下载: 全尺寸图片 幻灯片

案例 2. 针对4节点环形网络, 其中资源决策变量$\boldsymbol x_i\in{\bf{R}}^2$, 且节点目标函数分别为

其中

虚拟资源量$ \boldsymbol d $设定为$\boldsymbol d_1 \!=\! [2,1]^{\text{T}},\;\boldsymbol d_2 \!=\! [2,3]^{\text{T}}, \;\boldsymbol d_3 = $$ [2,4]^{\text{T}},\;\boldsymbol d_4 = [1,5]^{\text{T}}$. 算法初始值设定为$\boldsymbol x_1(0) = [2,0]^{\text{T}}$, $\boldsymbol x_2(0) \!=\! [1.5,0.5]^{\text{T}},\;\boldsymbol x_3(0) \!=\! [1,1]^{\text{T}},\;\boldsymbol x_4 = [4,6]^{\text{T}}$, 算法控制参数设定为$\alpha = 2,\;\beta = 75,\;\gamma = 30, \;\mu = $$ 0.2,\;\nu = 2$. 其仿真曲线分别如图4和图5所示, 星点表示各节点的最优解位置. 通过仿真案例2可以发现, 本文所设计算法对高维资源分配问题亦是有效的. 为了展示各算法在二维资源分配问题上的收敛速度, 现将本文所设计算法与文献[43]进行对比, 收敛速度图展示在图6, 从图中可以看出, 指数收敛算法 (9)比渐近收敛算法 (8)收敛稍快, 该结果与一维资源问题图3所得结论一致.

图 4  案例2中算法 (8)的各节点仿真结果轨迹

Fig. 4  The trajectories of each agent by algorithm (8) in case 2

下载: 全尺寸图片 幻灯片

图 5  案例2中算法 (9)各节点仿真结果轨迹

Fig. 5  The trajectories of each agent by algorithm (9) in case 2

下载: 全尺寸图片 幻灯片

图 6  基于案例2数据的算法性能对比

Fig. 6  The comparison between the existence algorithms and ours based on case 2

下载: 全尺寸图片 幻灯片

案例 3. 针对IEEE-39总线系统, 包含10个发电单元, 各发电单元的通信链路及其权重如下图7所示. 各发电单元的发电成本为$f_i(x_i) = a_{i,1}x_i^2+ a_{i2}x_i$且有${\boldsymbol{a}}_1=[0.096$$0.072,\;0.105,\;0.082,\;0.078,\; 0.90,\; 0.085,\; 0.092,\; 0.155,$$\;0.082]^{\text{T}}\in\mathbf{R}^{10}$ 和 $\boldsymbol a_2 = [1.22,\;\; 3.41,\;\;2.53,\;\;4.02,\; 2.90,$$2.72,\;1.22,\;4.20,\;\;3.53,\;\;3.02]^{\text{T}}\;\in \; \mathbf{R}^{10},$ 其中 $\boldsymbol a_1 = $$ [ a_{1,1},\cdots,a_{n,1}],\boldsymbol a_2 = [a_{1,1},\cdots,a_{n,1}]\in \mathbf{R}^{10}$. 各节点的局部不等式约束和虚拟分配量分别为$ x_i\in[0,300] $和$ d_i = 250 $. 以零初始状态出发, 在控制参数 $\alpha = 2.75,\; $$ \beta = 380,\;\gamma = 30$作用下, 算法 (8)和 (9)各自对应的状态轨迹如图8和图9所示. 此外两种算法的误差$ e $轨迹如图10所示, 其中误差$ e $定义同案例1. 从图10可知, 上述系统最终收敛到最优分配方案$\boldsymbol x^* =$[243.6524, 300, 216.5298, 268.1785, 289.1107, 251.5626, 275.1839, 238.0504, 143.4557, 274.2760] .

图 7  案例3中各发电单元的通信链路

Fig. 7  The communication links of each generator in case 3

下载: 全尺寸图片 幻灯片

图 8  案例3中算法 (8)的各节点仿真结果轨迹

Fig. 8  The trajectories of each agent by algorithm (8) in case 3

下载: 全尺寸图片 幻灯片

图 9  案例3中算法 (9)各节点仿真结果轨迹

Fig. 9  The trajectories of each agent by algorithm (9) in case 3

下载: 全尺寸图片 幻灯片

图 10  案例3中算法 (8)和 (9)的收敛速度比较

Fig. 10  The convergence rate comparison between algorithms (8) and (9) in case 3

下载: 全尺寸图片 幻灯片

4.   总结

本文针对有向平衡网络下的带约束分布式资源分配问题, 基于固定时间映射算子和基于多智能体一致性的梯度下降法, 设计了两种基于二阶多智能体系统的连续时间算法. 所设计算法对节点初始状态无约束, 并且控制参数全为常数. 通过对一维和二维资源分配问题数值仿真, 本文所提算法的有效性也得以保证. 然而在实际系统中, 通信网络通常会存在信息交互延迟, 而信息延迟给算法分析带来了很大的困难. 此外, 非平衡有向网络在实际系统中也是存在的, 所以在后续研究中, 我们将考虑由通信延迟等造成的有向非平衡网络下二阶多智能体资源分配问题.