CN101369002A - 利用实测轨迹和轨迹灵敏度校正发电机仿真参数的方法 - Google Patents

Abstract

本发明公开了一种利用实测轨迹和轨迹灵敏度校正发电机仿真参数的方法，该方法主要是针对发电机中某些参数难以通过出厂试验获取，且参数辨识难以摆脱电力系统仿真模型和参数的限制，从而难以在线获取发电机仿真参数的现状，通过评价实测轨迹与仿真轨迹的误差大小，获取发电机各参数对仿真轨迹的灵敏度，并对轨迹灵敏度进行排序，从而确定对实测轨迹影响较大的仿真参数集合，进而采用最优化方法对可校正的参数进行优化校正，从而获得使实测和仿真轨迹相吻合的发电机仿真参数。该方法避免了传统的发电机参数辨识方法对模型的依赖性，计算简单，具有在线应用潜力，也可为电力系统中其它元件仿真参数的校正中应用。

Description

利用实测轨迹和轨迹灵敏度校正发电机仿真参数的方法

技术领域：

本发明涉及一种利用实测轨迹和轨迹灵敏度校正发电机仿真参数的方法。

背景技术

电力系统数值仿真是研究电力系统动态行为的重要工具，对电力系统安全分析和运行具有十分重要的作用。精确的仿真模型和参数是数值仿真结果可信的保证。发电机是电力系统的重要设备，其参数的精确与否直接影响到数值仿真的可信度，从而影响了电力系统的安全分析水平。

获得准确的发电机仿真参数是电力系统安全分析的工作重点之一，虽然发电机出厂之前都经过实验确定了大部分的仿真参数，但是还有部分参数不易通过试验获得，且长期运行的发电机参数可能发生一定的改变；基于发电机模型线性化的参数辨识则受制于模型的影响，难以在线应用。迄今未见本发明方法的文献报导和应用。

发明内容

本发明针对本领域技术人员一直渴望解决发电机仿真参数难以在线获取的问题，提出一种利用实测轨迹和轨迹灵敏度的发电机仿真参数校正方法，通过校正发电机仿真参数使得实测轨迹与仿真轨迹相吻合，从而获得精确的发电机仿真参数。

本发明的技术方案是：一种利用实测轨迹和轨迹灵敏度校正发电机仿真参数的方法，其特征是：

1)给出实测轨迹与仿真轨迹的误差评价指标：当电力系统数值仿真的结果与系统实际量测轨迹有误差时，用能量误差指标、第一摆幅值误差、第一摆周期误差、第二摆幅值误差和第二摆周期误差指标评价系统电压相角的仿真与实测轨迹的误差，

能量误差指标(Error Energy)

$\mathrm{E\; E}=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(y_{\mathrm{sim\; u}}\left(i\right)-y_{\mathrm{meas}}\left(i\right))}^2}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(y_{\mathrm{meas}}\left(i\right)-y_{\mathrm{stab}})}^2}---\left(1\right)$

式(1)中y_simu(i)是仿真变量序列，y_meas(i)是实测变量序列，y_stab是实测变量的稳态值，N是实测变量与仿真变量个数，

对第一摆幅值误差(First Swing Magnitude Error)及第二摆幅值误差(Second Swing Magnitude Error)

$\mathrm{FSME}=\frac{{\mathrm{FMag}}_{\mathrm{simu}}-{\mathrm{FMag}}_{\mathrm{measur}}}{{\mathrm{FMag}}_{\mathrm{measur}}}---\left(2\right)$

$\mathrm{SSME}=\frac{{\mathrm{SMag}}_{\mathrm{simu}}-{\mathrm{SMag}}_{\mathrm{measur}}}{{\mathrm{SMag}}_{\mathrm{measur}}}---\left(3\right)$

式(2)、(3)中FMag_simu是仿真值第一摆的摆动幅值，FMag_measur是实测值第一摆的摆动幅值，SMag_simu是仿真第二摆的摆动幅值，SMag_measur是实测第二摆的摆动幅值，

第一摆摆动的周期误差(First Swing Period Error)及第二摆的周期误差(Second Swing Period Error)

$\mathrm{FSPE}=\frac{{\mathrm{FPer}}_{\mathrm{simu}}-{\mathrm{FPer}}_{\mathrm{measur}}}{{\mathrm{FPer}}_{\mathrm{measur}}}---\left(4\right)$

$\mathrm{SSPE}=\frac{{\mathrm{SPer}}_{\mathrm{simu}}-{\mathrm{SPer}}_{\mathrm{measur}}}{{\mathrm{SPer}}_{\mathrm{measur}}}---\left(5\right)$

式(4)(5)中FPer_simu是仿真的第一摆的摆动周期，FPer_measur是实测第一摆的摆动周期，SPer_simu是仿真的第二摆的摆动周期，SPer_measur是实测第二摆的摆动周期；

2)计算发电机各参数对轨迹的灵敏度，并对轨迹灵敏度进行排序：其一用基于数学模型的解析法计算轨迹灵敏度，轨迹关于参数的灵敏度是反映系统中某一参数发生微小变化时动态轨迹的变化程度，电力系统模型可用一组微分——代数方程组表示

