CN102609598B - 一种对大型电力系统进行电磁暂态仿真的方法 - Google Patents

Abstract

本发明公开了一种对大型电力系统进行电磁暂态仿真的方法。首先对外部系统进行等值，具体为：按照第M个子外部无源网络到第一个子外部无源网络的顺序，依次获取各子导纳有理函数式，最终获得第一个子外部无源网络的子导纳有理函数式，作为外部系统的端口导纳有理函数式；通过外部系统的端口导纳有理函数式获取外部系统的等值电路；通过外部系统的等值电路对外部电源进行处理，获取等效电源。在获得实际系统的等值系统之后，对该等值系统进行电磁暂态仿真分析，通过仿真结果检测设备及元件的谐波和电磁暂态特性，判断元件、设备及局部网络的谐波和暂态分量是否合理，如果否，更换电气设备、调整更合理参数及采取相应措施；如果是，流程结束。

Description

一种对大型电力系统进行电磁暂态仿真的方法

技术领域

本发明涉及电力系统电磁暂态领域，特别涉及一种对大型电力系统进行电磁暂态仿真的方法。

背景技术

对实际电力系统进行电磁暂态仿真是研究电力系统故障或者操作过后可能出现的暂态过电压和过电流是否满足要求，以及相应电力元件、设备及局部网络的谐波和电磁暂态分量是否合理，从而可以选择及调整更合适的设备和参数，及采取相应的预防措施等，是电力系统规划、设计、运行、分析及控制的一项基本且重要的工作内容。若对实际电力系统进行全面的电磁暂态仿真分析，往往耗时很长。而对电力系统进行暂态分析，一般也只是对于其中的某一部分相关变量的暂态变化过程感兴趣；因此，为了提高电磁暂态仿真过程的计算效率，可以只对感兴趣的部分网络(本发明称之为内部系统)，进行详细的仿真，而对于其余部分网络(本发明称之为外部系统)，进行简化等值而在保证研究系统仿真精度的前提下，使得大系统变成小系统而大大加快其仿真过程，并使得基于仿真的设备运行的合理性、参数调整、及预防措施的获得也大大加快，而备受关注^[1]。

电力系统电磁暂态等值研究因其作用明显，使其成为一个备受关注的研究热点。文献^[1-4]所介绍的ward等值法利用戴维南定理进行等值，其优点是简单方便，易于实现，在对精度要求不高的情况下可以实现对外部系统的等值；其不足在于只能在伪稳态情况下对外部系统进行等值，当系统中包括有比较多的谐波，及系统发生不对称故障，使得谐波影响不可忽略时，其等值效果不理想。文献^[5-7]提出利用对外部系统阻抗函数的零极点匹配而求得相应外部系统的等值电路，文献^[8]在文献^[5-7]的基础上提出了利用混合矩阵并结合QZ方法获得外部系统阻抗函数精度更高的零极点，然后再通过零极点匹配以期获得精度更高的等值电路；但是由于比例缩放环节的限制，这种方法只能在一个较小频率范围内具有较好的等值效果；当频带较宽时，等值效果仍不够理想^[9，10]。

针对上述方法存在等值精度不足的问题，文献^[11-14]提出了采用基于导纳有理函数式的VF(Vector Fitting)等值方法，该方法为目前常用的等值方法，可使等值网络能够在较宽频段内对外部系统均具有较高的近似度，取得了较高的精度。

发明人在实现本发明的过程中，发现现有用于对外部系统进行等值的技术中至少存在以下的缺点和不足：

传统方法存在等值精度不足的问题，而VF等值方法需要对外部系统进行频率采样，并利用最小二乘法求解系数矩阵为稠密矩阵的超定线性方程组来获得有理函数式的相应系数，计算量较大；且所求得的极点可能出现不稳定极点，而为了提高极点的稳定性，该方法往往还需采用其它使等值系统稳定的措施^[15-19]。因此，VF方法虽然相比传统等值方法在精度上有很大提高，但是在求解时间上花费巨大，使得等值电路的求解时间较长，不能方便地对整个系统进行电磁仿真，从而不能快速地验证电力系统中设备的谐波和暂态分量是否合理，进而会产生不必要的经济损失。

参考文献

[1]陈水明，黄璐璐，谢海滨，余占清.互联系统电磁暂态交互作用研究，II：交流侧系统等值及稳态调试.高电压技术，2007，37(3)：537-547

[2]柳勇军.电力系统机电暂态和电磁暂态混合仿真技术的研究.北京：清华大学，2006.

[3]李承.电力系统电磁暂态与机电暂态程序的混合仿真研究.天津：天津大学，2007.

[4]孙宏斌；张伯明；相年德；Ward型等值的非线性误差分析与应用.电力系统自动化，2008，32(13)：11-15.

[5]Watson N R，Arrillaga J.Frequency-dependent A.C.System Equivalents forHarmonic Studies and Transient Convertor Simulation.IEEE Trans.on PowerDelivery，1988，3(3)：1196-1203

[6]Hingorani N G，Burbery M F.Simulation of AC System Impedance in HVDCSystem Studies.IEEE Trans.on Power APPARATUS And SYSTEMS.1970，PAS-89(5)：820-828

[7]Morched A S，Ottevangers J H，Martf L.Multi-Port Frequency DependentNetwork Equivalents for the EMTP.IEEE Trans.on Power Delivery，1993，8(3)：1328-1335

[8]Ibrahima A I，Salama M M A.Frequency dependent network equivalents forelectromagnetic transient studies.Electrical Power and Energy Systems，1999，21(6)：395-404.

[9]Andrea Ubolli，Gustavsen B.Comparison of Methods for Rational Approximationof Simulated Time-Domain Responses：ARMA，ZD-VF，and TD-VF.IEEE Trans.onPower Delivery，2011，26(1)：279-288

[10]Neville Watson，Jos Arrillaga.Power systems electromagnetic transientssimulation.London：The Institution of Electrical Engineers，2003.

[11]Gustavsen B，Semlyen A.Simulation of Transmission Line Transients usingVector Fitting and Modal Decomposition.IEEE Trans.on Power Delivery，1998，13(2)：605-614

[12]Gustavsen B，Semlyen A.Rational approximation of frequency domain responsesby vector fitting.IEEE Trans.on Power Delivery，1999，14(3)：1052-1061

[13]Gustavsen B.Computer Code for Rational Approximation of FrequencyDependent Adittance Matrices.IEEE Trans.on Power Delivery，2002，17(4)：1093-1098

[14]Abner Ramirez.Vector Fitting-Based Calculation of Frequency-DependentNetwork Equivalents by Frequency Partitioning and Model-Order Reduction.IEEETrans.on Power Delivery，2009，24(1)：410-415

[15]Gustavsen B，Semlyen A.Enforcing passivity for Admittance MatricesApproximated by Rational Functions.IEEE Trans.on Power Systems，2001，16(1)：97-104

[16]Coelho C P，Phillips J，Silveira L M.A convex programming approach forgenerating guaranteed passive approximations to tabulated frequency-data.IEEETrans.Computer-Aided Design Integr.Circuits Syst.2004，23(2)：293-301

[17]Saraswat D，Achar R，Nakhla M S.Enforcing passivity for rational functionbased macromodels of tabulated data.Electrical Performance of ElectronicPackaging，Oct.27-29，2003：pp.295-298

[18]Lamecki A，Mrozowski M.Equivalent SPICE circuits with guaranteed passivityfrom nonpassive models.IEEE Trans.Microw.Theory Tech，2007，55(3)：526-32

[19]Gustavsen B.Fast Passivity enforcement for Pole-Residue Models byPerturbation of Residue Matrix Eigenvalues.IEEE Trans.on Power Delivery，2008，23(4)：2278-2285

发明内容

