惯性是系统稳定运行的基础保障^[1-2]，同步机惯性是系统的固有属性，它能够在系统受到大扰动后延缓转子电频率的偏移速率，为后续调频手段的启动争取时间^[3-4]. 惯量是描述惯性大小的物理量，是衡量电力系统频率稳定的重要参数，在频率变化的安全范围内，可以告知操作人员系统能承受的冲击功率^[5-7]. 新能源发电设备大量接入电力系统发电侧，部分传统的发电设备被替代，导致系统发电侧惯量逐渐下降，发生切机、短路和负荷变化等扰动事件的频率变化程度增加^[8-9]. 系统发电侧的惯量评估工作非常重要，准确刻画发电侧惯量是频率特性分析与调频策略制定的基础，对维护电力系统的稳定性具有重要意义.

学者针对系统发电侧的惯量评估开展了诸多研究. Cao等^[10]对发电侧惯量进行静态计算，认为发电侧惯量等于相同容量下所包含各发电机惯量的加权平均值；Guo等^[11]基于加权最小二乘法对发电侧惯量进行静态计算. 上述研究计算发电侧惯量的前提是已知单个发电机的惯量，否则无法计算发电侧惯量. 基于发电侧动态响应，直接对发电侧惯量参数进行辨识的方法随后兴起. 该辨识方法主要分为1）基于暂态响应信号的惯量辨识方法，2）基于小扰动信号的惯量辨识方法. 方法1）通常先采用如切负荷或切机的操作激发系统有功频率暂态特性，再借助发电机转子运动方程完成辨识估计. Inoue等^[12-14]结合发电机转子运动方程，利用频率测量数据和功率变化直接估计发电侧惯量. Wall等^[15-16]提出的在线算法估计了系统扰动发生的时间及系统经过扰动事件后的惯量，该算法须连续在线处理有功功率和频率变化率的滑动窗数据. Del Giudice等^[17]提出基于拓展卡尔曼滤波器的发电侧惯量估计方法，该方法必须在系统扰动发生的瞬间及时启动，否则会在惯量估计过程中引入较大误差. You等^[18]基于机电波理论，推导出机电波传播速度与惯量成反比关系，通过计算系统各节点机电波的传播速度，给出美国东部电网的惯量分布图. 方法1）存在局限：大扰动的扰动事件会改变系统的运行工作点，影响系统安全稳定运行. 在方法2）方面，陈志杰等^[19-21]根据电力系统运行时的频率变化情况，利用机器学习方法实现发电侧的惯量评估. 方法2）对历史数据的需求量大，须在已知的发电机惯量基础上迭代计算，较为复杂.

本研究利用同步相量测量单元(phase measurement unit, PMU)测量的机电振荡响应，提出基于小扰动下机电振荡参数的发电侧惯量评估方法. 1）基于各发电机与等值系统的惯量响应过程机理，推导在惯量响应结束时刻与各发电机惯量有关的不平衡功率分配公式，为求解系统中的发电侧惯量提供条件. 2）分析小扰动下发电机的机电振荡过程，得到系统特征方程的解与发电机功角振荡角频率的关系，利用韦达定理求取系统中各台发电侧的惯量，介绍单机系统发电侧惯量的计算方法. 3）在单机与多机系统仿真算例中验证所提评估方法的正确性.

忽略阻尼项，系统中第i台同步发电机的转子运动方程^[22]为

(1)$ \frac{2 H_i}{\omega_0} \frac{\mathrm{~d}^2 \delta_i}{\mathrm{~d} t^2}=2 H_i \frac{\mathrm{~d} \omega_i}{\mathrm{~d} t}=2 H_i \frac{\mathrm{~d} f_i}{\mathrm{~d} t}=P_{\mathrm{m} i}-P_{\mathrm{e} i}. $

式中：$ H_{{i}} $为第i台同步发电机惯量，$ \omega_{0} $为正常运行时同步发电机转子角速度，$ \delta_{i} $为第i台同步发电机的功角，t为时间，$ \omega_{i} $为第i台同步发电机转子角速度，$ f_{i} $为第i台同步发电机频率，$ P_{\mathrm{m} i} $为第i台发电机机械功率，$ P_{\mathrm{e} {i}} $为第i台发电机电磁功率. 假设在发电机稳态运行时段$ t \sim t+\Delta t $内，系统出现微小阶跃负荷扰动$ \Delta P_{{\mathrm{L}}} $，式中的机械功率维持不变，其余变量在稳定运行点附近波动，增量方程为

