近年来，经济调度在电力市场、金融市场、风力发电厂、智能电网和能源资源分配等领域引起了广泛的关注. 为了解决经济调度问题，许多学者提出大量的集中式算法，例如梯度搜索和拉格朗日迭代等^[1-2]. 集中式算法对于大型节点网络而言是不合适的. 因为集中式算法需要1个控制中心来收集整个网络的信息，会带来一些性能上的约束，比如较高的通信需求和通信消耗、灵活性差、计算负担重等问题.

为了克服上述集中式算法所带来的缺陷，许多学者提出大量基于离散系统的分布式优化算法^[3-18]，例如梯度算法^[3-6]、残差法^[7-9]、push-sum算法^[10-11]、二分法^[12]、一致性算法^[13-16]、交替乘子法^[17]、粒子群算法^[18]等. 近年来，更多的学者寻找更合适的连续系统算法来解决分布式优化问题^[19-25]. 从连续系统角度来看，通过微分包含理论等分析方法，可以得到更强的最优协调控制理论条件. 因为连续系统对应目标函数最陡的梯度方向. 通过罚函数法消掉不等式约束，Cherukuri等^[26]提出针对一阶连续多智能体系统的分布式协调算法. 受到上述文献的激发，He等^[27]设计集中式二阶算法.Wang等^[28]研究分布式2阶系统的动态方程. 与上述方法不同，Hoang等^[29]以鞍点动态方程出发，提出分布式控制算法用于解决最优资源分配问题。Zheng等^[30]以多智能体一致性为基础，提出分布式算法.

上述算法都是从无向图或有向平衡网络出发，无向网络算法中都需要构建一系列双随机矩阵，阻碍了这些算法的开发和实际应用，特别是在随时间变化的一般非平衡有向实际网络环境中，因为双随机矩阵的条件难以以分布式方式满足. 在现代电力系统的应用中，无向时变网络已经不能完全满足智能电网中复杂度更高以及不对等通信条件下的发展要求^[9].

本文针对有向网络下的电力系统经济分配问题，通过建立鞍点动态方程，设计基于连续系统的分布式1阶连续系统算法，且所设计算法只需每个节点各自的入度.

$G = (V,E,{A})$表示所研究问题对应的网络图，其中 $V{\rm{ = \{ 1,2,}}\cdots{\rm{,}}n{\rm{\} }}$为网络节点集合， $E \subseteq V \times V$为有向边集合. 对于任意的边 $(i,j) \in E$，节点 $i$记作是节点 $j$的出邻居或节点 $j$是节点 $i$的入邻居. 令 $N_i^{\rm{in}}: = \{ j|(i,j) \in E\} $和 $N_i^{\rm{out}}: = \{ j|(j,i) \in E\} $分别记作节点i的入邻居集和出邻居集. 本文的网络设置为不包含自回路，即 $i \notin N_i^{\rm{in}}.$令 $d_i^{\rm{in}} = \left| {N_i^{\rm{in}}} \right|$和 $d_i^{\rm{out}} = \left| {N_i^{\rm{out}}} \right|$分别为节点 $i$的入度与出度. 相应的，图 $G$所对应的伴随矩阵 ${A}$定义为 ${a_{ij}}: = {1}/({{1 + d_i^{\rm{in}}}}).$ 若有 $(i,j) \in E$，反之 ${a_{ij}}: = 0$. 此外，令 ${a_{ii}}: = {1}/({{1 + d_i^{\rm{in}}}})$. 特别地，针对无向图， ${A}$是对称矩阵. 节点 $i$到节点 $j$有通路即表示存在一组长度为 $l$的有向边 $(j,{i_{l - 1}}),({i_{l - 1}},{i_{l - 2}}),\cdots,({i_1},i)$. 若图 $G$中任意2个节点间都存在有向通路，则称图 $G$是强连通.

考虑带 $n$节点的经济调度问题，为每个节点 $i$设置一个变量 ${x_i}$，且每个节点拥有自己的局部代价函数 ${f_i}({x_i}):{\bf R} \to {\bf R}$. 本文解决的问题如下：

(1a) $\tag{1a}\min f({{x}}) = \sum\limits_{i = 1}^n {{f_i}({x_i})} ;$