$\left\{\begin{array}{c}g\\ x\left(t\right)=f(x\left(t\right),y\left(t\right),\mathrm{\θ})\\ 0=g(x\left(t\right),y\left(t\right),\mathrm{\θ})\end{array}\right.---\left(6\right)$

式(6)中，t为时间，x为系统状态变量，θ为系统元件的模型参数，变量x关于参数θ的灵敏度即为 $x_{\mathrm{\θ}}=\frac{\∂x}{\∂\mathrm{\θ}};$ 代数矢量y对参数θ的灵敏度即为 $y_{\mathrm{\θ}}=\frac{\∂y}{\∂\mathrm{\θ}},$ 用式(6)对参数θ进行求导可得下式

$\left\{\begin{array}{c}{x}_{\mathrm{\θ}}^g=\frac{\∂f}{\∂x}{x}_{\mathrm{\θ}}+\frac{\∂f}{\∂y}{y}_{\mathrm{\θ}}+\frac{\∂f}{\∂\mathrm{\θ}}\\ 0=\frac{\∂g}{\∂x}{x}_{\mathrm{\θ}}+\frac{\∂g}{\∂y}{y}_{\mathrm{\θ}}+\frac{\∂g}{\∂\mathrm{\θ}}\end{array}\right.---\left(7\right)$

从式(7)可以看出轨迹灵敏度x_θ，y_θ就是该方程的解，又因为式(7)依赖式(6)的解轨迹，通常式(6)(7)联立求解，即可得到轨迹关于参数的灵敏度；

其二用摄动法计算轨迹灵敏度，对参数θ作微小摄动Δθ，然后计算代数矢量相应的变化量Δx，Δy，即可利用下面公式近似计算x_θ，y_θ，

$x_{\mathrm{\θ}}=\frac{\∂x}{\∂\mathrm{\θ}}\≈\frac{\mathrm{\Δx}}{\mathrm{\Δ\θ}}---\left(8\right)$

$y_{\mathrm{\θ}}=\frac{\∂y}{\∂\mathrm{\θ}}\≈\frac{\mathrm{\Δy}}{\mathrm{\Δ\θ}}---\left(9\right)$

再根据灵敏度大小对参数排序，系统的实际观测轨迹的输出向量为Y＝[Y₁，Y₂，L，Y_n]^T，仿真轨迹的输出向量y＝[y₁，y₂，L，y_n]^T，系统仿真参数为θ＝(θ₁，θ₂，L，θ_m)^T，发电机仿真参数量纲并不相同，因此，各参数的轨迹灵敏度不具有可比较性，为使其灵敏度具有可比较性，将轨迹y关于参数的灵敏度定义为y_θθ，

轨迹对参数的灵敏度矩阵表示为：

$s(\mathrm{\θ},t)=\left[\begin{array}{ccc}s({\mathrm{\θ}}_1,t_1)& L& s({\mathrm{\θ}}_m,t_1)\\ s({\mathrm{\θ}}_1,t_2)& L& s({\mathrm{\θ}}_m,t_2)\\ M& O& M\\ s({\mathrm{\θ}}_1,t_n)& L& s({\mathrm{\θ}}_m,t_n)\end{array}\right],$ 其中 $s({\mathrm{\θ}}_i,t_j)=\frac{\∂y_j}{\∂{\mathrm{\θ}}_i}\mathrm{\θ}$

根据矩阵S的列向量的无穷大范数，对发电机参数排序，以此确定须校正的参数θ^*＝(θ₁ ^*，θ₂ ^*，L，θ_p ^*)^T，p<m；

3)选择轨迹灵敏度较大的参数组成待校正参数集合：参数独立性分析，根据高斯——牛顿法(Gauss——Newton Method)，将y(θ^*)在选定初值θ₀ ^*点泰勒展开，并忽略二阶以上的高阶项，

即： $y\left({\mathrm{\&theta;}}^{*}\right)=y\left({{\mathrm{\&theta;}}_0}^{*}\right)+{\left(\frac{\&PartialD;y}{\&PartialD;{\mathrm{\&theta;}}^{*T}}\right)}_{{{\mathrm{\&theta;}}_0}^{*}}({\mathrm{\&theta;}}^{*}-{{\mathrm{\&theta;}}_0}^{*})---\left(10\right)$

令 $\mathrm{\&Delta;}{\mathrm{\&theta;}}^{*}=\frac{{\mathrm{\&theta;}}^{*}-{{\mathrm{\&theta;}}_0}^{*}}{{{\mathrm{\&theta;}}_0}^{*}},$ $s={\left(\frac{\&PartialD;y}{\&PartialD;{\mathrm{\&theta;}}^{*T}}\right)}_{{{\mathrm{\&theta;}}_0}^{*}}{{\mathrm{\&theta;}}_0}^{*}$

则可得:y(θ^*)＝y(θ₀ ^*)+sΔθ^*     (11)

令Δy＝y(θ^*)-y(θ₀ ^*)，

则可得：Δy＝sgΔθ^*             (12)

只有当矩阵s^Tgs是满秩矩阵时，上式中Δθ^*可以唯一求解，这说明当矩阵s^Tgs是满秩矩阵时，可以通过轨迹误差校正参数θ^*＝(θ₁ ^*，θ₂ ^*，L，θ_p ^*)^T，当矩阵s^Tgs不是满秩矩阵时，矩阵s^Tgs中存在线性相关的列向量，这时需要把非线性相关的列向量找出，然后对这些非线性相关列向量所对应的参数进行校正；

找出矩阵s^Tgs中非线性相关列向量所对应参数：首先将矩阵s^Tgs进行特征值分解s^Tgs＝VΛV^-1，用式(13)确定矩阵s^Tgs的秩，由此确定该矩阵中存在多少个非线性相关的列向量，即有多少个参数校正，

$\mathrm{rank}(s^T\mathrm{gs},\mathrm{\&epsiv;})=\max \left\{i\right|\frac{\left|{\mathrm{\&sigma;}}_i\right|}{\left|{\mathrm{\&sigma;}}_1\right|}>\mathrm{\&epsiv;}\left|\right|s^T\mathrm{gs}\left|\right|)m---\left(13\right)$

式(13)中σ_i表示矩阵s^Tgs的特征值，m表示该矩阵的阶数，ε是一个估计误差精度的值，

然后，对矩阵V进行列选主元得到转换矩阵H，从而对矩阵s^Tgs中的列向量，按非相关性由强到弱排列，取其中前rank(s^Tgs，ε)个列向量所对应的参数即为校正的参数；

4)基于实测轨迹和仿真轨迹误差最小为目标，采用优化方法，对待校正的发电机仿真参数集合进行校正，直至精度满足要求：基于最小二乘法的参数校正，

