索江镭，胡志坚，刘宇凯，张子泳，王 尉

（1.武汉大学 电气工程学院，湖北 武汉 430072；2.西安交通大学 电力设备电气绝缘国家重点实验室，陕西 西安 710049；3.广州供电局有限公司，广东 广州 510260）

0 引言

近年来，随着“西电东送”战略的实施以及特高压交直流示范工程的建成，我国电网已初步形成了大规模、跨区域远距离的输电模式。这大幅提升了我国能源的优化调配能力，而同时也增加了系统发生区域低频振荡的危险，给电力系统造成了极大的安全隐患[1]。广域测量系统 WAMS（Wide Area Measurement System）的发展与应用给监测和控制系统低频振荡创造了有利的条件，基于WAMS的广域阻尼控制器 WADC（Wide Area Damping Controller）成为了解决系统区域振荡问题的有效方法之一[2]。目前我国包括清华大学、武汉大学、南方电网科学研究院、北京四方继保自动化股份有限公司等科研技术单位正致力于WADC的实际工程应用，这是实现未来智能电网“自愈性”的重要措施之一[3]。

目前，对于WADC的设计方法大多基于经典控制理论和现代控制理论，主要有极点配置法、留数法、鲁棒控制法、线性矩阵不等式法等[4]。这些方法需首先建立系统的详细模型，然后在系统稳定运行点进行线性化，再利用降阶技术得到系统状态空间模型，最后完成控制器设计。对于大型互联电力系统而言，其详细模型往往难以准确建立，并且随着电网规模的不断发展，大量新能源发电、FACTS、HVDC等先进技术的广泛应用也将增加系统建模的难度。而在实际运行中，系统参数的调整以及运行方式的变化也会对基于建模法得到的WADC产生不利影响[5]。因此，应寻求更为实用有效的模型处理方法，以满足WADC的实际工程应用。

系统辨识理论的发展给WADC模型的建立提供了有效的解决办法，由于其避免了建模法的诸多弊端，因此在近年来受到了广泛的关注[6]。文献[7]采用Prony算法对互联电力系统传递函数进行辨识，指出Prony算法具有计算速度快、辨识精度高的优点，但同时也存在对噪声敏感度大、对系统初始状态要求高的缺点。文献[8]采用了子空间辨识方法，该方法不需要进行迭代优化，仅依靠一些简单的线性算法，具有实现容易的优点，但由于该算法在辨识过程中需进行QR分解和奇异值分解，在数据量较大时，其计算量将剧增，这对于在线辨识极为不利。文献[9]采用了自回归滑动平均（ARMA）模型辨识方法，该方法具有对激励要求低、方法简便等优点，但在辨识过程中需通过最小二乘算法求解模型的自回归部分和滑动平均部分，若遇到系统输入、输出点较多时，将耗费大量时间用于计算协方差矩阵，这限制了该方法在电力系统的实际应用。

近年来，辨识领域出现了如辅助模型辨识思想、多新息辨识思想、耦合辨识思想等一系列新方法[10]，给电力系统状态空间辨识带来了新的思路。鉴于此，本文从理论上分析了采用传统最小二乘算法用于多输入多输出系统辨识时造成计算量大的根本原因，引入了耦合最小二乘算法解决上述问题，结合多新息思想推导了多新息耦合最小二乘算法的计算方法，实现了利用类噪声信号完成电力系统状态空间辨识。4机2区以及10机39节点系统验证了本文方法的有效性。

1 电力系统状态空间辨识方法

1.1 问题描述

电力系统在稳定运行点附近通常可用如下状态空间模型描述：

其中，AєRn×n、BєRn×r、CєRm×n为系统矩阵；x（t）єRn为状态变量；u（t）єRr为输入变量；y（t）єRm为输出变量。通过Laplace变换后得系统传递函数为：

其中，adj[]表示求矩阵的伴随矩阵；det[]表示求矩阵的行列式；I为单位矩阵；α（s）和 Q（s）分别为系统特征多项式和矩阵多项式，其表达式如式（3）所示。

将系统转化为传递函数模型，则有：

若系统存在随机干扰项 v（t）єRm，则有：

