CN102270845B

CN102270845B - 一种电力系统暂态稳定数字仿真计算模型及其算法

Info

Publication number: CN102270845B
Application number: CN201110218015.1A
Authority: CN
Inventors: 朱建国; 危雪; 刘淑军
Original assignee: China Three Gorges University CTGU
Current assignee: China Three Gorges University CTGU
Priority date: 2011-08-01
Filing date: 2011-08-01
Publication date: 2014-03-05
Anticipated expiration: 2031-08-01
Also published as: CN102270845A

Abstract

本发明提供了一种电力系统暂态稳定数字仿真计算模型及其算法，属于电力系统数字仿真领域。传统的机网接口算法将各发电机dq轴转换计算放在电网节点电压代数方程组的求解中进行，使得电网仿真模型成为时变的非线性代数方程组。本发明将各发电机dq轴转换计算放在各发电机组仿真模型中求解，使得电网仿真模型为时不变的线性代数方程组，大大降低了交替求解中电网仿真模型的计算量。该发明不仅使得采用交替求解的电力系统暂态稳定数字仿真计算得到简化，模型的物理意义更加明确，而且也显著降低了编程难度；由于交替求解中各发电机组仿真模型间的计算是解耦的，采用该发明可以构成一种实用的电力系统暂态稳定数字仿真并行计算模型。

Description

一种电力系统暂态稳定数字仿真计算模型及其算法

技术领域

本发明涉及一种电力系统暂态稳定数字仿真计算模型及其算法，属于电力系统数字仿真领域。

背景技术

在暂态稳定分析中，电力系统全系统数学模型的一般形式可写为：