(2)$ \begin{split} \frac{{2{H_i}}}{{{\omega _0}}}\frac{{{{\text{d}}^{\text{2}}}\Delta {\delta _i}}}{{{\text{d}}{t^{\text{2}}}}} =&\; 2{H_i}\frac{{{\text{d}}\Delta {\omega _i}}}{{{\text{d}}t}} = 2{H_i}\frac{{{\text{d}}\Delta {f_i}}}{{{\text{d}}t}} = \\ &\Delta {P_{{\text{m}}i}} - \Delta {P_{{\text{e}}i}}=- \Delta {P_{{\text{e}}i}} ;\\\frac{\mathrm{d}^2 \Delta \delta_i}{\mathrm{~d} t^2}= & \;\frac{\mathrm{d}^2 \delta_i(t+\Delta t)}{\mathrm{d} t^2}-\frac{\mathrm{d}^2 \delta_i(t)}{\mathrm{d} t^2}, \\ \frac{\mathrm{~d} \Delta \omega_i}{\mathrm{~d} t}=&\;\frac{\mathrm{d} \omega_i(t+\Delta t)}{\mathrm{d} t}-\frac{\mathrm{d} \omega_i(t)}{\mathrm{d} t}, \\ \frac{\mathrm{~d} \Delta f_i}{\mathrm{~d} t}=&\;\frac{\mathrm{d} f_i(t+\Delta t)}{\mathrm{d} t}-\frac{\mathrm{d} f_i(t)}{\mathrm{d} t}, \\ \Delta P_{\mathrm{m} i}=&\;P_{\mathrm{m} i}(t+\Delta t)-P_{\mathrm{m} i}(t)=0, \\ \Delta P_{\mathrm{e} i}=&\;P_{\mathrm{e} i}(t+\Delta t)-P_{\mathrm{e} i}(t), \quad \Delta P_{\mathrm{L}}=\sum_{i=1}^n \Delta P_{\mathrm{e} i} .\end{split} $

当不附加频率控制的新能源发电设备接入系统时，会造成系统等效惯量$ H_{\mathrm{sys}} $降低. 不计及新能源发电设备的惯量贡献，系统等效惯量$ H_{\mathrm{sys}} $的计算式为

(3)$ {H_{{\text{sys}}}} = {{\displaystyle\sum\limits_{i = 1}^n {({H_i} \cdot {S_i})} }}\Biggr/\left({{\displaystyle\sum\limits_{i = 1}^n {{S_i}} +\displaystyle\sum\limits_{j = 1}^q {{S_j}} }}\right) . $

式中：$ n $、$ q $分别为系统中同步发电机和新能源发电设备的数量，$ S_{{i}} $、$ S_{j} $分别为同步发电机$ {i} $和新能源发电设备$ \dot{j} $的容量. 当不附加频率控制的新能源发电设备接入发电侧时，发电侧的等效惯量由该侧同步发电机惯量聚合得到，发电侧等效惯量${H'_{{\mathrm{sys}}}}$的计算式为

(4)$ {H'_{{\mathrm{sys}}}} = {{\displaystyle\sum\limits_{i = 1}^{n'} {({H_i} \cdot {S_i})} }}\Biggr/\left({{\displaystyle\sum\limits_{i = 1}^{n'} {{S_i}+\displaystyle\sum\limits_{i = 1}^{q'} {{S_j}} } }}\right). $

式中：$n'$、$q'$分别为发电侧同步发电机和新能源发电设备的数量.