此时，将式（3）代入式（5），则有：

定义参数矩阵θ、α、输入信息向量φ（t）以及输出信息矩阵ψ（t）如下：

则式（6）可转化为如下形式：

其中，θ和α为待辨识参数矩阵。为方便计算，将参数矩阵θ和α合并为 ϑ，并定义信息矩阵Φ（t）如下：

其中，θi（i=1，2，…，m）为 θ 的第 i列；“⊗”表示求向量积；Im为单位矩阵，下标代表维数。

则式（7）可表示为：

上式即为本文所采用的待辨识模型，对其展开可得：

其中，ψi（t）为输出信息矩阵 ψ（t）的第 i行。

对于如式（9）所示的待辨识模型，通常可采用最小二乘法或者随机梯度法进行系统辨识，但值得注意的是，辨识模型中的信息矩阵Φ（t）包含了大量的零元。当系统输入输出点较多时，将造成辨识计算量非常大，难以满足在线辨识的要求。

1.2 最小二乘法

由于ϑ中包含了所有待辨识参数，因此通过下述递推算法则可得到第t时刻的系统辨识结果：

其中，L（t）和 P（t）分别为增益矩阵和协方差矩阵。而对于多输入多输出系统，需将上述辨识过程以行为单位分解成为若干个子系统再分别进行辨识，其过程如图1所示，图中 ym（t）、Φm（t）、Pm（t）分别为输出变量、信息矩阵和协方差矩阵第t时刻迭代的第m 行，和分别为第 t时刻和第t-1时刻迭代结果的第m行。

最小二乘法的基本原理是通过极小化如下最小二乘准则函数完成辨识过程：

图1 多输入多输出系统递推最小二乘法辨识示意图Fig．1 Schematic diagram of MIMO system identification by recursive least square algorithm

由图1可知，每个子系统的每次迭代过程都将计算高维协方差矩阵（总维数为（n+mnr）×（n+mnr），其中n、r、m分别为系统状态变量、输入变量和输出变量的维数）。而对于电力系统这样的大系统，虽然最小二乘法辨识精度较高，但在进行在线辨识时往往因其计算量较大而难以满足实际工程应用的要求。因此，减少辨识过程中的冗余计算是采用递推最小二乘法进行在线辨识的关键。

1.3 耦合最小二乘法

由图1可知，每个子系统的辨识过程是相对独立的，而这忽略了它们之间的联系，造成了大量的冗余计算。为消除这样的冗余计算过程，可在对第i个子系统进行辨识的时候，利用第i-1个子系统的辨识结果，辨识过程如图2所示。

由于待辨识系统中包含了共同的参数向量α，因此随着迭代的不断进行，α的估计值将不断收敛于真实值，即第i-1个子系统在第t时刻的估计值âi-1（t）将比第 i个子系统在第 t-1 时刻的估计值âi（t-1）更接近于真实值。由于待辨识系统中包含了共同的信息向量φ，因此利用第i-1个子系统的协方差矩阵 Pi-1（t）计算第 i个子系统的协方差矩阵 Pi（t）也将大幅减少计算量。耦合最小二乘法的具体计算公式如下所示：

1.4 多新息耦合最小二乘算法

多新息辨识理论是系统辨识领域新出现的重要分支，该理论将传统单新息修正技术拓展到多新息修正技术，通过这样的转变可提高辨识的精度[11]。在辨识过程中，通过不断对估计值进行修正，可逐渐逼近真值，这样的过程可由下式表示：

其中，e（t）为新息。 通过式（20）与式（11）的对照，有：

将上式中的 y（t）和 Φ（t）分别用式（22）和式（23）替代，则单新息 e（t）则转化为多新息形式 E（p，t），其中p为新息长度。

将式（22）和式（23）分别替代式（14）—（19）中的yi（t）和 Φi（t），则可得到多新息耦合最小二乘算法。

至此，上述辨识过程可总结如下。

b.收集输入数据 u（t）和输出数据 y（t），构成φ（t）、Y（p，t）以及 Γi（p，t）（i=1，2，…，m），其中 Γi（p，t）为Γ（p，t）的第 i行。