电力系统的实际观测输出向量为Y＝[Y₁，Y₂，L，Y_n]，仿真输出向量y＝[y₁，y₂，L，y_n]，初始仿真参数为θ₀ ^*，则根据最小二乘法的原理可得目标函数为：

J(θ^*)＝(Y-y)^T(Y-y)         (14)

将式(11)代入上式可得

J(Δθ^*)＝[Y-y(θ₀ ^*)-sΔθ^*]^T[(Y-y(θ₀ ^*)-sΔθ^*)]      (15)

对J求极小值，有

${\left(\frac{\&PartialD;J\left(\mathrm{\&Delta;}{\mathrm{\&theta;}}^{*}\right)}{\&PartialD;{\mathrm{\&theta;}}^{*T}}\right)}_{{{\mathrm{\&theta;}}_0}^{*}}=0$

解得：

$\mathrm{\&Delta;}{\mathrm{\&theta;}}^{*}={\left(s^{T}s\right)}^{-1}s(Y-y\left({\mathrm{\&theta;}}_0^{*}\right))---\left(16\right)$

参数的估计值为：

${\mathrm{\&theta;}}^{*}={\mathrm{\&theta;}}_0^{*}(1-\mathrm{\&Delta;}{\mathrm{\&theta;}}^{*})$

由于忽略了目标函数的高阶项，所以一般需要进行迭代求解，即

${\mathrm{\&theta;}}_{j+1}^{*}={\mathrm{\&theta;}}_j^{*}(1-\mathrm{\&Delta;}{\mathrm{\&theta;}}_j^{*})---\left(17\right)$

将上式迭代求解直到Δθ^*和

满足精度要求为止。

本发明利用实测轨迹和轨迹灵敏度校正发电机仿真参数的方法与传统的通过出厂实验获取发电机参数不同，本发明基于发电机联网运行中的受扰轨迹进行发电机仿真参数的校正，从而更能反映发电机运行时的参数特性；与传统的参数辨识相比，参数辨识方法必须基于发电机模型进行线性化，在给出输入输出的情况下辨识出发电机的仿真参数，因此，参数辨识对发电机模型依赖性较强，且采用线性系统的分析手段解决非线性系统的问题精度难以保证，而本发明避免了对系统模型的依赖性以及线性系统分析方法的局限性，仅基于轨迹和轨迹灵敏度，具有较高的应用价值。具有投资小，回收年限短，在利用可再生风力资源同时又起到配电网节能的作用。

附图说明

图1为发电机单机无穷大系统接线示意图。

图2为实测与原始仿真轨迹对比示意图。

图3能量误差随迭代变化情况示意图。

图4实测与原始仿真轨迹对比示意图。

图5实测与校正后仿真轨迹对比示意图。

具体实施方式

1 轨迹误差评价

当电力系统数值仿真的结果与系统实际量测轨迹有误差时，用下面5个指标评价系统电压相角的仿真与实测轨迹的误差。电压相角误差指标包括能量误差指标、第一摆幅值误差、第一摆周期误差、第二摆幅值误差和第二摆周期误差。

1.1 对相角全部变化过程进行比较和分析，定义为误差能量指标(ErrorEnergy)

$\mathrm{EE}=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_{\mathrm{simu}}\left(i\right)-y_{\mathrm{meas}}\left(i\right))}^2}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_{\mathrm{meas}}\left(i\right)-y_{\mathrm{stab}})}^2}---\left(1\right)$

式(1)中y_simu(i)是仿真变量序列，y_meas(i)是实测变量序列，y_stab是实测变量的稳态值，N是实测变量与仿真变量个数。

1.2 第一摆幅值误差(First Swing Magnitude Error)及第二摆幅值误差(Second Swing Magnitude Error)

$\mathrm{FSME}=\frac{{\mathrm{FMag}}_{\mathrm{simu}}-{\mathrm{FMag}}_{\mathrm{measur}}}{{\mathrm{FMag}}_{\mathrm{measur}}}---\left(2\right)$