本发明提供了一种对大型电力系统进行电磁暂态仿真的方法，在具有很高等值精度的前提下，本发明大幅度缩短了电磁暂态仿真时间，可方便快速地验证元件、设备及局部网络的谐波和暂态分量是否合理，进而有效地选择电气设备及采取相应措施，避免了不必要的经济损失，详见下文描述：

一种对大型电力系统进行电磁暂态仿真的方法，所述方法包括以下步骤：

(1)获取频域下的外部无源网络；

(2)对所述外部无源网络进行划分获取若干个子外部无源网络；

(3)按照第M个子外部无源网络到第一个子外部无源网络的顺序，依次获取各子导纳有理函数式，通过当前子外部无源网络与下一相邻子外部无源网络之间的端口，将当前子导纳有理函数式加入下一相邻子外部无源网络的导纳矩阵的自导和互导中，最终获得第一个子外部无源网络的子导纳有理函数式并作为外部系统的端口导纳有理函数式；

(4)通过所述外部系统的端口导纳有理函数式获取外部系统的等值电路；

(5)通过所述外部系统的等值电路对外部电源进行处理，获取等效电源；

(6)获得实际系统的等值系统，对所述等值系统进行电磁暂态仿真分析，通过仿真结果检测设备及元件的谐波和电磁暂态特性，判断元件、设备及局部网络的谐波和暂态分量是否合理，如果否，更换电气设备、调整更合理参数及采取相应措施；如果是，流程结束。

所述对所述外部无源网络进行划分获取若干个子外部无源网络具体为：

根据所述外部无源网络的结构，把所述外部无源网络划分成所述若干个子外部无源网络，按内部系统端口到外部系统末端的顺序依次编号，设共有M个子外部无源网络。

所述按照第M个子外部无源网络到第一个子外部无源网络的顺序，依次获取各子导纳有理函数式，通过当前子外部无源网络与下一相邻子外部无源网络之间的端口，将当前子导纳有理函数式加入下一相邻子外部无源网络的导纳矩阵的自导和互导中，最终获得第一个子外部无源网络的子导纳有理函数式并作为外部系统的端口导纳有理函数式，其具体过程为：

1)通过所述子外部无源网络获取混合矩阵；

2)通过所述混合矩阵获取子导纳有理函数式Y^r(s)；

3)通过所述子导纳有理函数式Y^r(s)依次求解最终获取第一个子外部无源网络的子导纳有理函数式Y¹(s)。

所述通过所述子外部无源网络获取混合矩阵具体为：

获得频域下的外部系统无源网络后，在子外部系统的某一连接端口处注入一个单位的电流，其余端口开路，得到式(1)：

Y^rU^r＝I^r _in put(1)

其中，Y^r为频域下的第r个子外部系统的节点导纳矩阵；U^r为第r个子外部系统和端口的各节点电压；I^r _in put为第r个子外部系统和端口的注入电流列向量，对应注入单位电流节点的元素为1，其余元素都为0；

若电压观测节点的电压用F^r表示，则有：

F^r＝(d^r)^TU^r    (2)

其中d为节点关联矩阵，即与电压观测节点相对应的位置元素为1，其余元素为0；

将式(1)与式(2)合并即可得如式(3)所示的混合矩阵：

$\left[\begin{array}{cc}{Y}^r& 0\\ -{\left(d^r\right)}^T& 1\end{array}\right]\left[\begin{array}{c}{U}^r\\ F^r\end{array}\right]=\left[\begin{array}{c}{I^r}_{\mathrm{input}}\\ 0\end{array}\right]---\left(3\right)$

利用克莱姆法则可求得F^r为：

$F^r=\frac{\det \left[\begin{array}{cc}{Y}^r& {I^r}_{\mathrm{imput}}\\ -{\left(d^r\right)}^T& 0\end{array}\right]}{\det \left[\begin{array}{cc}{Y}^r& 0\\ -{\left(d^r\right)}^T& 1\end{array}\right]}=\frac{\det H_1}{\det H}---\left(4\right)$

因式(4)中Y_r的各非零元素均为s的多项式，detH₁及detH均为s的多项式，其中，det为矩阵的行列式，T为矩阵转置。

所述通过所述混合矩阵获取子导纳有理函数式Y^r(s)具体为：

$Y^r\left(s\right)=\frac{I^r}{U^r}=\frac{1}{F^r}=\frac{\det H}{\det H_1}---\left(5\right)$

若注入电流节点和电压观测节点为同一节点，则Y^r(s)为自导纳；否则，Y^r(s)为两节点间的互导纳；

矩阵H和H₁由

求得，

$B_H^1={\left[\begin{array}{ccccc}{Y}_{\mathrm{1,1}}& ...& Y_{1,j}& ...& Y_{1,n}\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ Y_{i,1}& ...& Y_{i,j}& ...& Y_{i,n}\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ Y_{n,1}& ...& Y_{n,j}& ...& Y_{n,n}\end{array}\right]}_{(n+1)\×(n+1)}---\left(6\right)$

①对于式(7)所示的多项式：

a₀+a₁s+a₂s²+…+a_nsⁿ    (7)

记为式(8)的向量形式称为向量多项式：

[a₀，a₁，a₂，…，a_n](8)

②定义

为向量多项式运算中的加法与乘法，其运算规则分别为：

$[b_0,...,b_n]\⊕[c_0,...,c_m]=[b_0+c_0,...,b_n+c_n,c_{n+1}...,c_m],n\≤m---\left(9\right)$

$[b_0,...,b_n]\⊗[c_0,...,c_m]=[b_0{c}_0,(c_1{c}_0+b_0{c}_1),...,(b_{n-2}{c}_m+b_{n-1}{c}_{m-1}+b_n{c}_{m-2}),(b_{n-1}{c}_m+b_n{c}_{m-1}),b_n{c}_m]---\left(10\right)$

式(6)中的节点自导纳Y_ii与互导纳Y_ij如式(11)所示：

$\left\{\begin{array}{c}{Y}_{\mathrm{ii}}=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{1}{R_{\mathrm{ij}}+{\mathrm{sL}}_{\mathrm{ij}}}\right)+(G_{\mathrm{ii}}+{\mathrm{sC}}_{\mathrm{ii}})\\ Y_{\mathrm{ij}}=-\underset{j\≠i}{\overset{n}{\underset{j=1}{\mathrm{\Σ}}}}\frac{1}{R_{\mathrm{ij}}+{\mathrm{sL}}_{\mathrm{ij}}}\end{array}\right.---\left(11\right)$

式中R_ij、L_ij、C_ii、G_ii分别表示支路ij的电阻、电感及节点i的对地总导纳与总电容；n为系统节点个数，把Y_ii及Y_Ij采用向量多项式进行表示

$\left\{\begin{array}{c}{Y}_{\mathrm{ii}}=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{[R_{\mathrm{ij}},L_{\mathrm{ij}}]}+[G_{\mathrm{ii}},C_{\mathrm{ii}}]\\ Y_{\mathrm{ij}}=-\underset{j\≠i}{\overset{n}{\underset{j=1}{\mathrm{\Σ}}}}\frac{1}{[R_{\mathrm{ij}},L_{\mathrm{ij}}]}\end{array}\right.---\left(12\right)$

将

中的各个元素表示成式(12)的向量形式；然后对各个非零元素按照定义的运算规则进行下一步操作，以使得每个非零元素的分子、分母均为以s为变量的向量多项式；

将矩阵

简记为B¹，其行列式简记为det(B¹)；

变量说明：k为展开层数；矩阵为N_k×N_k的方阵；m_k为矩阵

列号(1≤m_k≤N_k)，

为矩阵

中第i行第m_k列元素，

为矩阵

中第1行第m_k列元素

对应的伴随矩阵；S_k、h_k为计算矩阵行列式的中间变量；

变量初始化：k＝1，m₀＝1，

N₁＝n，

①判断是否N_k＝0，如果是，

执行步骤⑥；如果否，执行步骤②；

②判断是否N_k＝1，如果是，

执行步骤⑥；如果否，执行步骤③；

③判断是否N_k＝2，如果是，

执行步骤⑥；如果否，执行步骤④；

④N_k＞2，按行展开求取行列式，即令m_k＝1，定义S_k＝0，h_k＝1，并进行如下操作：

1)判断是否满足m_k≤N_k且

