电力系统分析—潮流计算代码Python编程练习（基于极坐标形式的常规牛拉法）

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

一、前言

博主是一名电气专业在读学生，在学习《电力系统分析》这门课程的过程里尝试利用python编程实现潮流计算，在编程过程中遇到许多难题并发现代码中的不足之处，对代码也进行了多次改进，最终得到了一份比较满意的代码。写下此贴以记录代码改进过程。

---

提示：以下是本篇文章正文内容，下面案例可供参考

二、潮流计算是什么？

潮流计算是电力学名词，指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电力网中的分布。潮流计算是根据给定的电网结构、参数和发电机、负荷等元件的运行条件，确定电力系统各部分稳态运行状态参数的计算。通常给定的运行条件有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。待求的运行状态参量包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。
潮流计算是电力系统分析最基本的计算。事实上，潮流计算还是网损计算、静态安全分析、暂态稳定计算、小干扰静态稳定计算、短路计算、静态和动态等值计算的基础。因此对潮流计算收敛速度以及收敛精度的要求不言而喻。
本文潮流计算程序的实现是基于极坐标形式的牛顿—拉夫逊法，具体原理只在后文作简单介绍，有兴趣的读者可自行了解并查阅相关书籍（如《电力系统分析（下）》何仰赞第四版）

三、python实现过程

1.极坐标下相关方程

采用极坐标时，节点电压可表示为
$$\dot{U}=U_i\angle {\delta }_i=U_i(cos{\delta }_i+jsin{\delta }_i)$$
每一个PQ节点和PV节点有功功率不平衡量方程将写成
$$\Delta P_i=P_{is}-P_i=P_{is}-U_i{\Sigma }_{j=1}^n{U}_j(G_{ij}cos{\delta }_{ij}+B_{ij}sin{\delta }_{ij})=0$$
而对于每一个PQ节点还可以再列写一个无功功率不平衡量方程式
$$\Delta Q_i=Q_{is}-Q_i=Q_{is}-U_i{\Sigma }_{j=1}^n{U}_j(G_{ij}sin{\delta }_{ij}-B_{ij}cos{\delta }_{ij})=0$$
极坐标形式下的修正方程式为
$$\left[\begin{array}{c}\Delta P\\ \Delta Q\end{array}\right]=-\left[\begin{array}{cc}H& N\\ K& L\end{array}\right]\left[\begin{array}{c}\Delta \delta \\ U_{D2}^{-1}\Delta U\end{array}\right]$$
其中当$i\ne j$时，有
$$\begin{array}{c}{H}_{ij}=-U_i{U}_j(G_{ij}sin{\delta }_{ij}-B_{ij}cos{\delta }_{ij})\\ N_{ij}=-U_i{U}_j(G_{ij}cos{\delta }_{ij}+B_{ij}sin{\delta }_{ij})\\ K_{ij}=U_i{U}_j(G_{ij}cos{\delta }_{ij}+B_{ij}sin{\delta }_{ij})\\ L_{ij}=-U_i{U}_j(G_{ij}sin{\delta }_{ij}-B_{ij}cos{\delta }_{ij})\end{array}\}$$
当$i=j$时
$$\begin{array}{c}{H}_{ij}=U_i^2{B}_{ii}+Q\\ N_{ij}=-U_i{G}_{ii}-P_i\\ K_{ij}=U_i^2{G}_{ii}-P_i\\ L_{ij}=U_i^2{B}_{ii}-Q_i\end{array}\}$$
本文将基于以上主要原理方程式编写潮流计算程序。

2.具体实现过程

本次潮流计算程序设计是基于IEEE30节点系统实现的，文件采用IEEE标准数据格式，如下图所示

潮流计算程序流程大致可以分为以下几个部分：

1. 读取文件和设定初值
2. 生成节点导纳矩阵
3. 计算不平衡量
4. 判断是否符合收敛条件（如果符合则结束程序，反之则继续往下执行）
5. 计算雅可比矩阵
6. 计算修正方程式，得到修正量
7. 修正各节点电压和相角
8. 返回第三步
   其中最核心的部分在于计算雅可比矩阵，其计算正确性和计算速度对整个程序的收敛速度和收敛精度起决定性的影响。
   在最初编写雅可比计算程序时，笔者只是利用python的for循环和math库一个一个计算雅可比矩阵中的各个元素，如下所示