\{\begin{matrix} \overset{\cdot}{x} = f (x, U) \\ I (x, U) = YU \end{matrix} - - - (1)

其中，x=系统的状态变量；U=节点电压；I=节点注入电流；式（1）的第1个方程为描述发电机组或动态负荷的微分方程组；第2个方程为描述电网的代数方程组。

电力系统暂态稳定仿真计算的核心即：联立求解或交替求解方程1和2。

目前对机网接口的计算主要有以下4种方法：

1、发电机采用经典模型时的处理方法。可化为发电机等值导纳Y_G和发电机等值电流源I_G=Y_GE相并联的形式，将Y_G并入导纳阵，无操作时导纳阵不变，而每时步根据发电机转子角更新发电机注入网络的等值电流源即可求解网络方程，计算节点电压。

2、考虑凸积效应的直接解法。该方法将网络复数线性代数方程实、虚部分开，增阶化为xy同步坐标下实数线性代数方程，并将发电机方程由dq坐标化为xy坐标，，再和网络方程联立求解。这种解法对负荷非线性适应能力差，且发电机方程由dq坐标据转子角转化为xy坐标，引起导纳阵中发电机节点相对应的对角子块由于转子角的变化而为非定常元素，每一时步要重新计算因子表，耗费机时多且内存要增加一倍。

3、考虑凸积效应的迭代解法。该方法的特点是力求在复数域中求解线性代数方程来实现网络方程的求解，并使导纳阵元素在无操作时保持定常，不随发电机转子角而变化。但发电机凸积效应及转子角变化要通过修正发电机注入网络的电流源来计及，而电流源的计算还和t_n+1时刻的节点电压值有关。由于t_n+1时刻的节点电压正待计算，因此t_n+1时刻的电流源的值需通过迭代法计算而逼近准确值。该方法相对于直接解法有节省内存、因子表定常、计算速度快、便于适应非线性负荷模型等特点。该方法在一些实用暂态稳定分析程序中仍有应用。

4、考虑凸极效应的牛顿法。牛顿法是求解非线性代数方程组的优良方法，有良好的收敛性能。当发电机计及凸极效应，负荷计及非线性，系统中元件微分方程化为差分代数方程后，全网的代数方程联立，实质上是要求解一组非线性代数方程，故也可采用牛顿法求解。相对于直接解法和迭代解法，用牛顿法进行机网接口计算编程复杂，因为要计算雅可比矩阵元素，而雅可比矩阵元素随时间而变化，故计算机时也较多。但其最大优点是对非线性元件模型的适应性好，可将微分方程的差分代数方程和系统代数方程联立求解，无“交接误差”，故计算精度高、累计误差小，因而在暂态稳定分析中广泛应用。它常和隐式梯形积分法求解微分方程相结合。

电力系统暂态稳定数字仿真计算的核心是联立求解或交替求解（1）式中的方程1和方程2。在交替求解的传统方法中，如机网接口计算方法2、3所述，各发电机dq轴转换计算在电网代数方程组的求解中进行，使得电网仿真模型成为时变的非线性代数方程组，因而电网仿真模型计算量大且编程实现较复杂。而机网接口计算方法4所述牛顿法是联立求解（1）式中的方程1和方程2，虽然具有计算精度高等优点，但每时步的计算量也很大，上述方法均不易构成大规模电力系统暂态稳定实时数字仿真的并行计算模型。

发明内容

本发明所要解决的技术问题是提供一种电力系统暂态稳定数字仿真计算模型及其算法，将各发电机dq轴转换计算放在发电机组仿真模型中求解，使得电网仿真模型为时不变的线性代数方程组，大大降低了交替求解中电网仿真模型的计算量。该发明不仅使得采用交替求解的电力系统暂态稳定仿真计算得到简化，模型的物理意义更加明确，而且也显著降低了模型的编程难度和维护工作量，其更加重要的意义在于：采用该发明可以构成一种电力系统暂态稳定仿真并行计算模型。电网仿真模型子任务由PC0机计算，n个节点上的发电机组仿真模型子任务计算则平均分配给PC1～PCm机进行并行计算。电网仿真模型计算量的显著降低使得图3所示方案的实现成为可能。该方案的突出优点是：即便对于大规模电力系统的实时数字仿真，节点上的发电机组或动态负荷都可以使用详细模型。

本发明采取的技术方案是：一种电力系统暂态稳定数字仿真计算模型，所述电力系统暂态稳定数字仿真数学模型分割为电网数学模型和发电机组数学模型两部分（负荷动态模型的异步电机看作同步电机的特例），电网数学模型包括n节点的电网复数线性代数方程YU=I和发电机定子绕组电压方程，其中各发电机定子绕组电压方程中的变量均为同一dq轴上或xy同步坐标轴上描述的变量值，根据此数学模型形成电网仿真模型；发电机组数学模型包括dq轴转换方程和发电机组方程（不包括发电机定子绕组电压方程），根据此数学模型形成发电机组仿真模型。

所述电网仿真模型，是将n节点的电网复数线性代数方程YU=I，增阶化为同一dq轴上或xy同步坐标轴上的2n维实数线性代数方程Y_xyU_xy=I_xy，将权利要求1所述各发电机定子绕组电压方程和增阶网络方程Y_xyU_xy=I_xy一起联立，整理后可以形成如下方程：I_xy=AB_xyE′_xy，即为电网仿真模型，为时不变的实数线性代数方程，AB_xy为2n维线性时不变矩阵，I_xy为2n维节点注入电流向量，E′_xy为发电机定子绕组电压方程中的电势，为2n维向量。

所述发电机组仿真模型，是采用数值积分方法将发电机组的数学模型差分化，并将发电机组的差分化方程组中所包括的发电机定子绕组电压方程中的电势变量和机端电流变量，用电势E′_xy和节点注入电流I_xy来表示，即：将E′_xy和I_xy由同一dq轴或xy轴上根据dq轴转换方程转化至本机的dq轴上，形成各发电机组的仿真模型。

所述发电机组数学模型分割为同步电机数学模型，励磁系统数学模型，原动机及其调速系统数学模型三部分，其中同步电机数学模型包括转子绕组方程和dq轴转换方程，励磁系统数学模型包括励磁系统方程，原动机及其调速系统数学模型包括原动机及其调速系统方程和转子运动方程，采用数值积分方法将同步电机的转子绕组方程差分化，得到其差分化方程组，并根据权利要求3所述方法进行dq轴转换计算，形成同步电机的仿真模型；同样，采用数值积分方法将励磁系统、原动机及其调速系统的数学模型差分化，分别形成励磁系统仿真模型和原动机及其调速系统仿真模型，同步电机的仿真模型、励磁系统仿真模型和原动机及其调速系统仿真模型一起构成发电机组仿真模型。

一种电力系统暂态稳定数字仿真计算算法，由以下步骤构成：

步骤一：根据上一次计算得到的E′_xy的值，采用电网仿真模型计算发电机注入网络的电流I_xy的值；

步骤二：根据同步电机仿真模型本时步的状态量的值，及I_xy的值，及上一次励磁系统仿真模型计算得到的励磁电压E_f和原动机及其调速系统仿真模型计算得到的转子角δ的值，采用同步电机仿真模型计算下一时步的状态量的值E′_xy；

步骤三：根据E′_xy及I_xy的值可求取机端电压U_t的值，并根据励磁系统本时步状态量的值，采用励磁系统仿真模型计算励磁系统下一时步状态量的值；

步骤四：已知E′_xy及I_xy的值求取机端电磁力矩T_e的值，并根据原动机及其调速系统本时步状态量的值，采用原动机及其调速系统仿真模型计算其下一时步状态量的值；

步骤五：根据本次计算得到的发电机组状态量的值与上次计算得到的发电机组状态量的值之差的绝对值的大小判断计算是否收敛；

步骤六：若不收敛，则以本次计算得到的E′_xy的值，重复步骤一到步骤六的计算和判断；

步骤七：若收敛，则令t=t+Δt，相应各状态量t时刻的值替换为t+Δt时刻的值，重复步骤一到步骤七进行下一时步的计算。

一种电力系统暂态稳定数字仿真计算模型及其算法的有益效果如下：

发电机定子绕组与电网在电路上是直接相连的，电网的数学模型在包括发电机定子绕组电压方程的情况下不包括dq轴转换计算意味着发电机定子绕组的电压方程是在电网的数学模型所在的同步坐标轴上列写的。因此，电网仿真模型的计算量大大降低。而定子绕组电压方程中的电势变量和电流变量也出现在发电机转子绕组电压方程中，因此dq轴转换计算放在各发电机组的差分化方程组中进行。这不仅使得采用交替求解的电力系统暂态稳定仿真计算得到简化，模型的物理意义更加明确，也大大降低了编程难度。其更加重要的意义在于：由于交替求解中各发电机组模型间的计算是解耦的，采用该发明可以构成一种实用的电力系统暂态稳定仿真并行计算模型。如图3所示，电网仿真模型子任务由PC0机计算，n个节点上的发电机组仿真模型子任务计算则平均分配给PC1～PCm机进行并行计算。电网仿真模型计算量的显著降低使得图3所示方案的实现成为可能。该方案的突出优点是：即便对于大规模电力系统的实时数字仿真，节点上的发电机组或动态负荷都可以使用详细模型。

发电机组的计算机仿真模型按照其物理子系统建立，分别为同步电机仿真模型，励磁系统仿真模型，原动机及其调速系统仿真模型，可以分别建立这些物理子系统的仿真模型库，从而可以采用“搭积木”的方式构成不同的发电机组仿真模型，模型的编制和维护工作量都较大程度降低。

附图说明

下面结合附图和实施例对本发明作进一步说明：

图1是本发明的暂态稳定数字仿真全系统数学模型构成联系示意图；

图2是本发明的暂态稳定数字仿真计算模型算法的流程图；

图3是采用本发明的计算模型及其算法可能实现的电力系统暂态稳定实时数字仿真PC机群方案图。

具体实施方式

通过电网仿真模型的形成、发电机组仿真模型的形成、以及仿真模型的交替求解计算三部分来说明本发明的具体求解过程。

一、电网仿真模型的形成

电网数学模型包括发电机定子绕组电压方程和n节点的电网复数线性代数方程YU=I，其中各发电机定子绕组电压方程中的变量均为同一dq轴上或xy同步坐标轴上描述的变量值，根据此数学模型形成电网仿真模型。

当发电机计及暂态凸极效应，发电机定子绕组电压方程不能用简单的复数关系来表示，必须对d轴、q轴等值绕组分别列方程。当发电机采用四阶（或三阶）实用模型时，定子绕组电压方程为（三阶模型时，下列方程中E′_d=0，X′_q为X_q）

\{\begin{matrix} U_{d} = E_{d}^{'} + X_{q}^{'} I_{q} - r_{a} I_{d} \\ U_{q} = E_{q}^{'} - X_{d}^{'} I_{d} - r_{a} I_{q} \end{matrix} - - - (2)

当发电机采用五阶、六阶模型时

\{\begin{matrix} U_{d} = E_{d}^{''} + X_{q}^{''} I_{q} - r_{a} I_{d} \\ U_{q} = E_{q}^{''} - X_{d}^{''} I_{d} - r_{a} I_{q} \end{matrix} - - - (3)

由于式（2）和式（3）有相同的形式，故以下以发电机四阶实用模型为例，来说明电网仿真模型的形成。

式（2）写成矩阵形式：

[\begin{matrix} U_{d} \\ U_{q} \end{matrix}] = [\begin{matrix} E_{d}^{'} \\ E_{q}^{'} \end{matrix}] - [\begin{matrix} r_{a} & - X_{q}^{'} \\ X_{d}^{'} & r_{a} \end{matrix}] [\begin{matrix} I_{d} \\ I_{q} \end{matrix}] - - - (4)

这里设已将（4）式由dq坐标化为xy同步坐标（注：xy同步坐标也可为同一dq轴坐标，如第n台发电机的dq轴坐标dqn，以后不再注明），即发电机定子绕组电压方程中的变量为xy同步坐标轴上描述的变量值：

[\begin{matrix} U_{x} \\ U_{y} \end{matrix}] = [\begin{matrix} E_{x}^{'} \\ E_{y}^{'} \end{matrix}] - [\begin{matrix} r_{a} & - X_{q}^{'} \\ X_{d}^{'} & r_{a} \end{matrix}] [\begin{matrix} I_{x} \\ I_{y} \end{matrix}] - - - (5)

对于电网节点电压方程YU=I，Y为导纳矩阵，U为节点电压向量，I为节点注入电流向量。为了便于机网接口，将n节点的电网节点电压方程YU=I增阶化为2n维的实数线性代数方程：

[\begin{matrix} [\begin{matrix} I_{x 1} \\ I_{y 1} \end{matrix}] \\ \cdot \\ \cdot \\ \cdot \\ [\begin{matrix} I_{xi} \\ I_{yi} \end{matrix}] \\ \cdot \\ \cdot \\ \cdot \\ [\begin{matrix} I_{xn} \\ I_{yn} \end{matrix}] \end{matrix}] = [\begin{matrix} [\begin{matrix} G_{11} & {- B}_{11} \\ B_{11} & G_{11} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{1 i} & - B_{1 i} \\ B_{1 i} & G_{1 i} \end{matrix}] & \cdot \cdot & [\begin{matrix} G_{1 n} & {- B}_{1 n} \\ B_{1 n} & G_{1 n} \end{matrix}] \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ [\begin{matrix} G_{i 1} & {- B}_{i 1} \\ B_{i 1} & G_{i 1} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{ii} & {- B}_{ii} \\ B_{ii} & G_{ii} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{1 i} & {- B}_{1 i} \\ B_{1 i} & G_{1 i} \end{matrix}] \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ [\begin{matrix} G_{n 1} & {- B}_{n 1} \\ B_{n 1} & G_{n 1} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{ni} & {- B}_{ni} \\ B_{ni} & G_{ni} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{nn} & {- B}_{nn} \\ B_{nn} & G_{nn} \end{matrix}] \end{matrix}] [\begin{matrix} [\begin{matrix} U_{x 1} \\ U_{y 1} \end{matrix}] \\ \cdot \\ \cdot \\ \cdot \\ [\begin{matrix} U_{xi} \\ U_{yi} \end{matrix}] \\ \cdot \\ \cdot \\ \cdot \\ [\begin{matrix} U_{xn} \\ U_{yn} \end{matrix}] \end{matrix}] - - - (6)

式中，Gi_j+jB_ij=Y_ij为Y阵中i行j列元素。其中I_xi+jI_yi=I和U_xi+jU_yi=U分别为I、U中第i个元素

为不失一般性，设式（5）所描写的发电机接于网络第i个节点，则（5）式中U_x和U_y即为（6）式中U_xi和U_yi，将式（5）代入式（6）中第i个节点方程，消去U_xi和U_yi；对各个发电机节点均作相同处理后，可得如下网络方程。

[\begin{matrix} [\begin{matrix} I_{x 1} \\ I_{y 1} \end{matrix}] \\ \cdot \\ \cdot \\ \cdot \\ [\begin{matrix} I_{xi} \\ I_{yi} \end{matrix}] \\ \cdot \\ \cdot \\ \cdot \\ [\begin{matrix} I_{xn} \\ I_{yn} \end{matrix}] \end{matrix}] = [\begin{matrix} [\begin{matrix} G_{11} & {- B}_{11} \\ B_{11} & G_{11} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{li} & {- B}_{li} \\ B_{1 i} & G_{1 i} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{1 n} & {- B}_{1 n} \\ B_{1 n} & G_{1 n} \end{matrix}] \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ [\begin{matrix} G_{i 1} & {- B}_{i 1} \\ B_{i 1} & G_{i 1} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{ii} & {- B}_{ii} \\ B_{ii} & G_{ii} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{1 i} & {- B}_{1 i} \\ B_{1 i} & G_{1 i} \end{matrix}] \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ [\begin{matrix} G_{n 1} & {- B}_{n 1} \\ B_{n 1} & G_{n 1} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{ni} & {- B}_{ni} \\ B_{ni} & G_{ni} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{nn} & {- B}_{nn} \\ B_{nn} & G_{nn} \end{matrix}] \end{matrix}] [\begin{matrix} [\begin{matrix} E_{x 1}^{'} \\ E_{y 1}^{'} \end{matrix}] - [\begin{matrix} r_{a 1} & {- X}_{q 1}^{'} \\ X_{d 1}^{'} & r_{a 1} \end{matrix}] [\begin{matrix} I_{x 1} \\ I_{y 1} \end{matrix}] \\ \cdot \\ \cdot \\ \cdot \\ [\begin{matrix} E_{xi}^{'} \\ E_{yi}^{'} \end{matrix}] - [\begin{matrix} r_{a 1} & {- X}_{qi}^{'} \\ X_{d 1}^{'} & r_{a 1} \end{matrix}] [\begin{matrix} I_{xi} \\ I_{yi} \end{matrix}] \\ \cdot \\ \cdot \\ \cdot \\ [\begin{matrix} E_{xn}^{'} \\ E_{yn}^{'} \end{matrix}] - [\begin{matrix} r_{an} & {- X}_{qn}^{'} \\ X_{dn}^{'} & r_{an} \end{matrix}] [\begin{matrix} I_{xn} \\ I_{yn} \end{matrix}] \end{matrix}] - - - (7)

上式可化为：

其中：

\begin{matrix} [\begin{matrix} G_{Fji} & {- B}_{Fji} \\ B_{Fji} & G_{Fji} \end{matrix}] = [\begin{matrix} G_{ji} & {- B}_{ji} \\ B_{ji} & G_{ji} \end{matrix}] [\begin{matrix} r_{ai} & {- X}_{qi}^{'} \\ X_{di}^{'} & r_{ai} \end{matrix}], (i = 1 ~ n, j = 1 ~ n, j &NotEqual; i) \\ [\begin{matrix} G_{Fii} & {- B}_{Fii} \\ B_{Fii} & G_{Fii} \end{matrix}] [\begin{matrix} G_{ii} & {- B}_{ii} \\ B_{ii} & G_{ii} \end{matrix}] [\begin{matrix} r_{ai} & {- X}_{qi}^{'} \\ X_{di}^{'} & r_{ai} \end{matrix}] + [\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}], (j = i) \end{matrix}

（8）式可简写为：

A_xyI_xy=B_xyE_xy （9）

其中，

A_{xy} = [\begin{matrix} [\begin{matrix} G_{F 11} & {- B}_{F 11} \\ B_{F 11} & G_{F 11} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{F 1 i} & {- B}_{F 1 i} \\ B_{F 1 i} & G_{F 1 i} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{F 1 n} & {- B}_{F 1 n} \\ B_{F 1 n} & G_{F 1 n} \end{matrix}] \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ [\begin{matrix} G_{Fi 1} & {- B}_{Fi 1} \\ B_{Fi 1} & G_{Fi 1} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{Fii} & {- B}_{Fii} \\ B_{Fii} & G_{Fii} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{F 1 i} & {- B}_{F 1 i} \\ B_{F 1 i} & G_{F 1 i} \end{matrix}] \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ [\begin{matrix} G_{Fn 1} & {- B}_{Fn 1} \\ B_{Fn 1} & G_{Fn 1} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{Fni} & {- B}_{Fni} \\ B_{Fni} & G_{Fni} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{Fnn} & {- B}_{Fnn} \\ B_{Fnn} & G_{Fnn} \end{matrix}] \end{matrix}] I_{xy} = [\begin{matrix} [\begin{matrix} I_{x 1} \\ I_{y 1} \end{matrix}] \\ \cdot \\ \cdot \\ \cdot \\ [\begin{matrix} I_{xi} \\ I_{yi} \end{matrix}] \\ \cdot \\ \cdot \\ \cdot \\ [\begin{matrix} I_{xn} \\ I_{yn} \end{matrix}] \end{matrix}]

B_{xy} = [\begin{matrix} [\begin{matrix} G_{11} & {- B}_{11} \\ B_{11} & G_{11} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{1 i} & {- B}_{1 i} \\ B_{1 i} & G_{1 i} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{1 n} & {- B}_{1 n} \\ B_{1 n} & G_{1 n} \end{matrix}] \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ · & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ [\begin{matrix} G_{i 1} & {- B}_{i 1} \\ B_{i 1} & G_{i 1} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{ii} & {- B}_{ii} \\ B_{ii} & G_{ii} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{1 i} & {- B}_{1 i} \\ B_{1 i} & G_{1 i} \end{matrix}] \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & · & \cdot \\ [\begin{matrix} G_{n 1} & {- B}_{n 1} \\ B_{n 1} & G_{n 1} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{ni} & {- B}_{ni} \\ B_{ni} & G_{ni} \end{matrix}] & \cdot \cdot \cdot & [\begin{matrix} G_{nn} & {- B}_{nn} \\ B_{nn} & G_{nn} \end{matrix}] \end{matrix}] E_{xy} = [\begin{matrix} [\begin{matrix} E_{x 1}^{'} \\ E_{y 1}^{'} \end{matrix}] \\ \cdot \\ \cdot \\ \cdot \\ [\begin{matrix} E_{xi}^{'} \\ E_{yi}^{'} \end{matrix}] \\ \cdot \\ \cdot \\ \cdot \\ [\begin{matrix} E_{xn}^{'} \\ E_{yn}^{'} \end{matrix}] \end{matrix}]

Ixi为Ixy中第2i-1个元素，Iyi为Ixy中第2i个元素；E′xi为Exy中第2i-1个元素，E′yi为Exy中第2i个元素。

显然，由（9）式可得：

I_xy=AB_xyE_xy （10）

其中，AB_xy=A_xy ^-1B_xy，为2n维线性时不变矩阵。

二、发电机组仿真模型的形成

将发电机组数学模型分割为同步电机数学模型，励磁系统数学模型，原动机及其调速系统数学模型三部分，采用数值积分方法将同步电机的转子绕组方程差分化，得到其差分化方程组，并对定子绕组电压方程中也包含的电势和电流变量进行dq轴转换计算，形成同步电机的仿真模型；同样，采用数值积分方法将励磁系统、原动机及其调速系统的数学模型差分化，分别形成励磁系统仿真模型和原动机及其调速系统仿真模型，同步电机的仿真模型、励磁系统仿真模型和原动机及其调速系统仿真模型一起构成发电机组仿真模型。

当系统较大，且考虑各种调节器动态时，系统最大时间常数和最小时间常数之比可能很大，而呈现很强的“刚性”，隐式梯形法具有A稳定性，本文以微分方程采用隐式梯形法求解，发电机采用4阶实用模型、励磁系统如（15）式所示、原动机及其调速系统如（20）式所示为例，来说明发电机组模型的计算。

为不失一般性，令（11）至（24）式表示第i台发电机组模型的计算。其他发电机组模型的计算类似。

发电机转子绕组数学模型：

\{\begin{matrix} T_{d 0}^{'} {pE}_{q}^{'} = E_{f} - E_{q}^{'} - (X_{d} - X_{d}^{'}) i_{d} \\ T_{q 0}^{'} {pE}_{d}^{'} = - E_{d}^{'} + (X_{q} - X_{q}^{'}) i_{q} \end{matrix} - - - (11)

式中，Ed′、Eq′分别为d、q轴暂态电势，Xd、Xq分别为d、q轴同步电抗，Xd′、Xq′分别为d轴、q轴暂态电抗，Td0′、Tq0′分别为d、q轴开路暂态时间常数，Ef为励磁电压。

采用隐式梯形积分法差分化：

\begin{matrix} [\begin{matrix} T_{d 0}^{'} + Dt & 0 \\ 0 & T_{q 0}^{'} + Dt \end{matrix}] [\begin{matrix} E_{q}^{'} (t + Δt) \\ E_{d}^{'} (t + Δt) \end{matrix}] \\ = [\begin{matrix} T_{d 0}^{'} + Dt & 0 \\ 0 & T_{q 0}^{'} + Dt \end{matrix}] [\begin{matrix} E_{q}^{'} (t) \\ E_{d}^{'} (t) \end{matrix}] + [\begin{matrix} 2 Dt (X_{d}^{'} - X_{q}) & 0 \\ 0 & 2 Dt (X_{q} - X_{q}^{'}) \end{matrix}] [\begin{matrix} i_{d} \\ i_{q} \end{matrix}] + [\begin{matrix} 2 Dt \\ 0 \end{matrix}] E_{f} \end{matrix} - - - (12)

式中Dt=Δt，上式化为：

\begin{matrix} [\begin{matrix} E_{d}^{'} (t + Δt) \\ E_{q}^{'} (t + Δt) \end{matrix}] \\ = {[\begin{matrix} {2 T}_{q 0}^{'} + Dt & 0 \\ 0 & {2 T}_{d 0}^{'} + Dt \end{matrix}]}^{- 1} [\begin{matrix} {2 T}_{q 0}^{'} & 0 \\ 0 & {2 T}_{d 0}^{'} - Dt \end{matrix}] [\begin{matrix} E_{d}^{'} (t) \\ E_{q}^{'} (t) \end{matrix}] \\ + {[\begin{matrix} {2 T}_{q}^{'} + Dt & 0 \\ 0 & {2 T}_{d 0}^{'} + Dt \end{matrix}]}^{- 1} [\begin{matrix} 0 & 2 Dt (X_{q} - X_{q}^{'}) \\ 2 Dt (X_{d}^{'} - X_{d}) & 0 \end{matrix}] [\begin{matrix} i_{d} \\ i_{q} \end{matrix}] \\ + {[\begin{matrix} {2 T}_{q 0}^{'} + Dt & 0 \\ 0 & {2 T}_{d 0}^{'} + Dt \end{matrix}]}^{- 1} [\begin{matrix} 0 \\ 2 Dt \end{matrix}] E_{f} \end{matrix} - - - (13)

为了和网络方程接口，需将本机的dq轴坐标化为xy轴同步坐标，对（13）式两侧右乘坐标变换阵T，则式（13）化为：

\begin{matrix} [\begin{matrix} E_{x}^{'} (t + Δt) \\ E_{y}^{'} (t + Δt) \end{matrix}] \\ = T {[\begin{matrix} {2 T}_{q 0}^{'} + Dt & 0 \\ 0 & {2 T}_{d 0}^{'} + Dt \end{matrix}]}^{- 1} [\begin{matrix} {2 T}_{q 0}^{'} & 0 \\ 0 & {2 T}_{d 0}^{'} - Dt \end{matrix}] T^{- 1} [\begin{matrix} E_{x}^{'} (t) \\ E_{y}^{'} (t) \end{matrix}] \\ + T {[\begin{matrix} {2 T}_{q}^{'} + Dt & 0 \\ 0 & {2 T}_{d 0}^{'} + Dt \end{matrix}]}^{- 1} [\begin{matrix} 0 & 2 Dt (X_{q} - X_{q}^{'}) \\ 2 Dt (X_{d}^{'} - X_{d}) & 0 \end{matrix}] T^{- 1} [\begin{matrix} i_{x} \\ i_{y} \end{matrix}] \\ + T {[\begin{matrix} {2 T}_{q 0}^{'} + Dt & 0 \\ 0 & {2 T}_{d 0}^{'} + Dt \end{matrix}]}^{- 1} [\begin{matrix} 0 \\ 2 Dt \end{matrix}] E_{f} \end{matrix} - - - (14)

其中，

T = [\begin{matrix} \sin δ & \cos δ \\ - \cos & \sin δ \end{matrix}],

δ为第i台发电机组的转子角，即发电机组i的转子相对于发电机组n的转子的角度。

励磁系统数学模型：

\{\begin{matrix} (U_{Ref} - U_{t}) - U_{FP} = U_{1} \\ U_{1} * \frac{K_{A}}{1 + T_{A} * p} = U_{A} \\ (U_{A} - S_{E} * E_{qe}) \frac{1}{K_{E} + T_{E} * p} = E_{qe} \\ E_{qe} * \frac{K_{F} * p}{1 + T_{F} * p} = U_{FP} \end{matrix} - - - (15)

式中，U_Ref、U_t分别为参考电压、发电机端电压，K_E、S_E、T_E分别为励磁机类型常数、励磁机饱和系数及时间常数，K_A、T_A分别为调节器放大倍数及时间常数，K_F、T_F分别为反馈环节放大倍数及时间常数。消去代数量，得：

\{\begin{matrix} (U_{Ref} - U_{t} - U_{FP}) * \frac{K_{A}}{1 + T_{A} * p} = U_{A} \\ (U_{A} - S_{E} * E_{qe}) \frac{1}{K_{E} + T_{E} * p} = E_{qe} \\ E_{qe} * \frac{K_{F} * p}{1 + T_{F} * p} = U_{FP} \end{matrix} - - - (16)

采用隐式梯形积分法差分化：

\begin{matrix} [\begin{matrix} Dt + {2 T}_{A} & K_{A} Dt & 0 \\ - Dt & 0 & {2 T}_{E} + S_{E} Dt + K_{E} Dt \\ 0 & - Dt - {2 T}_{F} & {2 K}_{F} \end{matrix}] [\begin{matrix} U_{A} (t + Δt) \\ U_{FP} (t + Δt) \\ E_{qe} (t + Δt) \end{matrix}] \\ = [\begin{matrix} {2 T}_{A} - Dt & - K_{A} Dt & 0 \\ Dt & 0 & {2 T}_{E} - S_{E} Dt- K_{E} Dt \\ 0 & Dt - 2 T & {2 K}_{F} \end{matrix}] [\begin{matrix} U_{A} (t) \\ U_{FP} (t) \\ E_{qe} (t) \end{matrix}] - [\begin{matrix} {2 K}_{A} Dt \\ 0 \\ 0 \end{matrix}] U_{t} + [\begin{matrix} {2 K}_{A} Dt \\ 0 \\ 0 \end{matrix}] U_{Ref} \end{matrix} - - - (17)

上式化为：

\begin{matrix} [\begin{matrix} U_{A} (t + Δt) \\ U_{FP} (t + Δt) \\ E_{pe} (t + Δt) \end{matrix}] \\ = {[\begin{matrix} Dt + {2 T}_{A} & K_{A} Dt & 0 \\ - Dt & 0 & {2 T}_{E} + S_{E} Dt + K_{E} Dt \\ 0 & - Dt - {2 T}_{F} & {2 K}_{F} \end{matrix}]}^{- 1} [\begin{matrix} {2 T}_{A} - Dt & {- K}_{A} Dt & 0 \\ Dt & 0 & {2 T}_{E} - S_{E} Dt - K_{E} Dt \\ 0 & Dt - 2 T & {2 K}_{F} \end{matrix}] [\begin{matrix} U_{A} (t) \\ U_{FP} (t) \\ E_{qe} (t) \end{matrix}] \\ + {[\begin{matrix} Dt + {2 T}_{A} & K_{A} Dt & 0 \\ - Dt & 0 & {2 T}_{E} + S_{E} Dt + K_{E} Dt \\ 0 & - Dt - {2 T}_{F} & {2 K}_{F} \end{matrix}]}^{- 1} [\begin{matrix} {2 K}_{A} Dt \\ 0 \\ 0 \end{matrix}] U_{t} \\ + {[\begin{matrix} Dt + {2 T}_{A} & K_{A} Dt & 0 \\ - Dt & 0 & {2 T}_{E} + S_{E} Dt + K_{E} Dt \\ 0 & - Dt - {2 T}_{F} & {2 K}_{F} \end{matrix}]}^{- 1} [\begin{matrix} {2 K}_{A} Dt \\ 0 \\ 0 \end{matrix}] U_{Ref} \end{matrix} - - - (18)

其中，机端电压U_t可根据（5）式和（19）式求取：

U_{t} = \sqrt{U_{x} U_{x} + U_{y} U_{y}} - - - (19)

原动机及调速系统数学模型（包括转子运动方程）：

\{\begin{matrix} K_{δ} (ω_{0} - ω) = η \\ η - ξ = ρ \\ ξ = K_{i} Δu + \frac{K_{β} T_{β} p}{1 - T_{β} p} Δu \\ \frac{ρ}{T_{S} p} = Δu \\ Δu = u - u_{0} \\ u * \frac{1 - T_{ω} * p}{1 - 0.5 * T_{ω} * p} = \frac{T_{m}}{K_{mH}} \\ T_{J} * \frac{dω}{dt} = T_{m} - T_{e} - D * ω \\ \frac{dδ}{dt} = ω - ω_{n} \end{matrix} - - - (20)

式中，K_δ、T_S、K_i、K_β、T_i、T_ω、ω₀、ω_n、T_J、D分别为离心飞摆放大倍数、接力器时间常数、硬反馈放大倍数、软反馈放大倍数、软反馈时间常数、水锤效应时间常数、参考角速度、机n的角速度、转子惯性时间常数、转子机械阻尼系数。

消去代数量得：

\{\begin{matrix} \frac{K_{β} T_{β} p}{1 + T_{β} p} (u - u_{0}) = K_{δ} ω_{0} - K_{δ} ω -ρ- K_{i} (u - u_{0}) \\ \frac{ρ}{T_{S} p} = u - u_{0} \\ u * \frac{1 - T_{ω} * p}{1 - 0.5 * T_{ω} * p} = \frac{T_{m}}{K_{mH}} \\ T_{J} * \frac{dω}{dt} = T_{m} - T_{e} - D * ω \\ \frac{dδ}{dt} = ω - ω_{n} \end{matrix} - - - (21)

采用隐式梯形积分法差分化：

\begin{matrix} [\begin{matrix} K_{δ} Dt + {2 K}_{δ} T_{β} & {2 K}_{β} T_{β} + K_{i} Dt + {2 T}_{β} T_{i} & 0 & {2 T}_{β} + Dt \\ 0 & {2 T}_{S} & 0 & - Dt \\ 0 & K_{mH} Dt - {2 K}_{mH} T_{ω} & - Dt - T_{ω} & 0 \\ {2 T}_{J} & 0 & - Dt & 0 \end{matrix}] [\begin{matrix} ω (t + Δt) \\ μ (t + Δt) \\ T_{m} (t + Δt) \\ ρ (t + Δt) \end{matrix}] \\ = [\begin{matrix} {- K}_{δ} Dt + {2 K}_{δ} T_{β} & {2 K}_{β} T_{β} - K_{i} Dt + {2 T}_{β} T_{i} & 0 & {2 T}_{β} - Dt \\ 0 & {2 T}_{S} & 0 & Dt \\ 0 & {- K}_{mH} Dt - {2 K}_{mH} T_{ω} & Dt - T_{ω} & 0 \\ {2 T}_{J} & 0 & Dt & 0 \end{matrix}] [\begin{matrix} ω (t) \\ μ (t) \\ T_{m} (t) \\ ρ (t) \end{matrix}] + [\begin{matrix} 0 \\ 0 \\ 0 \\ - 2 Dt \end{matrix}] T_{e} + [\begin{matrix} {2 K}_{δ} Dt (ω_{0} + μ_{0}) \\ 0 \\ 0 \\ - 2 DtD \end{matrix}] \end{matrix} - - - (22)

上式化为：

其中，作用在转轴上的电磁力矩Te：

T_e=E′_xI_x+E′_yI_y+(X′_q-X′_d)I_xI_y （24）

下一时步的转子角：

δ(t+Δt)=δ(t)+0.5Dt(ω(t+Δt)+ω(t))-0.5(ω_n(t+Δt)+ω_n(t)) （25）

这里假定第n台发电机的dq轴即dq_n为同步坐标轴xy，所有发电机定子绕组电压方程的dq轴变量均转换到dq_n轴上。这里δ为第i台发电机组的转子角，即发电机组i的转子相对于发电机组n的转子的角度。这里，ω_n(t+Δt)及ω_n(t)为第n台发电机组仿真模型计算中的状态变量。

三、仿真模型的交替求解计算

这里仍然以微分方程采用隐式梯形法求解，发电机采用4阶实用模型、励磁系统如（15）式所示、原动机及其调速系统如（20）式所示为例，暂态稳定仿真中某一时步的计算为例，说明本发明中仿真模型的交替求解计算过程。

（1）、根据上一次计算得到的E′_xy的值，采用电网仿真模型（10）式计算发电机注入网络的电流I_xy的值。

（2）、根据同步电机仿真模型本时步的状态量的值及上一次计算得到的E_f、δ的值及I_xy的值，采用同步电机仿真模型（14）式计算下一时步的状态量的值E′_xy。

（3）、采用（5）式和（19）式求取的U_t的值，并根据励磁系统本时步状态量的值，采用励磁系统仿真模型（18）式计算下一时步状态量的值。

（4）、采用（24）式求取的T_e的值，并根据原动机及其调速系统本时步状态量的值，采用原动机及其调速系统仿真模型（23）和（25）式计算下一时步状态量的值。

（5）、根据本次计算得到的机组状态量的值与上次计算得到的机组状态量的值之差的绝对值的大小判断计算是否收敛。

（6）、若不收敛，则以本次计算得到的E′_xy的值，重复（1）～（6）的计算和判断。

（7）、若收敛，则令t=t+Δt，相应各状态量t时刻的值替换为t+Δt时刻的值，重复（1）～（7）进行下一时步的计算。

Claims

1.一种电力系统暂态稳定数字仿真计算模型，其特征在于：电力系统暂态稳定数字仿真数学模型分割为电网数学模型和发电机组数学模型两部分，电网数学模型包括n节点的电网复数线性代数方程YU=I和发电机定子绕组电压方程，其中各发电机定子绕组电压方程中的变量均为同一dq轴上或xy同步坐标轴上描述的变量值，根据此数学模型形成电网仿真模型；发电机组数学模型包括dq轴转换方程和发电机组方程，根据此数学模型形成发电机组仿真模型；

所述发电机组仿真模型，是采用数值积分方法将发电机组的数学模型差分化，并将发电机组的差分化方程组中所包括的发电机定子绕组电压方程中的电势变量和机端电流变量，用电势E′_xy和节点注入电流I_xy来表示，即：将E′_xy和I_xy由同一dq轴或xy轴上根据dq轴转换方程转化至本机的dq轴上，形成各发电机组的仿真模型；

所述发电机组数学模型分割为同步电机数学模型，励磁系统数学模型，原动机及其调速系统数学模型三部分，其中同步电机数学模型包括转子绕组方程和dq轴转换方程，励磁系统数学模型包括励磁系统方程，原动机及其调速系统数学模型包括原动机及其调速系统方程和转子运动方程，采用数值积分方法将同步电机的转子绕组方程差分化，得到其差分化方程组，并进行dq轴转换计算，形成同步电机的仿真模型；同样，采用数值积分方法将励磁系统、原动机及其调速系统的数学模型差分化，分别形成励磁系统仿真模型和原动机及其调速系统仿真模型，同步电机的仿真模型、励磁系统仿真模型和原动机及其调速系统仿真模型一起构成发电机组仿真模型。

2.根据权利要求1所述一种电力系统暂态稳定数字仿真计算模型，其特征在于：所述电网仿真模型，是将n节点的电网复数线性代数方程YU=I，增阶化为同一dq轴上或xy同步坐标轴上的2n维实数线性代数方程Y_xyU_xy=I_xy，将权利要求1所述各发电机定子绕组电压方程和增阶网络方程Y_xyU_xy=I_xy一起联立，整理后可以形成如下方程：I_xy=AB_xyE′_xy，即为电网仿真模型，为时不变的实数线性代数方程，AB_xy为2n维线性时不变矩阵，I_xy为2n维节点注入电流向量，E′_xy为发电机定子绕组电压方程中的电势，为2n维向量。

3.一种采用如权利要求1所述的一种电力系统暂态稳定数字仿真计算模型的仿真计算方法，其特征由以下步骤构成：

步骤三：根据E′_xy及I_xy的值求取机端电压U_t的值，并根据励磁系统本时步状态量的值，采用励磁系统仿真模型计算励磁系统下一时步状态量的值；

步骤六：若不收敛，则以本次计算得到的E′_xy的值，重复步骤一到步骤六的计算和判断；步骤七：若收敛，则令t=t+Δt，相应各状态量t时刻的值替换为t+Δt时刻的值，重复步骤一到步骤七进行下一时步的计算。