摘要—本文提出了一种基于几何的空间分支策略（ spatial branching strategy），用于解决集成电力-燃气系统中的圆锥非凸方程（ conic non-convex equations）。所提出的策略嵌入在空间分支定界算法中，以求解最优电力-燃气流问题的精确解。多个测试案例的数值结果验证了该方法的有效性和优越性。

关键词—分析几何学，多能源系统，非凸优化，空间分支。

I. INTRODUCTION

随着电力和燃气系统之间的相互依赖性不断增加，越来越多的关注被放在了最优电力-燃气流（OEGF）问题上[1]，在该问题中，最优电力流和最优燃气流模型是耦合的[2]。求解OEGF问题的精确解对于集成系统的安全性至关重要。然而，由于电力流和燃气流模型的非凸性，求解具有保证精度的OEGF问题具有挑战性。

二阶圆锥松弛（Second-order conic relaxation）是一种常用的方法，用于使OEGF问题具有凸性[3]，[4]，[5]，[6]，[7]。在电力侧（特别是配电系统）上，圆锥松弛支路流模型已被广泛应用[3]。至于燃气侧，稳态[4]，[5]和动态[6]燃气流模型中的非凸二次项也可以通过二阶锥规划（SOCP）进行松弛。然而，圆锥松弛问题的解可能不可行（紧约束）。即使通过实施非紧解的圆锥-凹程序（CCP）[7]可以恢复解，经过紧化的解也可能是次优解。

目前，空间分支定界（spatial branch-and-bound，S-B&B）是求解非凸模型的最佳可用算法之一[1]，并嵌入到多种通用求解器中，如SCIP和Gurobi。然而，求解非凸OEGF问题的计算时间过长，即使提高了精度，对于实际应用来说也是不可接受的。主要原因是，上述求解器的传统S-B&B算法仅为特定类型的非凸性设计，即双线性等式约束（$z=x\cdot y$，其中$(x,y,z)\in {\mathbb{R}}^3$）。然而，对于圆锥非凸性，形式为$z=\sqrt{x^2+y^2}$，其中$(x,y,z)\in {\mathbb{R}}^3$，并且$z\ge 0$，这在OEGF模型中更为常见，求解器必须将每个圆锥方程重新表述为三个双线性方程，导致非凸性的规模增加了三倍[1]。

为此，提出了一种新的空间分支策略，专门用于解决OEGF问题中的圆锥非凸性，具有更高的精度和效率。

II. MATHEMATICAL MODELING

鉴于本文的主要贡献集中在解决方法上，建模部分进行了简化。这里直接采用了[3]中提出的OEGF模型（其原始的非凸形式），并将支路流模型与准动态燃气流模型结合。

OEGF模型的广义矩阵形式称为$\mathcal{P}$，其表达式如下：

$$\begin{array}{c}\mathcal{P}:\min {\mathbf{c}}^T\mathbf{x}\text{s.t.}\mathbf{A}\mathbf{x}=\mathbf{e},\mathbf{B}\mathbf{x}\le \mathbf{f},{\left[\mathbf{x}\right]}_i\in \mathbb{Z}\\ {\left[\mathbf{x}\right]}_{r_n}=\sqrt{{\left[\mathbf{x}\right]}_{p_n}^2+{\left[\mathbf{x}\right]}_{q_n}^2},\forall n\in \Omega \end{array}\text{\hspace{1em}(1)(2)}$$

其中，$\mathbf{x}$ 是一个变量向量（其中部分为整数），$\mathbf{c},\mathbf{A},\mathbf{e},\mathbf{B},\mathbf{f}$是系数。表达式（1）表示OEGF模型的线性部分，而（2）表示非凸圆锥等式约束，包括视在功率流方程（见[3]中的（3）式）和Weymouth方程（见[3]中的（5）式）。此处引入了一个三元组索引，即$(p_n,q_n,r_n)$，对于每个$n\in \Omega$，其中$\Omega$是表示（2）中出现的变量的索引集。

注意，$\mathcal{P}$是为日调度集成电力-燃气系统而开发的，其中燃气流的动态不显著[2]。然而，所提议的S-B&B算法有可能被应用于短期运行调度，其中燃气流的瞬态由有限差分方法建模，且包含圆锥非凸性[6]。

III. PROPOSED SOLUTION METHODOLOGY

空间分支策略（spatial branching strategy）是使传统B&B算法能够处理非凸约束的核心技术。具体而言，我们研究了一个 relaxed problem $\mathcal{P}$，该问题可以通过将（2）中的“=”替换为“≥”来获得。放松后的 $\mathcal{P}$ 在此后称为 $\mathcal{RP}$。由于 $\mathcal{RP}$ 是一个混合整数SOCP模型，它可以通过传统的B&B（没有“空间”）算法，已经在求解器中嵌入（例如CPLEX），来求解。需要注意的是，求解器并不意识到非凸约束的存在，因此决策树某些节点的解可能是局部最优的，但违反（2）。为此，空间分支在这些节点上被采用，以创建两个子节点，这些子节点具有更严格的可行区域，从而将B&B更新为S-B&B。在本文中，提出了一种专门用于解决（2）中非凸性的空间分支策略，因此可以大大缓解传统双线性特定空间分支策略[1]带来的冗余问题。