当系统中出现扰动事件时，同步发电机组参与电网调频按时间顺序分为3个阶段：惯量响应阶段、一次调频阶段和二次调频阶段. 加入扰动类型为小扰动事件，一、二次调频一般不会发生动作，同步电机的频率变化特性主要由电机的惯量决定，意味着频率变化过程主要体现在惯量响应阶段，因此仅阐述同步电机的惯量响应过程. 当系统中有扰动事件发生时，会造成机械功率与电磁功率不平衡，此时同步电机的转子通过减速（加速）释放（吸收）旋转动能向系统中注入（吸收）有功功率，阻碍系统频率波动，这一阶段称为惯量响应阶段^[23]. 由式(1)可知，频率变化与功角变化密切相关，本研究主要关注惯量响应过程中的功角变化过程，利用功角振荡曲线评估发电侧的惯量.

如图1所示，多机系统中网络节点数为$ n+m+ 1 $，其中发电机内电势节点个数为n($ n \geqslant 2 $)，接地节点为0号节点，其他节点个数为m. 以此系统为例，推导多机系统的发电侧惯量. 以节点0(接地点)为参考节点，节点电压方程^[24]为

(5)$ \left[ \begin{gathered} {{\dot {\boldsymbol{I}}}_n} \\ {\boldsymbol{0}} \\ \end{gathered} \right] = {{\boldsymbol{Y}}_{\mathrm{s}}}\left[ \begin{gathered} {{\dot {\boldsymbol{E}}}_n} \\ {{\dot {\boldsymbol{V}}}_m} \\ \end{gathered} \right] = \left[ {\begin{array}{*{20}{c}} {{Y_{nn}}}&{{Y_{nm}}} \\ {{Y_{mn}}}&{{Y_{mm}}} \end{array}} \right]\left[ \begin{gathered} {{\dot {\boldsymbol{E}}}_n} \\ {{\dot {\boldsymbol{V}}}_m} \\ \end{gathered} \right] . $

式中：$ {\dot {\boldsymbol{I}}_n} $为发电机向系统注入的电流向量；$ {\dot {\boldsymbol{E}}_n} $、$ {\dot {\boldsymbol{V}}_m} $分别为发电机内电势节点和网络节点的电压列向量；$ {\boldsymbol{Y}}_{{\mathrm{s}}} $为系统的增广导纳矩阵，将$ {\boldsymbol{Y}}_{{\mathrm{s}}} $适当分解成4部分，即$ Y_{nn} $、$ Y_{n m} $、$ Y_{mn} $、$ Y_{m m} $. 将式(1)展开，消去$ {\dot {\boldsymbol{V}}_m} $，得到简化矩阵Y：

(6)$ {\dot {\boldsymbol{I}}_n} = {\boldsymbol{Y}}{\dot {\boldsymbol{E}}_n} = ({Y_{nn}} - {Y_{nm}}Y_{mm}^{ - 1}{Y_{mn}}){\dot {\boldsymbol{E}}_n} . $

矩阵Y的元素$ G_{ii} $、B_ii分别为发电机i的内电势节点的自电导和自电纳，$ G_{i j} $、$ B_{ij} $分别为发电机i、j之间的转移电导和电纳. 系统中第i台发电机电磁功率的表达式^[25]为

(7)$ P_{{\mathrm{e}i}}=E_i^2 G_{i i}+\sum_{j=1, j \neq i}^n E_i E_j\left(G_{i j} \cos \delta_{i j}+B_{i j} \sin \delta_{i j}\right). $

式中：$ \delta_{i j}=\delta_{i}-\delta_{j} $. 小信号线性化增量方程为

(8)$ \Delta P_{\mathrm{e} i}=P_{\mathrm{e} i}(t+\Delta t)-P_{\mathrm{e} i}(t)=\sum_{j=1, j \neq i}^n D_{i j} \Delta \delta_{i j}. $

式中：$\Delta \delta_{i j}=\delta_{i j}(t+\Delta t)-\delta_{i j}(t)$，$ D_{i j}=E_i E_j \left(B_{i j} \cos \delta_{i j}(t)-\right. \left.G_{i j} \sin \delta_{i j}(t)\right)$. 把式(2)代入式(8)，得到第i台发电机小信号增量方程为

(9)$ \frac{2 H_i}{\omega_0} \frac{\mathrm{~d}^2 \Delta \delta_i}{\mathrm{~d} t^2}+\sum_{j=1, j \neq i}^n D_{i j} \Delta \delta_{i j}=0. $