$\mathrm{SSME}=\frac{{\mathrm{SMag}}_{\mathrm{simu}}-{\mathrm{SMag}}_{\mathrm{measur}}}{{\mathrm{SMag}}_{\mathrm{measur}}}---\left(3\right)$

式(2)(3)中FMag_simu是仿真值第一摆的摆动幅值，FMag_measur是实测值第一摆的摆动幅值，SMag_simu是仿真第二摆的摆动幅值，SMag_measur是实测第二摆的摆动幅值。

第一摆摆动的周期误差(First Swing Period Error)及第二摆的周期误差(Second Swing Period Error)

$\mathrm{FSPE}=\frac{{\mathrm{FPer}}_{\mathrm{simu}}-{\mathrm{FPer}}_{\mathrm{measur}}}{{\mathrm{FPer}}_{\mathrm{measur}}}---\left(4\right)$

$\mathrm{SSPE}=\frac{{\mathrm{SPer}}_{\mathrm{simu}}-{\mathrm{SPer}}_{\mathrm{measur}}}{{\mathrm{SPer}}_{\mathrm{measur}}}---\left(5\right)$

式(4)(5)中FPer_simu是仿真的第一摆的摆动周期，FPer_measur是实测第一摆的摆动周期，SPer_simu是仿真的第二摆的摆动周期，SPer_measur是实测第二摆的摆动周期。

2 轨迹关于参数的灵敏度计算

2.1 基于数学模型的解析法计算轨迹灵敏度，轨迹关于参数的灵敏度是反映系统中某一参数发生微小变化时动态轨迹的变化程度。电力系统模型可用一组微分——代数方程组表示

