不同风机接入电力系统的潮流计算方法比较

第一篇：不同风机接入电力系统的潮流计算方法比较

不同风机接入电力系统的潮流计算方法比较

中文摘要：各类风机由于其结构特点的不同决定了其接入电力系统的潮流计算方法的差异，通过对普通异步风机、永磁直驱同步风机、双馈异步风机接入电力系统的潮流计算方法讨论和比较，可以更加清楚的掌握对各类风机的使用。

关键词：普通异步风机永磁直驱同步风机双馈异步风机电力系统潮流计算 正文：

近年来“随着人们对可持续发展战略的认识”风力发电在全球获得了迅猛发展，截止2005年，全球风电装机容量已达,59332MW，排名前几位分别为德国18428MW，西班牙10027MW，美国9149MW，中国以1246MW排第8位/。风电场建设包括风资源测量评估、风场选址、发电机选型、风场内部微观选址、风电场升压站建设、接入电网线路建设等一系列问题、其中“风电场选址与接入系统设计都是规划阶段的重要问题” 其重要性不言而喻、对含风电场电力系统进行潮流计算“是规划阶段的基础工作”是风电场升压站、接入系统线路方案确定的考虑因素之一，另外“在风电系统动态仿真中”状态变量的初值也由潮流计算获得，因此"对含风电场电力系统进行潮流计算是非常有意义的。

风电场的潮流计算主要是对风电并网后根据给定的网络结构和运行条件确定整个网络的电气状态，主要是对风电并网后的电网中各个节点的电压幅值和相角、网络中功率的分布及功率损耗等，并进行越界检查，以了解和评价风电场并网后的运行情况。常用于评估风电机组并网后对电网稳态运行的影响，同时也为分析风电场并网后分析对电网影响等其他理论研究工作的基础，具有重要的意义。国内外对风电潮流的研究有着几十年的历史，风电场的潮流计算主要包括

含普通异步电机的风电场潮流计算、含永磁直驱风机的潮流计算和含双馈异步电机的风电场潮流计算。风电场的潮流计算关键是正确建立风电机组的数学模型。电力系统中节点分为ＰＱ节点、ＰＶ节点和平衡节点。一般异步发电机本身没有励磁调节装置，不具有调整节点电压的能力，因此不能像常规的同步发电机一样将它视为电压幅值恒定的ＰＶ节点，它只能依靠无功补偿装置才能保持风电场出口电压恒定。同样风电场中的异步发电机向系统注入有功功率的同时还要从系统吸收一定的无功功率，吸收的无功大小与发电机的机端电压、发出的有功功率以及滑差有着密切相关，因此不能把它处理为恒功率的ＰＱ节点。

一、异步电机计算潮流的方法 首先给出了风电场的简化ＰＱ模型、简化ＲＸ模型的潮流计算模型，根据两种模型的缺点，提出了一种新的风电场潮流计算模型即为扩展潮流模型。扩展潮流模型将异步发电机直接纳入电力网络，同时将异步发电机等值电路模型的中间激磁支路的端点作为一个虚拟ＰＱ节点进行计算，建立约束方程；建立异步发电机转子侧和定子侧之间的转矩平衡方程。模型建立详细潮流约束方程，推导雅可比矩阵元素，建立扩展风电潮流算法，最后通过算例分析验证算法的有效性异步风力发电机潮流建模—牛顿拉夫逊法潮流计算，牛顿法是解非线性方程的有效的方法。它把非线性方程的求解变成反复的求解线性方程，逐步接近非线性方程的解的过程，通常称为逐次线性化过程。而且牛顿一拉夫逊法求解潮流计算具有平方项收敛的速度，能够使潮流快速收敛。极坐标的牛顿拉夫逊法求解潮流问题的步骤如下：（１）计算有功功率和无功功率的不平衡量 假设系统中有n号节点，第1-m号节点为PQ节点，第m+1~n-1号为PV节点，n号为平衡节点。PQ节点不平衡量为： PV节点的功率不平衡量为：（２）计算雅可比矩阵 修正方程如下： 当i≠j时： 当i=j时：（３）潮流方程的求解化成下面的方程 根据上述步骤进行潮流计算，直至潮流熟练，输出结果。风电机组容量等值计算，等值后的风电机组容量为： 简化结构异步电机潮流算法： 对于常规的含普通异步电机的风电场潮流计算时，一般只考虑风力发电机和所连接的电网，忽略了风力机部分以及风力机与发电机之间的联系。同时，在进行潮流计算时，只考虑风机发出的功率全部输送到与风机连接电网中，不考虑风机由于浆矩角的变化对输出功率的影响。在进行风电场的潮流计算时，由于风机的机械部分的变化跟不上电气部分的变化速度，一般将风电场节点作为ＰＱ节点进行处理，但ＰＱ节点的无功功率与风机自身的参数以及风场节点电压有关。忽略空间和时间谐波、忽略励磁饱和、忽略铁损，将转子电阻和转子电抗折算到定子侧，图３．１（ａ）给出异步电机稳态Ⅱ等值电路，其中是定子和转子电阻，是定子和转子电抗，是激磁回路电抗，Ｓ为异步发电机的滑差，是对地并联电容。由于激磁电抗远大于定子电抗，且定子电阻较小。因此把励磁支路外移，合并定子和转子支路。采用简化等值电路时，型简化 等值电路如图３．１（ｂ）所示。大型风电场中的普通异步发电机以超同步转速情况下发电运行，发电机将风力机提供的机械功率转化为电磁功率输出，同时从电网或无功补偿装置中吸 收无功功率来维持励磁电流的大小。图３—１（ｂ）所示简化等值电路中，定子电抗与转子电抗合并式。定义电纳容性为正，电抗感性为负，将激磁支路式。的电抗与对地电容并联等值电抗为 假定异步风机流向电网功率为正，则将异步风机视为发电机处理，定义风机的输出有功功率如下：、发电机滑差以及无功功率 无功功率对节点电压求导得： 将风电场节点作为ＰＱ节点，根据风机出力曲线确定不同风速下的风机的出力，根据潮流建模公式代入潮流计算，采用算例分析，计算出潮流结果。简化RX模型： 采用ＲＸ模型进行潮流计算的基本思想是两个迭代步骤：（１）通过常规潮流计算计算出发电机的端电压；（２）通过异步发电机的滑差迭代计算出发电机的滑差。采用简化ＲＸ模型避免了潮流计算的两步迭代，能够节省潮流计算的计算时间。图３．１（ａ）为普通异步发电机的稳态等值电路图，计算出普通异步发电机的电磁功率。由异步发电机原理知道，异步发电机发出的电磁功率式为 由此可见异步发电机的输送到电网的电磁功率Ｂ随着滑差ｓ的变化而变化，同时风力机的转速、叶尖速比、风能利用系数以及机械功率也与滑差ｓ有关。根据系统的功率平衡的关系，风力机的机械功率等于输送到电网的有功功率相等。由于初始的功率不相等，通过滑差ｓ的迭代使两个功率最终达到平衡。采用牛顿拉夫逊计算潮流，引入风力机的机械功率与送到电网的有功功率的差和滑差修正量，修正方程为：。与转子电流、滑差ｓ等有关，其表达 程序框图如下：

二、双馈异步风机的潮流计算方法