为了简化表述，略作标记法上的处理，我们专注于一个通用的不可行解$(\hat{x},\hat{y},\hat{z})$，其满足：

$$\delta =\hat{z}-\sqrt{{\hat{x}}^2+{\hat{y}}^2}>0,l_x\le \hat{x}\le u_x,l_y\le \hat{y}\le u_y\text{\hspace{1em}(3)}$$

其中，$l_x,u_x,l_y,u_y$是当前节点上$x$和$y$的下界和上界。可行区域的示意图如图1(a)中的绿色区域所示，$P_0$是不可行解$(\hat{x},\hat{y},\hat{z})$。为了通过分支消除$P_0$，我们研究了蓝色圆锥曲面上最紧的三角形，其顶点由图1(b)中的$P_1,P_3$和$O$给出（注意，边$P_1{P}_3$是一个圆弧，circular arc）。基于球坐标系，$P_1$和$P_3$的径向距离$(\rho )$、极角$(\theta )$和方位角$(\varphi )$分别表示为 $({\rho }_1,{\theta }_1,{\phi }_1)$，$({\rho }_3,{\theta }_3,{\phi }_3)$

ρ 1 = ρ 3 = 2 [ max ⁡ ( l x 2 , u x 2 ) + max ⁡ ( l y 2 , u y 2 ) ] θ 1 = θ 3 = arctan ⁡ ( 1 ) = π / 4 φ 1 = min ⁡ w y ∈ { l y , u y } , w x ∈ { l x , u x } arctan ⁡ ( w y / w x ) φ 3 = max ⁡ w y ∈ { l y , u y } , w x ∈ { l x , u x } arctan ⁡ ( w y / w x ) \begin{align*} \rho _{1} =& \rho _{3} = \sqrt{2 \left[ \max (l_{x}^{2}, u_{x}^{2}) + \max (l_{y}^{2}, u_{y}^{2}) \right]} \tag{4} \\ \theta _{1} =& \theta _{3} = \arctan (1) = \pi / 4 \tag{5} \\ \varphi _{1} =& \min _{ w_{y} \in \left\lbrace l_{y}, u_{y} \right\rbrace, w_{x} \in \left\lbrace l_{x}, u_{x} \right\rbrace } \arctan \left(w_{y} / w_{x}\right) \tag{6} \\ \varphi _{3} =& \max _{ w_{y} \in \left\lbrace l_{y}, u_{y} \right\rbrace, w_{x} \in \left\lbrace l_{x}, u_{x} \right\rbrace } \arctan \left(w_{y} / w_{x}\right) \tag{7} \end{align*} ρ1​=θ1​=φ1​=φ3​=​ρ3​=2[max(lx2​,ux2​)+max(ly2​,uy2​)] ​θ3​=arctan(1)=π/4wy​∈{ly​,uy​},wx​∈{lx​,ux​}min​arctan(wy​/wx​)wy​∈{ly​,uy​},wx​∈{lx​,ux​}max​arctan(wy​/wx​)​(4)(5)(6)(7)​

接下来，在$P_1$和$P_3$之间的弧上确定中点$P_2$，这样就可以将两个三角形切割平面$△O{P}_1{P}_2$和$△O{P}_2{P}_3$分别施加到下/上子节点。为了确保$P_0$在两个节点都被切割掉，$P_2$的一个备选选择为：

$$\begin{array}{c}{\rho }_2={\rho }_1={\rho }_3,{\theta }_2=\pi /4,{\phi }_2=\arctan (\hat{y}/\hat{x})\end{array}\text{\hspace{1em}(8)}$$

并在图1(b)中标注。注意，边$P_1{P}_2$和$P_2{P}_3$是直线，而不是弧线。

这里，$P_1,P_2,P_3$的球坐标$(\rho ,\theta ,\phi )$被转换为相应的笛卡尔坐标$(x,y,z)$，其转换关系为：