非零，如果是，执行步骤⑤；如果否，执行步骤3)；

2)k＝k-1，计算

h_k＝-h_k，m_k＝m_k+1，转步骤1)；

3)判断是否满足m_k≤N_k且为零，如果是，h_k＝-h_k，m_k＝m_k+1，转步骤1)；如果否，转步骤4)；

4)判断是否k＝1，如果是，则det(B1)＝S_k，执行步骤⑥；如果否，k＝k-1，

m_k＝m_k+1，转步骤1)；

⑤求取

的伴随矩阵，k＝k+1；去掉

的第m_k-1列和第一行即得到N_k＝N_k-1-1；转步骤①；

⑥判断是否k＝1，如果是，则计算矩阵H行列式，detH＝detB¹，流程结束；如果否，则转步骤④中的2)；

将detH、detH₁代人式(5)可得Y(s)：

$Y\left(s\right)=\det H\⊗(1/\det H_1)---\left(13\right)$

形式为：

$Y\left(s\right)=\frac{N\left(s\right)}{D\left(s\right)}=\frac{a_0+a_{1}s+a_2{s}^2+...+a_n{s}^m}{1+b_1{s}^1+b_2{s}^2+...+b_n{s}^n}---\left(14\right).$

所述通过所述子导纳有理函数式Y^r(s)依次求解最终获取第一个子外部无源网络的子导纳有理函数式Y¹(s)具体为：

(I)若已求得第r(2≤r≤M)个子网络端口q₁、q₂看入的自导纳有理函数式及互导纳有理函数式

则与之相邻的第r-1个小网络的导纳矩阵为：

$\left[\begin{array}{ccccc}.& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ ...& Y_{q_1}^{\′}\left(s\right)+Y_{q_1}^r\left(s\right)& ...& Y_{q_{1,}{q}_2}^{\′}\left(s\right)-Y_{q_1,q_2}^r\left(s\right)& ...\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ ...& Y_{q_1,q_2}^{\′}\left(s\right)-Y_{q_1,q_2}^r\left(s\right)& ...& Y_{q_2}^{\′}\left(s\right)+Y_{q_2}^r\left(s\right)& ...\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\end{array}\right]$ 式中

分别表示第r-1个小网络中q₁、q₂节点的自导及互导；

(II)第1子系统保留有理函数式的全部极点，从第2个子系统到第M个子系统，依次降低有理式中的高频极点的个数，并增加常数项的值；

按照由第M个子外部无源网络到第一个子外部无源网络的顺序依次求取，直到求得第一个子系统的等值自导及互导的有理式Y¹(s)，即获取到了整个外部系统的端口导纳有理式Y(s)，即Y(s)＝Y¹(s)。

所述通过所述外部系统的端口导纳有理函数式获取外部系统的等值电路具体为：

①Y(s)没有多重实根极点；

$Y\left(s\right)=\frac{N\left(s\right)}{D\left(s\right)}=\frac{a_0+a_{1}s+a_2{s}^2+...+a_n{s}^m}{\underset{i=1}{\overset{q}{\mathrm{\Π}}}(s-p_i)\underset{k=1}{\overset{r}{\mathrm{\Π}}}\left[(s-p_k^{\′}-p_k^{\′\′})(s-p_k^{\′}+p_k^{\′\′})\right]}---\left(15\right)$

式中，p_i为D(s)的实根，p′_k±jp″_k为D(s)的共轭复根，进行部分分式展开而得到式(16)：

$Y\left(s\right)=d+\mathrm{se}+\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\frac{c_i}{s-p_i}+\underset{k=1}{\overset{r}{\mathrm{\Σ}}}(\frac{c_k^{\′}+{\mathrm{jc}}_k^{\′\′}}{s-(p_k^{\′}+{\mathrm{jp}}_k^{\′\′})}+\frac{c_k^{\′}-{\mathrm{jc}}_k^{\′\′}}{s-(p_k^{\′}-{\mathrm{jp}}_k^{\′\′})})---\left(16\right)$

式中，d、e由多项式除法确定；系数c_i、c′_k、c″_k分别由式(17)、(18)获得：

c_i＝[(s-p_i)Y(s)]s＝p_i    (17)

$\left\{\begin{array}{c}{c}_k^{\′}={\left[(s-p_k^{\′}-{\mathrm{jp}}_k^{\′\′})Y\left(s\right)\right]}_{s=p_k^{\′}+{\mathrm{jp}}_k^{\′\′}}\\ c_k^{\′\′}={\left[(s-p_k^{\′}+{\mathrm{jp}}_k^{\′\′})Y\left(s\right)\right]}_{s=p_k^{\′}-{\mathrm{jp}}_k^{\′\′}}\end{array}\right.---\left(18\right)$

从而，按照式(19)至(25)所示的对应关系，即可求出式(16)各项相对应的简单支路，再进行并联，获得相应的等值电路；

常数项和一次项对应R₀和C₀：

C₀＝e，R₀＝1/d    (19)

实根对应RL支路：

R₁＝-p_i/c_i，L₁＝1/c_i    (20)

共轭复根：

$\frac{c_k^{\′}+{\mathrm{jc}}_k^{\′\′}}{s-(p_k^{\′}+{\mathrm{jp}}_k^{\′\′})}+\frac{c_k^{\′}-{\mathrm{jc}}_k^{\′\′}}{s-(p_k^{\′}-{\mathrm{jp}}_k^{\′\′})}---\left(21\right)$

与以下RLC支路相对应：

L＝1/(2c′_k)   (22)

R＝(-2p′_k+2(c′_kp′_k+c″_kp″_k)L)L    (23)

1/C＝(p′_k ²+p″_k ²+2(c′_kp′_k+c″_kp″_k)R)L    (24)

G＝-2(c′_kp′_k+c″_kp″_k)CL    (25)

②Y(s)含有多重实根极点；

将多重实根极点看作单重实根极点，只须获得与单重实根极点相对应的等值电路。

所述通过所述外部系统的等值电路对外部电源进行处理，获取等效电源具体为：

1)计算各端口与地之间以及各端口之间的等值电路的阻抗，分别记为Z_ii和Z_ij，j≠i；

2)等值前外部系统工频下各端口处电压、电流的幅值及相角为：

根据戴维南等值原理分别获取单端口外部系统的等效电源和多端口外部系统的等效电源；

单端口外部系统的等效电源 ${\stackrel{\·}{E}}_i={\stackrel{\·}{U}}_i-{\stackrel{\·}{I}}_i{Z}_{\mathrm{ii}}---\left(26\right)$

多端口外部系统的等效电源 ${\stackrel{\·}{E}}_i={\stackrel{\·}{U}}_i-{\stackrel{\·}{I}}_i{Z}_{\mathrm{ii}}+\underset{j\≠i}{\overset{m}{\underset{j=1}{\mathrm{\Σ}}}}\frac{{\stackrel{\·}{U}}_i-{\stackrel{\·}{U}}_j}{Z_{\mathrm{ij}}}{Z}_{\mathrm{ii}}(i=\mathrm{1,2}...,m)---\left(27\right)$

