一篇论文的整体思路和复现记录(一，数学篇)

一篇论文的整体思路和复现记录

论文：
L. Yang, J. Luo, Y. Xu, Z. Zhang and Z. Dong, “A Distributed Dual Consensus ADMM Based on Partition for DC-DOPF With Carbon Emission Trading,” in IEEE Transactions on Industrial Informatics, vol. 16, no. 3, pp. 1858-1872, March 2020, doi: 10.1109/TII.2019.2937513.

一个分布式对偶共识ADMM，基于直流动态最优潮流的分块，带有碳排放交易的ADMM

文章是怎么解决问题的呢？

摘要

文章给出了一个分布式的交替方向乘子法，用以解决在碳排放的情况下直流动态最优潮流的问题。一般来说，基于ADMM的分布式方法在相邻的子系统中公开边界节点和边界支路的信息。与这些方法相反，提出的方法（DC-ADMM-P）给出了一个新的策略，通过共识ADMM来解决DC-ADMM-CET的对偶问题，而仅仅需要公开相邻子系统之间的边界支路的公开信息。另外，DC-ADMM-P的收敛表现也因为减少了对偶乘子（dual multipliers）的数量并且采用了一个乘子提升的改进步长从而得到提升。
DC-ADMM-P在6到1062节点案例都得到了测试，并且与其他的分布式和去中心化方法进行了对比。仿真结果证实，DC-ADMM-P在解决DC-DOPF问题上的高效率，并且带有复杂的非线性因子，从而能够被看作凸的分布式函数。与此同时，也通过减少对偶乘子并且使用一个新的乘子更新策略实现了收敛表现的提升。

注意的点

从中可以看出文章有几个重要的点：
1.共识ADMM来解决DC-ADMM-P问题的对偶问题；
2.在求解问题的过程中，仅仅需要公开相邻子系统之间边界支路的公开信息，而不是像原来那样需要相邻子系统之间边界节点和边界支路的信息，从而更好的保护了不同子系统信息的私密性。
3.收敛表现因为对偶乘子数量的减少得到提升。

公式

一些符号及其含义：

符号	含义
$S^B$	所有节点集合
$S^G$	发电机组节点集合
$S^{TP}$	时间集合

目标函数和约束条件

目标函数： 交易总成本最低：
$$\begin{array}{r}{F}_{TC}=F_{TH}+F_{CET}\end{array}\text{\hspace{1em}(1)}$$
（1）热电机组的燃料成本：
$$\begin{array}{r}{F}_{TH}=\sum _{i\in S^G}\sum _{t\in S^{TP}}{f}_i(P_{i,t}^G)\end{array}\text{\hspace{1em}(2)}$$
机组$i$的生产成本与出力之间的关系用二次函数的形式刻画：
$$\begin{array}{r}{f}_i(P_{i,t}^G)={\alpha }_i+{\beta }_i{P}_{i,t}^G+{\gamma }_i(P_{i.t}^G{)}^2\end{array}\text{\hspace{1em}(3)}$$
（2）碳排放成本
$$\begin{array}{r}{F}_{CET}={\pi }_b\Delta E_b-{\pi }_s\Delta E_s\end{array}\text{\hspace{1em}(4)}$$