$\left\{\begin{array}{c}g\\ x\left(t\right)=f(x\left(t\right),y\left(t\right),\mathrm{\&theta;})\\ 0=g(x\left(t\right),y\left(t\right),\mathrm{\&theta;})\end{array}\right.---\left(6\right)$

式(6)中，t为时间，x为系统状态变量，如发电机的转角、转速和磁链等；y为系统的代数变量，如负荷节点电压的幅值和相角；θ为系统中各元件的模型参数，如发电机的同步电抗、阻尼系数、转动惯性时间常数等。变量x关于参数θ的灵敏度即为 $x_{\mathrm{\&theta;}}=\frac{\&PartialD;x}{\&PartialD;\mathrm{\&theta;}};$ 代数矢量y对参数θ的灵敏度即为 $y_{\mathrm{\&theta;}}=\frac{\&PartialD;y}{\&PartialD;\mathrm{\&theta;}}.$ 用式(6)对参数θ进行求导可得下式

$\left\{\begin{array}{c}{x}_{\mathrm{\&theta;}}^g=\frac{\&PartialD;f}{\&PartialD;x}{x}_{\mathrm{\&theta;}}+\frac{\&PartialD;f}{\&PartialD;y}{y}_{\mathrm{\&theta;}}+\frac{\&PartialD;f}{\&PartialD;\mathrm{\&theta;}}\\ 0=\frac{\&PartialD;g}{\&PartialD;x}{x}_{\mathrm{\&theta;}}+\frac{\&PartialD;g}{\&PartialD;y}{y}_{\mathrm{\&theta;}}+\frac{\&PartialD;g}{\&PartialD;\mathrm{\&theta;}}\end{array}\right.---\left(7\right)$

从式(7)可以看出轨迹灵敏度x_θ，y_θ就是该方程的解，又因为式(7)依赖式(6)的解轨迹。通常式(6)(7)联立求解，即可得到轨迹关于参数的灵敏度。

2.2 摄动法计算轨迹灵敏度，在复杂电力系统中式(7)计算困难，解析法很难应用。因此采用摄动法计算轨迹灵敏度，该方法无需对系统进行线性化处理，不涉及系统的物理本质和结构特点。即对参数θ作微小摄动Δθ，然后计算代数矢量相应的变化量Δx，Δy，即可利用下面公式近似计算x_θ，y_θ。

$x_{\mathrm{\&theta;}}=\frac{\&PartialD;x}{\&PartialD;\mathrm{\&theta;}}\&ap;\frac{\mathrm{\&Delta;x}}{\mathrm{\&Delta;\&theta;}}---\left(8\right)$

$y_{\mathrm{\&theta;}}=\frac{\&PartialD;y}{\&PartialD;\mathrm{\&theta;}}\&ap;\frac{\mathrm{\&Delta;y}}{\mathrm{\&Delta;\&theta;}}---\left(9\right)$

根据灵敏度大小对参数排序，假设系统的实际观测轨迹的输出向量为Y＝[Y₁，Y₂，L，Y_n]^T，仿真轨迹的输出向量y＝[y₁，y₂，L，y_n]^T，系统仿真参数为θ＝(θ₁，θ₂，L，θ_m)^T。发电机仿真参数量纲并不相同，因此各参数的轨迹灵敏度不具有可比较性。为使其灵敏度具有可比较性，将轨迹y关于参数的灵敏度定义为y_θθ。

轨迹对参数的灵敏度矩阵表示为：

$s(\mathrm{\&theta;},t)=\left[\begin{array}{ccc}s({\mathrm{\&theta;}}_1,t_1)& L& s({\mathrm{\&theta;}}_m,t_1)\\ s({\mathrm{\&theta;}}_1,t_2)& L& s({\mathrm{\&theta;}}_m,t_2)\\ M& O& M\\ s({\mathrm{\&theta;}}_1,t_n)& L& s({\mathrm{\&theta;}}_m,t_n)\end{array}\right],$ 其中 $s({\mathrm{\&theta;}}_i,t_j)=\frac{\&PartialD;y_j}{\&PartialD;{\mathrm{\&theta;}}_i}\mathrm{\&theta;}$

根据矩阵S的列向量的无穷大范数，对发电机参数排序，以此确定须校正的参数θ^*＝(θ₁ ^*，θ₂ ^*，L，θ_p ^*)^T，p<m。

2.3 待校正参数集的确定方法

在电力系统中，需要不断通过系统实际受扰轨迹对发电机仿真参数进行校正，来获得更加准确的发电机仿真参数。在电力系统数值仿真中，灵敏度小的参数对可观测轨迹的影响小，所以在基于量测轨迹对仿真模型参数校正时忽略灵敏度小的参数所存在的误差。

3 参数独立性分析：

根据高斯——牛顿法(Gauss——Newton Method)，将y(θ^*)在选定初值θ₀ ^*点泰勒展开，并忽略二阶以上的高阶项，

即：

$y\left({\mathrm{\&theta;}}^{*}\right)=y\left({{\mathrm{\&theta;}}_0}^{*}\right)+{\left(\frac{\&PartialD;y}{\&PartialD;{\mathrm{\&theta;}}^{*T}}\right)}_{{{\mathrm{\&theta;}}_0}^{*}}({\mathrm{\&theta;}}^{*}-{{\mathrm{\&theta;}}_0}^{*})---\left(10\right)$

令 $\mathrm{\&Delta;}{\mathrm{\&theta;}}^{*}=\frac{{\mathrm{\&theta;}}^{*}-{{\mathrm{\&theta;}}_0}^{*}}{{{\mathrm{\&theta;}}_0}^{*}},$ $s={\left(\frac{\&PartialD;y}{\&PartialD;{\mathrm{\&theta;}}^{*T}}\right)}_{{{\mathrm{\&theta;}}_0}^{*}}{{\mathrm{\&theta;}}_0}^{*}$

则可得：

y(θ^*)＝y(θ₀ ^*)+sΔθ^*        (11)

令Δy＝y(θ^*)-y(θ₀ ^*)，，

则可得：

Δy＝sgΔθ^*           (12)

只有当矩阵s^Tgs是满秩矩阵时，上式中Δθ^*可以唯一求解。这说明当矩阵s^Tgs是满秩矩阵时，可以通过轨迹误差校正参数θ^*＝(θ₁ ^*，θ₂ ^*，L，θ_p ^*)^T。

当矩阵s^Tgs不是满秩矩阵时，矩阵s^Tgs中存在线性相关的列向量，这时需要把非线性相关的列向量找出，然后对这些非线性相关列向量所对应的参数进行校正。

用下面方法可以找出矩阵s^Tgs中非线性相关列向量所对应参数：

首先将矩阵s^Tgs进行特征值分解s^Tgs＝VΛV^-1，用式(13)确定矩阵s^Tgs的秩，由此确定该矩阵中存在多少个非线性相关的列向量，即有多少个参数可以校正。

$\mathrm{rank}(s^T\mathrm{gs},\mathrm{\&epsiv;})=\max \{i|\frac{\left|{\mathrm{\&sigma;}}_i\right|}{\left|{\mathrm{\&sigma;}}_1\right|}>\mathrm{\&epsiv;}\left|\right|s^T\mathrm{gs}\left|\right|m\}---\left(13\right)$

式(13)中σ_i表示矩阵s^Tgs的特征值，m表示该矩阵的阶数，ε是一个估计误差精度的值。

然后，对矩阵V进行列选主元得到转换矩阵H。从而可以对矩阵s^Tgs中的列向量，按非相关性由强到弱排列，取其中前rank(s^Tgs，ε)个列向量所对应的参数即为可以校正的参数。

4 基于最小二乘法的参数校正

电力系统的实际观测输出向量为Y＝[Y₁，Y₂，L，Y_n]，仿真输出向量y＝[y₁，y₂，L，y_n]，初始仿真参数为θ₀ ^*，则根据最小二乘法的原理可得目标函数为：

J(θ^*)＝(Y-y)^T(Y-y)            (14)

将式(11)代入上式可得

J(Δθ^*)＝[Y-y(θ₀ ^*)-sΔθ^*]^T[(Y-y(θ₀ ^*)-sΔθ^*)]      (15)

对J求极小值，有

${\left(\frac{\&PartialD;J\left(\mathrm{\&Delta;}{\mathrm{\&theta;}}^{*}\right)}{\&PartialD;{\mathrm{\&theta;}}^{*T}}\right)}_{{{\mathrm{\&theta;}}_0}^{*}}=0$

可以解得：

$\mathrm{\&Delta;}{\mathrm{\&theta;}}^{*}={\left(s^{T}s\right)}^{-1}s(Y-y\left({\mathrm{\&theta;}}_0^{*}\right))---\left(16\right)$

参数的估计值为：

${\mathrm{\&theta;}}^{*}={\mathrm{\&theta;}}_0^{*}(1-\mathrm{\&Delta;}{\mathrm{\&theta;}}^{*})$

由于忽略了目标函数的高阶项，所以一般需要进行迭代求解，即

${\mathrm{\&theta;}}_{j+1}^{*}={\mathrm{\&theta;}}_j^{*}(1-\mathrm{\&Delta;}{\mathrm{\&theta;}}_j^{*})---\left(17\right)$

将上式迭代求解直到Δθ^*和

满足精度要求为止。

参照图1，在发电机单机无穷大系统中，在3号母线上发生600ms三相短路故障，通过PMU记录下1号母线的电压相角轨迹。使用电力系统分析综合程序(PSASP)对该故障进行仿真，发电机采用6阶模型，其参数如表1所示。基于本发明所提出方法，根据1母线的电压相角量测轨迹对发电机仿真参数进行校正。

表1 发电机仿真参数

Xd、Xq：直轴、交轴同步电抗

Xd’、Xq’：直轴、交轴暂态电抗

Xd″、Xq”：直轴、交轴次暂态电抗

T’d0：定子开路直轴暂态时间常数

T″d0：定子开路直轴次暂态时间常数

T’q0：定子开路交轴暂态时间常数

T″q0：定子开路交轴次暂态时间常数

Tj：转动惯性时间常数

D：阻尼系数

参照图2，可知母线1的电压相角仿真轨迹与实测轨迹存在很大偏差，说明发电机仿真参数与实际参数存在误差。因此，需要对发电机参数进行校正。

先进行参数灵敏度大小分析，把轨迹对参数的灵敏度矩阵记为S。各参数按其对应在矩阵S中列向量的无穷大范数由大到小排序：Tj、Xq、Xd”、Xq”、D、Td0’、Td0”、Xd’、Xq’、Xd、Tq0’、Tq0”。

因为排序在阻尼参数D之后的参数的灵敏度比较小，所以忽略这些参数存在的误差。因此，只对Tj、Xq、Xd”、Xq”、D这5个参数进行校正，将其记为θ^*＝[Tj，Xq，Xd＂，Xq＂，D]^T。

参数可校正性分析，对选取参数θ^*进行可校正性分析，将参数θ^*所对应的矩阵S^Tgs进行特征值分解S^TS＝VΛV^-1，其中

$\mathrm{\&Lambda;}=\left[\begin{array}{ccccc}0.02& 0& 0& 0& 0\\ 0& 0.166& 0& 0& 0\\ 0& 0& 0.4139& 0& 0\\ 0& 0& 0& 0.5879& 0\\ 0& 0& 0& 0& 1.9863\end{array}\right]$  

根据式(9)判断矩阵s^Tgs的秩为5，然后对矩阵V从右向左列选主元，可得到矩阵H

$H=\left[\begin{array}{ccccc}0& 0& 0& 0& 1\\ 0& 0& 0& 1& 0\\ 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 1& 0& 0& 0\end{array}\right]$

由H^Tθ^*可得到参数θ^*的非相关性排序，由强到弱为Tj、Xq、Xd”、Xq”、D。根据以上分析，选择Tj、Xq、Xd”、Xq”、D进行校正。

参照图3，参数的校正，用最小二乘法对参数进行校正，其误差收敛情况和参数的迭代过程见2表。

表2 迭代过程中各参数值收敛过程

 

迭代次数	能量误差	Tj	Xq	Xd”	Xq”	D
							0	13.65％	14	1.92	0.04	0.0602	0.06
1	4.06％	7.9299	1.2376	0.0748	0.0432	0.0461
							2	3.87％	8.0243	1.5691	0.0646	0.046	0.031
3	0.6％	8.7258	1.2502	0.0639	0.0469	0.0405
							4	0.39％	8.7895	1.2112	0.0641	0.0476	0.0411
5	1.23％	9.1586	1.2031	0.0668	0.0454	0.0459
							6	0.64％	8.6656	1.1699	0.0666	0.0491	0.0421
7	0.77％	8.5828	1.1934	0.0679	0.0507	0.0399
							8	0.49％	8.6296	1.1251	0.0682	0.052	0.0433

参照图4和5，校正结果分析，母线1电压相角实测与仿真误差指标见表3，发电机仿真参数校正前后比较见表4。

表3 母线1电压相角实测与仿真误差指标

表4 发电机仿真参数校正前后比较

Claims (1)

1.一种利用实测轨迹和轨迹灵敏度校正发电机仿真参数的方法，其特征是：

1)给出实测轨迹与仿真轨迹的误差评价指标：当电力系统数值仿真的结果与系统实际量测轨迹有误差时，用能量误差指标、第一摆幅值误差、第一摆周期误差、第二摆幅值误差和第二摆周期误差指标评价系统电压相角的仿真与实测轨迹的误差，

能量误差指标(Error Energy)

$\mathrm{EE}=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_{\mathrm{simu}}\left(i\right)-y_{\mathrm{meas}}\left(i\right))}^2}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_{\mathrm{meas}}\left(i\right)-y_{\mathrm{stab}})}^2}---\left(1\right)$