本发明提供的技术方案的有益效果是：本方法首先从外部系统的混合矩阵出发，设计有效方法把外部系统端口的导纳混合矩阵函数行列式直接转成有理式函数；通过把大系统分解为系列小系统，并逐个快速地获得各个子系统的有理函数式而获得外部系统的完整有理函数式；通过有理函数式获取外部系统的等值电路和等效电源；并对等值电路、等效电源和内部系统所组成的实际系统的等值系统进行电磁仿真；通过试验验证，本仿真方法所采取的等值方法明显缩短了获得外部等值系统的计算时间，具有很高的等值精度，同时也大大缩短了对实际大电力系统进行电磁暂态仿真时间；基于此，可方便快速地验证电力系统中元件、设备甚至局部网络的谐波和暂态分量是否合理及是否存在隐患，而可以有效地选择电气设备及采取有效措施，避免了大电力系统不必要的经济损失，而提高实际大电力系统运行的安全性及供电可靠性。

附图说明

图1为本发明提供的一种对大型电力系统进行电磁暂态仿真的方法的流程图；

图2为本发明提供的外部系统等值的示意图；

图3为本发明提供的外部大系统端口导纳有理函数式分布求解示意图；

图4为本发明提供的相邻子外部无源网络的等值关系图；

图5为本发明提供的注入电流节点和电压观测节点的示意图；

图6为本发明提供的RLC等值电路模型示意图；

图7为本发明提供的单端口含电源等值电路模型示意图；

图8为本发明提供的多端口含电源等值电路模型示意图；

图9(a)和图9(b)为本发明提供的IEEE118节点系统3号节点电压曲线对比图；

图10为本发明提供的IEEE118节点外部系统划分方案示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

为了缩短电磁暂态仿真时间，方便快速地验证电力系统中元件、设备及局部网络的谐波和暂态分量是否合理，进而有效地选择电气设备及采取相应措施，避免经济损失，本发明实施例提供了一种对大型电力系统进行电磁暂态仿真的方法，参见图1，详见下文描述：

101：获取频域下的外部无源网络；

其中，参见图2，将外部系统中电源置零，具体为电压源接地和电流源开路；若外部系统中包含有用分布参数表示的长线路，可用多个π型等值电路的串联来表示。对外部系统中的各元件进行拉普拉斯变换，从而得到频域下的外部无源网络。

102：对外部无源网络进行划分获取若干个子外部无源网络；

对于实际大电力系统要被等值的外部无源网络一般比较大，往往包括有众多的电源，节点数小则数百个，多则数千个，计算量过大。因此，本方法提出了首先将外部无源网络进行划分，分步求解各网络的端口导纳有理函数式，从而使得计算量得到明显的降低。

其中，参见图3，该步骤具体为：根据外部无源网络的结构，把外部无源网络划分成若干个子外部无源网络，按内部系统端口到外部系统末端的顺序依次编号，设共有M个子外部无源网络，划分时，尽量使得这些子外部无源网络之间的端口数较少。M的数量根据外部无源网络的结构和实际应用中的需要确定，优选每个子外部无源网络中的节点个数较少，具体实现时，本发明实施例对此不做限制。

103：按照第M个子外部无源网络到第一个子外部无源网络的顺序，依次获取各子导纳有理函数式，通过当前子外部无源网络与下一相邻子外部无源网络之间的端口，将当前子导纳有理函数式加入下一相邻子外部无源网络的导纳矩阵的自导和互导中，最终获得第一个子外部无源网络的子导纳有理函数式并作为外部系统的端口导纳有理函数式；

参见图4，例如：当前子外部无源网络为第r(2≤r≤M)个子外部无源网络，下一相邻子外部无源网络为第r-1个子外部无源网络，通过端口q₁、q₂将第r个子导纳有理函数式加入第r-1个子外部无源网络的导纳矩阵的自导

和互导

中。

在求取端口节点间自导时，该端口节点同时作为注入电流节点和观测节点；在求取端口两节点间互导时，任取二者之一作为注入电流节点，则另一节点为电压观测节点。

以求取第r(2≤r≤M)个子外部无源网络的端口导纳有理函数式为例，具体步骤如下：

1)通过子外部无源网络获取混合矩阵；

获得频域下的外部系统无源网络后，在子外部系统的某一连接端口处(例如：图2中的端口i)注入一个单位的电流，其余端口开路，得到式(1)：

Y^rU^r＝I^r _in put(1)

其中，Y^r为频域下的第r个子外部系统的节点导纳矩阵；U^r为第r个子外部系统和端口的各节点电压；I^r _in put为第r个子外部系统和端口的注入电流列向量，对应注入单位电流节点的元素为1，其余元素都为0。

若电压观测节点(如图2中端口j)的电压用F^r表示，则有：

F^r＝(d^r)^TU^r    (2)

其中d为节点关联矩阵，即与电压观测节点相对应的位置元素为1，其余元素为0。

参见图5，例如1、2和7节点作为内部系统，3、4、5和6节点作为外部系统，则可以看出，2、7节点为端口节点。将2节点作为注入电流节点，7节点作为电压观测节点，则有：

$\left[\begin{array}{cccccc}{Y}_{\mathrm{2,2}}& Y_{\mathrm{2,3}}& Y_{\mathrm{2,4}}& Y_{2,5}& Y_{\mathrm{2,6}}& Y_{\mathrm{2,7}}\\ Y_{\mathrm{3,2}}& Y_{\mathrm{3,3}}& Y_{\mathrm{3,4}}& Y_{\mathrm{3,5}}& Y_{\mathrm{3,6}}& Y_{\mathrm{3,7}}\\ Y_{\mathrm{4,2}}& Y_{\mathrm{4,3}}& Y_{\mathrm{4,4}}& Y_{\mathrm{4,5}}& Y_{\mathrm{4,6}}& Y_{\mathrm{4,7}}\\ Y_{\mathrm{5,2}}& Y_{\mathrm{5,3}}& Y_{\mathrm{5,4}}& Y_{\mathrm{5,5}}& Y_{\mathrm{5,6}}& Y_{\mathrm{5,7}}\\ Y_{\mathrm{6,2}}& Y_{\mathrm{6,3}}& Y_{\mathrm{6,4}}& Y_{\mathrm{6,5}}& Y_{6,6}& Y_{\mathrm{6,7}}\\ Y_{7,2}& Y_{\mathrm{7,3}}& Y_{\mathrm{7,4}}& Y_{7,5}& Y_{7,6}& Y_{\mathrm{7,7}}\end{array}\right]\left[\begin{array}{c}{U}_2\\ U_3\\ U_4\\ U_5\\ U_6\\ U_7\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right]$ $F=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{U}_2\\ U_3\\ U_4\\ U_5\\ U_6\\ U_7\end{array}\right]$

将式(1)与式(2)合并即可得如式(3)所示的混合矩阵：

$\left[\begin{array}{cc}{Y}^r& 0\\ -{\left(d^r\right)}^T& 1\end{array}\right]\left[\begin{array}{c}{U}^r\\ F^r\end{array}\right]=\left[\begin{array}{c}{I^r}_{\mathrm{input}}\\ 0\end{array}\right]---\left(3\right)$

利用克莱姆法则可求得F^r为：

$F^r=\frac{\det \left[\begin{array}{cc}{Y}^r& {I^r}_{\mathrm{imput}}\\ -{\left(d^r\right)}^T& 0\end{array}\right]}{\det \left[\begin{array}{cc}{Y}^r& 0\\ -{\left(d^r\right)}^T& 1\end{array}\right]}=\frac{\det H_1}{\det H}---\left(4\right)$

因式(4)中Y^r的各非零元素均为s的多项式，故求得的detH₁及detH均为s的多项式，其中，det为矩阵的行列式，T为矩阵转置。

2)通过第r个子外部系统的混合矩阵获取第r个子外部系统的子导纳有理函数式Y^r(s)；