由于$ \delta_{12}+\delta_{23}+\cdots+\delta_{(n-1) n}+\delta_{n 1}=0 $，可知式(9)不是n个独立的二阶微分方程，而是$ n-1 $个独立方程组. 将第$ {i} $($ i \neq n $)个方程减去第n个方程得到

(10)$ \begin{split}& \frac{\mathrm{d}^2 \Delta \delta_{i n}}{\mathrm{~d} t^2}+\frac{\omega_0}{2 H_i} \sum_{j=1, j \neq i}^n D_{i j} \Delta \delta_{i n}- \\&\quad \frac{\omega_0}{2 H_i} \sum_{j=1, j \neq i}^n D_{i j} \Delta \delta_{j n}+\frac{\omega_0}{2 H_n} \sum_{j=1}^{n-1} D_{n j} \Delta \delta_{j n}=0.\end{split} $

定义$ \boldsymbol{X}\;=\;\left[\;\Delta \delta_{1 n},\; \Delta \delta_{2 n},\; \cdots,\; \Delta \delta_{(n-1) n},\; \Delta \omega_{1 n},\; \Delta \omega_{1 n},\; \cdots,\; \Delta \omega_{(n-1) n}\;\right]^{\mathrm{T}}$，其中$ \Delta \omega_{i n}=\Delta \dot{\delta}_{i n}$，由式(9)得到小信号线性化状态方程为

(11)$ \begin{split} &\quad\quad\quad\quad\quad\quad\dot{\boldsymbol{X}}=\boldsymbol{A X},\\&{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {{{\mathbf{0}}_{(n - 1) \times (n - 1)}}}&{{{\boldsymbol{E}}_{(n - 1) \times (n - 1)}}} \\ {\boldsymbol{M}}&{{{\mathbf{0}}_{(n - 1) \times (n - 1)}}} \end{array}} \right] .\end{split}$

式中：$ {{\boldsymbol{E}}_{(n - 1) \times (n - 1)}} $为单位矩阵，

$ \boldsymbol{M}=\left[ \begin{array}{ccccc}\dfrac{\omega_0}{2 H_1} \displaystyle\sum_{j=2}^n D_{1 j}+\dfrac{\omega_0}{2 H_n} D_{n 1} & -\dfrac{\omega_0}{2 H_1} D_{12}+\dfrac{\omega_0}{2 H_n} D_{n 2} & -\dfrac{\omega_0}{2 H_1} D_{13}+\dfrac{\omega_0}{2 H_n} D_{n 3} & \ldots & -\dfrac{\omega_0}{2 H_1} D_{1(n-1)}+\dfrac{\omega_0}{2 H_n} D_{n(n-1)} \\-\dfrac{\omega_0}{2 H_2} D_{21}+\dfrac{\omega_0}{2 H_n} D_{n 1} & \dfrac{\omega_0}{2 H_2} \displaystyle\sum_{j=1, j \neq 2}^n D_{2 j}+\dfrac{\omega_0}{2 H_n} D_{n 2} & -\dfrac{\omega_0}{2 H_2} D_{23}+\dfrac{\omega_0}{2 H_n} D_{n 3} & \ldots & -\dfrac{\omega_0}{2 H_2} D_{2(n-1)}+\dfrac{\omega_0}{2 H_n} D_{n(n-1)} \\-\dfrac{\omega_0}{2 H_3} D_{31}+\dfrac{\omega_0}{2 H_n} D_{n 1} & -\dfrac{\omega_0}{2 H_3} D_{32}+\dfrac{\omega_0}{2 H_n} D_{n 2} & \dfrac{\omega_0}{2 H_3} \displaystyle\sum_{j=1, j \neq 3}^n D_{3 j}+\dfrac{\omega_0}{2 H_n} D_{n 3} & \ldots & -\dfrac{\omega_0}{2 H_3} D_{3(n-1)}+\dfrac{\omega_0}{2 H_n} D_{n(n-1)} \\\vdots & \vdots & \vdots & {\text{ }} & \vdots \\-\dfrac{\omega_0}{2 H_{n-1}} D_{(n-1) 1}+\dfrac{\omega_0}{2 H_n} D_{n 1} & -\dfrac{\omega_0}{2 H_{n-1}} D_{(n-1) 2}+\dfrac{\omega_0}{2 H_n} D_{n 2} & -\dfrac{\omega_0}{2 H_{n-1}} D_{(n-1) 3}+\dfrac{\omega_0}{2 H_n} D_{n 3} & \cdots & \dfrac{\omega_0}{2 H_{n-1}} \displaystyle\sum_{j=1, j \neq n-1}^n D_{(n-1) j}+\dfrac{\omega_0}{2 H_n} D_{n(n-1)}\end{array} \right] . $

