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

一、前言

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

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

二、潮流计算是什么？

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

三、python实现过程

1.极坐标下相关方程

采用极坐标时，节点电压可表示为
U˙=Ui∠δi=Ui(cosδi+jsinδi)\dot{U}=U_i\angle\delta_i=U_i(cos\delta_i+jsin\delta_i)U˙=Ui∠δi=Ui(cosδi+jsinδi)
每一个PQ节点和PV节点有功功率不平衡量方程将写成
ΔPi=Pis−Pi=Pis−UiΣj=1nUj(Gijcosδij+Bijsinδij)=0\Delta P_i=P_{is}-P_i=P_{is}-U_i\Sigma^n_{j=1}U_j(G_{ij}cos\delta_{ij}+B_{ij}sin\delta_{ij})=0ΔPi=Pis−Pi=Pis−UiΣj=1nUj(Gijcosδij+Bijsinδij)=0
而对于每一个PQ节点还可以再列写一个无功功率不平衡量方程式
ΔQi=Qis−Qi=Qis−UiΣj=1nUj(Gijsinδij−Bijcosδij)=0\Delta Q_i=Q_{is}-Q_i=Q_{is}-U_i\Sigma^n_{j=1}U_j(G_{ij}sin\delta_{ij}-B_{ij}cos\delta_{ij})=0ΔQi=Qis−Qi=Qis−UiΣj=1nUj(Gijsinδij−Bijcosδij)=0
极坐标形式下的修正方程式为
[ΔPΔQ]=−[HNKL][ΔδUD2−1ΔU]\left[ \begin{matrix} \Delta P \\ \Delta Q \end{matrix} \right] = -\left[ \begin{matrix} H & N \\ K & L \end{matrix} \right]\left[ \begin{matrix} \Delta \delta \\ U^{-1}_{D2}\Delta U \end{matrix} \right] [ΔPΔQ]=−[HKNL][ΔδUD2−1ΔU]
其中当i≠ji\not=ji=j时，有
Hij=−UiUj(Gijsinδij−Bijcosδij)Nij=−UiUj(Gijcosδij+Bijsinδij)Kij=UiUj(Gijcosδij+Bijsinδij)Lij=−UiUj(Gijsinδij−Bijcosδij)}\left. \begin{matrix} H_{ij}=-U_iU_j(G_{ij}sin\delta_{ij}-B_{ij}cos\delta_{ij}) \\ N_{ij}=-U_iU_j(G_{ij}cos\delta_{ij}+B_{ij}sin\delta_{ij}) \\ K_{ij}=U_iU_j(G_{ij}cos\delta_{ij}+B_{ij}sin\delta_{ij}) \\ L_{ij}=-U_iU_j(G_{ij}sin\delta_{ij}-B_{ij}cos\delta_{ij}) \end{matrix} \right\} Hij=−UiUj(Gijsinδij−Bijcosδij)Nij=−UiUj(Gijcosδij+Bijsinδij)Kij=UiUj(Gijcosδij+Bijsinδij)Lij=−UiUj(Gijsinδij−Bijcosδij)⎭⎪⎪⎬⎪⎪⎫
当i=ji=ji=j时
Hij=Ui2Bii+QNij=−UiGii−PiKij=Ui2Gii−PiLij=Ui2Bii−Qi}\left. \begin{matrix} H_{ij}=U_i^2B_{ii}+Q \\ N_{ij}=-U_iG_{ii}-P_i \\ K_{ij}=U_i^2G_{ii}-P_i \\ L_{ij}=U_i^2B_{ii}-Q_i \end{matrix} \right\} Hij=Ui2Bii+QNij=−UiGii−PiKij=Ui2Gii−PiLij=Ui2Bii−Qi⎭⎪⎪⎬⎪⎪⎫
本文将基于以上主要原理方程式编写潮流计算程序。

2.具体实现过程

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

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

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

for i in range(total_count): #求Hfor 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_Psigma=0for 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]]=sigmadelta_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_busself.Z_bus = Z_busself.branch_resistance = branch_resistanceself.branch_reactance = branch_reactanceself.line_charging = line_chargingself.ratio = ratio

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

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

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