式(1)中y_simu(i)是仿真变量序列，y_meas(i)是实测变量序列，y_stab是实测变量的稳态值，N是实测变量与仿真变量个数，

对第一摆幅值误差(First Swing Magnitude Error)及第二摆幅值误差(Second Swing Magnitude Error)

$\mathrm{FSME}=\frac{{\mathrm{FMag}}_{\mathrm{simu}}-{\mathrm{FMag}}_{\mathrm{measur}}}{{\mathrm{FMag}}_{\mathrm{measur}}}---\left(2\right)$

$\mathrm{SSME}=\frac{{\mathrm{SMag}}_{\mathrm{simu}}-{\mathrm{SMag}}_{\mathrm{measur}}}{{\mathrm{SMag}}_{\mathrm{measur}}}---\left(3\right)$

式(2)、(3)中FMag_simu是仿真值第一摆的摆动幅值，FMag_measur是实测值第一摆的摆动幅值，SMag_simu是仿真第二摆的摆动幅值，SMag_measur是实测第二摆的摆动幅值，

第一摆摆动的周期误差(First Swing Period Error)及第二摆的周期误差(Second Swing Period Error)

$\mathrm{FSPE}=\frac{{\mathrm{FPer}}_{\mathrm{simu}}-{\mathrm{FPer}}_{\mathrm{measur}}}{{\mathrm{FPer}}_{\mathrm{measur}}}---\left(4\right)$