式(11)的特征值方程为

(12)$ \left| {\lambda _{(n - 1) \times (n - 1)}^2{{\boldsymbol{E}}_{(n - 1) \times (n - 1)}} - {\boldsymbol{M}}} \right| = 0. $

由式(12)可以得到$ n-1 $对共轭复特征根，特征根的实部为0，虚部为$ \Delta \delta_{{in }} $的固有振荡频率，则该系统有$ n-1 $个固有振荡角频率，记为$ w_{{\mathrm{z}} i n}(i \neq n) $. 记$ \lambda_{i}^{2}(i \neq n) $为系统状态方程每对共轭复特征根的平方，则$ \lambda_{{i}}^{2} $与$ w_{{{\mathrm{z}}in}} $关系为

(13)$ \lambda_i^2=-w_{{\mathrm{z}} i n}^2(i \neq n). $

系统的特征方程表示为

(14)$ \left(\lambda-\lambda_1^2\right)\left(\lambda-\lambda_2^2\right) \cdots\left(\lambda-\lambda_{n-2}^2\right)\left(\lambda-\lambda_{n-1}^2\right)=0. $

将式(14)展开得到

(15)$ a_{n-1} \lambda^{2(n-1)}+a_{n-2} \lambda^{2(n-2)}+\cdots+a_1 \lambda^{2 \cdot 1}+a_0=0. $

式中：${a_0},{a_1}, \cdots ,{a_{n - 2}},{a_{(n - 1)}} \in {\bf{R}}$. 对式(15)运用韦达定理得到

(16)$ \begin{split}\lambda_1^2+\lambda_2^2+\cdots+\lambda_{n-1}^2 & =-\frac{a_{n-2}}{a_{n-1}}= \\& -\frac{w_0}{2} \sum_{i=1}^{n-1}\left( \sum_{j=1, j \neq i}^n \frac{D_{i j}}{H_i}+\frac{D_{n i}}{H_n} \right) .\end{split} $

其中$ w_{0} $、$ D_{{ij}} $由系统稳态运行的参数求得，通过$ n-1 $个固有振荡频率可以得到关于H_i的$ n-1 $个方程. 系统中有n个H_i，须找到发电机惯量满足的另一个方程. 对于式(2)，将系统中所有的同步发电机转子运动方程式相加得到

(17)$ 2 \sum_{i=1}^n H_i \frac{\mathrm{~d} \Delta \omega_i}{\mathrm{~d} t}=-\sum_{i=1}^n \Delta P_{\mathrm{e} i}=-\Delta P_{\mathrm{L}}. $

定义$ \omega_{{\mathrm{sys}}} $为系统转子角速度的等效动态过程或平均动态过程^[26]. 在惯量响应阶段，各发电机的$ {\mathrm{d} \Delta \omega_i}/ {\mathrm{d} t} $逐渐与$ {\mathrm{d} \Delta \omega_{{\mathrm{sys}}}}/{\mathrm{d} t}$相等，$ \omega_{{\mathrm{sys}}}$的增量形式为

(18)$ \begin{split}&\Delta \omega_{ {{\mathrm{sys}} }}={\frac{1}{H_{\text {sys }}}\displaystyle\sum_{i=1}^n H_i \Delta \omega_i}, \\&\sum_{i=1}^n H_i \Delta \omega_i=\sum_{i=1}^n H_i \omega_i(t+\Delta t)-\sum_{i=1}^n H_i \omega_i(t).\end{split} $

将式(17)代入式(18)，得到

(19)$ \dfrac{\mathrm{d} \Delta \omega_{{\mathrm{sys}}}}{\mathrm{d} t}=\dfrac{-\Delta P_{\mathrm{L}}}{2H_{{\mathrm{sys}}}}. $