{ x { ⋅ } = ρ { ⋅ } cos ⁡ φ { ⋅ } sin ⁡ θ { ⋅ } ; y { ⋅ } = ρ { ⋅ } sin ⁡ φ { ⋅ } sin ⁡ θ { ⋅ } ; z { ⋅ } = ρ { ⋅ } cos ⁡ θ { ⋅ } . ∀ { ⋅ } ∈ { 1 , 2 , 3 } \begin{equation*} {\begin{cases}x_{\lbrace \cdot \rbrace } = \rho _{\lbrace \cdot \rbrace } \cos \varphi _{\lbrace \cdot \rbrace } \sin \theta _{\lbrace \cdot \rbrace }; \\ y_{\lbrace \cdot \rbrace } = \rho _{\lbrace \cdot \rbrace } \, \sin \varphi _{\lbrace \cdot \rbrace } \sin \theta _{\lbrace \cdot \rbrace }; \\ z_{\lbrace \cdot \rbrace } = \rho _{\lbrace \cdot \rbrace } \cos \theta _{\lbrace \cdot \rbrace }. \\ \end{cases}} \forall \lbrace \cdot \rbrace \in \lbrace 1, 2, 3\rbrace \tag{9} \end{equation*} ⎩ ⎨ ⎧​x{⋅}​=ρ{⋅}​cosφ{⋅}​sinθ{⋅}​;y{⋅}​=ρ{⋅}​sinφ{⋅}​sinθ{⋅}​;z{⋅}​=ρ{⋅}​cosθ{⋅}​.​∀{⋅}∈{1,2,3}​(9)​

然后，$△O{P}_1{P}_2$和$△O{P}_2{P}_3$的切割平面，如图1 $(\mathrm{c})$，(d)所示，以行列式（线性约束）形式表示为：

$$\begin{array}{r}{f}_{\Delta O{P}_1{P}_2}(x,y,z)=\left|\begin{array}{ccc}x& y& z\\ x_1-0& y_1-0& z_1-0\\ x_2-0& y_2-0& z_2-0\end{array}\right|\le 0\\ f_{\Delta O{P}_2{P}_3}(x,y,z)=\left|\begin{array}{ccc}x& y& z\\ x_2-0& y_2-0& z_2-0\\ x_3-0& y_3-0& z_3-0\end{array}\right|\le 0\end{array}\text{\hspace{1em}(10)(11)}$$

这进一步推导出空间分支策略为：

$$[f_{△O{P}_1{P}_2}(x,y,z)\le 0]\vee [f_{△O{P}_2{P}_3}(x,y,z)\le 0]\text{\hspace{1em}(12)}$$

其中，(12)式的左右两边可以分别施加在当前节点上，创建下/上子节点。

命题1：不可行解$P_0$可以通过在两个子节点中应用（12）来消除。
证明：由于（12）左右两边的证明过程相似，本文以向下分支（10）为例。对于一个不可行解$(\hat{x},\hat{y},\hat{z})$，假设其满足矛盾条件：
$$f_{△O{P}_1{P}_2}(\hat{x},\hat{y},\hat{z})\le 0\text{\hspace{1em}(13)}$$
将（8）和（9）代入（13）式，我们得到：
$$\begin{array}{rl}{f}_{\Delta O{P}_1{P}_2}(\hat{x},\hat{y},\hat{z})=& \frac{1}{2}{\rho }_1{\rho }_2(\sin {\phi }_1-\sin {\phi }_2)\hat{x}\\ & -\frac{1}{2}{\rho }_1{\rho }_2(\cos {\phi }_1-\cos {\phi }_2)\hat{y}\\ & +\frac{1}{2}{\rho }_1{\rho }_2(\cos {\phi }_1\sin {\phi }_2\\ & -\sin {\phi }_1\cos {\phi }_2)\hat{z}\le 0\end{array}\text{\hspace{1em}(14)}$$
给定${\rho }_1={\rho }_2={\rho }_3$和${\phi }_2=\arctan (\hat{y}/\hat{x})$，（14）式可以进一步简化为：
$$\hat{x}\sin {\phi }_1-\hat{y}\cos {\phi }_1\le \frac{\hat{z}}{\sqrt{{\hat{x}}^2+{\hat{y}}^2}}(\hat{x}\sin {\phi }_1-\hat{y}\cos {\phi }_1)\text{\hspace{1em}(15)}$$
根据（6）式的定义，${\phi }_1\le \arctan (\hat{y}/\hat{x})$，这意味着$\hat{x}\sin {\phi }_1-\hat{y}\cos {\phi }_1\le 0$。因此，（15）式简化为：
$$\sqrt{{\hat{x}}^2+{\hat{y}}^2}\ge \hat{z}\text{\hspace{1em}(16)}$$
注意，式（16）违反了（3）中的前提条件，因此（13）式无效，从而完成证明。□

通过CPLEX实现所提出的空间分支策略的算法如算法1所示。该算法基于混合整数SOCP的解集$RP$，并与内置的分支与切割（B&C）算法一起工作。由于CPLEX无法识别（2）中的SOC方程，因此需要通过回调函数（callbacks）定制B&C算法，以应对这种情况[9]。每当一个解被找到并且对$RP$可行时，现有回调将被调用以判断该解是否对原始问题$P$可行，如果不可行，则拒绝该解。在CPLEX做出分支决策之前，它将调用分支回调请求建议。在这里，所提出的空间分支策略被用来指导CPLEX分支到那些对$RP$可行但违反了原问题约束的节点。该算法的详细框架在作者之前的研究中有所描述[1]，并且示例代码可以在[8]中找到。