其中， $Y^r\left(s\right)=\frac{I^r}{U^r}=\frac{1}{F^r}=\frac{\det H}{\det H_1}---\left(5\right)$

若注入电流节点和电压观测节点为同一节点，则Y^r(s)为自导纳；否则，Y^r(s)为两节点间的互导纳。

当H(s)、H₁(s)的各个元素均为分式多项式函数时，其问题的难点在于如何合理地存储多项式以及对分式多项式函数的计算。本方法提出了基于向量技术的求解方法，即用两个一维向量分别存储分式中分子、分母的系数，并按多项式的运算规则对系数进行计算，从而可排除多项式中变量的影响，以计算detH为例，其过程如下：

若外部系统的节点数为n，则式(5)中的H结构为：

$H={\left[\begin{array}{cccccc}{Y}_{\mathrm{1,1}}& ...& Y_{1,j}& ...& Y_{1,n}& 0\\ .& & .& & .& .\\ .& & .& & .& .\\ .& & .& & .& .\\ Y_{i,1}& ...& Y_{i,j}& ...& Y_{i,n}& 0\\ .& & .& & .& .\\ .& & .& & .& .\\ .& & .& & .& .\\ Y_{n,1}& ...& Y_{n,j}& ...& Y_{n,n}& 0\\ 0& 0& -1& 0& 0& 1\end{array}\right]}_{(n+1)\×(n+1)}---\left(6\right)$

不妨令：

$B_H^1={\left[\begin{array}{ccccccc}{Y}_{\mathrm{1,1}}& ...& Y_{1,i}& ...& Y_{1,j}& ...& Y_{1,n}\\ .& & .& & .& & .\\ .& & .& & .& & .\\ .& & .& & .& & .\\ Y_{i,1}& ...& Y_{i,i}& ...& Y_{i,j}& ...& Y_{i,n}\\ .& & .& & .& & .\\ .& & .& & .& & .\\ .& & .& & .& & .\\ Y_{j,1}& ...& Y_{j,i}& ...& Y_{j,j}& ...& Y_{j,n}\\ .& & .& & .& & .\\ .& & .& & .& & .\\ .& & .& & .& & .\\ Y_{n,1}& ...& Y_{n,i}& ...& Y_{n,j}& ...& Y_{n,n}\end{array}\right]}_{(n\&times;n)}(i<j)---\left(7\right)$

因此，若求得的多项式，也就求得detH的多项式，为了实现对

多项式的求取，本方法引入如下定义及操作运算规则：

①对于式(8)所示的多项式：

a₀+a₁s+a₂s²+…a_nsⁿ    (8)

把其记为式(9)的向量形式，本方法称其为向量多项式：

[a₀，a₁，a₂，…，a_n](9)

②定义

为向量多项式运算中的加法与乘法，其运算规则分别为(设n≤m)：

$[b_0,...,b_n]\&CirclePlus;[c_0,...,c_m]=[b_0+c_0,...,b_n+c_n,c_{n+1}...,c_m]---\left(10\right)$

$[b_0,...,b_n]\&CircleTimes;[c_0,...,c_m]=[b_0{c}_0,(c_1{c}_0+b_0{c}_1),...,(b_{n-2}{c}_m+b_{n-1}{c}_{m-1}+b_n{c}_{m-2}),(b_{n-1}{c}_m+b_n{c}_{m-1}),b_n{c}_m]---\left(11\right)$

式(7)中的节点自导纳Y_ii与互导纳Y_ij如式(12)所示：

$\left\{\begin{array}{c}{Y}_{\mathrm{ii}}=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\frac{1}{R_{\mathrm{ij}}+{\mathrm{sL}}_{\mathrm{ij}}}\right)+(G_{\mathrm{ii}}+{\mathrm{sC}}_{\mathrm{ii}})\\ Y_{\mathrm{ij}}=-\underset{j\&NotEqual;i}{\overset{n}{\underset{j=1}{\mathrm{\&Sigma;}}}}\frac{1}{R_{\mathrm{ij}}+{\mathrm{sL}}_{\mathrm{ij}}}\end{array}\right.---\left(12\right)$

式中R_ij、L_ij、C_ii、G_ii分别表示支路ij的电阻、电感及节点i的对地总导纳与总电容(总电容等于本节点对地电容与本节点相连的各条支路的对地电容一半之和)；n为系统节点个数。把Y_ii及Y_ij采用向量多项式进行表示，则可以改写为如下的向量形式(以下的向量的含义均与此类似)：

$\left\{\begin{array}{c}{Y}_{\mathrm{ii}}=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{1}{[R_{\mathrm{ij}},L_{\mathrm{ij}}]}+[G_{\mathrm{ii}},C_{\mathrm{ii}}]\\ Y_{\mathrm{ij}}=-\underset{j\&NotEqual;i}{\overset{n}{\underset{j=1}{\mathrm{\&Sigma;}}}}\frac{1}{[R_{\mathrm{ij}},L_{\mathrm{ij}}]}\end{array}\right.---\left(13\right)$

对

矩阵中各元素进行如下预处理：首先把

中的各个元素表示成式(13)的向量形式；然后对各个非零元素按照定义的向量多项式运算法则进行下一步操作，以使得每个非零元素的分子、分母均为以s为变量的向量多项式。

下面进一步给出按行展开求取

的递归计算方法。为了表示方便，将矩阵

简记为B¹，其行列式简记为det(B¹)：(以下

分别表示向量多项式运算中的加法与乘法；“+”、“×”分别表示数值与向量多项式以及数值间运算中的加法与乘法)。

变量说明：k为展开层数；矩阵

为N_k×N_k的方阵；m_k为矩阵

列号(1≤m_k≤N_k)，

为矩阵

中第i行第m_k列元素，

为矩阵

中第1行第m_k列元素对应的伴随矩阵；S_k、h_k为计算矩阵行列式的中间变量。

变量初始化：k＝1，m₀＝1，

N₁＝n，

①判断是否N_k＝0，如果是，

执行步骤⑥；如果否，执行步骤②；

②判断是否N_k＝1，如果是，

执行步骤⑥；如果否，执行步骤③；

③判断是否N_k＝2，如果是，

执行步骤⑥；如果否，执行步骤④；

④N_k＞2，按行展开求取行列式，即令m_k＝1，定义S_k＝0，h_k＝1，并进行如下操作：

1)判断是否满足m_k≤N_k且

非零，如果是，执行步骤⑤；如果否，执行步骤3)；

2)k＝k-1，计算

h_k＝-h_k，m_k＝m_k+1，转步骤1)；

3)判断是否满足m_k≤N_k且

为零，如果是，h_k＝-h_k，m_k＝m_k+1，转步骤1)；如果否，转步骤4)；

4)判断是否k＝1，如果是，则det(B¹)＝S_k，执行步骤⑥；如果否，

k＝k-1，

m_k＝m_k+1，转步骤1)。

⑤求取

的伴随矩阵，k＝k+1；去掉

的第m_k-1列和第一行即得到

N_k＝N_k-1-1；转步骤①。

⑥判断是否k＝1，如果是，则计算矩阵H行列式，detH＝detB¹，流程结束；如果否，则转步骤④中的2)。

通过上述方法即可计算出多项式矩阵的函数行列式。对于阶数为n的矩阵H，求取其函数行列式的多项式函数表达式的计算复杂性为o(n！)。由于实际电力系统网络的稀疏性很高，一般可以达到90％以上，因此，单纯行列式的计算复杂性大约为o(0.1ⁿn！)。但行列式计算的每一步均涉及到多项式运算，多项式运算的计算复杂性与参与运算的多项式阶数有关，而多项式阶数又随系统储能元件和节点数的增大而增大。因此，当n比较小，例如小于15时，计算速度很快；但当n很大时，计算量显然过大，故进一步采用分步等值的措施，以减少计算复杂性，加快计算速度。