约束：
（1）节点功率平衡约束
$$\begin{array}{r}{P}_{i,t}^G-P_{i,t}^D=\sum _{j\in S^B}{B}_{ij}{\theta }_{j,t}\end{array}\text{\hspace{1em}(5)}$$
(2)传输线的容量约束
$$\begin{array}{r}|\frac{{\theta }_{i,t}-{\theta }_{j,t}}{x_{ij}}|\le F_{ij}^{max}\forall t\in S^{TP}\forall ij\in S^L\forall i,j\in S^B\end{array}\text{\hspace{1em}(6)}$$
(3)参考节点
$$\begin{array}{r}\theta =0\forall i\in S^{ref}\forall t\in S^{TP}\end{array}\text{\hspace{1em}(7)}$$
(4)除发电机节点以外的节点功率约束
$$\begin{array}{r}{P}_{i,t}^G=0\forall i\in S^B-S^G\forall t\in S^{TP}\end{array}\text{\hspace{1em}(8)}$$
(5)发电机的容量控制
$$\begin{array}{r}\underline{P_i^G}\le P_{i,t}^G\le \overline{P_i^G}\forall i\in S^G\forall t\in S^{TP}\end{array}\text{\hspace{1em}(9)}$$
(6)发电机的爬坡上升约束
$$\begin{array}{r}{P}_{i,t}^G-P_{i,t-1}^G\le P_{i,t-1}^{G,up}\forall i\in S^G\forall t\in S^{TP}\end{array}\text{\hspace{1em}(10)}$$
(7)发电机的爬坡下降约束
$$\begin{array}{r}{P}_{i,t-1}^G-P_{i,t}^G\le P_{i,t-1}^{G,down}\forall i\in S^G\forall t\in S^{TP}\end{array}\text{\hspace{1em}(11)}$$
(8)碳排放容量约束
$$\begin{array}{r}\sum _{i\in S^G}\sum _{t\in S^{TP}}{E}_i(P_{i,t}^G)\le E_0+\Delta E_b-\Delta E_s\end{array}\text{\hspace{1em}(12)}$$
其中：出力对应的碳排放可用二次函数表示：$E_i(P_{i,t}^G)=a_i+b_i{P}_{i,t}^G+c_i(P_{i,t}^G{)}^2$
(9)可供购买的最大排放限额
$$\begin{array}{r}0\le \Delta E_b\le \Delta E_b^{max}\end{array}\text{\hspace{1em}(13)}$$
(10)可供销售的最大排放限额
$$\begin{array}{r}0\le \Delta E_s\le \Delta E_x^{max}\end{array}\text{\hspace{1em}(14)}$$

化简步骤

在确定数学模型后，文中令$i\in S^B$,
$P^D\in [P_1^D,P_2^D,\dots ,P_N^D]$,
$P_i^G=[P_{i,1}^G,P_{i,2}^G,\dots ,P_{i,T}^G]$，
${\theta }_i=[{\theta }_{i1},{\theta }_{i2},\dots ,{\theta }_{iT}]$,
$P_i^D=[P_{i1}^D,P_{i2}^D,\dots ,P_{iT}^D]$,
$x_0=[\Delta E_b,\Delta E_s]$.
带入前述（5）式：
$$\begin{array}{r}{P}_i^G-P_i^D=\sum _{j\in S^B}{B}_{ij}{\theta }_j\forall i\in S^B\end{array}\text{\hspace{1em}(15)}$$
$$\begin{array}{r}-F_{ij}^{max}\mathbf{e}\le \frac{{\theta }_i-{\theta }_j}{x_{ij}}\le F_{ij}^{max}\mathbf{e}\forall ij\in S^L,\forall i,j\in S^B\end{array}\text{\hspace{1em}(16)}$$
另外，设：$b_{PF}=P^D$, $\widetilde{A_0}=0$,
$$\begin{array}{c}\widetilde{A_i}=\left[\begin{array}{cc}0& -B_{1i}I\\ \dots & \\ I& -B_{ii}I\\ \dots & \\ 0& -B_{Ni}I\end{array}\right]\end{array}$$（此处由于KATEX不支持\multicolumn，无法将省略号居中放置\，凑合看看就好）
由（15）得到：
∑ i ∈ S B ∪ { 0 } A ~ i x i = b P F (17) \begin{aligned} \sum_{i\in S^B\cup\{0\}}\tilde{A}_ix_i=b_{PF} \end{aligned}\tag{17} i∈SB∪{0}∑​A~i​xi​=bPF​​(17)
相似地，设${\tilde{M}}_0=0$，由（16）可得：
∑ i ∈ S B ∪ { 0 } M ~ i x i ≤ b B R (18) \begin{aligned} \sum_{i\in S^B\cup\{0\}}\tilde{M}_ix_i\leq b_{BR} \end{aligned}\tag{18} i∈SB∪{0}∑​M~i​xi​≤bBR​​(18)

另外，对应碳排放约束，有${\tilde{h}}_0(x_0)=[-1,1]x_0=\Delta E_s-\Delta E_b$， 对（12）式可做变换：
∑ i ∈ S B ∪ { 0 } h ~ i ( x ) ≤ b E (19) \begin{aligned} \sum_{i\in S^B\cup\{0\}} \tilde{h}_i(x)\leq b_{E} \end{aligned}\tag{19} i∈SB∪{0}∑​h~i​(x)≤bE​​(19)
其中，$b_E=E_0$, ${\tilde{h}}_i(x_i)={\sum }_{t\in S^{TP}}{E}_i(P_{i,t}^G)（i\in S^G）$以及${\tilde{h}}_i(x_i)=0（i\in S^B-S^G）$