当各台发电机转子的$ {\mathrm{d} \Delta w_i}/{\mathrm{d} t}$与$ {\mathrm{d} \Delta \omega_{{\mathrm{sys}}}}/{\mathrm{d} t}$相等时，结合式(17)和式(19)得到

(20)$ \Delta P_{\mathrm{e} i}^{\prime}={H_i} \Delta P_{\mathrm{L}}\Biggr/{\displaystyle\sum_{i=1}^n H_i}. $

定义惯量比为$ c_i(i=1,2, \cdots, n)$，表达式为

(21)$ c_i={H_i}\Biggr/{\displaystyle\sum_{i=1}^n H_i}={\Delta P_{\mathrm{e} i}^{\prime}}/{\Delta P_{\mathrm{L}}} . $

联立求解式(16)和式(21)，得到

(22)$ \left. \begin{gathered} {H_1} = {{{c_1}{\omega _0}\displaystyle\sum\limits_{i = 1}^{n - 1} {\left(\displaystyle\sum\limits_{j = 1,j \ne i}^n {\frac{{{D_{ij}}}}{{{c_i}}}} +\frac{{{D_{ni}}}}{{{c_n}}}\right)} }}\Biggr/\left({2{\displaystyle\sum\limits_{i = 1}^{n - 1} {\omega _{zin}^2} }}\right), \\ {H_2} = {{{c_2}{\omega _0}\displaystyle\sum\limits_{i = 1}^{n - 1} {\left(\displaystyle\sum\limits_{j = 1,j \ne i}^n {\frac{{{D_{ij}}}}{{{c_i}}}} +\frac{{{D_{ni}}}}{{{c_n}}}\right)} }}\Biggr/\left({2{\displaystyle\sum\limits_{i = 1}^{n - 1} {\omega _{zin}^2} }}\right), \\ {\text{ }} \quad\quad\quad\quad\quad\quad\quad\quad\vdots \\ {H_n} = {{{c_n}{\omega _0}\displaystyle\sum\limits_{i = 1}^{n - 1} {\left(\displaystyle\sum\limits_{j = 1,j \ne i}^n {\frac{{{D_{ij}}}}{{{c_i}}}} +\frac{{{D_{ni}}}}{{{c_n}}}\right)} }}\Biggr/\left({2{\displaystyle\sum\limits_{i = 1}^{n - 1} {\omega _{zin}^2} }}\right). \\ \end{gathered} \right\} $

式(22)是多机系统的各发电机惯量计算式. 计算各发电机惯量，关键是根据系统中各台发电机的有功功率曲线和功角振荡曲线求出惯量比c_i和固有振荡频率$ \omega_{{\mathrm{z}}in} $.

简单描述单机无穷大系统发电侧的惯量求取如下. 在单机无穷大系统，式(1)的增量方程^[25]为

(23)$ \begin{split} &\frac{{2{H_{\text{1}}}}}{{{\omega _0}}}\frac{{{{\text{d}}^{\text{2}}}\Delta {\delta _{\text{1}}}}}{{{\text{d}}{t^{\text{2}}}}}+{P_{\max }}\cos {\delta _{\text{1}}}(t)\Delta {\delta _{\text{1}}} = {\text{0}}, \\& \frac{\mathrm{d}^2 \Delta \delta_1}{\mathrm{~d} t^2}=\frac{\mathrm{d}^2 \delta_1(t+\Delta t)}{\mathrm{d} t^2}-\frac{\mathrm{d}^2 \delta_1(t)}{\mathrm{d} t^2} .\end{split} $

式中：$ P_{\text {max }} $为同步发电机的最大电磁功率. 求式(23)二阶微分方程特征根的虚部${\omega _{z1}}$=$\sqrt {{{{\omega _0}{P_{\max }}\cos {\delta _{\text{1}}}(t)}}/{(2{{H_1}})}} $，${\omega _{z1}}$为发电机的固有振荡频率. ${\omega _0}$、$ P_{\text {max }} $、$ \cos \delta_{10} $均为系统稳态运行时的参数，可以直接求取. 发电侧的惯量为