for i in range(total_count): #求H
        for j in range(total_count):
            if i != j:
                H[i][j]=-U[P_row[i]-1]*U[P_row[j]-1]*(NAM[P_row[i]-1][P_row[j]-1].real*sin(delta[P_row[i]-1]-delta[P_row[j]-1])
                                                      -NAM[P_row[i]-1][P_row[j]-1].imag*cos(delta[P_row[i]-1]-delta[P_row[j]-1]))
            else:
                H[i][j]=U[P_row[i]-1]**2*NAM[P_row[i]-1][P_row[i]-1].imag+Q_i.get(P_row[i],0)

包括在计算不平衡量时也是利用了类似的方法

for i in range(total_count): #求delta_P
        sigma=0
        for j in range(num1):
            sigma+=U[P_row[i]-1]*U[j]*(NAM[P_row[i]-1][j].real*cos(delta[P_row[i]-1]-delta[j])+NAM[P_row[i]-1][j].imag*sin(delta[P_row[i]-1]-delta[j]))
        P_i[P_row[i]]=sigma
        delta_P[i]=(Init_LoadMW[i]/100-sigma)

但当我运行时，代码花了将近半分钟的时间才跑完全程，且结果也与原文件给出的最终电压和最终相角差异较大（部分节点得出的结果相差甚至超过0.5）.
事实上，这种for循环+行列单独索引的方法虽然操作简单，但是运行效率极低，几乎每迭代一次都需要几百毫秒（这在实际工程是绝对不满足要求的，因为电网计算往往是实时计算甚至超实时计算），同时大规模地采用这种方法不仅极大地降低了程序运行速度，而且会使整个代码看起来十分混乱，还有就是对于像节点导纳矩阵、雅可比矩阵等的概念必须有一个较为通透的了解，如果在某个细节上的理解存在错误，那么很可能导致最终的程序出现收敛次数过多，甚至无法收敛的问题，因此我们开始考虑如何用更为简洁且高效的方法来实现雅可比矩阵的计算。
首先，我们分析了潮流计算的特点以及可能用到的变量，并考虑的代码的整洁性和可读性，我们基于python面向对象编程引入了三个类: net、bus以及branch，分别用来存储潮流计算的函数（包括计算雅可比矩阵，计算不平衡量等）、节点的相关属性以及支路的相关属性，其中bus和branch如下所示

class bus:
    def __init__(self, num, Type, LoadMW_In, LoadMW_Out, LoadMVAR_In,
                 LoadMVAR_Out, U, delta, shunt_conductance, shunt_susceptance):
        self.num = num  #编号
        self.U = U  #电压幅值
        self.Type = Type #节点类型
        self.delta = delta  #电压相角
        self.LoadMW_In = LoadMW_In  #注入有功功率
        self.LoadMW_Out = LoadMW_Out  #输出有功功率
        self.LoadMVAR_In = LoadMVAR_In  #注入无功功率
        self.LoadMVAR_Out = LoadMVAR_Out  #输出无功功率
        self.shunt_conductance = shunt_conductance  #接地电阻
        self.shunt_susceptance = shunt_susceptance  #接地电抗

    def print_bus(self):
        print('编号', self.num)
        print('电压幅值', self.U)
        print('电压相角', self.delta)
        print('注入有功功率', self.LoadMW_In)
        print('输出有功功率', self.LoadMW_Out)
        print('注入无功功率', self.LoadMVAR_In)
        print('输出无功功率', self.LoadMVAR_Out)
        print('接地电阻', self.shunt_conductance)
        print('接地电抗', self.shunt_susceptance)

class branch:
    def __init__(self, Tap_bus, Z_bus, branch_resistance, branch_reactance,
                 line_charging, ratio):
        self.Tap_bus = Tap_bus
        self.Z_bus = Z_bus
        self.branch_resistance = branch_resistance
        self.branch_reactance = branch_reactance
        self.line_charging = line_charging
        self.ratio = ratio

为避免篇幅过长，net的具体代码便不在这展示了。

接下来主要重点讲述一下我们改进后的雅可比矩阵计算实现方法

首先我们以有功功率计算公式为例
$$P_i=U_i{\Sigma }_{j=1}^n{U}_j(G_{ij}cos{\delta }_{ij}+B_{ij}sin{\delta }_{ij})$$
观察表达式可以发现我们一共要用到3个量：节点电压幅值$U$，节点电压相角$\delta$以及节点导纳矩阵相应位置的元素实部和虚部