s.t.

(1b) $\tag{1b}\sum\limits_{i = 1}^n {{x_i}} = X,$

(1c) $\tag{1c}{x_i} \in {X_i} = [{\underline x _i},{\bar x_i}].$

式中： $X$表示整体的网络资源需求， ${\underline x _i}{\text{、}}{\bar x_i}$分别表示各自节点的生产能力约束. 定义如下拉格朗日函数：

(2) $\begin{split} L({{x}},\mu ,{{{\lambda}} _m},{{{\lambda}} _M}) =& f({{x}}) + \mu ({\bf{1}}_n^{\rm{T}}{{x}} - X)+ \\ & \sum\limits_{i = 1}^n {{\lambda _{{m_i}}}({{\underline x }_i} - {x_i})} + \sum\limits_{i = 1}^n {{\lambda _{{M_i}}}({x_i} - {{\bar x}_i})} . \end{split} $

式中： $\mu {\text{、}}{\lambda _{{m_i}}}{\text{、}}{\lambda _{{M_i}}}$分别为针对等式约束、节点生产能力下限和生产能力上限的拉格朗日乘子。定义 ${{{\lambda}} _m} = {[{\lambda _{{m_1}}},{\lambda _{{m_2}}},\cdots,{\lambda _{{m_n}}}]^{\rm{T}}},{{{\lambda}} _M} = {[{\lambda _{{M_1}}},{\lambda _{{M_2}}},\cdots,{\lambda _{{M_n}}}]^{\rm{T}}},$针对式（2），下面的引理很重要.

引理 1^[29] 假设 ${f_i}$为凸函数。 $({{{x}}^ * },{\mu ^ * },{{\lambda}} _m^ * ,{{\lambda }}_M^ * )$为函数 $L$的鞍点，即 ${{{x}}^ * }$为问题（1）的最优解当且仅当下面的条件都满足：

(3a) $\tag{3a}\nabla f({{{x}}^ * }) + {\bf 1}_n{\mu ^ * } - {{\lambda}} _m^ * + {{\lambda}} _M^ * = 0;$

(3b) $\tag{3b}{\bf{1}}_n^{\rm{T}}{{{x}}^ * } = X,\quad x_i^ * \in {X_i},\quad\forall i \in V;$

(3c) $\tag{3c}\lambda _{{m_i}}^ * ({\underline x _i} - x_i^ * ) = 0,\quad\forall i \in V;$

(3d) $\tag{3d}\lambda _{{M_i}}^ * (x_i^ * - {\bar x_i}) = 0,\quad\forall i \in V.$

式（3c）表明： $\lambda _{{m_i}}^ * \ne 0 \Rightarrow x_i^ * = {\underline x _i},$和 $x_i^ * \ne {\underline x _i} \Rightarrow \lambda _{{m_i}}^ * = 0.$ 同理， $\lambda _{{M_i}}^ * $可得类似结论.

假设 1 为了解决有向网络下的经济调度问题，作出如下假设.

1） 网络节点图 $G$是强连通的.

2） 对于任意的节点 $i$，代价函数 ${f_i}({x_i})$都是 ${\varpi _i}$强凸函数，即 $(x - y)({f_i}(x) - {f_i}(y)) \geqslant {\varpi _i}{\left| {x - y} \right|^2},$ $\forall x,y \in {X_i}.$

上述假设在分布式优化算法中是合理的，更多的相关知识介绍可以参考文献[12].

为了实现以分布式的方式获取最优解,对于每个节点 $i$，分配一个拉格朗日乘子 ${\lambda _i}$. 当所有的局部拉格朗日乘子一致收敛于 ${\mu ^ * }$时，可以获得最优解.Hoang等^[29]设计算法如下：

(4a) $\tag{4a}{\dot y_i} = \sum\limits_{i = 1}^n {{a_{ij}}({\lambda _i} - {\lambda _j})} ,$

(4b) $\tag{4b}{\dot \lambda _i} = ({x_i} - {r_i}) - \sum\limits_{i = 1}^n {{a_{ij}}({\lambda _i} - {\lambda _j})} - {y_i},$