对（7）-（11）式，由于每个机组单独解耦，定义 X i = { x i ∣ ( 7 ) − ( 11 ) , i ∈ S B } \Chi_i=\{x_i|(7)-(11), i\in S^B\} Xi​={xi​∣(7)−(11),i∈SB}，即为对于每个节点（bus）均满足（7）到（11）的约束条件；
对应热电机组燃料成本：$g_i(x_i)=0（i\in S^B-S^G）$，$g_i(x_i)={\sum }_{t\in S^{TP}}{f}_i(P_{i,t}^G)（i\in S^G）$
对剩余的（13）-（14）式，为碳交易约束，对应：
$g_0(x_0)=[{\pi }_b,-{\pi }_b]x_0$, 得到：定义 X 0 = { x 0 ∣ ( 13 ) − ( 14 ) } \Chi_0=\{x_0|(13)-(14)\} X0​={x0​∣(13)−(14)}

在以上简化在以上简化的基础上，令${\psi }_0(x_0)=[{\tilde{h}}_0(x_0)-b_E;{\tilde{M}}_0]i\in S^B$

在以上简化的基础上，

1. 令$i=0$时，${\psi }_0(x_0)=[{\tilde{h}}_0(x_0)-b_E;{\tilde{M}}_0],A_0(x_0)={\tilde{A}}_0$;
2. $i\in S^B$时，$A_i(x_i)=[\widetilde{A_i}{x}_i-b_{PF}^i],{\psi }_i(x_i)=[{\tilde{h}}_i(x_i);{\tilde{M}}_i{x}_i-b_{BR}^i]$
3. 另外：${\sum }_{i\in S^B}{b}_{BR}^i=b_{BR},{\sum }_{i\in S^B}{b}_{BF}^i=b_{BF}$

综合上述简化，DC-DOPF-CRT模型显示如下：
Minimize ∑ i ∈ S B ∩ { 0 } g i ( x i ) subject to ∑ i ∈ S B ∩ { 0 } ψ i ( x i ) ≤ 0 ∑ i ∈ S B ∩ { 0 } A i ( x i ) = 0 x i ∈ X i i ∈ { 0 } ∩ S B (20) \begin{aligned} \text{Minimize} \quad &\sum_{i\in S^B \cap\{0\}}g_i(x_i) \\ \text{subject to} \quad & \sum_{i\in S^B \cap\{0\}} \psi_i({x_i})\leq0\\ & \sum_{i\in S^B \cap\{0\}} A_i({x_i})=0\\ & x_i\in\Chi_i\quad i\in\{0\}\cap S^B \end{aligned}\tag{20} Minimizesubject to​i∈SB∩{0}∑​gi​(xi​)i∈SB∩{0}∑​ψi​(xi​)≤0i∈SB∩{0}∑​Ai​(xi​)=0xi​∈Xi​i∈{0}∩SB​(20)

对偶ADMM求解

首先原来整个 S B ∪ { 0 } S^B\cup\{0\} SB∪{0}集合分块。设 ∣ π ∣ = n + 1 , π = { A 0 , A 1 , … , A n } |\pi|=n+1, \pi = \{\mathcal{A}_0, \mathcal{A}_1, \dots, \mathcal{A}_n\} ∣π∣=n+1,π={A0​,A1​,…,An​}分别是集合 S B ∪ { 0 } S^B\cup\{0\} SB∪{0}的一个分块，即满足：