c.分别由式（12）和式（13）计算第1个子系统的增益矩阵 L1（t）和协方差矩阵 P1（t），并由式（14）刷新估计值。

d.对于第 i个子系统，分别由式（15）和式（16）计算增益矩阵 Li（t）和协方差矩阵 Pi（t），并由式（17）刷新估计值，直至第 m 个子系统。

2 电力系统状态空间辨识相关问题

2.1 激励信号的选取

就系统辨识问题而言，激励信号与系统的可辨识性息息相关，在选取激励信号时，通常期望其能够激发系统所有的振荡模态，并且具有一定的持续性。常用的激励信号有方波信号、白噪声信号或者如下激励信号[12]：

图2 耦合递推最小二乘法辨识示意图Fig.2 Schematic diagram of coupling recursive least square algorithm

其中，0≤t≤2T0，T0=T／2，T 为持续激励信号的周期；P0为激励信号频谱的幅值；[f1，fh]为持续信号频带范围。上述激励信号的优点是在频带范围内的不同频率点都具有相同的激励强度，其波形如图3所示（图中Pd为标幺值）。

图3 激励信号波形Fig.3 Waveform of exciting signal

由于电力系统是一个复杂庞大的系统，并且其电源结构、网架结构逐年发生着新的变化，若在每次辨识时都在不同节点处额外注入如上文所述的特殊激励信号，这无疑增加了辨识过程的复杂性，也不利于实际工程的应用。因此，本文将利用电力系统中存在的因负荷变化、机组投切等随机小扰动信号作为系统的输入信号进行系统辨识。

2.2 数据采样的处理

对于电力系统状态空间在线辨识问题，信号采样周期与采样窗口长度的选取对于辨识速度和精度都会造成一定的影响。本文在进行数据采样的时候，主要根据以下3条原则：

a.根据奈奎斯特采样定理，对于周期为T的激励信号，采样周期Ts应满足Ts＜T／2，但采样周期不宜过小，否则将增加不必要的计算量；

b.应保证采样数据窗长度ts＞T，通过增加数据窗长度ts，可一定程度上提高辨识精度；

c.采样数据点总数N=ts／Ts不宜过大，否则将降低辨识速度。

2.3 输入输出位置的选取

由于本文对电力系统状态空间进行辨识的主要目的是用于WADC的设计，因此在选择系统输入输出信号的时候，可采用相关模式的能控／能观指标作为选取依据。但当WADC采用多输入多输出控制结构时，该方法不宜采用。因为在计算能控／能观指标时需首先得到系统矩阵A、B、C，这需要收集全部待选信号（如全网发电机功角、联络线有功等），得到一个高维多输入多输出系统，再通过辨识方法才能得到系统矩阵，这大幅增加了计算负担，不利于实际应用。为解决此问题，本文通过参与因子法作为选取的依据[13]。

其中，φji和ψji分别为左、右特征向量第j行第i列元素，反映了第j个状态变量与第i个振荡模式之间的相关性，是一个能控和能观的综合性指标。由于矩阵A反映系统固有特性，与控制器输入输出点的选择无关，因此可任意选择一组观测信号，进行单输入单输出的系统辨识，得到系统矩阵A，然后可通过左、右特征向量方便计算相关因子。

2.4 电力系统在线辨识流程与阻尼控制结构

通过上述分析，电力系统在线辨识流程可用图4表示。其中数据预处理环节主要是通过去趋势、零均值化、规格化以及低通滤波等处理过程，将辨识信号中存在的信号偏移、直流分量、干扰噪声等对辨识结果不利的部分去除，以达到提高辨识精度的目的。而对于辨识模型定阶问题，本文将采用贝叶斯准则BIC（Bayesian Information Criterion）进行处理。

图4 电力系统在线辨识流程图Fig.4 Flowchart of online power system identification

通过上述辨识过程，可得到电力系统的状态空间模型，然后根据如图5所示的WADC设计结构，可实现WADC的在线设计。

图5 WADC在线设计示意图Fig.5 Schematic diagram of WADC online design

3 仿真算例

3.1 4机2区测试系统

为验证本文算法的有效性，首先采用4机2区系统进行仿真试验，如图6所示，模型具体参数见文献[14]。