(4c) $\tag{4c}{\dot \lambda _{{m_i}}} = {k_{{m_i}}}{\lambda _{{m_i}}}({\underline x _i} - {x_i}),$

(4d) $\tag{4d}{\dot \lambda _{{M_i}}} = {k_{{M_i}}}{\lambda _{{M_i}}}({x_i} - {\bar x_i}),$

(4e) $\tag{4e}{\dot x_i} = - {k_i}(\nabla {f_i}({x_i}) + \lambda _i^{} - {\lambda _{{m_i}}} + {\lambda _{{M_i}}}).$

式中： ${k_i}{\text{、}}{k_{{m_i}}}{\text{、}}{k_{{M_i}}}$为控制参数； ${r_i}$表示节点 $i$的虚拟分配方案， $\displaystyle\sum\nolimits_{i = 1}^n {{r_i}} = X$； $\displaystyle\sum\nolimits_{i = 1}^n {{y_i}} = 0$. 设定 ${\lambda _{{m_i}}} > 0,$则由指数函数的性质可知， ${\lambda _{{m_i}}}(t) \geqslant 0,\forall i \in V$. 同理可得， ${\lambda _{{M_i}}}(t) \geqslant $0. ${ y}{\text{、}}{\lambda}{\text{、}}{{\lambda} _m}{\text{、}}{{\lambda} _M}{\text{、}}{ x}$分别为 ${y_i}{\text{、}}{\lambda _i}{\text{、}} $ ${\lambda _{{m_i}}}{\text{、}}{\lambda _{{M_i}}}{\text{、}}{x_i}$的向量形式.

1.3节中，针对无向网络下的经济调度问题，Hoang等^[29]设计基于多智能体连续系统和不动点理论的分布式算法。但该算法无法解决有向网络下的经济调度问题，因为有向图的拉普拉斯矩阵 ${{L}}$不再具有性质 ${\bf{1}}_n^{\rm{T}}{{L}} = {\bf{0}}_n^{\rm{T}}.$由文献[7]可知，存在正向量 ${{\xi}} = {[{\xi _1},{\xi _2},\cdots,{\xi _n}]^{\rm{T}}},$满足 ${{{\xi}} ^{\rm{T}}}{{L}} = {\bf{0}}_n^{\rm{T}}.$对式（4b）左乘 ${{{\xi}} ^{\rm{T}}}$，可得 $\sum\nolimits_{i = 1}^n {{\xi _i}({x_i} - {r_i})} = 0$，即经济调度问题（1）中的负荷等式约束无法得到满足. 条件 ${\bf{1}}_n^{\rm{T}}{{L}} = {\bf{0}}_n^{\rm{T}}$主要是保证了无向图中各发电机组权重一致. 为了实现有向网络下的发电机组权重平衡，对算法（4）作出如下修改：

(5a) $\tag{5a}{\dot {{z}}_i} = - \sum\limits_{i = 1}^n {{a_{ij}}({{{z}}_i} - {{{z}}_j})} ,$

(5b) $\tag{5b}{\dot y_i} = \alpha \beta \sum\limits_{i = 1}^n {{a_{ij}}({\lambda _i} - {\lambda _j})} ,$

(5c) $\tag{5c}\begin{array}{l} {{\dot \lambda }_i} = \alpha {(z_i^i)^{ - 1}}({x_i} - {r_i}) - \theta {y_i} - \beta \displaystyle\sum\limits_{i = 1}^n {{a_{ij}}({\lambda _i} - {\lambda _j})} , \\ \end{array} $

(5d) $\tag{5d}{\dot \lambda _{{m_i}}} = {k_{{m_i}}}{\lambda _{{m_i}}}({\underline x _i} - {x_i}),$

(5e) $\tag{5e}{\dot \lambda _{{M_i}}} = {k_{{M_i}}}{\lambda _{{M_i}}}({x_i} - {\bar x_i}),$

(5f) $\tag{5f}{\dot x_i} = - {k_i}(\nabla {f_i}({x_i}) + \lambda _i^{} - {\lambda _{{m_i}}} + {\lambda _{{M_i}}}),$