1. $\varnothing \notin \pi$, ⋃ A j ∈ π , j ∈ 0 , … , n A j = S B ∪ { 0 } \bigcup_{\mathcal{A}_j\in\pi, j\in 0,\dots, n}\mathcal{A}_j = S^B\cup\{0\} ⋃Aj​∈π,j∈0,…,n​Aj​=SB∪{0}
2. $(\forall A,B\in \pi )A\ne B\Rightarrow A\cap B=\varnothing$
   对这其中每一个分块均能够写出一个Lagrange对偶函数：
   L ( λ , ν ) = ∑ j = 0 n L j ( λ , ν ) L j ( λ , ν ) = inf ⁡ x i ∈ X i , i ∈ A j { ∑ i ∈ A j g i ( x i ) + ⟨ ∑ i ∈ A j ψ i ( x i ) , λ ⟩ + ⟨ ∑ i ∈ A j A i ( x i ) , ν ⟩ } (21) \begin{align*} L(\lambda, \nu) &= \sum^n_{j=0}L_j(\lambda, \nu)\\ L_j(\lambda, \nu) &= \inf_{x_i\in\Chi_i, i\in A_j}\{\sum_{i\in A_j}g_i(x_i)+\langle\sum_{i\in \mathcal{A}_j}\psi_i(x_i), \lambda\rangle+ \langle\sum_{i\in \mathcal{A}_j}\mathcal{A}_i(x_i), \nu\rangle\} \end{align*}\tag{21} L(λ,ν)Lj​(λ,ν)​=j=0∑n​Lj​(λ,ν)=xi​∈Xi​,i∈Aj​inf​{i∈Aj​∑​gi​(xi​)+⟨i∈Aj​∑​ψi​(xi​),λ⟩+⟨i∈Aj​∑​Ai​(xi​),ν⟩}​(21)
   其中，$\lambda ,\nu$为相应项分别对应相应约束的对偶乘子，且$\lambda \ge 0$。
   此时，令$y=⟨\lambda ,\nu ⟩$，那么（20）式子的对偶问题为：
   $$\begin{array}{r}\max L(y)\end{array}\text{\hspace{1em}(22)}$$

Minimized { − ∑ j = 0 n L j ( x j ) } subject to y − z j = 0 j = 0 , … , n (23) \begin{aligned} \text{Minimized} \quad & \{-\sum_{j=0}^n L_j(x_j)\}\\ \text{subject to} \quad & y-z_j= 0\quad j = 0,\dots,n \end{aligned}\tag{23} Minimizedsubject to​{−j=0∑n​Lj​(xj​)}y−zj​=0j=0,…,n​(23)

应用ADMM解决这一问题：
$$\begin{array}{r}{y}^{k+1}:=\underset{y}{\mathrm{argmin}}\{⟨\sum _{j=0}^n{p}_j^k,y⟩+\frac{\rho }{2}\sum _{j=0}^n{\Vert y-z_j^k\Vert }^2\}\end{array}\text{\hspace{1em}(24)}$$
z j k + 1 : = a r g m i n z j { − L j ( z j ) − ⟨ p j k , z j ⟩ + ρ 2 ∑ j = 0 n ∥ y k + 1 − z j ∥ 2 } (25) \begin{align} z^{k+1}_j:=\mathop{argmin}\limits_{z_j}\{-L_j(z_j)-\langle p_j^k, z_j\rangle+\frac{\rho}{2}\sum_{j=0}^n{\Vert y^{k+1}-z_j\Vert}^2 \} \end{align}\tag{25} zjk+1​:=zj​argmin​{−Lj​(zj​)−⟨pjk​,zj​⟩+2ρ​j=0∑n​∥yk+1−zj​∥2}​(25)
$$\begin{array}{r}{p}_j^{k+1}:=p_j^k+\rho (y^{k+1}-z_j^{k+1})\end{array}\text{\hspace{1em}(26)}$$