$\mathrm{SSPE}=\frac{{\mathrm{SPer}}_{\mathrm{simu}}-{\mathrm{SPer}}_{\mathrm{measur}}}{{\mathrm{SPer}}_{\mathrm{measur}}}---\left(5\right)$

式(4)(5)中FPer_simu是仿真的第一摆的摆动周期，FPer_measur是实测第一摆的摆动周期，SPer_simu是仿真的第二摆的摆动周期，SPer_measur是实测第二摆的摆动周期；

2)计算发电机各参数对轨迹的灵敏度，并对轨迹灵敏度进行排序：其一用基于数学模型的解析法计算轨迹灵敏度，轨迹关于参数的灵敏度是反映系统中某一参数发生微小变化时动态轨迹的变化程度，电力系统模型可用一组微分——代数方程组表示

$\left\{\begin{array}{c}g\\ x\left(t\right)=f(x\left(t\right),y\left(t\right),\mathrm{\&theta;})\\ 0=g(x\left(t\right),y\left(t\right),\mathrm{\&theta;})\end{array}\right.---\left(6\right)$

式(6)中，t为时间，x为系统状态变量，θ为系统元件的模型参数，变量x关于参数θ的灵敏度即为 $x_{\mathrm{\&theta;}}=\frac{\&PartialD;x}{\&PartialD;\mathrm{\&theta;}};$ 代数矢量y对参数θ的灵敏度即为 $y_{\mathrm{\&theta;}}=\frac{\&PartialD;y}{\&PartialD;\mathrm{\&theta;}},$ 用式(6)对参数θ进行求导可得下式

$\left\{\begin{array}{c}{x}_{\mathrm{\&theta;}}^g=\frac{\&PartialD;f}{\&PartialD;x}{x}_{\mathrm{\&theta;}}+\frac{\&PartialD;f}{\&PartialD;y}{y}_{\mathrm{\&theta;}}+\frac{\&PartialD;f}{\&PartialD;\mathrm{\&theta;}}\\ 0=\frac{\&PartialD;g}{\&PartialD;x}{x}_{\mathrm{\&theta;}}+\frac{\&PartialD;g}{\&PartialD;y}{y}_{\mathrm{\&theta;}}+\frac{\&PartialD;g}{\&PartialD;\mathrm{\&theta;}}\end{array}\right.---\left(7\right)$

从式(7)可以看出轨迹灵敏度x_θ，y_θ就是该方程的解，又因为式(7)依赖式(6)的解轨迹，通常式(6)(7)联立求解，即可得到轨迹关于参数的灵敏度；

其二用摄动法计算轨迹灵敏度，对参数θ作微小摄动Δθ，然后计算代数矢量相应的变化量Δx，Δy，即可利用下面公式近似计算x_θ，y_θ，

$x_{\mathrm{\&theta;}}=\frac{\&PartialD;x}{\&PartialD;\mathrm{\&theta;}}\&ap;\frac{\mathrm{\&Delta;x}}{\mathrm{\&Delta;\&theta;}}---\left(8\right)$

$y_{\mathrm{\&theta;}}=\frac{\&PartialD;y}{\&PartialD;\mathrm{\&theta;}}\&ap;\frac{\mathrm{\&Delta;y}}{\mathrm{\&Delta;\&theta;}}---\left(9\right)$

再根据灵敏度大小对参数排序，系统的实际观测轨迹的输出向量为Y＝[Y₁，Y₂，L，Y_n]^T，仿真轨迹的输出向量y＝[y₁，y₂，L，y_n]^T，系统仿真参数为θ＝(θ₁，θ₂，L，θ_m)^T，发电机仿真参数量纲并不相同，因此，各参数的轨迹灵敏度不具有可比较性，为使其灵敏度具有可比较性，将轨迹y关于参数的灵敏度定义为y_θθ，

轨迹对参数的灵敏度矩阵表示为：

$s(\mathrm{\&theta;},t)=\left[\begin{array}{ccc}s({\mathrm{\&theta;}}_1,t_1)& L& s({\mathrm{\&theta;}}_m,t_1)\\ s({\mathrm{\&theta;}}_1,t_2)& L& s({\mathrm{\&theta;}}_m,t_2)\\ M& O& M\\ s({\mathrm{\&theta;}}_1,t_n)& L& s({\mathrm{\&theta;}}_m,t_n)\end{array}\right],$ 其中 $s({\mathrm{\&theta;}}_i,t_j)=\frac{\&PartialD;y_j}{\&PartialD;{\mathrm{\&theta;}}_i}\mathrm{\&theta;}$

根据矩阵S的列向量的无穷大范数，对发电机参数排序，以此确定须校正的参数θ^*＝(θ₁ ^*，θ₂ ^*，L，θ_p ^*)^T，p<m；

3)选择轨迹灵敏度较大的参数组成待校正参数集合：参数独立性分析，根据高斯——牛顿法(Gauss——Newton Method)，将y(θ^*)在选定初值θ₀ ^*点泰勒展开，并忽略二阶以上的高阶项，

即： $y\left({\mathrm{\&theta;}}^{*}\right)=y\left({{\mathrm{\&theta;}}_0}^{*}\right)+{\left(\frac{\&PartialD;y}{\&PartialD;{\mathrm{\&theta;}}^{*T}}\right)}_{{{\mathrm{\&theta;}}_0}^{*}}({\mathrm{\&theta;}}^{*}-{{\mathrm{\&theta;}}_0}^{*})---\left(10\right)$

令 $\mathrm{\&Delta;}{\mathrm{\&theta;}}^{*}=\frac{{\mathrm{\&theta;}}^{*}-{{\mathrm{\&theta;}}_0}^{*}}{{{\mathrm{\&theta;}}_0}^{*}},$ $s={\left(\frac{\&PartialD;y}{\&PartialD;{\mathrm{\&theta;}}^{*T}}\right)}_{{{\mathrm{\&theta;}}_0}^{*}}{{\mathrm{\&theta;}}_0}^{*}$

则可得：y(θ^*)＝y(θ₀ ^*)+sΔθ^*    (11)

令Δy＝y(θ^*)-y(θ₀ ^*)，

则可得：Δy＝sgΔθ^*            (12)

只有当矩阵s^Tgs是满秩矩阵时，上式中Δθ^*可以唯一求解，这说明当矩阵s^Tgs是满秩矩阵时，可以通过轨迹误差校正参数θ^*＝(θ₁ ^*，θ₂ ^*，L，θ_p ^*)^T，当矩阵s^Tgs不是满秩矩阵时，矩阵s^Tgs中存在线性相关的列向量，这时需要把非线性相关的列向量找出，然后对这些非线性相关列向量所对应的参数进行校正；

找出矩阵s^Tgs中非线性相关列向量所对应参数：首先将矩阵s^Tgs进行特征值分解s^Tgs＝VΛV^-1，用式(13)确定矩阵s^Tgs的秩，由此确定该矩阵中存在多少个非线性相关的列向量，即有多少个参数可校正，

$\mathrm{rank}(s^T\mathrm{gs},\mathrm{\&epsiv;})=\max \left\{i|\frac{\left|{\mathrm{\&sigma;}}_i\right|}{\left|{\mathrm{\&sigma;}}_1\right|}>\mathrm{\&epsiv;}\left|\right|s^T\mathrm{gs}\left|\right|m\right\}---\left(13\right)$

式(13)中σ_i表示矩阵s^Tgs的特征值，m表示该矩阵的阶数，ε是一个估计误差精度的值，

然后，对矩阵V进行列选主元得到转换矩阵H，从而对矩阵s^Tgs中的列向量，按非相关性由强到弱排列，取其中前rank(s^Tgs，ε)个列向量所对应的参数即为可校正的参数；

4)基于实测轨迹和仿真轨迹误差最小为目标，采用优化方法，对待校正的发电机仿真参数集合进行校正，直至精度满足要求：基于最小二乘法的参数校正，