其中所添加的平衡变量 ${z_i}$用来平衡有向网络的权重，且初始值为 ${{ z}_i} = k{[0,\cdots,0,{1_i},0,\cdots,0]^{\rm{T}}}.$ 算法（5）的收敛判据为： $\Delta {{\lambda}} = \displaystyle\sum\nolimits_{i,j} {\left| {{\lambda _i} - {\lambda _j}} \right|} < \varepsilon $且 ${{\lambda}} $持续不变.

算法的运行流程如下。

1）初始化控制参数 $\alpha {\text{、}}\beta {\text{、}}{k_i}{\text{、}}{k_{{m_i}}}{\text{、}}{k_{{M_i}}}$和虚拟分配负荷 ${r_i}$，变量 ${\lambda _i}{\text{、}}{y_i}{\text{、}}{x_i}$零初始化和变量 ${\lambda _{{m_i}}}{\text{、}}{\lambda _{{M_i}}}$置1， ${{{z}}_i} = {[0,\cdots,0,{1_i},0,\cdots,0]^{\rm{T}}}.$

2）通过式（5a）计算平衡变量 ${{z}}.$

3）通过式（5b）计算辅助变量 ${{y}}.$

4）通过式（5c）~（5e）获取拉格朗日乘子 ${{\lambda}} {\text{、}}{{{\lambda}} _m}{\text{、}}{{{\lambda}} _M}.$

5）通过映射（5f）获取功率分配值 ${{x}}.$

6）判断算法是否收敛至最优值。若是，则结束；若否，则转至2）步继续计算.

在给出收敛性证明之前，先给出以下引理.

引理2^[23] 　在假设1满足的条件下，存在一个正交矩阵 ${ Q} = [{ q},{ R}] \in {{\bf R}^{n \times n}}$，其中 ${ q} = {{\bf{1}}_n}/{{\sqrt n }}$，有

(6a) $\tag{6a}{{L}} = \left[ {\begin{array}{*{20}{c}} { q}&{ R} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\bf 0}&{{\bf{0}}_{n - 1}^{\rm{T}}} \\ {{{\bf{0}}_{n - 1}}}&{{{{\varLambda}} _1}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{ q}^{\rm{T}}}} \\ {{{ R}^{\rm{T}}}} \end{array}} \right],$

(6b) $\tag{6b}{ R}{({{{\varLambda}} _1})^{ - 1}}{{ R}^{\rm{T}}}{{L}} = {{LR}}{({{{\varLambda}} _1})^{ - 1}}{{ R}^{\rm{T}}} = {{{K}}_n}.$

式中：

${{{K}}_n} = {{{I}}_n} - {{\bf{1}}_n}{\bf{1}}_n^{\rm{T}}/{n},\;{{{\varLambda}} _1} = {\rm{diag}}\;[ {\rho _2}({{L}}),\;{\rho _3}({{L}}),\cdots, {\rho _n}({{L}})]. $

引理3^[25] 　假设本文讨论有向网络 $G$是强连通的，则对应的拉普拉斯矩阵 ${{L}}$满足如下性质.

1） 矩阵 ${{L}}$存在一个零特征根和对应的正左特征向量 ${{\xi}} = {[{\xi _1},{\xi _2},\cdots,{\xi _n}]^{\rm{T}}},$满足 $\displaystyle\sum\nolimits_{i = 1}^n {{\xi _i}} = 1$，即 ${\xi ^{\rm{T}}}{{L}} = {\bf{0}}_n^{\rm{T}}.$

2） $\exp\; ( - {{L}}t)$是一个对角元素非负的对角矩阵.

3） $B = \mathop {\lim }\limits_{t \to \infty } \exp\; ( - {{L}}t) = {{\bf{1}}_n}{{\xi} ^{\rm{T}}}.$

定理1 令假设1成立且控制参数满足

$\quad\quad\quad\dfrac{\alpha }{{2(1 + \alpha ){\varpi _{\min }}}} < \delta < \theta {\xi _{\min }},\;\beta > \dfrac{{{\alpha ^2}{\theta ^2}}}{{2{\rho _2}({{L}})\delta }},$

则算法（5）收敛到最优解，即

(7) $ \mathop {\lim }\limits_{t \to \infty } {x_i}(t) = x_i^ * ,\quad\forall i \in V.$