图6 4机2区测试系统Fig.6 4-machine 2-area test system

对上述系统进行特征值分析可知，系统存在3组振荡模式，其中频率为1.14 Hz和1.16 Hz的为区内模式，频率为0.61 Hz的为区间模式。区内模式的阻尼通过本地电力系统稳定器PSS（Power System Stabilizer）能得到有效提升，而对于区间模式的阻尼，可通过WADC进行附加控制加以改善。

为模拟实际系统中因负荷的随机投切产生的小扰动类噪声信号，本文向节点7和节点9注入如图7所示的经低通滤波的高斯白噪声信号（图中幅值为标幺值），其幅值为负荷水平的1%，系统各发电机功角与区域联络线有功时域响应如图8所示（图中有功功率为标幺值）。

图7 随机小扰动信号Fig.7 Stochastic small disturbance signal

图8 各发电机功角以及区域联络线有功的时域响应Fig.8 Time-domain response of generator power angle and areal tie-line active power

根据上文所述原则，选取采样周期Ts为0.05 s，数据窗长度ts为20 s。并将发电机G1励磁电压和区域联络线有功作为系统矩阵A辨识的观测信号。由于该辨识过程目的主要用于参与因子的计算，对于辨识精度要求要低于系统完整状态空间的辨识，因此选用耦合最小二乘法进行初步辨识，并通过贝叶斯准则进行模型定阶，得到一个9阶的初步辨识结果，然后计算各发电机功角的参与因子。

表1为初步辨识结果与实际系统特征值的比较，由于系统主导振荡模式的频率和阻尼比能够反映系统矩阵A的辨识精度，因此该辨识结果表明通过耦合最小二乘算法得到的辨识结果能够基本表征系统矩阵A的主要特征。图9为各发电机功角在模式1下参与因子的辨识结果与实际理论计算结果对比。由图9可知，虽然辨识系统与原系统的参与因子计算存在一定的误差，但各发电机之间的相对关系并没有发生本质变化，即同步发电机G1、G3、G4的参与因子相对于同步发电机G2要大，因此本文选择同步发电机 G3、G4与同步发电机 G1的功角差 Δδ13和 Δδ14作为待辨识系统的反馈信号，将同步发电机G1和G3的发电机励磁电压作为待辨识系统的控制信号。进行离线辨识时，可选取相对较大的新息长度。而若采用本文方法进行在线辨识时，则新息长度的选取应在辨识误差和耗时上进行折中，根据实际系统进行调整试验。值得注意的是：表2中的采样长度也是影响辨识精度与速度的关键因素，其选取不宜过大或过小，应根据2.2节原则b与c进行选取，以保证辨识精度的同时尽量避免耗时较长。

表1 实际模型和辨识模型特征值分析对比Table 1 Comparison of eigenvalues between actual model and identification model

图9 辨识系统与实际系统的参与因子对比Fig.9 Comparison of participation factors between identification system and actual system

为直观说明本文辨识方法的有效性，将辨识模型与实际系统在相同持续激励下的输出响应进行比较。图10和图11分别表示的是在新息长度为1、3和5这3种情况下，同步发电机G3、G4与同步发电机G1的功角差 Δδ13和 Δδ14（均为标幺值）的辨识结果与实际测量结果的对比。图中数据都是从第6秒开始，这是因为辨识算法需经过一段时间才能收敛。由图10和图11可知，当新息长度的增加时，辨识结果与测量结果更为接近。

为进一步比较不同新息长度对辨识结果的影响，本文将测试系统的2个输出通道的误差取平均值进行比较，如表2所示。

上述仿真结果为本文采用的辨识算法在主频为2.80 GHz的CPU、内存为1.85 GB的PC上通过MATLAB 7.10.0版本进行的试验数据。由表2可知，采用传统方法，即新息长度为1时，辨识误差较大。随着新息长度的提升，可降低辨识的相对误差，但当新息长度增加到一定程度时，辨识误差的减少幅度并不显著，而耗时则大幅提升。因此，若采用本文方法

图10 不同新息长度下G3与G1功角差的辨识值与测量值的对比Fig.10 Comparison between identified and measured power angle difference between G1and G3for different innovation lengths