首先我们将U=[U1,U2,U3…U30]U = \left[ U_1 ,U_2,U_3 …U_{30} \right]U=[U1,U2,U3…U30]也就是一个30x1的列向量(以本文所采用的30节点系统为例)变成1x30的行向量 UT=[U1U2U3⋮U30]U^T=\left[\begin{matrix} U_1 \\ U_2 \\ U_3 \\ \vdots \\ U_{30}\end{matrix}\right]UT=⎣⎢⎢⎢⎢⎢⎡U1U2U3⋮U30⎦⎥⎥⎥⎥⎥⎤
然后利用相关函数将其复制转化成30x30的矩阵，即
Ucopy=[U1U2…U30U1U2…U30U1U2…U30⋮⋮⋮⋮U1U2…U30]U_{copy}=\left[\begin{matrix} U_1 & U_2 & \dots & U_{30} \\ U_1 & U_2 & \dots & U_{30} \\ U_1 & U_2 & \dots & U_{30} \\ \vdots & \vdots & \vdots & \vdots \\ U_1 & U_2 & \dots & U_{30}\end{matrix}\right]Ucopy=⎣⎢⎢⎢⎢⎢⎡U1U1U1⋮U1U2U2U2⋮U2………⋮…U30U30U30⋮U30⎦⎥⎥⎥⎥⎥⎤
之后将这个30x30的矩阵每一行都乘以对应行号的电压，也就是第一行每个元素都乘以U1U_1U1，第二行每个元素都乘以U2U_2U2，以此类推，也就相当于等式中的Ui∗UjU_i*U_jUi∗Uj，得到矩阵如下
UiUj=[U1U1U1U2…U1U30U2U1U2U2…U2U30U3U1U3U2…U3U30⋮⋮⋮⋮U30U1U30U2…U30U30]U_iU_j=\left[\begin{matrix} U_1U_1 & U_1U_2 & \dots & U_1U_{30} \\ U_2U_1 & U_2U_2 & \dots & U_2U_{30} \\ U_3U_1 & U_3U_2 & \dots & U_3U_{30} \\ \vdots & \vdots & \vdots & \vdots \\ U_{30}U_1 & U_{30}U_2 & \dots & U_{30}U_{30}\end{matrix}\right] UiUj=⎣⎢⎢⎢⎢⎢⎡U1U1U2U1U3U1⋮U30U1U1U2U2U2U3U2⋮U30U2………⋮…U1U30U2U30U3U30⋮U30U30⎦⎥⎥⎥⎥⎥⎤
这一过程可以通过上述行向量和列向量，利用Numpy库中的matmul函数直接实现。

然后将δ=[δ1,δ2,…,δ30]δ=[\delta_1 , \delta_2 ,\dots , \delta_{30}]δ=[δ1,δ2,…,δ30]也做类似的复制处理，得到
δcopy=[δ1δ2…δ30δ1δ2…δ30δ1δ2…δ30⋮⋮⋮⋮δ1δ2…δ30]\delta_{copy}=\left[\begin{matrix} \delta_1 & \delta_2 & \dots & \delta_{30} \\ \delta_1 & \delta_2 & \dots & \delta_{30} \\ \delta_1 & \delta_2 & \dots & \delta_{30} \\ \vdots & \vdots & \vdots & \vdots \\ \delta_1 & \delta_2 & \dots & \delta_{30}\end{matrix}\right] δcopy=⎣⎢⎢⎢⎢⎢⎡δ1δ1δ1⋮δ1δ2δ2δ2⋮δ2………⋮…δ30δ30δ30⋮δ30⎦⎥⎥⎥⎥⎥⎤
并将复制处理后的30x30格式的δδδ矩阵减去原来30*1格式的δδδ列向量，即
δij=[δ1−δ1δ2−δ1…δ30−δ1δ1−δ2δ2−δ2…δ30−δ2δ1−δ3δ2−δ3…δ30−δ3⋮⋮⋮⋮δ1−δ30δ2−δ30…δ30−δ30]\delta_{ij}=\left[\begin{matrix} \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 \\ \vdots & \vdots & \vdots & \vdots \\ \delta_1-\delta_{30} & \delta_2-\delta_{30} & \dots & \delta_{30}-\delta_{30}\end{matrix}\right] δij=⎣⎢⎢⎢⎢⎢⎡δ1−δ1δ1−δ2δ1−δ3⋮δ1−δ30δ2−δ1δ2−δ2δ2−δ3⋮δ2−δ30………⋮…δ30−δ1δ30−δ2δ30−δ3⋮δ30−δ30⎦⎥⎥⎥⎥⎥⎤
这样就得到了δijδ_{ij}δ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为节点导纳矩阵