本节在对双馈风机进行潮流建模时，将双馈风电机组的风电场作为恒定功率因数的ＰＱ节点进行处理，功率因数值设为０。９８。然后根据双馈异步发电机的稳态等值电路以及双馈风机的发电系统示意图分析各个功率之间的关系，及转矩平衡关系，建立约束方程，推导雅可比矩阵元素，建立潮流计算模型。图４—１（ａ）为双馈异步发电机的稳态等值电路，图４．１（ｂ）为双馈异步电机系统发电示意图。在图４．１（ａ），所为定子电压，复功率，以为转子电压为转子端到中间节点的为定子端到中间节点的复功率，其他参数同普通异步电机等值电路的参数意义相同。在图４—１（ｂ）中，风电场的节点设为i节点，为网侧变流器出口电压，分别为风机输入到电网的有功功率和无功功率，为变压器的阻抗。系统中的普通节点的潮流计算约束方程按照节点类型建立约束方程，推导雅可比矩阵元素。风电场节点的潮流计算变量的约束方程根据双馈风电机组的功率平衡关系以及风力机与发电机之间的转矩平衡条件来建立，并推导出对应的雅可比矩阵元素。令，以下定义类似。节点i流向节点Ｗ的复功率： 节点w流向定子侧的复功率为： 节点w流向转子侧的复功率为： 转子端流向节点w的复功率为： 网侧变流器端口流出的复功率为： 节点i流向变流器的复功率为：（１）风机与外网之间的连接 根据基尔霍夫电流定律，流入节点电流之和为零。分别注入ｉ节点的有功功率和无功功率。节点ｉ对应于电网有功功率平衡约束方程为： 对应变量的雅可比矩阵的元素为： 节点i对应于电网的无功功率平衡约束方程为： 对应变量的雅可比矩阵元素为：（２）双馈电机的等值电路 在双馈异步电机稳态等值电路中，虚拟内节点w连接定子支路和转子支路，以及激磁支路，图４．２为虚拟节点的电路结构。相比较普通异步电机，转子电压值不为零。根据节点功率平衡，建立内节点w有功功率平衡约束方程： 推导对应变量的雅可比矩阵元素为： 虚拟节点w的无功平衡约束方程为： 对应变量的雅可比矩阵的元素为：（３）双馈风机的外部电路 对于风电场出口节点来讲，不仅连接外电网和双馈异步发电机的定子侧，还通过变压器连接着网侧变流器，变压器的阻抗为为风电场出口处的结构，根据基尔霍夫电流原理，对ｉ节点注入功率之和为０。建立ｉ节点与发电机之间功率平衡约束方程。对i节点与发电机部分建立有功功率平衡约束方程： 对应的雅可比矩阵元素为： 节点i对与发电机部分建立无功功率平衡约束方程： 对应变量的雅可比矩阵的元素为： 在稳态运行时，忽略变换器的开关损耗。根据能量守恒原理，背靠背变流器输出有功功率之和等于零，即向两侧流出（或注入）的有功功率不变。网侧变流器输出无功为给定值一般取零，以防止ＰＶ控制方式下，电网需求无功过小，网侧变流器输出无功倒流至异步电机。对变流器建立约束方程，并推导雅可比 矩阵元素。背靠背变流器有功功率的约束方程如下： 对应变量的雅可比矩阵的元素为： 由于变流器之间不考虑无功功率的传输，网侧变流器采用单位功率因数的控制方式运行，的值设为０。对应的网侧变流器的无功功率约束方程为： 对应变量的雅可比矩阵为：（４）转矩平衡简化模型 双馈异步发电机的电磁功率为： 转矩平衡方程为： 双馈风机的机械功率采用最大功率跟踪计算，其表达式为： 和普通异步电机相比，在双馈异步电机的转矩平衡方程时，考虑了转子电压对功率平衡的影响。对应变量的雅可比矩阵元素为： 根据潮流计算约束方程以及推导雅可比矩阵元素，建立潮流计算模型： 双馈风电机组潮流计算的程序框图：

三、永磁直驱同步风机的潮流计算

牛顿法是解非线性方程的有效的方法。它把非线性方程的求解变成反复的求解线性方程，逐步接近非线性方程的解的过程，通常称为逐次线性化过程。这是牛顿法的核心。用牛顿法解三相潮流问题的步骤如下：（1）、计算功率不平衡方程

式中表示了一个有n+1个母线的系统功率不平衡矩阵，其中有m个PQ母线，n-m个PV母线，1个平衡节点。

PQ节点的功率不平衡量为该节点的功率给定值与当前电压计算出来的实际功率差，可表示为：

其中i=1,2,3….，p=a，b，c。

而对ＰＶ节点来说，节点电压幅值是给定的，不再作为变量。同时，该点无法预先给定无功功率。这样，该点的无功不平衡量也就失去了约束作用。因此，在迭代过程中无须计算与ＰＶ节点有关的无功功率方程式。

只有当迭代结束后，即各节点的电压向量求得后，再求PV节点应当维持的无功功率。（2）计算雅可比矩阵

（４）含风电场的电力系统三相潮流的求解最终能化成求解下列方程

风电场在牛顿法中的处理

风电场等各种分布式电源可以建立成ＰＱ节点，ＰＶ节点，ＰＩ节点和Ｐ－Ｑ（Ｖ）节点这四种节点类型。对ＰＱ类型的分布式电源只需将它们简单处理成功率值是负的负荷即可。本节主要分析其他三种类型的分布式电源在程序中的处理。（１）Ｐ恒定，Ｖ恒定的ＰＶ节点

ＰＶ节点可以直接代入牛顿法中处理。若ＰＶ母线与系统通过ｎ（ｎ＝ｌ，２，３…）相线路连接，则母线上各节点注入功率为母线总注入功率的ｎ分之一。在每次迭代后，可以求出节点的电压相角和无功功率。若计算出的节点无功越限，则将其转换成对应的ＰＱ节点，Ｑ值等于该分布式电源能发出的最大无功值。如果在后续迭代中，又出现该节点电压越界，重新将其转换成ＰＶ节点。（２）Ｐ恒定，电流幅值Ｉ恒定的ＰＩ节点

ＰＩ节点不可以直接代入牛顿法中处理，所以在每次迭代前须做一定的处理。若ＰＩ母线与系统通过ｎ（ｎ＝ｌ，２，３…）相线路连接，则母线上各节点注入功率为母线总注入功率的ｎ分之一。相应的无功功率可以由上一次迭代得到的电压，给定的电流幅值和有功功率计算得出：

其中，为第ｋ＋１次迭代的分布式电源的无功功率值；

五分别为第ｋ次迭代得到的电压的实部和虚部；Ｉ为恒定的分布式电源的电流相量的幅值；尸为恒定的有功功率值。

因此在进行潮流计算时，第抖１次迭代前可以把ＰＩ节点的无功注入量求出，在第ｋ＋ｌ迭代过程中便可将ＰＩ节点处理成有功和无功输出分别为Ｐ和的ＰＱ节点。在每次迭代后，可以求出节点的电压相角和无功功率。ＰＩ型的分布式电源也有无功输出的限制，但从式（３－６）可以看出，的标幺值一般在１．０附近，Ｐ和，是

。值的只是该ＰＩ两个必需维持的值，所以影响最后计算出来的节点的给定有功功率和电流幅值，即Ｐ和，若给定得合理，则计算得出的无功功率不会越限。

（３）Ｐ恒定，Ｖ不定，Ｑ受Ｐ、Ｖ限定的Ｐ—Ｑ（Ｖ）节点Ｐ－Ｑ（Ｖ）节点不可以直接代入牛顿法中处理，所以在每次迭代前须做一定的处理。Ｐ－Ｑ（Ｖ）节点给定的输出有功功率只为异步电机的输出有功功率，节点电压Ｕ在每次迭代后都得到修正，节点的注入无功功率Ｑ计算公式如下：

其中，ｓ为异步电机的转差率；为发电机定子电抗与转子电抗之和：为励磁电抗；ｒ为转子电阻；Ｑ’为异步电机的吸收无功；为异步电机的功率因数：

仍为并联电容器后节点的功率因数；一般要求在０．９以上；鳞为并联电容器需要补偿的无功；刀为投入的并联电容器组数；

为每组电容器补偿的无功；Ｑ’为电容器组实际补偿的无功；Ｑ为参与潮流迭代的节点注入无功。这里提出的动态调整并联电容器组接入组数的方法，相比于将风机视为功率因数恒定的考虑，更加符合实际异步风机运行的情况，因而，基于该模型的分析计算结果更为精确。Ｐ－Ｑ（Ｖ）母线与系统通过ｎ（ｎ＝ｌ，２，３…）相线路连接，则母线上各节点注入功率为母线总注入功率的ｎ分之一。在进行潮流计算时，第ｋ次迭代后可以把Ｐ-Ｑ（Ｖ）节点的无功吸收量求出，在第ｋ＋ｌ迭代过程中便可将Ｐ-Ｑ（Ｖ）节点处理成有功和无功输出分别为Ｐ和的ＰＱ节点。