首先我们将$$U=\left[U_1,U_2,U_3\dots U_{30}\right]$$也就是一个30x1的列向量(以本文所采用的30节点系统为例)变成1x30的行向量 $$U^T=\left[\begin{array}{c}{U}_1\\ U_2\\ U_3\\ ⋮\\ U_{30}\end{array}\right]$$
然后利用相关函数将其复制转化成30x30的矩阵，即
$$U_{copy}=\left[\begin{array}{cccc}{U}_1& U_2& \dots & U_{30}\\ U_1& U_2& \dots & U_{30}\\ U_1& U_2& \dots & U_{30}\\ ⋮& ⋮& ⋮& ⋮\\ U_1& U_2& \dots & U_{30}\end{array}\right]$$
之后将这个30x30的矩阵每一行都乘以对应行号的电压，也就是第一行每个元素都乘以$U_1$，第二行每个元素都乘以$U_2$，以此类推，也就相当于等式中的$U_i\ast U_j$，得到矩阵如下
$$U_i{U}_j=\left[\begin{array}{cccc}{U}_1{U}_1& U_1{U}_2& \dots & U_1{U}_{30}\\ U_2{U}_1& U_2{U}_2& \dots & U_2{U}_{30}\\ U_3{U}_1& U_3{U}_2& \dots & U_3{U}_{30}\\ ⋮& ⋮& ⋮& ⋮\\ U_{30}{U}_1& U_{30}{U}_2& \dots & U_{30}{U}_{30}\end{array}\right]$$
这一过程可以通过上述行向量和列向量，利用Numpy库中的matmul函数直接实现。

然后将$\delta =[{\delta }_1,{\delta }_2,\dots ,{\delta }_{30}]$也做类似的复制处理，得到
$${\delta }_{copy}=\left[\begin{array}{cccc}{\delta }_1& {\delta }_2& \dots & {\delta }_{30}\\ {\delta }_1& {\delta }_2& \dots & {\delta }_{30}\\ {\delta }_1& {\delta }_2& \dots & {\delta }_{30}\\ ⋮& ⋮& ⋮& ⋮\\ {\delta }_1& {\delta }_2& \dots & {\delta }_{30}\end{array}\right]$$
并将复制处理后的30x30格式的$\delta$矩阵减去原来30*1格式的$\delta$列向量，即
$${\delta }_{ij}=\left[\begin{array}{cccc}{\delta }_1-{\delta }_1& {\delta }_2-{\delta }_1& \dots & {\delta }_{30}-{\delta }_1\\ {\delta }_1-{\delta }_2& {\delta }_2-{\delta }_2& \dots & {\delta }_{30}-{\delta }_2\\ {\delta }_1-{\delta }_3& {\delta }_2-{\delta }_3& \dots & {\delta }_{30}-{\delta }_3\\ ⋮& ⋮& ⋮& ⋮\\ {\delta }_1-{\delta }_{30}& {\delta }_2-{\delta }_{30}& \dots & {\delta }_{30}-{\delta }_{30}\end{array}\right]$$
这样就得到了${\delta }_{ij}$矩阵。最后利用python中的复数操作方法取出节点导纳矩阵中的实数部分和虚数部分，如下所示

# 用对非平衡节点P_i的定义计算所有节点的P_i(包括平衡节点)，方便后期使用索引
U_ij = np.matmul(self.U, self.U.T)  # 对应位置的U_i*U_j
delta_ij = self.delta - self.delta.T  # 对应位置的delta_i-delta_j
G_ij, B_ij = self.NAM.real, self.NAM.imag  # 对应位置的G_ij,B_ij,NAM为节点导纳矩阵