将detH、detH₁代人式(5)可得第r个子外部系统的子导纳有理函数式Y^r(s)：

$Y^r\left(s\right)=\det H\&CircleTimes;(1/\det H_1)---\left(14\right)$

此时Y^r(s)的分母及分子均是以向量形式表示的多项式，按上述向量多项式的定义，将它们转化为一般意义下的多项式。其形式为：

$Y^r\left(s\right)=\frac{N\left(s\right)}{D\left(s\right)}=\frac{a_0+a_{1}s+a_2{s}^2+...+a_n{s}^m}{1+b_1{s}^1+b_2{s}^2+...+b_n{s}^n}---\left(15\right)$

式(15)即为子外部系统注入电流节点与电压观测节点之间(或端口之间)的导纳有理函数式。

3)通过第r个子外部系统的子导纳有理函数式Y^r(s)依次求解最终获取第一个子外部无源网络的子导纳有理函数式Y¹(s)。

在由已求得的子导纳有理函数式Y^r(s)求解其相邻子网络的导纳有理函数式Y^r-1(s)时，本方法用到了以下两个措施：

(I)若已求得第r(2≤r≤M)个子网络端口q₁、q₂看入的自导纳有理函数式

及互导纳有理函数式则与之相邻的第r-1个小网络的导纳矩阵为：

$\left[\begin{array}{ccccc}.& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ ...& Y_{q_1}^{\&prime;}\left(s\right)+Y_{q_1}^r\left(s\right)& ...& Y_{q_{1,}{q}_2}^{\&prime;}\left(s\right)-Y_{q_1,q_2}^r\left(s\right)& ...\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ ...& Y_{q_1,q_2}^{\&prime;}\left(s\right)-Y_{q_1,q_2}^r\left(s\right)& ...& Y_{q_2}^{\&prime;}\left(s\right)+Y_{q_2}^r\left(s\right)& ...\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\end{array}\right]$ 式中

分别表示第r-1个小网络中q₁、q₂节点的自导及互导。

(II)由于外部系统中包含元件较多，其总体导纳有理函数式的阶数会很高，相应地有理式求解的计算量会比较大。本方法采取舍去离内部系统边界节点较远的子系统导纳有理函数式中高频极点的对应项，并增加常数项的值来弥补由此产生工频下导纳值的差别来达到降低有理式阶数的目的(即舍去式(17)中极点p_i、p′_k±jp″_k较大的项，并增加d值来保证工频时Y(s)值不变)。具体做法是第1子系统保留有理函数式的全部极点，从第2个子系统到第M个子系统，依次降低有理式的极点个数，从而使得求解各个子系统有理式的计算量得到降低。

每一子网络的等值端口自导及互导均按上述过程进行计算，按照由第M个子外部无源网络到第一个子外部无源网络的顺序依次求取，直到求得第一个(也就是与内部网络直接相连)子系统的等值自导及互导的有理式Y¹(s)，从而也就求得了整个外部系统的端口导纳有理式Y(s)，即Y(s)＝Y¹(s)。

104：通过外部系统的端口导纳有理函数式获取外部系统的等值电路；

对于式(15)所示的导纳有理函数式，其等值电路的求解本方法分两种情况加以讨论：①Y(s)没有多重实根极点；②Y(s)有多重实根极点。

①Y(s)没有多重实根极点

意味着式(15)的分母D(s)中至多包含q个单实根和r个共轭复根(q+2r＝n)。与Y(s)相对应的简单网络的求解过程为：

$Y\left(s\right)=\frac{N\left(s\right)}{D\left(s\right)}=\frac{a_0+a_{1}s+a_2{s}^2+...+a_n{s}^m}{\underset{i=1}{\overset{q}{\mathrm{\&Pi;}}}(s-p_i)\underset{k=1}{\overset{r}{\mathrm{\&Pi;}}}\left[(s-p_k^{\&prime;}-p_k^{\&prime;\&prime;})(s-p_k^{\&prime;}+p_k^{\&prime;\&prime;})\right]}---\left(16\right)$

式中，p_i为D(s)的实根，p′_k±jp″_k为D(s)的共轭复根。根据留数定理，将式(16)进行部分分式展开而得到式(17)：

$Y\left(s\right)=d+\mathrm{se}+\underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}\frac{c_i}{s-p_i}+\underset{k=1}{\overset{r}{\mathrm{\&Sigma;}}}(\frac{c_k^{\&prime;}+{\mathrm{jc}}_k^{\&prime;\&prime;}}{s-(p_k^{\&prime;}+{\mathrm{jp}}_k^{\&prime;\&prime;})}+\frac{c_k^{\&prime;}-{\mathrm{jc}}_k^{\&prime;\&prime;}}{s-(p_k^{\&prime;}-{\mathrm{jp}}_k^{\&prime;\&prime;})})---\left(17\right)$

式中，d、e由多项式除法确定；系数c_i、c′_k、c″_k分别由式(18)、(19)获得：

ci＝[(s-p_i)Y(s)]s＝p_i    (18)

$\left\{\begin{array}{c}{c}_k^{\&prime;}={\left[(s-p_k^{\&prime;}-{\mathrm{jp}}_k^{\&prime;\&prime;})Y\left(s\right)\right]}_{s=p_k^{\&prime;}+{\mathrm{jp}}_k^{\&prime;\&prime;}}\\ c_k^{\&prime;\&prime;}={\left[(s-p_k^{\&prime;}+{\mathrm{jp}}_k^{\&prime;\&prime;})Y\left(s\right)\right]}_{s=p_k^{\&prime;}-{\mathrm{jp}}_k^{\&prime;\&prime;}}\end{array}\right.---\left(19\right)$

从而，按照式(20)至(26)所示的对应关系，即可求出式(17)各项相对应的简单支路，再进行并联，即可获得相应的等值网络，如图6所示。

常数项和一次项对应R₀和C₀：

C₀＝e，R₀＝1/d    (20)

实根对应RL支路：

R₁＝-p_i/c_i，L₁＝1/c_i    (21)

共轭复根：

$\frac{c_k^{\&prime;}+{\mathrm{jc}}_k^{\&prime;\&prime;}}{s-(p_k^{\&prime;}+{\mathrm{jp}}_k^{\&prime;\&prime;})}+\frac{c_k^{\&prime;}-{\mathrm{jc}}_k^{\&prime;\&prime;}}{s-(p_k^{\&prime;}-{\mathrm{jp}}_k^{\&prime;\&prime;})}---\left(22\right)$

与以下RLC支路相对应：

L＝1/(2c′_k)(23)

R＝(-2p′_k+2(c_kp′_k+c″_kp″_k)L)L    (24)

1/C＝(p′_k ²+p″_k ²+2(c′_kp′_k+c″_kp″_k)R)L(25)

G＝-2(c′_kp′_k+c″_kp″_k)CL    (26)

②Y(s)含有多重实根极点

将多重实根极点看作单重实根极点，只须获得与单重实根极点相对应的等值电路。

105：通过外部系统的等值电路对外部电源进行处理，获取等效电源；

前面均是把外部系统当作无源网络进行处理的，实际应用时，还需考虑外部电源对内部系统的影响。由于外部电源只产生工频下的激励，对内部系统的影响也只有在系统频率为工频时有影响。故所求的等效电源只需保证工频时各端口处的电压电流在等值前后保持不变即可。参见图7和图8，可按如下步骤求取外部系统的等效电源：

1)计算各端口与地之间以及各端口之间的等值电路的阻抗，分别记为Z_ii(端口i与地之间)和Z_ij(端口i与j之间，j≠i)；

2)等值前外部系统工频下各端口处电压、电流的幅值及相角为：