电力系统的实际观测输出向量为Y＝[Y₁，Y₂，L，Y_n]，仿真输出向量y＝[y₁，y₂，L，y_n]，初始仿真参数为

，则根据最小二乘法的原理可得目标函数为：

J(θ^*)＝(Y-y)^T(Y-y)                                           (14)

将式(11)代入上式可得

J(Δθ^*)＝[Y-y(θ₀ ^*)-sΔθ^*]^T[(Y-y(θ₀ ^*)-sΔθ^*)]  (15)

对J求极小值，有

${\left(\frac{\&PartialD;J\left(\mathrm{\&Delta;}{\mathrm{\&theta;}}^{*}\right)}{\&PartialD;{\mathrm{\&theta;}}^{*T}}\right)}_{{{\mathrm{\&theta;}}_0}^{*}}=0$

解得：

$\mathrm{\&Delta;}{\mathrm{\&theta;}}^{*}={\left(s^{T}s\right)}^{-1}s(Y-y\left({\mathrm{\&theta;}}_0^{*}\right))---\left(16\right)$

参数的估计值为：

${\mathrm{\&theta;}}^{*}={\mathrm{\&theta;}}_0^{*}(1-\mathrm{\&Delta;}{\mathrm{\&theta;}}^{*})$

由于忽略了目标函数的高阶项，所以一般需要进行迭代求解，即

${\mathrm{\&theta;}}_{j+1}^{*}={\mathrm{\&theta;}}_j^{*}(1-\mathrm{\&Delta;}{\mathrm{\&theta;}}_j^{*})---\left(17\right)$

将上式迭代求解直到Δθ^*和

满足精度要求为止。

Images