把准备好的变量相乘后相加，得到的矩阵中每个元素的表达式都是$V_i{V}_j\ast (G_{ij}cos{\delta }_{ij}+B_{ij}sin{\delta }_{ij})$,其中的$i$和$j$分别对应元素所在的行和列，最后用求和操作即得到一个30x1的矩阵，每行元素对应每个节点的有功功率。而无功功率的操作基本与有功一致，不再赘述。
基于这种想法，我们再次观察前文所述的修正方程式以及雅可比矩阵中各子块元素的表达式
仔细观察可以发现，常规牛拉法的雅可比矩阵有四个分块矩阵，四个分块矩阵中元素的表达式也都用到3个量:$U$，$\delta$以及节点导纳矩阵相应位置的元素。再观察一下我们还能发现一个特点: 在求取雅可比矩阵子块元素时，$i=j$和$i\ne j$的表达式其实是可以合并的，就拿H子块来说，当$i=j$时，如果认为$Q_i=0$，那么无论是$i=j$还是$i\ne j$，二者表达式完全一致，这就意味着我们完全可以先按照之前处理Pi和Qi那样用类似的方法去处理雅可比矩阵的生成，生成之后再单独索引提取出对角线元素并给其加上对应的$P_i$和$Q_i$即可。这样一来不仅简化了代码，而且用到了Numpy的高级索引和内置函数，效率相比于python内置函数提高了上百倍，具体实现过程如下（以H子块为例）

 def NL_Jacob(self):
        #生成雅可比矩阵
        U = np.vstack([self.U[self.P_row - 1], self.U[self.Q_row - 1]])
        U = np.repeat(U.T, len(self.P_row) + len(self.Q_row), axis=0)
        UU = U * np.vstack([self.U[self.P_row - 1], self.U[self.Q_row - 1]])

        delta = np.vstack(
            [self.delta[self.P_row - 1], self.delta[self.Q_row - 1]])
        delta = np.repeat(delta.T, len(self.P_row) + len(self.Q_row), axis=0)
        delta = -(delta - np.vstack(
            [self.delta[self.P_row - 1], self.delta[self.Q_row - 1]]))
        
        self.H = -UU[:len(self.P_row), :len(self.P_row)] * (
            self.NAM[self.P_row - 1][:, self.P_row - 1].real *
            np.sin(delta[:len(self.P_row), :len(self.P_row)]) -
            self.NAM[self.P_row - 1][:, self.P_row - 1].imag *
            np.cos(delta[:len(self.P_row), :len(self.P_row)]))
        self.H[np.arange(len(self.P_row)),
               np.arange(len(self.P_row))] += self.Q_i[self.P_row - 1]

其中P_row是一个python列表，记录储存PQ节点和PV节点在文件中所处位置，而Q_row则是用来储存PQ节点在文件中所处位置。（因为当时想的是在极坐标情况下每一个PQ节点和PV节点都要计算有功不平衡量，而PQ节点还要额外计算无功不平衡量，因此直接这样命名了，现在想来好像有点草率了😂）
需要注意的是，虽然基本原理和$P_i$和$Q_i$的计算无异，但由于整个矩阵是53x53的矩阵（同样以本文所采用的30节点系统为例），所以我们事先生成好了53x53的$U\ast U$的电压矩阵和53*53的${\delta }_{ij}$矩阵，之后利用索引分别提取其中的部分分块来生成对应的分块矩阵（H,N,K,L），同时在计算完一个分块矩阵的元素后再单独提取出其中的对角线元素，加上对应的$P_i$和$Q_i$才算完成。
至此我们便完成了基于Python的潮流计算程序，实测程序运行结果如下

各节点电压幅值为:
[[1.06       1.045      1.0214564  1.01264102 1.01       1.01059014
  1.00257545 1.01       1.05122907 1.04534075 1.082      1.05585053
  1.071      1.04110288 1.03679673 1.0439088  1.03985247 1.02766251
  1.0253907  1.02960022 1.03285831 1.03336256 1.02655487 1.02128763
  1.01691309 0.99922814 1.02273634 1.00706292 1.00288636 0.99140546]]
 
各节点电压相角为:
[[  0.          -5.38170753  -7.52177212  -9.27052684 -14.16163305
  -11.07717343 -12.8706939  -11.82113871 -14.41204162 -16.07882848
  -14.41204162 -15.52840747 -15.52840747 -16.39520983 -16.46241103
  -16.0267584  -16.27755164 -17.02264801 -17.16368539 -16.94979429
  -16.52688949 -16.51439444 -16.79921577 -16.90215113 -16.42226392
  -16.84227857 -15.86550164 -11.70519499 -17.09670324 -17.98046186]]
 
程序运行时间为0.021942138671875s

与文件给出的各节点电压幅值和相角对比图如下


细心的读者可以发现其实精度仍然不算完美，这是因为笔者水平有限，只能写出最基本的潮流计算程序，如果需要更精确的结果，则应该还需要考虑负荷静特性等问题（也许？大概？😂）,但收敛速度已经达到了较高水平。

四、 总结

潮流计算程序的编写花费了很大的精力，但个人认为这样一次编写程序的经历能极大提高自己对电力系统潮流计算的理解，而且也是一个很好的锻炼机会，培养自己动手实践的能力，以及在编写过程中遇到的诸如程序本身，各个概念的理解等问题，这些都是极为宝贵的经验。在编写过程中也得到了学院师兄以及各位老师同学的很多帮助，在此致以衷心的感谢。也希望这篇文章能给正忙于编写潮流程序相关大作业的读者同学们一点启发。
最后，人生苦短，我用python。👏