根据戴维南等值原理分别获取单端口外部系统的等效电源和多端口外部系统的等效电源。

单端口外部系统的等效电源 ${\stackrel{\&CenterDot;}{E}}_i={\stackrel{\&CenterDot;}{U}}_i-{\stackrel{\&CenterDot;}{I}}_i{Z}_{\mathrm{ii}}---\left(27\right)$

多端口外部系统的等效电源 ${\stackrel{\&CenterDot;}{E}}_i={\stackrel{\&CenterDot;}{U}}_i-{\stackrel{\&CenterDot;}{I}}_i{Z}_{\mathrm{ii}}+\underset{j\&NotEqual;i}{\overset{m}{\underset{j=1}{\mathrm{\&Sigma;}}}}\frac{{\stackrel{\&CenterDot;}{U}}_i-{\stackrel{\&CenterDot;}{U}}_j}{Z_{\mathrm{ij}}}{Z}_{\mathrm{ii}}(i=\mathrm{1,2}...,m)---\left(28\right)$

至此，经过上述5个步骤，对于一个实际外部大系统，即可求得其等值系统。若外部系统的总节点数为n，采用分步求解后，每次待等值网络节点数为np(n_p≤15)，相应的求解其端口导纳有理函数式所需计算量为A，则整个外部系统的端口导纳有理函数式计算复杂度近似为

其中

表示上取整运算。按步骤103所示的求解单个子网络的端口导纳有理式方法，在子网络的节点数比较小时(如小于15)相应的端口导纳有理式计算时间就很小，因此A是一个相对比较小的值，故外部系统端口导纳有理函数式分步求解的总体计算复杂性与外部系统节点数n成近似线性关系，相应地其计算量及所需的时间均大幅降低，且均可被工程所接受。

106：获得实际系统的等值系统，对等值系统进行电磁暂态仿真分析，通过仿真结果检测设备及元件的谐波和电磁暂态特性，判断元件、设备及局部网络的谐波和暂态分量是否合理，如果否，更换电气设备、调整更合理参数及采取相应措施，如果是，流程结束。

其中，实际系统的等值系统由等值电路、等效电源和内部系统所组成的。基于本方法所得到的实际大电力系统的等值系统，在拥有很高等值精度的前提下，大大缩短了电磁暂态仿真时间，从而可以方便快捷地检测电力系统中元件、设备及局部网络的谐波和暂态分量，进而避免了不必要的经济损失。

下面以具体的实施例来说明本发明实施例提供的一种对大型电力系统进行电磁暂态仿真的方法的可行性，参见图9和图10，详见下文描述：

以IEEE-118电力节点系统(考虑发电机内阻)为例验证有效性。若节点1的发电机出口处在0.305s时发生A相接地故障，以其相邻线路1-2、1-3及节点1、2、3为内部系统，2号节点、3号节点为边界节点，其余为外部系统。

对于IEEE118系统，外部系统电源置零后，由于发电机存在内阻，其为115节点的网络。利用本方法和VF方法分别对外部系统进行等值，得到3号节点的电压波形及未等值时的电压波形曲线，参见图9(a)和图9(b)。可以看出本方法与VF方法均能与未等值时电压波形很好匹配，具有较高的精确度。但VF方法需要经过采样和有理式的求取两个阶段，耗时很长，VF方法总计耗时23.563s。利用本方法，依次从末端小系统向端口节点等值，即从子系统11向子系统1逐步计算，对于离研究系统较远的网络，如子系统2到子系统11，忽略其高频部分，以降低计算量，本方法总的计算时间为12.645s，相当于VF方法时间的一半。

当仿真时间长度为1s，简化后的系统及未简化系统的仿真步长均为50μs，未简化系统的仿真计算时间为26.859s，而对简化后系统的仿真计算时间缩短为7.235s，仿真时间大大减小，使得可以及时的发现电力故障并对其进行检修，进而避免了不必要的经济损失。

综上所述，本发明实施例提供了一种对大型电力系统进行电磁暂态仿真的方法，本方法采取的对外部系统的等值方法首先从外部系统的混合矩阵出发，设计有效方法把外部系统端口的导纳混合矩阵函数行列式直接转成有理式函数；把外部系统合理地划分为系列小系统，并逐个地得到各个小系统的有理式而得到外部系统的完整有理式；通过有理函数式获取外部系统的等值电路和等效电源；并对等值电路、等效电源和内部系统所组成的实际系统的等值系统进行电磁仿真；通过试验验证，本方法所采取的对外部系统的等值大幅度地缩短了外部系统等值网络的计算时间，所得到的等值系统拥有很高的等值精度；基于本发明实施例得到的实际大电力系统等值系统，大幅度地缩短了对实际大电力系统电磁暂态仿真时间，从而可以方便快捷地检测电力系统中元件、设备及局部网络的谐波和暂态分量，而大大提高了仿真效率，及时快速地发现元件、设备及网络中存在的问题，从而避免了不必要的经济损失。

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种对大型电力系统进行电磁暂态仿真的方法，其特征在于，所述方法包括以下步骤：

(1)获取频域下的外部无源网络；

(2)对所述外部无源网络进行划分获取若干个子外部无源网络；

(3)设共有M个子外部无源网络，按照第M个子外部无源网络到第一个子外部无源网络的顺序，依次获取各子导纳有理函数式，通过当前子外部无源网络与下一相邻子外部无源网络之间的端口，将当前子导纳有理函数式加入下一相邻子外部无源网络的导纳矩阵的自导和互导中，最终获得第一个子外部无源网络的子导纳有理函数式并作为外部系统的端口导纳有理函数式；

(4)通过所述外部系统的端口导纳有理函数式获取外部系统的等值电路；

(5)通过所述外部系统的等值电路对外部电源进行处理，获取等效电源；

(6)获得实际系统的等值系统，对所述等值系统进行电磁暂态仿真分析，通过仿真结果检测设备及元件的谐波和电磁暂态特性，判断元件、设备及局部网络的谐波和暂态分量是否合理，如果否，更换电气设备、调整更合理参数及采取相应措施；如果是，流程结束；

其中，步骤（3）具体过程为：

1）通过所述子外部无源网络获取混合矩阵；

2）通过所述混合矩阵获取子导纳有理函数式Y^r(s)；

3）通过所述子导纳有理函数式Y^r(s)依次求解最终获取第一个子外部无源网络的子导纳有理函数式Y¹(s)；

其中，r表示外部无源网络划分成的M个子网络中的任意一个，即2≤r≤M；其中，所述通过所述子外部无源网络获取混合矩阵具体为：

获得频域下的外部系统无源网络后，在子外部系统的某一连接端口处注入一个单位的电流，其余端口开路，得到式（1）：

Y^rU^r=I^r _input          （1）

其中，Y^r为频域下的第r个子外部系统的节点导纳矩阵；U^r为第r个子外部系统和端口的各节点电压；I^r _input为第r个子外部系统和端口的注入电流列向量，对应注入单位电流节点的元素为1，其余元素都为0；

若电压观测节点的电压用F^r表示，则有：

F^r=(d^r)^TU^r          （2）

其中d为节点关联矩阵，即与电压观测节点相对应的位置元素为1，其余元素为0；d^r表示第r个外部无源子网络的关联矩阵；

将式（1）与式（2）合并即可得如式（3）所示的混合矩阵：

$\left[\begin{array}{cc}{Y}^r& 0\\ -{\left(d^r\right)}^T& 1\end{array}\right]\left[\begin{array}{c}{U}^r\\ F^r\end{array}\right]=\left[\begin{array}{c}{I^r}_{\mathrm{input}}\\ 0\end{array}\right]---\left(3\right)$

利用克莱姆法则可求得F^r为：