把准备好的变量相乘后相加，得到的矩阵中每个元素的表达式都是ViVj∗(Gijcosδij+Bijsinδij)V_iV_j*(G_{ij}cos\delta_{ij}+B_{ij}sin\delta_{ij})ViVj∗(Gijcosδij+Bijsinδij),其中的iii和jjj分别对应元素所在的行和列，最后用求和操作即得到一个30x1的矩阵，每行元素对应每个节点的有功功率。而无功功率的操作基本与有功一致，不再赘述。
基于这种想法，我们再次观察前文所述的修正方程式以及雅可比矩阵中各子块元素的表达式
仔细观察可以发现，常规牛拉法的雅可比矩阵有四个分块矩阵，四个分块矩阵中元素的表达式也都用到3个量:UUU，δδδ以及节点导纳矩阵相应位置的元素。再观察一下我们还能发现一个特点: 在求取雅可比矩阵子块元素时，i=ji=ji=j和i≠ji\neq ji=j的表达式其实是可以合并的，就拿H子块来说，当i=ji=ji=j时，如果认为Qi=0Q_i=0Qi=0，那么无论是i=ji=ji=j还是i≠ji\neq ji=j，二者表达式完全一致，这就意味着我们完全可以先按照之前处理Pi和Qi那样用类似的方法去处理雅可比矩阵的生成，生成之后再单独索引提取出对角线元素并给其加上对应的PiP_iPi和QiQ_iQi即可。这样一来不仅简化了代码，而且用到了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节点还要额外计算无功不平衡量，因此直接这样命名了，现在想来好像有点草率了

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

电力系统matlab程序设计,电力系统分析潮流计算课程序设计及其MATLAB程序设计
目录一.程序设计目的 (2) 二.程序设计要求 (4) 三.13节点配网潮流计算 (4) 3.1主要流程............................................... ...
powerworld电力系统仿真，潮流计算，短路计算，电力系统分析。潮流计算对比，牛拉法，PQ分解法对比
powerworld电力系统仿真,潮流计算,短路计算,电力系统分析.潮流计算对比,牛拉法,PQ分解法对比编号:1710662437866344电气女博士
电力系统潮流计算及Matlab编程实现
目录 1. 潮流计算: 2. 潮流计算常用算法: 2.1 牛顿-拉夫逊算法 2.1.1 牛顿-拉夫逊法的基本原理 2.1.2 潮流计算的修正方程 2.1.3 节点电压用极坐标表示时的牛顿-拉夫逊法潮 ...
环形网络潮流计算matlab，利用matlab编程计算任意环形网络牛拉法潮流计算程序
环形网络潮流计算matlab,利用matlab编程计算任意环形网络牛拉法潮流计算程序,程序通用性强,通过修改参数可以得到任意节点和网络的环形网络牛拉法潮流计算. YID:696064261479453 ...
环形网络潮流计算matlab 利用matlab编程计算任意环形网络牛拉法潮流计算程序
环形网络潮流计算matlab 利用matlab编程计算任意环形网络牛拉法潮流计算程序,程序通用性强,通过修改参数可以得到任意节点和网络的环形网络牛拉法潮流计算. YID:856064261479453 ...
环形网络潮流计算matlab 利用matlab编程计算任意环形网络牛拉法潮流计算程序，程序通用性强
环形网络潮流计算matlab 利用matlab编程计算任意环形网络牛拉法潮流计算程序,程序通用性强,通过修改参数可以得到任意节点和网络的环形网络牛拉法潮流计算. 现有:6960642614794538 ...
环形网络潮流计算matlab 利用matlab编程计算任意环形网络牛拉法潮流计算程序，程序通用性强，通过修改参数
环形网络潮流计算matlab 利用matlab编程计算任意环形网络牛拉法潮流计算程序,程序通用性强,通过修改参数可以得到任意节点和网络的环形网络牛拉法潮流计算. 现有:6960642614794538 ...
MATLAB牛拉法计算潮流,matlab潮流计算
<matlab潮流计算>由会员分享,可在线阅读,更多相关<matlab潮流计算(14页珍藏版)>请在装配图网上搜索. 1.附录1使用牛顿拉夫逊法进行潮流计算的Matlab程序代 ...
用matlab计算潮流牛拉法,matlab潮流计算
<matlab潮流计算>由会员分享,可在线阅读,更多相关<matlab潮流计算(14页珍藏版)>请在人人文库网上搜索. 1.附录 1使用牛顿拉夫逊法进行潮流计算的Matlab ...
用matlab计算潮流牛拉法,Matlab牛拉法潮流计算程序
Matlab牛拉法潮流计算程序 V(k)=sqrt(e(k)^2+f(k)^2); %计算各节点电压的模值 sida(k)=atan(f(k)./e(k))*180./pi; %计算各节点电压的角度 ...

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

目录