更新时候，求解（24），得到求解$y_{k+1}$更新值的计算式：
$$\begin{array}{r}{y}^{k+1}=\frac{1}{n+1}\sum _{n+1}^n{z}_j^k+-\frac{1}{(n+1)\rho }\sum _{j=0}^n{p}_j^k\end{array}\text{\hspace{1em}(26)}$$
联立（21）带入（25）,得到求解$z_j^{k+1}$更新值的计算式：（固定$x_i\in {\mathrm{X}}_i,i\in {\mathcal{A}}_j$的情况）
z j = [ λ z j ν z j ] = [ m a x { 0 , λ y k + 1 + 1 ρ ( ∑ i ∈ A j ψ i ( x i ) + λ p j k ) } ν y k + 1 + 1 ρ ( ∑ i ∈ A j A i ( x i ) + λ p j k ) ] (29) \begin{align} z_j = \left[ \begin{matrix} \lambda_{zj}\\ \nu_{zj} \end{matrix} \right] = \left[ \begin{matrix} max\{0, \lambda_{y^{k+1}}+\frac{1}{\rho}(\sum_{i\in \mathcal{A}_j \psi_i(x_i)+\lambda_{p_j^k}})\}\\ \nu_{y^{k+1}}+\frac{1}{\rho}(\sum_{i\in \mathcal{A}_j}\mathcal{A}_i(x_i)+\lambda_{p_j^k}) \end{matrix} \right] \end{align}\tag{29} zj​=[λzj​νzj​​]=[max{0,λyk+1​+ρ1​(∑i∈Aj​ψi​(xi​)+λpjk​​​)}νyk+1​+ρ1​(∑i∈Aj​​Ai​(xi​)+λpjk​​)​]​(29)
在所有的子集${\mathcal{A}}_j$中，取最优的子集 { x i ∈ X i , i ∈ A j } \{x_i\in \Chi_i, i\in\mathcal{A}_j\} {xi​∈Xi​,i∈Aj​}:
{ x i ∗ , i ∈ A j } : = a r g m i n x i ∈ χ i , i ∈ A j { ∑ i ∈ A j g i ( x i ) + ρ 2 ∥ [ m a x { 0 , λ y k + 1 + 1 ρ ( ∑ i ∈ A j ψ i ( x i ) + λ p j k ) } ν y k + 1 + 1 ρ ( ∑ i ∈ A j A i ( x i ) + λ p j k ) ] ∥ 2 } (30) \begin{align} \{x_i^*,i\in\mathcal{A}_j \} :=\mathop{argmin}\limits_{x_i \in \chi_i, i\in\mathcal{A}_j}\left\{\sum_{i\in\mathcal{A}_j}g_i(x_i)+\frac{\rho}{2}{\left\Vert\left[ \begin{matrix} max\{0, \lambda_{y^{k+1}}+\frac{1}{\rho}(\sum_{i\in A_j \psi_i(x_i)+\lambda_{p_j^k}})\}\\ \nu_{y^{k+1}}+\frac{1}{\rho}(\sum_{i\in A_j}A_i(x_i)+\lambda_{p_j^k}) \end{matrix} \right]\right\Vert}^2\right\} \end{align}\tag{30} {xi∗​,i∈Aj​}:=xi​∈χi​,i∈Aj​argmin​⎩ ⎨ ⎧​i∈Aj​∑​gi​(xi​)+2ρ​ ​[max{0,λyk+1​+ρ1​(∑i∈Aj​ψi​(xi​)+λpjk​​​)}νyk+1​+ρ1​(∑i∈Aj​​Ai​(xi​)+λpjk​​)​] ​2⎭ ⎬ ⎫​​(30)
令 τ = m a x { 0 , λ y k + 1 + 1 ρ ( ∑ i ∈ A j ψ i ( x i ) + λ p j k ) } \tau = max\{0, \lambda_{y^{k+1}}+\frac{1}{\rho}(\sum_{i\in A_j}\psi_i(x_i)+\lambda_{p_j^k})\} τ=max{0,λyk+1​+ρ1​(∑i∈Aj​​ψi​(xi​)+λpjk​​)}
解最小值：
$$\begin{array}{rl}\min & \sum _{i\in {\mathcal{A}}_j}{g}_i(x_i)+\frac{\rho }{2}{\Vert \left[\begin{array}{c}\tau \\ {\nu }_{y^{k+1}}+\frac{1}{\rho }({\sum }_{i\in {\mathcal{A}}_j}{A}_i(x_i)+{\lambda }_{p_j^k})\end{array}\right]\Vert }^2\\ s.t.& \tau \ge 0\\ & \tau \ge {\lambda }_{y^{k+1}}+\frac{1}{\rho }(\sum _{i\in A_j}{\psi }_i(x_i)+{\lambda }_{p_j^k})\\ & x_i\in {\chi }_i,i\in {\mathcal{A}}_j\end{array}\text{\hspace{1em}(31)}$$
综合上述，可得算法1（DC-ADMM-P）的流程

• 初始化，有$z_j^0$, $p_j^0$，$M>0$, $\rho >0$, ${ϵ}^{feasible}>0$已知

  for $k=0,\dots ,M$
  - 根据前述（27）式获得$y^{k+1}$

  • for（并行求解）对应每一个子集$j=0,\dots ,n$
    • 根据（31）式获得 { x i ∗ , i ∈ A j } \{x_i^*, i\in\mathcal{A}_j\} {xi∗​,i∈Aj​}
    • 根据（29）式获得$z_j^{k+1}$
    • 根据（26）式获得$p_j^{k+1}$
  • end
  • if $ϵ:=\Vert y^{k+1}-y^k{\Vert }_{\infty }<{ϵ}^{feasible}$,退出循环

• end