通过对以上三种风机的潮流计算方法的分析，我们对风机接入电力系统有了更深入的了解，对比各风机的特点，可得到不同的结论。首先普通异步风机具有一般异步电机的特点，只是对PQ节点的处理方法稍有不同。通过对双馈风电机组的潮流建模建立及验证，计及双馈电机详细内部结构和各种稳态安全约束，采用优化算法求解不同风速下无功出力范围；基于电机有功损耗最小，求解无功最优分布。根据双馈风电机组的潮流计算，获得风电机组的稳态运行情况和发电机运行参数，验证了双馈风电机组潮流算法的有效性。在进行无功优化时，发出无功时，最大无功出力受转子绕组电流限制。吸收无功时，最小无功出力当定子绕组电流限制。调整网侧变流器参考无功出力，可以增加双馈电机无功出力范围，但是在恒ＰＱ运行方式下，可能引起无功环流，增加电机的损耗。随风速增加或者定子电压下降，双馈异步电机无功出力范围变窄。采用恒功率因数方式运行时，低风速下双馈电机无功出力得不到充分利用，而高风速、低定子电压下可能达不到给定功率因数要求。根据最小有功损耗优化计算出的无功功率，作为双馈异步电机稳态运行时的输出无功功率值，可提高双馈异步电机的运行效率。而虽然双馈异步电机占主导地位，但永磁直驱也解决了一些双馈电机解决不了的问题。参考文献：

第二篇：电力系统三种潮流计算方法的比较

电力系统三种潮流计算方法的比较

一、高斯-赛德尔迭代法：

以导纳矩阵为基础，并应用高斯--塞德尔迭代的算法是在电力系统中最早得到应用的潮流计算方法，目前高斯一塞德尔法已很少使用。将所求方程

f(x)



0

改写为

x(x)不能直接得出方程的根，给一个猜测值 x0得 x1(x0)又可取x1为猜测值，进一步得：

x2(x1)

迭代 反复猜测

xk1(xk)

*mxk则方程的根 xlik

优点：

1.原理简单，程序设计十分容易。

2.导纳矩阵是一个对称且高度稀疏的矩阵，因此占用内存非常节省。

3.就每次迭代所需的计算量而言，是各种潮流算法中最小的，并且和网络所包含的节点数成正比关系。缺点：

1.收敛速度很慢。

2.对病态条件系统，计算往往会发生收敛困难：如节点间相位角差很大的重负荷系统、包含有负电抗支路(如某些三绕组变压器或线路串联电容等)的系统、具有较长的辐射形线路的系统、长线路与短线路接在同一节点上，而且长短线路的长度比值又很大的系统。

3.平衡节点所在位置的不同选择，也会影响到收敛性能。

二、牛顿-拉夫逊法：

求解 f(x)0设

x



x 0

x

，则

f(x0x)0

按牛顿二项式展开：

f(x)f(x)x1f(x)(x)21f(x)(x)3000002!3!

当△x不大，则取线性化（仅取一次项）

f(x0)f(x0)x0

1则可得修正量 xf(x0)f(x0)

x1x0x

x对

f

（）



0

得： f(x0)xf(x0)作变量修正：

k 

1

x

k





x

k，求解修正方程

x牛顿法是数学中求解非线性方程式的典型方法，有较好的收敛性。自从20世纪60年代中期采用了最佳顺序消去法以后，牛顿法在收敛性、内存要求、计算速度方面都超过了其他方法，成为直到目前仍被广泛采用的方法。优点：

1.收敛速度快，若选择到一个较好的初值，算法将具有平方收敛特性，一般迭代4—5次便可以收敛到一个非常精确的解。而且其迭代次数与所计算网络的规模基本无关。

2.具有良好的收敛可靠性，对于前面提到的对以节点导纳矩阵为基础的高斯一塞德尔法呈病态的系统，牛顿法均能可靠地收敛。

3.牛顿法所需的内存量及每次迭代所需时间均较前述的高斯一塞德尔法为多，并与程序设计技巧有密切关系。缺点：

牛顿法的可靠收敛取决于有一个良好的启动初值。如果初值选择不当，算法有可能根本不收敛或收敛到一个无法运行的解点上。解决方法：

对于正常运行的系统，各节点电压一般均在额定值附近，偏移不会太大，并且各节点间的相位角差也不大，所以对各节点可以采用统一的电压初值(也称为“平直电压”)，“平直电压”法假定：

Ui01i00

或

ei01fi00(i1,2,...,n;is)

这样一般能得到满意的结果。但若系统因无功紧张或其它原因导致电压质量很差或有重载线路而节点间角差很大时，仍用上述初始电压就有可能出现问题。可以先用高斯一塞德尔法迭代1-2次；以此迭代结果作为牛顿法的初值，也可以先用直流法潮流求解一次以求得一个较好的角度初值，然后转入牛顿法迭代。

三、P-Q分解法：

电力系统中常用的PQ分解法派生于以极坐标表示的牛顿—拉夫逊法，其基本思想是把节点功率表示为电压向量的极坐标形式，以有功功率误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，把有功和无功分开进行迭代其主要特点是以一个（n-1）阶和一个m阶不变的、对称的系数矩阵B,B代替原来的（n+m-1）阶变化的、不对称的系数矩阵M，以此提高计算速度，降低对计算机贮存容量的要求。P-Q分解法在计算速度方面有显著的提高，迅速得到了推广。原理：

PH修正方程为：QKNV LV雅克比矩阵元素的表达如下： a)当i≠j时

HijViVj(GijsinijBijcosij)NijViVj(GijcosijBijsinij)

KijViVj(GijcosijBijsinij)LijViVj(GijsinijBijcosij)b)当i＝j时

HiiVi2BiiQi2NiiViGiiPi 2KiiViGiiPiLiiVi2BiiQi 对修正方程的第一个简化是：

PH0Q0LV V上式可分别写成以下两式

PH

QLVV

在一般情况下，线路两端电压的相角差是不大的(不超过100~200)，因此可以认为 cosij1,Gijsinij《Bij

因此可得：HijViVjBij

（i，j=1，2，„，n-1）

LijViVjBij

（i，j=1，2，„，m）HiiVi2BiiQiLiiVBiiQi2i

经一系列化简得P—Q分解法的修正方程式：

PB

QBVBVV原P—Q分解法的修正方程的简化形式为：



QBVV

PPQ分解法的修正方程式的特点：

1.以一个(n-1)阶和一个(m-1)阶系数矩阵B、B替代原有的系数矩阵J，提高了计算速度，降低了对贮存容量的要求。

2.以迭代过程中保持不变的系数矩阵B、B替代原有的系数矩阵J，显著的提高了计算速度。

3.以对称的系数矩阵B、B替代原有的系数矩阵J，使求逆等运算量和所需的储存容量都大为减少。P-Q分解法两个主要特点：

1.降阶在潮流计算的修正方程中利用了有功功率主要与节点电压相位有关,无功功率主要与节点电压幅值有关的特点,实现P-Q分解,使系数矩阵由原来的2N×2N阶降为N×N阶,N为系统的节点数(不包括缓冲节点)。

2.因子表固定化利用了线路两端电压相位差不大的假定,使修正方程系数矩阵元素变为常数,并且就是节点导纳的虚部。

由于以上两个特点,使快速分解法每一次迭代的计算量比牛顿法大大减少。P-Q分解法只具有一次收敛性,因此要求的迭代次数比牛顿法多,但总体上快速分解法的计算速度仍比牛顿法快。快速分解法只适用于高压网的潮流计算,对中、低压网,因线路电阻与电抗的比值大,线路两端电压相位差不大的假定已不成立,用快速分解法计算,会出现不收敛问题。

第三篇：电力系统潮流计算

南 京 理 工 大 学

《电力系统稳态分析》

课程报告

姓名

XX

学 号： 5*** 自动化学院 电气工程

基于牛顿-拉夫逊法的潮流计算例题编程报学院(系): 专

业: 题

目: 任课教师 硕士导师 告

杨伟 XX

2015年6月10号

基于牛顿-拉夫逊法的潮流计算例题编程报告