$F^r=\frac{\left[\begin{array}{cc}{Y}^r& {I^r}_{\mathrm{input}}\\ -{\left(d^r\right)}^T& 0\end{array}\right]}{\det \left[\begin{array}{cc}{Y}^r& 0\\ -{\left(d^r\right)}^T& 1\end{array}\right]}=\frac{{\det H}_1}{\det H}---\left(4\right)$

因式（4）中Y^r的各非零元素均为s的多项式，detH₁及detH均为s的多项式，其中，det为矩阵的行列式，T为矩阵转置；

其中，所述通过所述混合矩阵获取子导纳有理函数式Y^r(s)具体为：

$Y^r\left(s\right)=\frac{I^r}{U^r}=\frac{1}{F^r}=\frac{\det H}{\det H_1}---\left(5\right)$

若注入电流节点和电压观测节点为同一节点，则Y^r(s)为自导纳；否则，Y^r(s)为两节点间的互导纳；

矩阵H和H₁由

求得，

$B_H^1={\left[\begin{array}{ccccc}{Y}_{\mathrm{1,1}}& ...& Y_{1,j}& ...& Y_{1,n}\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ Y_{i,1}& ...& Y_{i,j}& ...& Y_{i,n}\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ Y_{n,j}& ...& Y_{n,j}& ...& Y_{n,n}\end{array}\right]}_{(n\&times;n)}---\left(6\right)$

①对于式（7）所示的多项式：

a₀+a₁s+a₂s²+…+a_nsⁿ          （7）

记为式（8）的向量形式称为向量多项式：

[a₀,a₁,a₂,…,a_n]          （8）

②定义⊕、

为向量多项式运算中的加法与乘法，其运算规则分别为：

[b₀,…,b_n]⊕[c₀,…,c_m]=[b₀+c₀,…,b_n+c_n,c_n+1…,c_m]，n≤m          （9）

$[b_0,...,b_n]\&CircleTimes;[c_0,...,c_m]=[b_0{c}_0,(b_1{c}_0+b_0{c}_1),...,(b_{n-2}{c}_m+b_{n-1}{c}_{m-1}+b_n{c}_{m-2}),(b_{n-1}{c}_m+b_n{c}_{m-1}),b_n{c}_m]---\left(10\right)$

式（6）中的节点自导纳Y_ii与互导纳Y_ij如式（11）所示：

$\left\{\begin{array}{c}{Y}_{\mathrm{ii}}=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\frac{1}{R_{\mathrm{ij}}+s{L}_{\mathrm{ij}}}\right)+(G_{\mathrm{ii}}+s{C}_{\mathrm{ii}})\\ Y_{\mathrm{ij}}=-\underset{j\&NotEqual;i}{\overset{n}{\underset{j=1}{\mathrm{\&Sigma;}}}}\frac{1}{R_{\mathrm{ij}}+s{L}_{\mathrm{ij}}}\end{array}\right.---\left(11\right)$

式中Ri_j、Li_j、C_ii、G_ii分别表示支路ij的电阻、电感及节点i的对地总导纳与总电容；n为系统节点个数，把Y_ii及Y_ij采用向量多项式进行表示

$\left\{\begin{array}{c}{Y}_{\mathrm{ii}}=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{1}{[R_{\mathrm{ij}},L_{\mathrm{ij}}]}+[G_{\mathrm{ii}},C_{\mathrm{ii}}]\\ Y_{\mathrm{ij}}=-\underset{j\&NotEqual;i}{\overset{n}{\underset{j=1}{\mathrm{\&Sigma;}}}}\frac{1}{{[R}_{\mathrm{ij}},L_{\mathrm{ij}}]}\end{array}\right.---\left(12\right)$

将

中的各个元素表示成式（12）的向量形式；然后对各个非零元素按照定义的运算规则进行下一步操作，以使得每个非零元素的分子、分母均为以s为变量的向量多项式；

将矩阵

简记为B¹，其行列式简记为det(B¹)；

变量说明：k为展开层数；矩阵

为N_k×N_k的方阵，N_k表示展开第k层后矩阵

的阶数；m_k为矩阵列号1≤m_k≤N_k，

为矩阵中第i行第m_k列元素，

为矩阵

中第1行第m_k列元素

对应的伴随矩阵；S_k、h_k为计算矩阵

行列式的中间变量；

变量初始化：k=1，m₀=1，

N₁=n，

1≤i,j≤n；

①判断是否N_k=0，如果是，

执行步骤⑥；如果否，执行步骤②；

②判断是否N_k=1,如果是，执行步骤⑥；如果否，执行步骤③；

③判断是否N_k=2，如果是，

执行步骤⑥；如果否，执行步骤④；

④N_k>2，按行展开求取行列式，即令m_k=1，定义S_k=0，h_k=1，并进行如下操作：

1）判断是否满足m_k≤N_k且

非零，如果是，执行步骤⑤；如果否，执行步骤3）；

2）k=k-1，计算

h_k=-h_k，m_k=m_k+1，转步骤1）；

3）判断是否满足m_k≤N_k且为零，如果是，h_k=-h_k，m_k=m_k+1，转步骤1）；如果否，转步骤4）；

4）判断是否k=1，如果是，则det(B¹)=S_k，执行步骤⑥；如果否，

k=k-1，

m_k=m_k+1，转步骤1）；

⑤求取的伴随矩阵，k=k+1；去掉

的第m_k-1列和第一行即得到

N_k=N_k-1-1；转步骤①；

⑥判断是否k=1，如果是，计算矩阵H行列式，detH=detB¹，流程结束；如果否，则转步骤④中的2）；

将detH、detH₁代人式（5）可得Y^r(s)：

$Y^r\left(s\right)=\det H\&CircleTimes;(1/\det H_1)---\left(13\right)$

形式为：

$Y^r\left(s\right)=\frac{N\left(s\right)}{D\left(s\right)}=\frac{a_0+a_{1}s+a_2{s}^2+...+a_n{s}^m}{1+b_1{s}^1+b_2{s}^2+...+b_n{s}^n}---\left(14\right)$

；

其中，所述通过所述子导纳有理函数式Y^r(s)依次求解最终获取第一个子外部无源网络的子导纳有理函数式Y¹(s)具体为：

(Ⅰ)若已求得第r,2≤r≤M个子网络端口q₁、q₂外部等效的自导纳有理函数式

及互导纳有理函数式

则与之相邻的第r-1个子网络的导纳矩阵为：

$\left[\begin{array}{ccccc}.& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ ...& Y_{q_1}^{\&prime;}\left(s\right)+Y_{q_1}^r\left(s\right)& ...& Y_{q_1,q_2}^{\&prime;}\left(s\right)-Y_{q_1,q_2}^r\left(s\right)& ...\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ ...& Y_{q_1,q_2}^{\&prime;}\left(s\right)-Y_{q_1,q_2}^r\left(s\right)& ...& Y_{q_2}^{\&prime;}\left(s\right)+Y_{q_2}^r\left(s\right)& ...\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\end{array}\right]$ 式中

分别表示第r-1个子网络中q₁、q₂节点的自导及互导；

(Ⅱ)第1子系统保留有理函数式的全部极点，从第2个子系统到第M个子系统，依次降低有理式中的高频极点的个数，并增加常数项的值；

按照由第M个子外部无源网络到第一个子外部无源网络的顺序依次求取，直到求得第一个子系统的等值自导及互导的有理式Y¹(s)，即获取到了整个外部系统的端口导纳有理式Y(s)，即Y(s)=Y¹(s)。

2.根据权利要求1所述的一种对大型电力系统进行电磁暂态仿真的方法，其特征在于，所述对所述外部无源网络进行划分获取若干个子外部无源网络具体为：

根据所述外部无源网络的结构，把所述外部无源网络划分成所述若干个子外部无源网络，按内部系统端口到外部系统末端的顺序依次编号。

Images