(24)$ {H_1} = {{{\omega _0}{P_{\max }}\cos {\delta _{10}}}}/(2{{\omega _{z1}^2}}). $

由式(24)可知，求取单机无穷大系统发电机的惯量，关键是根据功角振荡曲线求出固有振荡频率$ {\omega _{{\mathrm{z}}{\text{1}}}} $.

同步PMU广泛应用于监测电力系统的动态过程^[27-28]，能够测量系统中各台发电机的有功功率曲线和功角振荡曲线. 假设在$t = 0\;{\text{s}}$时系统中出现微小阶跃负荷扰动，在$ t=t_{{\mathrm{s}}} $时$ {\mathrm{d} \Delta \omega_i}/{\mathrm{d} t}={\mathrm{d} \Delta \omega_{{\mathrm{sys}}}}/{\mathrm{d} t} $，由于惯量响应时间一般在扰动发生后$ 0\sim 2\; {\mathrm{s}} $，则${t_{\mathrm{s}}} < 2\;{\text{s}}$或${t_{\mathrm{s}}} > 2\;{\text{s}}$. 当${t_{\mathrm{s}}} < 2\;{\text{s}}$时，在$t = 2\;{\text{s}}$时刻，有$ {\mathrm{d} \Delta \omega_i}/ {\mathrm{d} t}= {\mathrm{d} \Delta \omega_{{\mathrm{sys}}}}/{\mathrm{d} t}$；当${t_{\mathrm{s}}} > 2\;{\text{s}}$时，在$t = 2\;{\text{s}}$时刻，有$ {\mathrm{d} \Delta \omega_i}/ {\mathrm{d} t} \approx {\mathrm{d} \Delta \omega_{{\mathrm{sys}}}}/{\mathrm{d} t}$. 为了统一起见，假设在扰动发生后的第2 s，$ {\mathrm{d} \Delta \omega_i}/{\mathrm{d} t}={\mathrm{d} \Delta \omega_{{\mathrm{sys}}}}/{\mathrm{d} t}$. 在$ t=0^{+}$时，发电机承担的不平衡功率与同步功率系数有关^[29]，在$t = 2\;{\text{s}}$时发电机承担的不平衡功率与$ H_i$有关. 因此，在时间段$ 0 \sim 2\;{\text{s}}$，发电机瞬时有功功率改变量曲线($ \Delta P_{\mathrm{e} i}$关于时间的曲线)的振荡中心在时刻变化，在$t = 2\;{\text{s}}$时，振荡中心为$ \Delta P'_{\text{e} i} $. 求取$ \Delta P'_{\mathrm{e} i}$的方法：记最接近$t = 2\;{\text{s}}$之前的发电机瞬时有功功率改变量曲线波峰与功率波谷分别为$ \Delta P_{\text {e}i{\mathrm{t}}} $、$ \Delta P_{{\mathrm{e}}i{\mathrm{l}}} $，则$ \Delta P'_{{\mathrm{e}}i} $的表达式为

(25)$ \Delta P_{\mathrm{e} i}^{\prime}=\frac{\Delta P_{\mathrm{e} i {\mathrm{t}}}+\Delta P_{\mathrm{e} i {\mathrm{l}}}}{2}. $

依据式(21)得到惯量比$ c_{i}(i=1,2, \cdots, n) $. 将${\omega _{{\mathrm{z}}in}}$、${\omega _{{\text{z1}}}}$统称为发电机功角(或相对功角)振荡曲线的固有振荡角频率${\omega _{\mathrm{z}}}$. 记发电机功角(或相对功角)振荡曲线在$t = 0 \sim 2\;{\text{s}}$阶段相邻波峰波谷对应时刻分别为$ t_{1} $、$ t_{2} $，

(26)$ {\omega _{\mathrm{z}}} = \frac{2}{{{t_{\mathrm{z}}}}} = \frac{{2{\text{π }}}}{{2\left| {{t_1} - {t_2}} \right|}} = \frac{{\text{π }}}{{\left| {{t_1} - {t_2}} \right|}}. $