图11 不同新息长度下G4与G1功角差的辨识值与测量值的对比Fig.11 Comparison between identified and measured power angle difference between G1and G4for different innovation lengths

表2 新息长度对辨识误差和计算耗时的影响Table 2 Impact of innovation length on identification error and calculation time

3.2 New England测试系统

为说明本文算法对于不同系统的适用性，选取New England测试系统进行上述仿真试验，测试系统单线图如图12所示，测试系统相关参数见文献[15]。为模拟实际系统中存在的因负荷投切和参数波动产生的随机扰动，向各负荷点注入经低通滤波的高斯白噪声信号，此时测试系统46条传输线有功响应（标幺值）如图13所示。

图12 New England测试系统Fig．12 New England test system

图13 New England测试系统各线路有功响应Fig．13 Active power response of each line of New England test system

选取同步发电机G1的励磁电压和节点9-39的联络线有功作为待辨识系统的观测信号，通过贝叶斯准则得到13阶系统，对其进行特征值分析并与原系统进行对比，如图14所示。

从图14可知，辨识系统能基本反映原系统的振荡模式，由此可进行系统参与因子的计算。考虑到本文辨识目的主要用于WADC的设计，因此选取频率为0.61 Hz和0.84 Hz的2个区域模式作为选择依据，并将发电机功角参与因子较大者列入表3。

根据上述参与因子计算结果，选择同步发电机G5、G7、G10励磁侧电压作为待辨识系统输入信号，选择同步发电机功角差 Δδ5-10、Δδ7-10、Δδ9-10（同步发电机G10为平衡机）作为待辨识系统输出信号。选择新息长度p=5，数据窗长度为2000个采样点，在此情况下采用多新息耦合最小二乘算法进行系统辨识，其结果与实际系统的在相同持续激励下的输出响应如图15所示（纵轴均为标幺值）。

图14 辨识系统与原系统特征值对比Fig．14 Comparison of eigenvalues between identification system and original system

表3 发电机功角参与因子Table 3 Power angle participation factor of synchronous generators

图15 辨识模型与实际模型的输出响应对比Fig．15 Comparison of output response between identification model and actual model

由图15可知，系统辨识模型与实际模型输出响应基本吻合，这说明本文辨识结果精度较高。为进一步量化辨识误差与辨识时间，将结果列入表4，并与多新息最小二乘算法进行比较。经分析可知，本文方法在相同新息长度情况下，误差要略高于多新息最小二乘算法，但由于采用了耦合辨识思想，辨识时间要明显优于多新息最小二乘算法，这更符合WADC在线设计的要求。

表4 新息长度对辨识误差和计算耗时的影响Table 4 Impact of innovation length on identification error and calculation time

4 结论

a.本文从理论上分析了采用传统最小二乘算法用于多输入多输出系统辨识时将造成计算量大的根本原因在于，每次迭代过程都重复计算高维协方差矩阵，这使得最小二乘算法难以应用于WADC的在线设计。

b.本文推导了多新息耦合最小二乘算法的基本计算方法，并提出采用该方法对电力系统状态空间进行辨识，避免了采用传统最小二乘算法计算量大的问题，克服了单新息最小二乘算法误差高的问题，具有一定在线辨识的能力。

c.本文提出的多新息耦合最小二乘算法随着新息长度的增加，辨识精度也会进一步提升。但增加到一定程度后，辨识精度的提升并不明显，而辨识时间却大幅增加。因此在利用该方法进行在线辨识时，应根据实际情况进行调整，避免新息长度选取过大或者过小造成的辨识时间过长或者辨识精度过低的问题。

d.本文提出的多新息耦合最小二乘算法在相同新息长度上的辨识误差上要略高于多新息最小二乘算法，但在耗时上要明显优于后者。因此，多新息最小二乘算法更符合在线辨识的要求。

e.本文提出利用电力系统中因符合投切、参数波动等小幅扰动进行状态空间辨识，避免了传统辨识方法在每次辨识时需额外附加激励才能进行辨识的弊端，能方便在线完成系统状态空间模型的更新，具有较高的实用价值。