摘要：电力系统潮流计算的目的在于：确定电力系统的运行方式、检查系统中各元件是否过压或者过载、为电力系统继电保护的整定提供依据、为电力系统的稳定计算提供初值、为电力系统规划和经济运行提供分析的基础。潮流计算的计算机算法包含高斯—赛德尔迭代法、牛顿-拉夫逊法和P—Q分解法等，其中牛拉法计算原理较简单、计算过程也不复杂，而且由于人们引入泰勒级数和非线性代数方程等在算法里从而进一步提高了算法的收敛性和计算速度。同时基于MATLAB的计算机算法以双精度类型进行数据的存储和运算, 数据精确度高，能进行潮流计算中的各种矩阵运算，使得传统潮流计算方法更加优化。

一 研究内容

通过一道例题来认真分析牛顿-拉夫逊法的原理和方法（采用极坐标形式的牛拉法），同时掌握潮流计算计算机算法的相关知识，能看懂并初步使用MATLAB软件进行编程，培养自己电力系统潮流计算机算法编程能力。

例题如下：用牛顿-拉夫逊法计算下图所示系统的潮流分布，其中系统中5为平衡节点，节点5电压保持U=1.05为定值，其他四个节点分别为PQ节点，给定的注入功率如图所示。计算精度要求各节点电压修正量不大于10-6。

二 牛顿-拉夫逊法潮流计算 1 基本原理

牛顿法是取近似解x(k)之后，在这个基础上，找到比x(k)更接近的方程的根，一步步地迭代，找到尽可能接近方程根的近似根。牛顿迭代法其最大优点是在方程f(x)=0的单根附近时误差将呈平方减少，而且该法还可以用来求方程的重根、复根。电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点的电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成节点电压新的初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。2 基本步骤和设计流程图

形成了雅克比矩阵并建立了修正方程式，运用牛顿-拉夫逊法计算潮流的核心问题已经解决，已有可能列出基本计算步骤并编制流程图。由课本总结基本步骤如下：

1)形成节点导纳矩阵Y；

2)设各节点电压的初值，如果是直角坐标的话设电压的实部e和虚部f；如果是极坐标的话则设电压的幅值U和相角a；

3)将各个节点电压的初值代入公式求修正方程中的不平衡量以及修正方程的系数矩阵的雅克比矩阵；

4)解修正方程式，求各节点电压的变化量，即修正量； 5)计算各个节点电压的新值，即修正后的值；

6)利用新值从第（3）步开始进入下一次迭代，直至达到精度退出循环； 7)计算平衡节点的功率和线路功率，输出最后计算结果； ① 公式推导

② 流程图

三

matlab编程代码

clear;

% 如图所示1，2，3，4为PQ节点,5为平衡节点

y=0;

% 输入原始数据，求节点导纳矩阵

y(1,2)=1/(0.07+0.21j);

y(4,5)=0;y(1,3)=1/(0.06+0.18j);

y(1,4)=1/(0.05+0.10j);

y(1,5)=1/(0.04+0.12j);

y(2,3)=1/(0.05+0.10j);

y(2,5)=1/(0.08+0.24j);

y(3,4)=1/(0.06+0.18j);

for i=1:5

for j=i:5

y(j,i)=y(i,j);

end

end

Y=0;

% 求节点导纳矩阵中互导纳

for i=1:5

for j=1:5

if i~=j

Y(i,j)=-y(i,j);

end

end

end

% 求节点导纳矩阵中自导纳

for i=1:5

Y(i,i)=sum(y(i,:));

end

Y

% Y为导纳矩阵

G=real(Y);

B=imag(Y);% 输入原始节点的给定注入功率

S(1)=0.3+0.3j;

S(2)=-0.5-0.15j;

S(3)=-0.6-0.25j;

S(4)=-0.7-0.2j;

S(5)=0;

P=real(S);

Q=imag(S);

% 赋初值，U为节点电压的幅值，a为节点电压的相位角

U=ones(1,5);

U(5)=1.05;

a=zeros(1,5);

x1=ones(8,1);

x2=ones(8,1);

k=0;

while max(x2)>1e-6

for i=1:4

for j=1:4

H(i,j)=0;

N(i,j)=0;

M(i,j)=0;

L(i,j)=0;

oP(i)=0;

oQ(i)=0;

end

end

% 求有功、无功功率不平衡量

for i=1:4

for j=1:5

oP(i)=oP(i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));

oQ(i)=oQ(i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));

end

oP(i)=oP(i)+P(i);

oQ(i)=oQ(i)+Q(i);

end

x2=[oP,oQ]';

% x2为不平衡量列向量

% 求雅克比矩阵

% 当i~=j时，求H,N,M,L

for i=1:4

for j=1:4

if i~=j

H(i,j)=-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));

N(i,j)=-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));

L(i,j)=H(i,j);

M(i,j)=-N(i,j);

end

end

end

% 当i=j时，求H,N,M,L

for i=1:4

for j=1:5

if i~=j H(i,i)=H(i,i)+U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));N(i,i)=N(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));

M(i,i)=M(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));

L(i,i)=L(i,i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))

end

end

N(i,i)=N(i,i)-2*(U(i))^2*G(i,i);

L(i,i)=L(i,i)+2*(U(i))^2*B(i,i);

end

J=[H,N;M,L]

% J为雅克比矩阵

x1=-((inv(J))*x2);

% x1为所求△x的列向量

% 求节点电压新值，准备下一次迭代

for i=1:4

oa(i)=x1(i);

oU(i)=x1(i+4)*U(i);

end

for i=1:4

a(i)=a(i)+oa(i);

U(i)=U(i)+oU(i);

end

k=k+1;

end

k,U,a

% 求节点注入功率

i=5;

for j=1:5

P(i)=U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)))+P(i);

Q(i)=U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))+Q(i);

end

S(5)=P(5)+Q(5)*sqrt(-1);

S

% 求节点注入电流

I=Y*U'

四

运行结果

节点导纳矩阵

经过五次迭代后的雅克比矩阵

迭代次数以及节点电压的幅值和相角（弧度数）

节点注入功率和电流

五 结果分析

在这次学习和实际操作过程里：首先，对电力系统分析中潮流计算的部分特别是潮流计算的计算机算法中的牛顿-拉夫逊法进行深入的研读，弄明白了其原理、计算过程、公式推导以及设计流程。牛顿-拉夫逊法是求解非线性方程的迭代过程，其计算公式为FJX，式中J为所求函数的雅可比矩阵；X为需要求的修正值；F为不平衡的列向量。利用x(*)=x(k+1)+X(k+1)进行多次迭代，通过迭代判据得到所需要的精度值即准确值x(*)。六 结论

通过这个任务，自己在matlab编程，潮流计算，word文档的编辑功能等方面均有提高，但也暴漏出一些问题：理论知识储备不足，对matlab的性能和特点还不能有一个全面的把握，对word软件也不是很熟练，相信通过以后的学习能弥补这些不足，达到一个新的层次。

第四篇：电力系统潮流计算程序设计

电力系统潮流计算程序设计

姓名：韦应顺

学号：2011021052 电力工程学院

牛顿—拉夫逊潮流计算方法具有能够将非线性方程线性化的特点，而使用MATLAB语言是由于MATLAB语言的数学逻辑强，易编译。

【】【】1.MATLAB程序12

Function tisco %这是一个电力系统潮流计算的程序 n=input（‘n请输入节点数：n=’）； m=input(‘请输入支路数：m=’);ph=input(‘n请输入平衡母线的节点号：ph=’)； B1=input（‘n请输入支路信号：B1=’）;%它以矩阵形式存贮支路的情况，每行存贮一条支路 %第一列存贮支路的一个端点 %第二列存贮支路的另一个端点 %第三列存贮支路阻抗

%第四列存贮支路的对地导纳

%第五列存贮变压器的变比，注意支路为1 %第六列存贮支路的序号

B2=input(‘n请输入节点信息：B2=’)； %第一列为电源侧的功率 %第二列为负荷侧的功率 %第三列为该点的电压值

%第四列为该点的类型：1为PQ，2为PV节点，3为平衡节点 A=input（‘n请输入节点号及对地阻抗：A=’）； ip=input（‘n请输入修正值：ip=’）； %ip为修正值);Y=zeros（n）；