在Matlab/Simulink平台搭建单机无穷大系统，结构图如图2所示. 设置系统基准容量为100 MW，发电机内部阻抗仅由暂态电抗组成，参数如表1所示. 表中，参数无单位标记的均为标幺值. 在发电机侧附近加入5 MW阶跃负荷(有功功率为系统基准容量的5%)，得到发电机的功角曲线如图3所示. 由式(24)得到系统发电机的惯量H=4.57 s，误差为8.6%，与实际发电机惯量接近，验证了本研究所提评估方法的正确性.

在Matlab/Simulink平台建立双机互联系统，结构图如图4所示. 设置系统基准容量均为100 MW，各发电机内部阻抗仅由暂态电抗组成，参数如表2所示. 表中，参数无单位标记的均为标幺值. 在负载处接入5 MW阶跃负荷(有功功率为负荷功率的5%)，得到各发电机瞬时有功功率改变量曲线与发电机1、2相对功角曲线分别如图5、6所示，根据式(22)得到双机惯量结果如表3所示. 由表3可知，测量值与实际值接近，验证了所提评估方法的正确性.

在Matlab/Simulink平台搭建WSCC3机9节点系统，结构图如图7所示. 该系统有3台发电机、9条母线和3处负荷，系统详细信息参见文献[30]. 各发电机额定容量、惯量如表4所示. 在节点6处接入20 MW的阶跃负荷(有功功率约为负荷功率的5%)，各台发电机瞬时有功功率改变量曲线与各发电机相对功角曲线分别如图8、图9所示. 根据式(22)得到的惯量评估结果如表5所示. 可以看出，测量值与实际值接近，验证了所提评估方法的正确性.

为了验证所提评估方法在新能源接入发电侧时评估等效惯量的正确性，在2号发电机处以1台相同容量且不附加频率控制的风力发电机组代替其中一台同步发电机组. 2号发电侧的实际等效惯量由式(4)计算得到，利用所提出惯量评估的方法求得发电侧惯量，验证评估方法的正确性. 在节点6处加入20 MW的冲击负荷(有功功率约为负荷功率的5%)，如图10、图11所示分别为各台发电机瞬时有功功率改变量曲线、发电机相对功角曲线，如表6所示为各发电侧的惯量评估结果. 可以看出，测量值与实际值接近，验证了所提方法的正确性.

在Matlab/Simulink平台搭建10机39节点系统，结构图如图12所示. 该系统有10台发电机、39条母线和19处负荷，系统详细信息参见文献[31]. 在5号发电机机组添加不附加频率控制的新能源发电设备，5号发电机等效惯量由式(4)计算得到. 各发电机以系统额定容量(1 000 MW)为基准的惯量如表7所示. 使节点39处的负荷有功功率突增300 MW(约为系统总负荷功率的5%)，各发电机有功功率改变量$ \Delta P_{\mathrm{e}i}^{\prime} $如表8所示，2~9号发电机与1号发电机相对功角曲线振荡角频率${\omega _z}$如表9所示. 由惯量评估式(22)计算得到系统各发电机或等效发电机惯量如表10所示. 由表10可以看出，测量值与实际值接近，验证了所提评估方法的正确性.

本研究提出基于小扰动下机电振荡参数的发电侧惯量评估方法，在单机系统与多机系统仿真算例中验证了所提评估方法的正确性，为发电侧惯量优化配置提供了一定的理论基础. 1) 当系统中出现扰动时，各发电机先进入惯量响应阶段，各发电机的转子角速度变化率逐渐与系统等效角速度变化率相等，在惯量响应结束时刻，系统的不平衡功率分配与各发电机惯量有关. 2) 对于单机系统，通过单台发电机转子运动的增量方程可以得到发电侧惯量计算式，根据发电机功角振荡曲线求取固有振荡频率，得到单机无穷大系统的发电侧惯量. 3) 对于多机系统，通过多机系统的小信号状态方程与特征根的关系，结合惯量响应结束时刻与系统中各发电机惯量有关的不平衡功率分配公式得到发电侧惯量计算式. 根据发电机瞬时有功功率改变量曲线、相对功角振荡曲线得到惯量比与固有振荡频率，得到多机系统的发电侧惯量. 后续将深入研究计及新能源设备调频动态的发电侧惯量评估方法.