证：收敛性证明可参考文献[29,30].

针对经济调度问题（1），考虑有向网络情形，提出基于连续系统的分布式优化算法. 利用数值仿真案例来验证上述算法的有效性. 采用著名的IEEE-14标准测试系统，连接网络如图1所示. 节点的代价函数定义为 ${f_i}({x_i}) = {a_{i,1}}x_i^2 + {a_{i,2}}{x_i} + {a_{i,3}},$具体参数定义如表1所示^[31]. 设定总负荷约束为 $X = 300$，每个节点的虚拟负荷设定为 ${{r}} = {[60,60,60,60,60]^{\rm{T}}}.$

在该仿真案例中，设置系统参数如下： ${k_i} = 10, $ ${k_{{m_i}}} = {k_{{M_i}}} = 0.1,\alpha = 3,\beta = 2,\theta = 2.$仿真结果如图2~6所示. 从图2可知，问题（1）的最优解为 ${{ x}^ * } = [66.60, $ $70,47.55,55.48,60.36]^{\rm{T}}.$该结果与Yang等^[11]所提的集中式算法结果一致. 从引理1的平衡条件（3a）可知，算法（5c）中的λ_i要协同趋于最优值−7.328 8，从图3可以看出该结论. 因为要考虑经济调度问题中各节点生产能力约束问题，算法（5）中设定有辅助拉格朗日乘子 ${\lambda _{{m_i}}}{\text{、}}{\lambda _{{M_i}}}.$该变量的设定目标是将超出生产能力约束的节点功率分配值拉回约束范围内. 从图4可以看出，节点2的功率分配值达到节点2的生产上限，所以有 $\lambda _{{M_2}}^ * \ne 0$，表示节点2的最优分配值是边界上限值，该结果也可以从图2看出，节点2的最优分配值为70 MW. 平衡变量 $z_i^i$的作用是解决有向图下的权重平衡问题. 若算法（5a）中的初始值设定为 ${{{z}}_i} = [0,\cdots,0,{1_i}, $ $0,\cdots,0]^{\rm{T}},$则 $z_i^i$收敛至引理3中的正实数 ${\xi _i}$且满足 $\displaystyle\sum\nolimits_{i = 1}^n {{\xi _i}} = 1.$该结论可以从图5得出. 辅助变量 ${y_i}$的作用是平衡各节点真实分配值 ${x_i}$与虚拟设定负荷 ${r_i}$的差值. 以节点1为例，差值为6.6，所以 ${y_1}= $ ${{\alpha ({x_1} - {r_1})}}/({{{\xi _1}\theta }}) = $ ${{3 \times 6.6}}/({{0.29 \times 2}}) = 34.18$，该结果与图6一致. 算法（5）对形如问题（1）的经济调度问题是完全有效的. 为了验证平衡变量 $z_i^i$ 初始值对收敛性的影响，定义误差为 $e = \ln $ $ \left(\displaystyle\sum\nolimits_{i = 1}^n {{{({x_i} - x_i^ * )}^2}}\bigg/2 \right)$，针对不同的初始值，收敛速度的对比如图7所示. 随着 $k$的增大，算法收敛速度增大；反之，收敛速度下降. 参数 $k$对算法的收敛速度影响原理如下：随着 $k$的减小，对应的辅助变量 ${{y}}$变大，因而其从零增长至最优值的时间增加；反之收敛速度变快. 一般 $k$默认为1.

针对智能电网中的经济调度问题，基于鞍点动态方程理论，本文研究有向网络下的基于连续系统的分布式优化算法. 在带约束的分布式优化问题中，最优点的局部梯度不一定都是零，本文通过增加一辅助变量y_i来平衡各节点的局部最优梯度，可以实现控制参数的常量化. 针对非平衡有向网络的经济调度问题，由于此时网络的权值是不平衡的，本文通过引入变量 $z_i^i$来消除这种不平衡. 经仿真案例验证可知，所设计算法是非常有效的. 上述算法迭代过程中需要节点之间进行连续通信，在大规模网络中，这对网络通信消耗巨大. 在接下来的研究中，考虑到连续通信所带来的巨大资源消耗，将事件触发通信机制引入经济调度问题是未来的研究方向之一.