Y(p,q)=Y(p,q)-1./(B1(i3)*B1(i5);e=zeros（1，n）；

Y(p,q)=Y(p,q);f=zeros（1，n）；

no=2*ph=1； Y(q,q)=Y(q,q)+1./B1(i3)+B1(i4)/2；

End for i=1：n

G=real(Y);if A（i2）=0

B=imag(Y);p=A(i1)；

Y(p p)=1./A(i2);for i=1：n End e(i)=real(B2(i3));End f(i)=imag(B2(i3));For i=1:m S(i)=B2(i1)-B2(i2);p=B1(i1);V(i)=B2(i3);p=B1(i2);end Y(p,p)=Y(p,p)+1./(B1(i3)*B1(i5)^2+B1(i4)./2P=real(S);Q=imag(S);[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no);J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no);[De,Di]=hxf(J,D,F,ph,n,no);t=0;while

max(abs(De))>ip&max(abs(Dfi)>ip

t=t+1;

e=e+De；

f=f+Df;

[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)；

J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no)；

[De,Df]=hxf(J,Df,ph,n,no);end v=e+f*j;for i=1:n hh(i)=conj(Y(ph,i)*v(i));end S(ph)=sum(hh)*v(ph);B2(ph,1)=S(ph);V=abs(v);

jd=angle(v)*180/p;resulte1=[A(:,1),real(v),imag(v),V,jd,real(S’)，imag(S’),real(B2(:1)),imag(B2(:1)),real(B2(:2)),imag(B2(:,2))];for i=1:m

a(i)=conj((v(B1(i1))/B1(i5)-v(B1(i2))/B1(i3));

b(i)=v(B1(i1))*a(i)-j*B1(i4)*v(B1(i))^2/2;

c(i)=-v(B1(i2))*a(i)-j*B1(i4)*v(B1(i2))^2/2;end result2=[B1(:,6),B1(:,1),B1(:,2),real(b’),imag(b’),real(c’),imag(c’), real(b’+c’),imag(b’+c’)];printcut(result1,S,b,c,result2);type resultm function [C,D,Df]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)%该子程序是用来求取Df for i=1:n

If

i=ph

C(i)=0;

D(i)=0;

For j=i:n

C(i)=C(i)+G(i,j)*e(j)-B(i,j)*f(j);D(i)=D(i)+G(i,j)*f(j)+B(i,j)*e(j);end

P1=C(i)*e(i)+D(i)*f(i);Q1=C(i)*f(i)-D(i)*e(i);V1=e(i)^2+f(i)^2;If

B2(i4)=2 p=2*i-1;

Df(p)=P(i)-P1;p=p+1;else p=2*i-1;

Df(p)=P(i)-P1;p=p+1;

Df(p)=Q(i)-Q1;end end end Df=Df’;If ph=n Df(no=[];end

function [De,Df]=hxf(J,Df,ph,n,no)%该子函数是为求取De Df DX=JDf;DX1=DX;

x1=length(DX1);if ph=n DX(no)=0;DX(no+1)=0;

For i=(no+2):(x1+2)DX(i)=DX1(i-2);End Else

DX=[DX1,0,0];End k=0;

[x,y]=size(DX);For i=1:2:x K=k+1;

Df(k)=DX(i);De(k)=DX(i+1);End End case 2 Function for j=1:n J=jacci(Y,G,B,PQ,e,f,V,C,D,B2,n,ph,no)X1=G(i,j)*f(i)-B(i,j)*e(i);

X2=G(i,j)*e(i)+B(i,j)*f(i);%该子程序是用来求取jacci矩阵

for i=1:n X3=0;switch B2(i4)X4=0;case 3 P=2*i-1;continue q=2*j-1;case 1 J(p,q)=X1;for j=1:n m=p+1;if

J=&J=ph J(m,q)=X3;X1=G(i)*f(i)-B(i,j)*e(i);q=q+1;X2=G(i,j)*e(i)+B(i,j)*f(i);J(p,q)=X2;X3=-X2;J(m,q)=X4;X4=X1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);p=2*i-1;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=2*j-1;X3=0;J(p,q)=X1;X4=0;m=p+1;P=2*i-1;J(p,q)=X2;q=2*j-1;J(m,q)=X4;J(p,q)=X1;Else if j=&j=jph m=p+1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X3;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=q+1;X3= C(i)+G(i,j)*e(i)-B(i,j)*f(i);J(p,q)=X2;X4= C(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X4;P=2*i-1;end q=2*j-1;end J(p,q)=X1;end m=p+1;end J(m,q)=X3;if ph=n q=q+1;J(no:)=[];J(p,q)=X2;J(no:)=[];J(m,q)=X4;J(:,no)=[];End J(:,no)=[];End

2实例验证 【例题】设有一系统网络结线见图1，各支路阻抗和各节点功率均已以标幺值标示于图1中，其中节点2连接的是发电厂，设节点1电压保持U1＝1.06定值，试计算其中的潮流分布，请输入节点数：n=5 请输入支路数：m=7 请输入平衡母线的节点号：ph=l 请输入支路信息：

BI=[ l 2 0．02+0．06i O l 1；1 3 0．08+0．24i 0 1 2；2 3 0．06+0．18i 0 l 3： 2 4 0．06+0．18i O l 4： 2 5 0．04+0．12i 0 l 5： 3 4 0．01+0．03i 0 l 6： 4 5 0．08+0．24i O 1 7] 请输入节点信息：

B2=[ 0 0 1．06 3；0．2+0．20i 0 1 1；一O．45一O．15i 0 l l；一0．4-0．05i 0 l 1；一0．6—0．1i 0 1 l] 请输入节点号及对地阻抗： A=[l 0；2 0；3 0；4 0；5 O ] 请输入修正值：ip=0．000 0l

参考文献

[1]陈珩．电力系统稳定分析[M]．北京：中国电力出版社，2002：139—187．

[2]郑阿奇．MATLAB实用教程[M]．北京：电子工业出版社，2005：1-243．

[3] 束洪春，孙士云，等．云电送粤交商流混联系统全过 程动态电压研究[J】．中国电力，2008，4l(10)：l-4． SHU Hong—ch吼，SUN Shi-yun，et a1．Research on fun prc'cess dyn锄ic Voltage stabil时of hybrid AC／DC poWer tmnsmission System舶m Yu衄an proVince to G啪gdong province【J】．Electric Power，2008，4l(10)： l-4．

[4] 朱新立，汤涌，等．大电网安全分析的全过程动态仿 真技术[J】．电网技术，2008，32(22)：23—28． SONG Xin—Ii，TANG Yof唱，et a1． Full dyn锄ic simulation for the stabilhy a眦lysis of large power system【J】．Power System融IlrIolo影，2008，32(22)： 23．28．

[5]Roytelm锄I，Shallidehpour S M．A comprehcnsivc long teml dynaIIlic simulation for powcr system recoVery【J】． IEEE Transactions 0n Power Systems，1994，9(3)． [6] 石雩梅，汪志宏，等．发电机励磁系统数学模型及参 数对电网动态稳定性分析结果影响的研究[J】．继电 器，2007，35(21)：22-27．

SHI Xue．mei，WANG Zlli-hon舀et a1．Iksearch on the innuence of g锄e翰to璐baScd ∞de诅iled excitation system models柚d parameterS t0 power铲id dyn锄ic stabil时【J】．Relay，2007，35(2 1)：22-27．

[7] 方思立，朱方．快速励磁系统对系统稳定的影响[J】．中 国电机工程学报，1986，6(1)：20．28．

FANG Si．1i，ZHU Fang．The effbct of f弧t．respon∞

excitation system on the stability of power netwofk【J】． Proceedings ofthe CSEE，1986，6(1)：20-28．

[8] 刘取．电力系统稳定性及发电机励磁控制[M】．北京： 中国电力出版社，2007．

LIU Qu．Power system S诅bility锄d generator excitation control【M】．BeUing：ChiIla Electric Powef Press，2007． [9] Dallachy J L，Anderson T．EXperience with rcplacing ro诅ting exciters wim static exciters【J】．1k InStitution of Electrical Engineers，1 996．

[10] 陈利芳，陈天禄．浅谈自并励励磁系统在大容量机组 中的应用【J】．继电器，2007，35(1)：8l培4． CHEN Li-f抽岛CHEN Tian—lIL Application of 辩l仁exci组tion mode in large capacity髫memtor unit【J】． ReIay'2007，35(1)：81-84．

[11] 方思立，刘增煌，孟庆和．大型汽轮发电机自并励励 磁系统的应用条件【J]．中国电力，1994，27(12)：61．63． FANG Si．Ii，LIU Zeng-hu锄g，MENG Qin争hc．m application conditions of large turbine generator self-excitation system【J】．Electric Powef，1994，27(12)： 61．63．

[12]梁小冰，黄方能．利用EMTDC进行长持续时间过程 的仿真研究【J】．电网技术，2002，26(9)：55．57． LIANG Xiao-bing，HUANG Fan争眦ng．How to cany out simulalion of long dul‘ation processes by use of EMTDC【J】．Power System 11echnology，2002，26(9)： 55-57．

[13]王卉，陈楷，彭哲，等．数字仿真技术在电力系统中 的应用及常用的几种数字仿真工具【J】．继电器，2004，32(21)：7l一75．

wANG Hui，CHEN Kai，PENG zhe，et a1．Application of digital simulation眦hniques棚d severaJ simulation tools in power system[J】．Relay，2004，32(21)：71·75．

[14]IEEE Power Engmeering Socie哆．IEEE std 421．5．2005 IEEE玎ccOmmended practice for excitation system models for power system stabiI时studies【s】．

第五篇：电力系统潮流计算程序

电力系统潮流计算c语言程序，两行，大家可以看看，仔细研究，然后在这个基础上修改。谢谢

#include “stdafx.h” #include #include #include

#include“Complex.h” #include“wanjing.h” #include“gauss.h” using namespace std;

int _tmain(int argc, _TCHAR* argv[]){

int i;

//i作为整个程序的循环变量

int N=Bus::ScanfBusNo();//输入节点个数

int L=Line::ScanflineNo();//输入支路个数

if((L&&N)==0){return 0;} //如果找不到两个文件中的任意一个，退出

Line *line=new Line[L];//动态分配支路结构体

Line::ScanfLineData(line);//输入支路参数

Line::PrintfLineData(line,L);//输出支路参数

Bus *bus=new Bus[N];//动态分配结点结构体

for(int i=0;i

bus[i].Sdelta.real=0;

bus[i].Sdelta.image=0;}

Bus::ScanfBusData(bus);//输入节点参数

Bus::PrintfBusData(bus,N);//输出结点参数

Complex **X;X=new Complex *[N];for(i=0;i

Bus::JisuanNodeDnz(X,line,bus,L,N);//计算节点导纳矩阵

Bus::PrintfNodeDnz(X,N);//输出节点导纳矩阵

int NN=(N-1)*2;double **JacAug;JacAug=new double *[NN];for(i=0;i

double *x;x=new double[NN];int count=1;

LOOP:

Bus::JisuanNodeI(X,bus,N);//计算节点注入电流

Bus::JisuanNodeScal(X,bus,N);//计算节点功率

Bus::JisuanNodeScal(X,bus,N);//计算节点功率

Bus::JisuanNodeSdelta(bus,N);//计算节点功率差值

Bus::PrintfNodeScal(X,bus,N);//输出节点功率差值

int icon=wehcon1(bus,N);//whether converbence看迭代是否结束

if(icon==1){

cout<<“icon=”<

Bus::JisuanJacAug(JacAug,X,bus,N);//计算雅可比增广矩阵 // Bus::PrintfJacAug(JacAug,N);

gauss::gauss_slove(JacAug,x,NN);//解方程组求出电压差值

Bus::ReviseNodeV(bus,x,N);//修正节点电压

// Bus::PrintfNodeV(bus,N);

count++;

goto LOOP;}

else

{

for(i=0;i

{

int statemp,endtemp;

Complex aa,bb,cc,dd,B;

B.real=0;

B.image=-line[i].B;

statemp=line[i].start;

endtemp=line[i].end;

aa=Complex::productComplex(Complex::getconj(bus[statemp-1].V), B);

bb=Complex::subComplex

(Complex::getconj(bus[statemp-1].V), Complex::getconj(bus[endtemp-1].V));

cc=Complex::productComplex(bb , Complex::getconj(line[i].Y));

dd=Complex::CaddC(aa,cc);

line[i].stoe=Complex::productComplex(bus[statemp-1].V,dd);

aa=Complex::productComplex(Complex::getconj(bus[endtemp-1].V), B);

bb=Complex::subComplex

(Complex::getconj(bus[endtemp-1].V), Complex::getconj(bus[statemp-1].V));

cc=Complex::productComplex(bb , Complex::getconj(line[i].Y));

dd=Complex::CaddC(aa,cc);

line[i].etos=Complex::productComplex(bus[endtemp-1].V,dd);

}

cout<<“icon=”<

Bus::JisuanNodeScal(X,bus,N);//计算节点功率

for(i=0;i

{

bus[i].Scal.real = bus[i].Scal.real + bus[i].Load.real;//发电机功率=注入功率+负荷功率

bus[i].Scal.image= bus[i].Scal.image+ bus[i].Load.image;

bus[i].V=Complex::Rec2Polar(bus[i].V);

}

cout<<“====节点电压===============发电机发出功率======”<

for(i=0;i

{

cout<<“节点”<<(i+1)<<'t';

Complex::PrintfComplex(bus[i].V);

coutt(bus[i].Scal.real);

coutt(bus[i].Scal.image);

cout<

}

cout<<“======线路传输功率==========”<

for(i=0;i

{

int statemp,endtemp;

statemp=line[i].start;

endtemp=line[i].end;

cout<

Complex::PrintfComplex(Complex::ComDivRea(line[i].stoe,0.01));

Complex::PrintfComplex(Complex::ComDivRea(line[i].etos,0.01));

cout<

} }

return 0;}

#include “stdafx.h” #include #include #include

#include“Complex.h” #include“wanjing.h” #include“gauss.h” using namespace std;

int _tmain(int argc, _TCHAR* argv[]){

int i;

//i作为整个程序的循环变量

int N=Bus::ScanfBusNo();//输入节点个数

int L=Line::ScanflineNo();//输入支路个数

if((L&&N)==0){return 0;} //如果找不到两个文件中的任意一个，退出

Line *line=new Line[L];//动态分配支路结构体

Line::ScanfLineData(line);//输入支路参数

Line::PrintfLineData(line,L);//输出支路参数

Bus *bus=new Bus[N];//动态分配结点结构体

for(int i=0;i

bus[i].Sdelta.real=0;

bus[i].Sdelta.image=0;}

Bus::ScanfBusData(bus);//输入节点参数

Bus::PrintfBusData(bus,N);//输出结点参数

Complex **X;X=new Complex *[N];for(i=0;i

Bus::JisuanNodeDnz(X,line,bus,L,N);//计算节点导纳矩阵

Bus::PrintfNodeDnz(X,N);//输出节点导纳矩阵

int NN=(N-1)*2;double **JacAug;JacAug=new double *[NN];for(i=0;i

double *x;x=new double[NN];int count=1;

LOOP:

Bus::JisuanNodeI(X,bus,N);//计算节点注入电流

Bus::JisuanNodeScal(X,bus,N);//计算节点功率

Bus::JisuanNodeScal(X,bus,N);//计算节点功率

Bus::JisuanNodeSdelta(bus,N);//计算节点功率差值 Bus::PrintfNodeScal(X,bus,N);//输出节点功率差值

int icon=wehcon1(bus,N);//whether converbence看迭代是否结束

if(icon==1){

cout<<“icon=”<

Bus::JisuanJacAug(JacAug,X,bus,N);//计算雅可比增广矩阵

// Bus::PrintfJacAug(JacAug,N);

gauss::gauss_slove(JacAug,x,NN);//解方程组求出电压差值

Bus::ReviseNodeV(bus,x,N);//修正节点电压

// Bus::PrintfNodeV(bus,N);

count++;

goto LOOP;}

else

{

for(i=0;i

{

int statemp,endtemp;

Complex aa,bb,cc,dd,B;

B.real=0;

B.image=-line[i].B;

statemp=line[i].start;

endtemp=line[i].end;

aa=Complex::productComplex(Complex::getconj(bus[statemp-1].V), B);

bb=Complex::subComplex

(Complex::getconj(bus[statemp-1].V), Complex::getconj(bus[endtemp-1].V));

cc=Complex::productComplex(bb , Complex::getconj(line[i].Y));

dd=Complex::CaddC(aa,cc);

line[i].stoe=Complex::productComplex(bus[statemp-1].V,dd);

aa=Complex::productComplex(Complex::getconj(bus[endtemp-1].V), B);

bb=Complex::subComplex

(Complex::getconj(bus[endtemp-1].V), Complex::getconj(bus[statemp-1].V));

cc=Complex::productComplex(bb , Complex::getconj(line[i].Y));

dd=Complex::CaddC(aa,cc);

line[i].etos=Complex::productComplex(bus[endtemp-1].V,dd);

}

cout<<“icon=”<

Bus::JisuanNodeScal(X,bus,N);//计算节点功率

for(i=0;i

{

bus[i].Scal.real = bus[i].Scal.real + bus[i].Load.real;//发电机功率=注入功率+负荷功率

bus[i].Scal.image= bus[i].Scal.image+ bus[i].Load.image;

bus[i].V=Complex::Rec2Polar(bus[i].V);

}

cout<<“====节点电压===============发电机发出功率======”<

for(i=0;i

{

cout<<“节点”<<(i+1)<<'t';

Complex::PrintfComplex(bus[i].V);

coutt(bus[i].Scal.real);

coutt(bus[i].Scal.image);

cout<

}

cout<<“======线路传输功率==========”<

for(i=0;i

{

int statemp,endtemp;

statemp=line[i].start;

endtemp=line[i].end;

cout<

Complex::PrintfComplex(Complex::ComDivRea(line[i].stoe,0.01));

Complex::PrintfComplex(Complex::ComDivRea(line[i].etos,0.01));

cout<

} }

return 0;}

#include using namespace std;

class Complex//定义复数类 { public: double real;double image;int

RecPolar;//0表示直角坐标，1表示极坐标

static Complex CaddC(Complex c1,Complex c2);//求两个复数和

static Complex subComplex(Complex c1,Complex c2);//求两个复数差

static Complex productComplex(Complex c1,Complex c2);//求两个复数积

static Complex divideComplex(Complex c1,Complex c2);//求两个复数商

static Complex ComDivRea(Complex c1,double r2);//除数

static Complex getconj(Complex c1);//求一个复数共轭

static Complex getinverse(Complex c1);//取倒数

static double getComplexReal(Complex c1);//求一个复数实部

static double getCompleximage(Complex c1);//求一个复数虚部

static void PrintfComplex(Complex c1);//显示一个复数

static void PrintfmultiComplex(Complex C,int N);//显示多个复数

static void zeroComplex(Complex c1);//将复数复零

static Complex Rec2Polar(Complex c1);//取极坐标

Complex(){

RecPolar=0;} };

Complex Complex::Rec2Polar(Complex c1)//极坐标表示 { Complex Node;Node.real=sqrt(c1.real*c1.real+c1.image*c1.image);Node.image=atan2(c1.image,c1.real)*180/3.1415926;Node.RecPolar=1;return Node;}

Complex Complex::CaddC(Complex c1,Complex c2)//复数加法 {

Complex Node;

Node.real=c1.real+c2.real;

Node.image=c1.image+c2.image;

return Node;}

Complex Complex::subComplex(Complex c1,Complex c2)//复数减法 {

Complex Node;

Node.real=c1.real-c2.real;

Node.image=c1.image-c2.image;

return Node;}

Complex Complex::productComplex(Complex c1,Complex c2)//复数乘法 {

Complex Node;

Node.real=c1.real*c2.real-c1.image*c2.image;

Node.image=c1.image*c2.real+c2.image*c1.real;

return Node;} Complex Complex::divideComplex(Complex c1,Complex c2)//复数除法 {

Complex Node;

Node.real=(c1.real*c2.real+c1.image*c2.image)/(pow(c2.real,2)+pow(c2.image,2));Node.image=(c1.image*c2.real-c1.real*c2.image)/(pow(c2.real,2)+pow(c2.image,2));return Node;} Complex Complex::ComDivRea(Complex c1,double r1)//复数除数 { Complex Node;Node.real=c1.real/(r1);Node.image=c1.image/(r1);return Node;} Complex Complex::getconj(Complex c1)//取共轭 {

Complex Node;

Node.real=c1.real;Node.image=-c1.image;

return Node;}

Complex Complex::getinverse(Complex c1)//取倒数 { Complex Node;Node.real=1;Node.image=0;Node=(Complex::divideComplex(Node,c1));return Node;}

double Complex::getComplexReal(Complex c1)//取实部 {

return c1.real;}

double

Complex::getCompleximage(Complex c1)//取虚部 {

return c1.image;}

void

Complex::PrintfComplex(Complex c1)//按直角坐标输出 { if(c1.RecPolar==0){ cout.precision(6);

cout.width(8);

cout.setf(ios::right);

cout<

”;

cout.precision(6);

cout.width(8);

cout.setf(ios::left);

cout<

”;} else {

cout<

Complex::zeroComplex(Complex c1)//清零 { c1.real=0;c1.image=0;}

class gauss { public: static void gauss_slove(double **a,double *x,int NN);static void gauss_output();};

void gauss::gauss_slove(double **a,double *x,int NN){

int n,i,j,k,*pivrow,**pivcol;double eps,pivot,sum,aik,al;

n=NN;pivrow=new int[n];pivcol=new int *[n];

for(i=0;i

pivot= fabs(a[k][k]);

pivrow[k]=k;//行

pivcol[k][0]=k;pivcol[k][1]=k;//列n*2矩阵

for(i=k;i

{

for(j=k;j

{

if(pivot

{

pivot=fabs(a[i][j]);

pivrow[k]=i;//行

pivcol[k][1]=j;//列

}

}

}

if(pivot

{

cout<<“error”<

getchar();

exit(0);

}

if(pivrow[k]!=k)//行变换

{

for(j=k;j<(n+1);j++)

{

al=a[pivrow[k]][j];

a[pivrow[k]][j]=a[k][j];

a[k][j]=al;

}

}

if(pivcol[k][1]!=k)//列变换

{

for(i=0;i

{

al=a[i][pivcol[k][1]];

a[i][pivcol[k][1]]=a[i][k];

a[i][k]=al;

}

}

if(k!=(n-1))//将矩阵化为上三角形

式

{

for(i=(k+1);i

{

aik=a[i][k];

for(j=k;j<(n+1);j++)

{

a[i][j]-=aik*a[k][j]/a[k][k];

}

}

} } x[n-1]=a[n-1][n]/a[n-1][n-1];//解方程

for(i=(n-2);i>=0;i--){

sum=0;

for(j=(i+1);j

{

sum +=a[i][j]*x[j];0.182709

0.016894-0.0310701

-0.0402051 0.156702

-0.0355909-0.0668055

-0.00703229-0.0886481

-0.0129814-0.0390805

-0.0135062-0.1023

-0.0460568

-0.0342827

-0.00382402-0.102896

-0.0184062

}

x[i]=(a[i][n]-sum)/a[i][i];} for(k=(n-2);k>=0;k--){

al=x[pivcol[k][1]];

x[pivcol[k][1]]=x[pivcol[k][0]];

x[pivcol[k][0]]=al;}

cout<<“节点电压修正量”<

cout<

}

====节点功率计算值==== 0.935261

-0.159048 0.573909

0.0789973-0.00289889

-0.00796623-0.0791247

-0.0168362-0.436255

-0.0580392 0.0359139

-0.0106592-0.229118

-0.0885419-0.136179

-0.148207 0.0446243

0.0111298-0.0223764

-0.00695775-0.0237482

-0.198318

-5.24266e-015

-0.0354071

-0.0925078

-1.05629e-015

-0.0391348

0.014529

0.00158644

-0.0258771

-0.109514

icon=1进行第2次迭代 节点电压修正量

=================-0.00164889-0.000540034-0.00261067-0.00532027-0.00235315-0.00600971-0.00189677-0.00643874-0.0023631-0.00650659-0.00170949-0.0074907-0.00164545-0.00485415-0.00493977-0.0119042-0.00331285-0.0175611-0.00207908

-0.00347744-0.0869347-9.48909e-015-0.0110778-0.0538236-7.53784e-016-0.0168097 7.049e-005-0.00146487-0.00458276 0.00251645

-0.00336375-0.00530645-0.0147816-0.000326161-0.00640487-0.00251701-0.0169829-0.00175286-0.0174333-0.0239063

-0.0119192-0.076014

-0.0160104-0.441997

-0.0750285 0.000250012

3.72542e-005-0.228052

-0.108844-0.100078

-0.105634 0.000410707

0.000378067-0.057497

-0.0195879 0.200039

0.0582563-0.00307326-0.0163809-0.00232773-0.0175806 8.74293e-005-0.0192018 0.000558996-0.0197776-0.000247851-0.0193784-0.00115346-0.0185848-0.00127275-0.0186244-0.00010108-0.0188966 0.000553585-0.0200901-3.76315e-005-0.0208303 0.00308341-0.0219386-0.00195916-0.0205356-0.00184757-0.0076401 0.00197593-0.0245534 0.00434657-0.027534

====节点功率计算值==== 0.98623

-0.134163 0.583136

0.166278-0.111173

0.199792

-0.0621041

-0.0821379

-0.0350785

-0.0902383

-0.0320461

-0.0951562

-0.0220362

-0.175458

4.72557e-015

-0.0320661

-0.0871134

-7.03489e-017

-0.0350769

0.000273455

1.51804e-005

-0.0240417

-0.10604

icon=1进行第3次迭代 节点电压修正量

=================-2.67079e-005-2.30128e-006-2.20543e-005-6.00686e-005-2.33043e-005-6.85601e-005-3.22294e-005-2.61107e-005-2.80198e-005-6.6167e-005-2.34528e-005

-0.0739846 0.0227868-0.0158709-0.0248173-0.0179447-0.0578368-0.00890719-0.0337091-0.00693706-0.111601 1.21429e-014-0.0159145-0.0667319 9.24355e-016-0.0228592 7.10354e-005-6.6188e-006-0.00889343-0.0184098

-5.66132e-005-4.4646e-005-1.74668e-005-4.50947e-005-0.000181763-3.81763e-006-0.000286581-6.68993e-005-1.28441e-005-5.17172e-005-0.000223284-4.54717e-005-2.47586e-005 4.32335e-007-0.000258494 1.82635e-005-0.000272051-6.95195e-006-0.000251969 1.11318e-005-0.000279418 5.74737e-005-0.000307368 6.86998e-005-0.000320274 5.38112e-005-0.00031447 3.59531e-005-0.00030494 3.37607e-005-0.000307449 5.26532e-005-0.000310721 6.92761e-005-0.000350373 5.60942e-005-0.00040977 0.000123641-0.000440259 1.36149e-005-0.000426973-1.70227e-005-9.37794e-005 0.000113675-0.000544011 0.000176034-0.000636202

====节点功率计算值====

0.986878

-0.133979 0.583

0.167193-0.024

-0.012-0.076

-0.016-0.442

-0.0748606

1.43501e-008

1.07366e-008-0.228

-0.109

-0.0999999

-0.104049 4.51318e-008

8.98835e-008-0.0579999

-0.0199999 0.2

0.0591018-0.112

-0.0749997 0.2

0.0242519-0.062

-0.016-0.082

-0.025-0.035

-0.018

-0.0900001

-0.058-0.032

-0.00899997-0.095

-0.0339999-0.022

-0.00699998-0.175

-0.112

-6.07156e-015

-1.19217e-014-0.032

-0.016-0.087

-0.0669999

7.03078e-017

-9.23979e-016-0.035

-0.0229999

1.09492e-007

4.45699e-008 1.54958e-009

-2.01531e-010-0.024

-0.00899994-0.106

-0.0189996

icon=0,迭代结束。

====节点电压===============发电机发出功率======

节点1

1.05

0。

98.6878-13.3979

节点2

1.045

-1.846。

29.4193

节点3

1.02384-3.83352。

0

节

点25 1.01216-9.68486。

0

0 0 节点4

1.01637-4.55698。

0

节

点26 0.994393

-10.1089。

0 0

0 节点5

1.01

-6.48617。

节 点27 1.02012-9.42025。

0

11.5139 0 节点6

1.01332-5.38073。

0

节

点28 1.00992-5.86244。

0

0 0 节点7

1.00489-6.38368。

0

节

点29 1.00022-10.6579。

0

0 节点8 19.5951 节点9 0 节点10 0 节点11 5.91018 节点12 0 节点13 2.42519 节点14 0 节点15 0 节点16 0 节点17 0 节点18 0 节点19 0 节点20 0 节点21 0 节点22 0 节点23 0 节点24 0 1.01

-5.62974。

1.03905-6.78143。

1.03595-8.69362。

-4.5962。

1.04711-7.80323。

1.05

-6.34392。

1.03242-8.7401。

1.02788-8.86784。

1.03458-8.45044。

1.03051-8.83678。

1.01845-9.5141。

1.01604-9.70326。

1.02022-9.50938。

1.0237-9.17478。

1.02432-9.17024。

1.01802-9.36719。

1.01339-9.68362。

0 20

节 点30 0.988705

-11.5464。

0

0 0

======

线路传输功率========== 2to1

-57.7373

5.41674i

58.3454

0

-15.1827i

3to1

-39.659

-7.75964i

40.3424

1.78481i

4to2

-30.87

-9.74186i

31.4153

0

3.58352i

4to3

-37.0772

-7.78596i

37.259

6.55964i

5to2

-44.3717

-9.78456i

45.2968

0

4.84242i

6to2

-38.4766

-8.22625i

39.3252

0

2.87667i

6to4

-34.946

1.92384i

35.0885

0

-3.28202i

7to5

-0.16304

-6.41767i

0.171702

0

2.2985i

7to6

-22.637

-4.48233i

22.7745

0

1.44238i

8to6

-11.8939

-5.48098i

11.913

0

3.70557i

6to9

12.3737

-12.3826i

-12.3737

0

13.0033i

6to10

10.9107

-3.80907i

-10.9107

0

4.53223i

11to9

5.91018i

0

-5.08963i

10to9

-32.652

-2.3712i

32.652

0

3.46974i

4to12

23.5411

-11.5375i

-23.5411

0

13.2407i

13to12

2.42519i

1.05

-1.90978i 1.66484i 14to12

-7.9019

-2.06732i

7.97894

30to29

-3.6702

-0.542564i

3.70398

2.22749i 0.606393i 15to12

-18.254

-5.74885i

18.4835

28to8

-1.89152

-3.79982i

1.89395

6.20089i-4.9239i 16to12-7.53872

-2.90237i

7.59633

28to6

-14.7868

-2.82565i

14.8234

3.02352i 0.294601i 15to14-1.69544

-0.461488i

1.70189

请按任意键继续...0.467323i 17to16-4.03014 1.10238i 18to15-6.08074 1.46028i 19to18-2.87549 0.478389i 20to19

6.6418-2.93222i 20to10

-8.8418 3.85077i 17to10-4.96987 4.76656i 21to10-16.1562 9.42843i 22to10-7.87782 4.21401i 22to21

1.34443-2.01837i 23to15-5.59369 2.25006i 24to22-6.48186 2.08163i 24to23-2.38596 0.579814i 25to24-0.167617 0.281364i 26to25

-3.5 2.3674i 27to25

3.39433-2.08638i 28to27

16.1446 3.13006i 29to27-6.10398 1.67047i 30to27-6.92979-1.07089i-1.37839i-0.467767i

2.96679i-3.66679i-4.72911i-9.18162i-4.10132i

2.01969i-2.17981i-2.00141i-0.56401i

-0.28102i-2.29999i

2.11848i-2.10093i-1.50639i

-1.3574i

4.03872

6.12096

2.88074

-6.62452

8.9242

4.98423

16.2709

7.93248

-1.34378

5.62846

6.53339

2.39369

0.167814

3.54513

-3.37751

-16.1446

6.19083

7.09313

高等电力系统分析 IEEE30节点潮流程序

班级：电研114班

姓名：王大伟

学号：2201100151