CN103618492A - 一种基于时频变换的同步发电机参数辨识方法 - Google Patents

Abstract

本发明提供一种基于时频变换的同步发电机参数辨识方法，包括：S1，将发电机模型中的时域内发电机d、q轴微分方程进行拉氏变换得到频域内发电机d、q轴方程，再经过拉氏反变换得到时域内发电机d、q轴积分方程；S2，对时域内发电机d、q轴积分方程里的状态变量进行替换；S3，将时域内发电机d、q轴积分方程离散为矩阵与积分方程系数列向量相乘等于可测列向量或者可测列向量积分的形式，求解得到积分方程系数列向量；S4，根据积分方程系数列向量与发电机参数之间的关系求解得到发电机模型参数。将发电机微分方程转换为积分方程，有利于提高数值计算的稳定性，解决辨识结果不稳定的问题。

Description

一种基于时频变换的同步发电机参数辨识方法

技术领域

本发明涉及电力系统分析技术领域，具体涉及一种基于时频变换的同步发电机参数辨识方法。 

背景技术

同步发电机的模型和参数的准确性关系到电网的安全稳定运行，目前，电力系统分析计算仿真所用的发电机参数是厂家所提供的理论值或设计值，或不得已采用模型，由于不能计及涡流、磁滞、饱和等实际运行工况的影响，所得结果常常与实际工况不符，严重影响了系统分析计算的准确度和可信度，因此，发电机参数的辨识工作是智能调度技术的重要组成部分。近年来，随着数学技术的进步和辨识理论的发展，以及相角测量技术(PMU)、全球定位技术(CPS)等新技术在电力系统的广泛应用，建立在这些新技术平台上的发电机参数辨识取得了很大的发展。 

发电机参数辨识方法主要可归结为时域辨识法和频域辨识法两类。 

时域辨识法基于励磁阶跃或负荷小扰动的动态响应来进行辨识，属于在线辨识法，能够自然计及工况的影响，不用附加过多的假设条件。现有算法主要有最小二乘法、卡尔曼滤波法、进化策略法，基因法、小波分析法、神经网络法以及粒子群优化法等。这些辨识算法在基于各种仿真数据辨识时往往可以得到满意的结果，然而，在针对实测数据时，常常会遇到不同试验甚至重复的试验辨识所得参数相差很大的情况。这是由于在线试验扰动较小，次暂态过程的可观测性较差，加上环境噪声和功角测量精度的影响，导致上述算法很难获得准确的辨识结果。 

频域分析法可分为直流衰减法、静态频域法(SSFR)和动态频域法(OLFR)。直流衰减法由于响应持续过程短，影响频率特性的量测精度，目前较少使用。静态频域法需要大功率的变频电源作为信号源，不便于试验的开展且试验结果无法反映饱和效应。动态频域法得到的参数能够反映实际运行工况，更适用于动态稳定的研究。静态频域法和动态频域法的结合是近年来的发展趋势，即由SSFR给出参数的初值再根据一定运行条件下的OLFR法来修正参数的方法可改善算法的稳定性，且具有一定的滤波能力，不过对输人扰动信号的波形、幅值大小及其相关性要求严格。此外，频域响应分析建立在线性系统的基础上，不能反映同步电机参数的非线性特点。 

此外，通过抛载法和数值计算法可以直接获取发电机参数。抛载法优点是易于实现，且计算方法相对简单，但无法计及励磁电压对暂态过程的影响，不适合主流的自并励机组。将 抛载法和加权最小二乘等算法有效的时域辨识法有望解决励磁电压的影响，然而q轴抛载时需要准确的功角信息，运行工况难以调整，以往研究表明试验所测d轴参数较准确，而q轴所得参数误差较大。数值计算法主要包括磁路磁导法和有限元法，前者需要进行大量的等值，精度不高。有限元方法计算量非常庞大大，且对端部和槽楔等处的电阻、电抗比较困难，计算精度也有待进一步提高。 

发明内容

本发明针对现有技术的不足，提出一种基于时频变换的同步发电机参数辨识方法，包括： 

步骤S1，将发电机模型中的时域内发电机d、q轴微分方程进行拉氏变换得到频域内发电机d、q轴方程，再进行拉氏反变换得到时域内发电机d、q轴积分方程； 

步骤S2，根据所述发电机模型中发电机定子d、q轴电压与状态变量的关系式，对所述时域内发电机d、q轴积分方程里的状态变量进行替换； 

步骤S3，假设发电机定子d、q轴电压和电流为由实测采样点构成的列向量，将所述时域内发电机d、q轴积分方程离散为矩阵与积分方程系数列向量相乘等于可测列向量或者可测列向量及其积分的形式，求解得到所述积分方程系数列向量； 

步骤S4，根据所述积分方程系数列向量与发电机参数之间的关系通过求解方程组得到发电机模型参数。 

本发明提供的第一优选实施例中：所述方法基于机组投运前必做的甩负荷试验或其它在线大扰动试验。 

本发明提供的第二优选实施例中：所述步骤S1中在对时域内发电机模型进行拉氏变换时，需要包含所述状态变量的初值； 

所述步骤S1中得到所述频域内发电机d、q轴方程后对所述频域内发电机d、q轴方程经过处理后再进行拉氏反变换，处理的过程包括： 

消去所述频域内发电机d轴或q轴方程的中间状态变量，使频域内发电机d轴或q轴方程中分别只保留一个状态变量，再在所述频域内发电机d轴或q轴方程两侧分别除以本方程最高次幂的拉氏算子。 

本发明提供的第三优选实施例中：所述步骤后S2中利用所述发电机模型中发电机定子d、q轴电压与状态变量的关系式将所述状态变量用可测电气量替换，所述可测电气量包括机端电压或机端电流； 

所述步骤S3中所述时域内发电机d、q轴积分方程离散为矩阵形式为AX=Y； 

其中A为由各可测列向量及其各种积分形式组成的矩阵，X为待辨识的发电机参数构成 的方程系数列向量，其各分量是发电机实用参数的代数组合，Y为可测列向量或某些可测列向量及其积分的代数和。 

本发明提供的第四优选实施例中：所述发电机模型为： 

$\left\{\begin{array}{c}{U}_d={\mathrm{\ωE}}_d^{\′\′}+{\mathrm{\ωX}}_q^{\′\′}{i}_q-R_a{i}_d\\ U_q={\mathrm{\ωE}}_q^{\′\′}-{\mathrm{\ωX}}_d^{\′\′}{i}_d-R_a{i}_q\\ T_{d0}^{\′}{\mathrm{pE}}_q^{\′}=E_f-E_q^{\′}-\frac{X_d-X_d^{\′}}{X_d^{\′}-X_d^{\′\′}}(E_q^{\′}-E_q^{\′\′})\\ T_{d0}^{\′\′}{\mathrm{pE}}_q^{\′\′}=E_q^{\′}-E_q^{\′\′}-(X_d^{\′}-X_d^{\′\′})i_d\\ T_{q0}^{\′}{\mathrm{pE}}_d^{\′}={-E}_d^{\′}-\frac{X_q-X_q^{\′}}{X_q^{\′}-X_q^{\′\′}}(E_d^{\′}-E_d^{\′\′})\\ T_{q0}^{\′\′}{\mathrm{pE}}_d^{\′\′}=E_d^{\′}-E_d^{\′\′}+(X_q^{\′}-X_q^{\′\′})i_q\end{array}\right.---\left(1\right);$

其中，E″_d，E′_d，E″_q，E′_q，δ，ω为六个状态变量，通常称为六阶实用模型；X_d、X_q、X′_d、X′_q、X″_d、X″_q、T′_d0、T′_q0、T″_d0、T″_q0、R_a为发电机参数，X_d为直轴同步电抗；X_q为交轴同步电抗；X′_d为直轴暂态电抗；X′_q交轴暂态电抗；X″_d直轴次暂态电抗；X″_q交轴次暂态电抗；T′_d0、T′_q0、T″_d0、T″_q0、为时间常数；R_a为电枢电阻；U_d、U_q分别为发电机定子d、q轴电压；i_d、i_q分别为发电机定子d、q轴电流；E_f为励磁电压；p为微分算子

所述方程（1）中d轴的状态方程为： 

$\left\{\begin{array}{c}{T}_{d0}^{\′\′}{\mathrm{pE}}_q^{\′\′}=E_q^{\′}-E_q^{\′\′}-(X_d^{\′}-X_d^{\′\′})i_d\\ T_{d0}^{\′}{\mathrm{pE}}_q^{\′}=E_f-E_q^{\′}-\frac{X_d-X_d^{\′}}{X_d^{\′}-X_d^{\′\′}}(E_q^{\′}-E_q^{\′\′})\end{array}\right.---\left(2\right);$

所述方程（1）中q轴的状态方程为： 

$\left\{\begin{array}{c}{T}_{q0}^{\′\′}{\mathrm{pE}}_d^{\′\′}=E_d^{\′}-E_d^{\′\′}+(X_q^{\′}-X_q^{\′\′})i_q\\ T_{q0}^{\′}{\mathrm{pE}}_d^{\′}=-E_d^{\′}-\frac{X_q-X_q^{\′}}{X_q^{\′}-X_q^{\′\′}}(E_d^{\′}-E_d^{\′\′})\end{array}\right.---\left(3\right).$

本发明提供的第五优选实施例中：所述步骤S1中对所述发电机模型中的所述时域内发电机d、q轴微分方程进行拉氏变换得到所述频域内发电机d、q轴方程的方法包括： 

步骤S101，对于所述方程（2），令E″_q＝u,E′_q＝v,E_f＝w,T″_do＝λ,T′_do＝μ，k₁＝X′_d-X″_d， 

i_d＝i；对于所述方程（3），令E″_d＝u,E′_d＝v,T″_qo＝λ,T′_qo＝μ，k₁＝-(X_q′-X_q″)，  $k_2=\frac{X_q-X_q^{\′}}{X_q^{\′}-X_q^{\′\′}},$ i_q＝i； 

步骤S102，将方程（2）和方程（3）改写成： 

$\left\{\begin{array}{c}\mathrm{\λpu}=v-u-k_{1}i\\ \mathrm{\μpv}=w-v-k_2(v-u)\end{array}\right.---\left(4\right);$

其中，对于所述方程（3），w＝0； 

步骤S103，对所述方程（4）两端进行拉式变换可得： 

$\left\{\begin{array}{c}\mathrm{\λ}(\mathrm{sU}-u_0)=V-U-k_{1}I\\ \mathrm{\μ}(\mathrm{sV}-v_0)=W-V-k_2(V-U)\end{array}\right.---\left(5\right);$

其中，U为u的拉氏变换，V为由v的拉氏变换，W为w的拉氏变换，I为i的拉氏变换，u₀为u的初值，v₀为v的初值； 

所述步骤S1中得到所述频域内发电机d、q轴方程后对所述频域内发电机d、q轴方程经过处理后再进行拉氏反变换，对所述方程（5）进行处理后得到方程（6）： 

$\mathrm{\μ\λU}+\frac{(\mathrm{\μ}+\mathrm{\λ}+k_2\mathrm{\λ})U}{s}+\frac{{\mathrm{\μk}}_{1}I}{s}-\frac{\mathrm{\μ\λ}{u}_0}{s}+\frac{(k_1+k_1{k}_2)I}{s^2}-\frac{{\mathrm{\μV}}_0+{\mathrm{\λu}}_0+k_2{\mathrm{\λu}}_0}{s^2}=\frac{W-U}{s^2}---\left(6\right);$

对所述方程（6）进行拉氏反变换得到所述时域内发电机d、q轴积分方程（7）： 

$\mathrm{\μ\λu}+(\mathrm{\μ}+\mathrm{\λ}+k_2\mathrm{\λ}){\∫}_{0\_}^t\mathrm{udt}+{\mathrm{\μk}}_1{\∫}_{0\_}^t\mathrm{idt}+(k_1+k_1{k}_2){\∫}_{0\_}^t{\∫}_{0\_}^t\mathrm{idtdt}-({\mathrm{\μV}}_0+{\mathrm{\λu}}_0+k_2{\mathrm{\λu}}_0)t$

${-\mathrm{\μ\λu}}_0={\∫}_{0\_}^t{\∫}_{0\_}^t[w-u]\mathrm{dtdt}---\left(7\right).$

本发明提供的第六优选实施例中：所述步骤S2中根据所述发电机模型中发电机定子d轴和q轴电压与状态变量的关系式分别为： 

d轴：E″_q＝U_q+X″_di_d；q轴：E″_d＝U_d-X″_qi_q； 

对所述时域内发电机d轴和q轴积分方程里的状态变量进行替换后得到： 

d轴： $\begin{array}{c}\mathrm{\μ\λ}{U}_q+(\mathrm{\μ}+(1+k_2)\mathrm{\λ}){\∫}_{0\_}^t{U}_q\mathrm{dt}+\mathrm{\μ\λ}{X}_d^{\′\′}{i}_d+(\mathrm{\μ}{X}_d^{\′}+(1+k_2)\mathrm{\λ}{X}_d^{\′\′}){\∫}_{0\_}^t{i}_d\mathrm{dt}\\ +X_d{\∫}_{0\_}^t{\∫}_{0\_}^t{i}_d\mathrm{dtdt}-(\mathrm{\μ}{V}_0+\mathrm{\λ}{u}_0+k_2\mathrm{\λ}{u}_0)t-\mathrm{\μ\λ}{u}_0={\∫}_{0\_}^t{\∫}_{0\_}^t[w-U_q]\mathrm{dtdt}\end{array}---\left(8\right);$

q轴： $\begin{array}{c}\mathrm{\μ\λ}{U}_d+(\mathrm{\μ}+\mathrm{\λ}+k_2\mathrm{\λ}){\∫}_{0\_}^t{U}_d{d}_t-\mathrm{\μ\λ}{X_q}^{\′\′}{i}_q-(\mathrm{\μ}{X_q}^{\′}+(1+k_2)\mathrm{\λ}{X_q}^{\′\′}){\∫}_{0\_}^t{i}_q{d}_t\\ -X_q{\∫}_{0\_}^t{\∫}_{0\_}^t{i}_q{d}_t{d}_t-(\mathrm{\μ}{V}_0+\mathrm{\λ}{u}_0+k_2\mathrm{\λ}{u}_0)t-\mathrm{\μ\λ}{u}_0=-{\∫}_{0\_}^t{\∫}_{0\_}^t{U}_d{d}_t{d}_t\end{array}---\left(9\right).$

本发明提供的第七优选实施例中：所述步骤S3中所述时域内发电机d、q轴积分方程离散为矩阵形式为AX=Y； 

将所述方程（8）或方程（9）写成矩阵的形式如方程（10）： 

AX=Y      (10)； 

其中对于d轴状态方程（8）： 

$A=[U_q,{\∫}_{0\_}^t{U}_q\mathrm{dt},i_d,{\∫}_{0\_}^t{i}_d\mathrm{dt},{\∫}_{0\_}^t\left({\∫}_{0\_}^t{i}_d\mathrm{dt}\right)\mathrm{dt},-t,-1],Y={\∫}_{0\_}^t\left({\∫}_{0\_}^t(w-U_q)\mathrm{dt}\right)\mathrm{dt},$

X＝[μλ,(μ+λ+k₂λ),μλX″_d,μX′_d+(1+k₂)λX″_d,X_d,(1+k₂)λu₀+μv₀,μλu₀]； 

对于q轴状态方程（9）： 

$A=[U_d,{\∫}_{0\_}^t{U}_d{d}_t,-i_q,-{\∫}_{0\_}^t{i}_q{d}_t,-{\∫}_{0\_}^t{\∫}_{0\_}^t{i}_q{d}_t{d}_t,-t,-1],Y=-{\∫}_{0\_}^t{\∫}_{0\_}^t{U}_d{d}_t{d}_t,$

X＝[μλ,(μ+λ+k₂λ),μλX_q″,μX_q′+(1+k₂)λX_q″,X_q,(1+k₂)λu₀+μv₀,μλu₀]。 

本发明提供的第八优选实施例中：所述步骤S4中： 

所述X与所述发电机d轴实用参数之间的关系为： 

$\left\{\begin{array}{c}{T}_{d0}^{\′}{T}_{d0}^{\′\′}=X\left(1\right)\\ T_{d0}^{\′}+T_{d0}^{\′\′}\frac{X_d-X_d^{\′\′}}{X_d^{\′}-X_d^{\′\′}}=X\left(2\right)\\ T_{d0}^{\′}{T}_{d0}^{\′\′}{X}_d^{\′\′}=X\left(3\right)\\ T_{d0}^{\′}{X}_d^{\′}+X_d^{\′\′}{T}_{d0}^{\′\′}\frac{X_d-X_d^{\′\′}}{X_d^{\′}-X_d^{\′\′}}=X\left(4\right)\\ X_d=X\left(5\right)\\ T_{d0}^{\′\′}\frac{X_d-X_d^{\′\′}}{X_d^{\′}-X_d^{\′\′}}{E}_{q\left(0\right)}^{\′\′}+T_{d0}^{\′}{E}_{q\left(0\right)}^{\′}=X\left(6\right)\\ T_{d0}^{\′}{T}_{d0}^{\′\′}{E}_{q\left(0\right)}^{\′\′}=X\left(7\right)\end{array}\right.---\left(11\right);$

所述X与所述发电机q轴实用参数之间的关系为： 

$\left\{\begin{array}{c}{T}_{q0}^{\′}{T}_{q0}^{\′\′}=X\left(1\right)\\ T_{q0}^{\′}+T_{q0}^{\′\′}\frac{X_q-X_q^{\′\′}}{X_q^{\′}-X_q^{\′\′}}=X\left(2\right)\\ T_{q0}^{\′}{T}_{q0}^{\′\′}{X}_q^{\′\′}=X\left(3\right)\\ T_{q0}^{\′}{X}_q^{\′}+X_q^{\′\′}{T}_{q0}^{\′\′}\frac{X_q-X_q^{\′\′}}{X_q^{\′}-X_q^{\′\′}}=X\left(4\right)\\ X_q=X\left(5\right)\\ T_{q0}^{\′\′}\frac{X_q-X_q^{\′\′}}{X_q^{\′}-X_q^{\′\′}}{E}_{d\left(0\right)}^{\′\′}+T_{q0}^{\′}{E}_{d\left(0\right)}^{\′}=X\left(6\right)\\ T_{q0}^{\′}{T}_{q0}^{\′\′}{E}_{d\left(0\right)}^{\′\′}=X\left(7\right)\end{array}\right.---\left(12\right).$

本发明提供的第九优选实施例中：所述当试验前为稳态工况时，将d轴状态变量初值u₀＝E″_q(0)＝U_q0+X″_di_d0，v₀＝E′_q(0)＝U_q0+X′_di_d0或q轴状态变量初值 u₀＝E″_d(0)＝U_d0-X″_qi_q0v₀＝E′_d(0)＝U_d0-X′_qi_q0代入方程（8）或（9）中从而消去两个初值u₀和v₀参数从而减少辨识向量X的个数。 

本发明提供的一种基于时频变换的同步发电机参数辨识方法的有益效果包括： 

1、本发明提供的一种基于时频变换的同步发电机参数辨识方法，将发电机微分方程转换为积分方程，有利于提高数值计算的稳定性，较好地克服环境噪音和数值计算噪声对时域辨识法的干扰，解决辨识结果不稳定的问题。整个辨识过程中无需迭代，计算速度很快，便于利用PMU数据实现在线辨识。 

2、辨识方法可以基于模型结构发生变化的大扰动试验，由于大扰动试验时阻尼绕组中能够产生足够大的感应电流，由此对机端电气量产生足够可观测的变化，提高机端电气量的信噪比，有利于辨识的开展。可以基于机组投运前必做的甩负荷试验，工程实施非常方便，具有很强的工程应用价值。 

具体实施方式

下面对本发明的具体实施方式作进一步详细说明。 

本发明提供一种基于时频变换的同步发电机参数辨识方法，包括： 

步骤S1，将发电机模型中的时域内发电机d、q轴微分方程进行拉氏变换得到频域内发电机d、q轴方程，再经过拉氏反变换得到时域内发电机d、q轴积分方程。 

步骤S2，根据发电机模型中发电机定子d、q轴电压与状态变量的关系式，对时域内发电机d、q轴积分方程里的状态变量进行替换。 

步骤S3，假设发电机定子d、q轴电压和电流为由实测采样点构成的列向量，将时域内发电机d、q轴积分方程离散为矩阵与积分方程系数列向量相乘等于可测列向量或者可测列向量及其积分的形式，求解得到该积分方程系数列向量。 

步骤S4，根据该积分方程系数列向量X与发电机参数之间的关系得到发电机模型参数。 

进一步的，步骤S1中在对时域内发电机模型进行拉氏变换时，需要包含状态变量的初值； 

步骤S1中得到所述频域内发电机d、q轴方程后对所述频域内发电机d、q轴方程经过处理后再进行拉氏反变换，处理的过程包括： 

消去频域内发电机d轴或q轴方程的中间状态变量，使该频域内发电机d轴或q轴方程中分别只保留一个状态变量，再在该频域内发电机d轴或q轴方程两侧分别除以本方程最高次幂的拉氏算子。 

步骤后S2中利用发电机模型中发电机定子d、q轴电压与状态变量的关系式将状态变量 用机端电压或机端电流等可测电气量替换。 

时域内发电机d、q轴积分方程可以离散为矩阵形式为AX=Y。 

其中A为由各可测列向量及其各种积分形式组成的矩阵，X为待辨识的参数构成的积分方程系数列向量，其各分量是发电机实用参数的代数组合，Y为可测列向量或某些可测列向量及其积分的代数和。 

步骤S1中将发电机微分方程转换为积分方程，有利于提高数值计算的稳定性，较好地克服环境噪音和数值计算噪声对时域辨识法的干扰，解决辨识结果不稳定的问题。整个辨识过程中无需迭代，计算速度很快，便于利用PMU数据实现在线辨识。 

此外，本发明提供的辨识方法可以基于模型结构发生变换的大扰动试验，由于大扰动试验时阻尼绕组中能够产生足够大的感应电流，由此对机端电气量产生足够可观测的变化，提高机端电气量的信噪比，有利于辨识的开展。可以基于机组投运前必做的甩负荷试验，工程实施非常方便，具有很强的工程应用价值。 

进一步的，本发明提供一种基于时频变换的同步发电机参数辨识方法的实施例中，辨识过程采用的发电机模型的实施例为BPA仿真软件模型，如方程（1）所示： 

$\left\{\begin{array}{c}{U}_d={\mathrm{\ωE}}_d^{\′\′}+{\mathrm{\ωX}}_q^{\′\′}{i}_q-R_a{i}_d\\ U_q={\mathrm{\ωE}}_q^{\′\′}-{\mathrm{\ωX}}_d^{\′\′}{i}_d-R_a{i}_q\\ T_{d0}^{\′}{\mathrm{pE}}_q^{\′}=E_f-E_q^{\′}-\frac{X_d-X_d^{\′}}{X_d^{\′}-X_d^{\′\′}}(E_q^{\′}-E_q^{\′\′})\\ T_{d0}^{\′\′}{\mathrm{pE}}_q^{\′\′}=E_q^{\′}-E_q^{\′\′}-(X_d^{\′}-X_d^{\′\′})i_d\\ T_{q0}^{\′}{\mathrm{pE}}_d^{\′}=-E_d^{\′}-\frac{X_q-X_q^{\′}}{X_q^{\′}+X_q^{\′\′}}(E_d^{\′}-E_d^{\′\′})\\ T_{q0}^{\′\′}{\mathrm{pE}}_d^{\′\′}=E_d^{\′}-E_d^{\′\′}+(X_q^{\′}-X_q^{\′\′})i_q\\ T_J\mathrm{p\ω}=T_m-[E_q^{\′\′}{i}_q+E_d^{\′\′}{i}_d-(X_d^{\′\′}-X_q^{\′\′})i_d{i}_q]-D(\mathrm{\ω}-1)\\ \mathrm{p\δ}=\mathrm{\ω}-1\end{array}\right.---\left(1\right)$

方程(1)中E″_d，E′_d，E″_q，E′_q，δ，ω为六个状态变量，通常称为六阶实用模型；X_d、X_q、X′_d、X′_q、X″_d、X″_q、T′_d0、T′_q0、T″_d0、T″_q0、R_a、T_J为发电机参数，X_d为直轴同步电抗；X_q为交轴同步电抗；X′_d为直轴暂态电抗；X′_q交轴暂态电抗；X″_d直轴次暂态电抗；X″_q交轴次暂态电抗；T′_d0、T′_q0、T″_d0、T″_q0、为时间常数；R_a为电枢电阻，通常很小可以忽略；T_J为转动惯量；U_d、U_q分别为发电机定子d、q轴电压；i_d、i_q分别为发电机定子d、q轴电流；E_f为励磁电压；T_m为原动机转矩；D为阻尼系数；p为微分算子

方程（1）中前两个方程为发电机定子d、q轴电压与状态变量的关系式，工程应用中认为定子绕组的暂态对电力稳定运行的影响不大，故发电机模型中忽略了定子绕组的暂态而将定子回路以代数方程描述。 

方程（1）中关于d轴的状态方程如（2）所示 

$\left\{\begin{array}{c}{T}_{d0}^{\′\′}{\mathrm{pE}}_q^{\′\′}=E_q^{\′}-E_q^{\′\′}-(X_d^{\′}-X_d^{\′\′})i_d\\ T_{d0}^{\′}{\mathrm{pE}}_q^{\′}=E_f-E_q^{\′}-\frac{X_d-X_d^{\′}}{X_d^{\′}-X_d^{\′\′}}(E_q^{\′}-E_q^{\′\′})\end{array}\right.---\left(2\right)$

方程（1）中关于q轴的状态方程为： 

$\left\{\begin{array}{c}{T}_{q0}^{\′\′}{\mathrm{pE}}_d^{\′\′}=E_d^{\′}-E_d^{\′\′}+(X_q^{\′}-X_q^{\′\′})i_q\\ T_{q0}^{\′}{\mathrm{pE}}_d^{\′}=-E_d^{\′}-\frac{X_q-X_q^{\′}}{X_q^{\′}-X_q^{\′\′}}(E_d^{\′}-E_d^{\′\′})\end{array}\right.---\left(3\right)$

方程（1）中最后两式为转子运动方程，是为了确定功角δ和转速ω。若能通过测量的方法获得准确的δ和ω，则（1）中的电气方程和机械方程可以实现解耦，参数辨识。本发明中所指的辨识参数是指其电气参数，不包含对T_J和D的辨识。由于可以通过齿盘测速信号和机端电压信号获得准确的发电机机功角δ和转速ω，因此本发明在进行参数辨识时只需用到方程（1）中的前六式即电气方程。 

步骤S1中将发电机模型中的时域内发电机d、q轴微分方程进行拉氏变换得到频域内发电机d、q轴方程，该时域内发电机d、q轴微分方程即为方程（2）和（3）。 

进行拉氏变换之前进行变量的替换： 

对于方程（2），令E″_q＝u,E′_q＝v,E_f＝w,T″_do＝λ,T′_do＝μ，k₁＝X′_d-X″_d，

i_d＝i。 

对于方程（3），令E″_d＝u,E′_d＝v,T″_qo＝λ,T′_qo＝μ，k₁＝-(X_q′-X_q″)，

i_q＝i。 

则方程（2）和方程（3）可以改写成： 

$\left\{\begin{array}{c}\mathrm{\λpu}=v-u-k_{1}i\\ \mathrm{\μpv}=w-v-k_2(v-u)\end{array}\right.---\left(4\right)$

其中，对于q轴微分方程（3），w＝0。 

对方程（4）两端进行拉式变换可得： 

$\left\{\begin{array}{c}\mathrm{\λ}(\mathrm{sU}-u_0)=V-U-k_{1}I\\ \mathrm{\μ}(\mathrm{sV}-v_0)=W-V-k_2(V-U)\end{array}\right.---\left(5\right)$

其中，该拉氏变换需要包含状态变量的初值，U为u的拉氏变换，V为由v的拉氏变换， W为w的拉氏变换，I为i的拉氏变换，u₀为u的初值，v₀为v的初值。 

消去（5）中的中间状态变量V并整理可得： 

μλs²U+(μ+λ+k₂λ)sU+μk₁sI-μλu₀s+(k₁+k₁k₂)I-μV₀-λu₀-k₂λu₀＝W-U  （6） 

将方程（6）两边同时除以本方程最高次幂的拉氏算子s²可得方程（7）。 

$\mathrm{\μ\λU}+\frac{(\mathrm{\μ}+\mathrm{\λ}+k_2\mathrm{\λ})U}{s}+\frac{{\mathrm{\μk}}_{1}I}{s}-\frac{\mathrm{\μ\λ}{u}_0}{s}+\frac{(k_1+k_1{k}_2)I}{s^2}-\frac{{\mathrm{\μV}}_0+{\mathrm{\λu}}_0+k_2{\mathrm{\λu}}_0}{s^2}=\frac{W-U}{s^2}---\left(7\right)$

对方程（7）进行拉氏反变换可得方程（8）。 

$\mathrm{\μ\λu}+(\mathrm{\μ}+\mathrm{\λ}+k_2\mathrm{\λ}){\∫}_{0\_}^t\mathrm{udt}+{\mathrm{\μk}}_1{\∫}_{0\_}^t\mathrm{idt}+(k_1+k_1{k}_2){\∫}_{0\_}^t{\∫}_{0\_}^t\mathrm{idtdt}-({\mathrm{\μV}}_0+{\mathrm{\λu}}_0+k_2{\mathrm{\λu}}_0)t$

                                                             （8） 

$-\mathrm{\μ\λ}{u}_0={\∫}_{0\_}^t{\∫}_{0\_}^t[w-u]\mathrm{dtdt}$

在方程（8）中，u为q轴或d轴次暂态电动势E″_q或E″_d，通常不可测，因此方程（8）并不能直接应用。 

以d轴状态方程为例： 

由U_q＝E″_q-X″_di_d-R_ai_q并略去定子电阻可得E″_q＝U_q+X″_di_d，将其代入方程（8）并整理可得方程（9）。 

$\mathrm{\μ\λ}{U}_q+(\mathrm{\μ}+(1+k_2)\mathrm{\λ}){\∫}_{0\_}^t{U}_q\mathrm{dt}+\mathrm{\μ\λ}{X}_d^{\′\′}{i}_d+(\mathrm{\μ}{X}_d^{\′}+(1+k_2)\mathrm{\λ}{X}_d^{\′\′}){\∫}_{0\_}^t{i}_d\mathrm{dt}$

$+X_d{\∫}_{0\_}^t{\∫}_{0\_}^t{i}_d\mathrm{dtdt}-({\mathrm{\μV}}_0+{\mathrm{\λu}}_0+k_2{\mathrm{\λu}}_0)t-{\mathrm{\μ\λu}}_0={\∫}_{0\_}^t{\∫}_{0\_}^t[w-U_q]\mathrm{dtdt}---\left(9\right)$

假设U_q和i_d等为由实测采样点构成的列向量，则将方程（9）写成矩阵的形式如方程（10）。 

AX=Y                  （10） 

其中 $A=[U_q,{\∫}_{0\_}^t{U}_q\mathrm{dt},i_d,{\∫}_{0\_}^t{i}_d\mathrm{dt},{\∫}_{0\_}^t\left({\∫}_{0\_}^t{i}_d\mathrm{dt}\right)\mathrm{dt},-t,-1],Y={\∫}_{0\_}^t\left({\∫}_{0\_}^t(w-U_q)\mathrm{dt}\right)\mathrm{dt},$

X＝[μλ,(μ+λ+k₂λ),μλX″_d,μX′_d+(1+k₂)λX″_d,X_d,(1+k₂)λu₀+μv₀,μλu₀]。 

求解方程（10）可得X，X中包含了5个d轴参数X_d、X′_d、X″_d、T′_d0、T″_d0和两个初值u₀和v₀共七个未知参数，根据X的七个分量和七个未知参数的关系可计算出d轴所有参数，X与发电机d轴实用参数之间的关系如方程（11）所示。 

$\left\{\begin{array}{c}{T}_{d0}^{\′}{T}_{d0}^{\′\′}=X\left(1\right)\\ T_{d0}^{\′}+T_{d0}^{\′\′}\frac{X_d-X_d^{\′\′}}{X_d^{\′}-X_d^{\′\′}}=X\left(2\right)\\ T_{d0}^{\′}{T}_{d0}^{\′\′}{X}_d^{\′\′}=X\left(3\right)\\ T_{d0}^{\′}{X}_d^{\′}+X_d^{\′\′}{T}_{d0}^{\′\′}\frac{X_d-X_d^{\′\′}}{X_d^{\′}-X_d^{\′\′}}=X\left(4\right)\\ X_d=X\left(5\right)\\ T_{d0}^{\′\′}\frac{X_d-X_d^{\′\′}}{X_d^{\′}-X_d^{\′\′}}{E}_{q\left(0\right)}^{\′\′}+T_{d0}^{\′}{E}_{q\left(0\right)}^{\′}=X\left(6\right)\\ T_{d0}^{\′}{T}_{d0}^{\′\′}{E}_{q\left(0\right)}^{\′\′}=X\left(7\right)\end{array}\right.---\left(11\right)$

同理，对于q轴状态方程，将E″_d＝U_d-X″_qi_q代入方程（8）中得到方程（12）。 

${\mathrm{\μ\λU}}_d+(\mathrm{\μ}+\mathrm{\λ}+k_2\mathrm{\λ}){\∫}_{0\_}^t{U}_d{d}_t-{{\mathrm{\μ\λX}}_q}^{\′\′}{i}_q-({{\mathrm{\μX}}_q}^{\′}+(1+k_2){{\mathrm{\λX}}_q}^{\′\′}){\∫}_{0\_}^t{i}_q{d}_t---\left(12\right)$

${-X}_q{\∫}_{0\_}^t{\∫}_{0\_}^t{i}_q{d}_t{d}_t-({\mathrm{\μV}}_0+{\mathrm{\λu}}_0+k_2{\mathrm{\λu}}_0)t-{\mathrm{\μ\λu}}_0=-{\∫}_{0\_}^t{\∫}_{0\_}^t{U}_d{d}_t{d}_t$

写成矩阵的形式中： 

$A=[U_d,{\∫}_{0\_}^t{U}_d{d}_t,-i_q,-{\∫}_{0\_}^t{i}_q{d}_t,-{\∫}_{0\_}^t{\∫}_{0\_}^t{i}_q{d}_t{d}_t,-t,-1],Y=-{\∫}_{0\_}^t{\∫}_{0\_}^t{U}_d{d}_t{d}_t,$

X＝[μλ,(μ+λ+k₂λ),μλX_q″,μX_q′+(1+k₂)λX_q″,X_q,(1+k₂)λu₀+μV₀,μλu₀]。 

X与发电机q轴实用参数之间的关系如方程（13）所示。 

$\left\{\begin{array}{c}{T}_{q0}^{\′}{T}_{q0}^{\′\′}=X\left(1\right)\\ T_{q0}^{\′}+T_{q0}^{\′\′}\frac{X_q-X_q^{\′\′}}{X_q^{\′}-X_q^{\′\′}}=X\left(2\right)\\ T_{q0}^{\′}{T}_{q0}^{\′\′}{X}_q^{\′\′}=X\left(3\right)\\ T_{q0}^{\′}{X}_q^{\′}+X_q^{\′\′}{T}_{q0}^{\′\′}\frac{X_q-X_q^{\′\′}}{X_q^{\′}-X_q^{\′\′}}=X\left(4\right)\\ X_q=X\left(5\right)\\ T_{q0}^{\′\′}\frac{X_q-X_q^{\′\′}}{X_q^{\′}-X_q^{\′\′}}{E}_{d\left(0\right)}^{\′\′}+T_{q0}^{\′}{E}_{d\left(0\right)}^{\′}=X\left(6\right)\\ T_{q0}^{\′}{T}_{q0}^{\′\′}{E}_{d\left(0\right)}^{\′\′}=X\left(7\right)\end{array}\right.---\left(13\right)$

当试验前为稳态工况时，还可以将d轴状态变量初值u₀＝E″_q(0)＝U_q0+X″_di_d0，v₀＝E′_q(0)＝U_q0+X′_di_d0或q轴状态变量初值u₀＝E″_d(0)＝U_d0-X″_qi_q0v₀＝E′_d(0)＝U_d0-X′_qi_q0代入方程（9）或（12）中从而消去两个初值u₀和v₀参数从而减少辨识向量X的个数。 

根据上述辨识方法并结合现场某台发电机甩50％负荷试验辨识所得的发电机参数如表1 所示。 

表1甩负荷辨识数据 

表1中/表示未给出设计值。 

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求范围当中。 

Claims (10)

1.一种基于时频变换的同步发电机参数辨识方法，其特征在于，所述方法包括： 

步骤S1，将发电机模型中的时域内发电机d、q轴微分方程进行拉氏变换得到频域内发电机d、q轴方程，再进行拉氏反变换得到时域内发电机d、q轴积分方程； 

步骤S2，根据所述发电机模型中发电机定子d、q轴电压与状态变量的关系式，对所述时域内发电机d、q轴积分方程里的状态变量进行替换； 

步骤S3，假设发电机定子d、q轴电压和电流为由实测采样点构成的列向量，将所述时域内发电机d、q轴积分方程离散为矩阵与积分方程系数列向量相乘等于可测列向量或者可测列向量积分的形式，求解得到积分方程系数列向量； 

步骤S4，根据所述积分方程系数列向量与发电机参数之间的关系求解得到发电机模型参数。 

2.如权利要求1所述的方法，其特征在于，所述方法基于机组投运前必做的甩负荷试验或其它在线大扰动试验。 

3.如权利要求1所述的方法，其特征在于， 

所述步骤S1中在对时域内发电机模型进行拉氏变换时，需要包含所述状态变量的初值； 

所述步骤S1中得到所述频域内发电机d、q轴方程后对所述频域内发电机d、q轴方程经过处理后再进行拉氏反变换，处理的过程包括： 

消去所述频域内发电机d轴或q轴方程的中间状态变量，使频域内发电机d轴或q轴方程中分别只保留一个状态变量，再在所述频域内发电机d轴或q轴方程两侧分别除以本方程最高次幂的拉氏算子。 

4.如权利要求1所述的方法，其特征在于，所述步骤后S2中利用所述发电机模型中发电机定子d、q轴电压与状态变量的关系式将所述状态变量用可测电气量替换，所述可测电气量包括机端电压或机端电流； 

所述步骤S3中所述时域内发电机d、q轴积分方程离散为矩阵形式为AX=Y； 

其中A为由各可测列向量及其各种积分形式组成的矩阵，X为待辨识的发电机参数构成的方程系数列向量，其各分量是发电机实用参数的代数组合，Y为可测列向量或某些可测列向量及其积分的代数和。 

5.如权利要求1所述的方法，其特征在于，所述发电机模型为： 

其中，E″_d，E′_d，E″_q，E′_q，δ，ω为六个状态变量，通常称为六阶实用模型；X_d、X_q、X′_d、X′_q、X″_d、X″_q、T′_d0、T′_q0、T″_d0、T″_q0、R_a为发电机参数，X_d为直轴同步电抗；X_q为交轴同步电抗；X′_d为直轴暂态电抗；X′_q交轴暂态电抗；X″_d直轴次暂态电抗；X″_q交轴次暂态电抗；T′_d0、T′_q0、T″_d0、T″_q0、为时间常数；R_a为电枢电阻；U_d、U_q分别为发电机定子d、q轴电压；i_d、i_q分别为发电机定子d、q轴电流；E_f为励磁电压；p为微分算子

所述方程（1）中d轴的状态方程为： 

所述方程（1）中q轴的状态方程为： 

6.如权利要求5所述的方法，其特征在于， 

所述步骤S1中对所述发电机模型中的所述时域内发电机d、q轴微分方程进行拉氏变换得到所述频域内发电机d、q轴方程的方法包括： 

步骤S101，对于所述方程（2），令E″_q＝u,E′_q＝v,E_f＝w,T″_do＝λ,T′_do＝μ，k₁＝X′_d-X″_d， 

i_d＝i；对于所述方程（3），令E″_d＝u,E′_d＝v,T″_qo＝λ,T′_qo＝μ，k₁＝-(X_q′-X_q″)， i_q＝i； 

步骤S102，将方程（2）和方程（3）改写成： 

其中，对于所述方程（3），w＝0； 

步骤S103，对所述方程（4）两端进行拉式变换可得： 

其中，U为u的拉氏变换，V为由v的拉氏变换，W为w的拉氏变换，I为i的拉氏变换，u₀为u的初值，v₀为v的初值； 

所述步骤S1中得到所述频域内发电机d、q轴方程后对所述频域内发电机d、q轴方程经过处理后再进行拉氏反变换，对所述方程（5）进行处理后得到方程（6）： 

对所述方程（6）进行拉氏反变换得到所述时域内发电机d、q轴积分方程（7）： 

                                                            （7） 

。

7.如权利要求6所述的方法，其特征在于，所述步骤S2中根据所述发电机模型中发电机定子d轴和q轴电压与状态变量的关系式分别为： 

d轴：E″_q＝U_q+X″_di_d；q轴：E″_d＝U_d-X″_qi_q； 

对所述时域内发电机d轴和q轴积分方程里的状态变量进行替换后得到： 

d轴：

q轴：

8.如权利要求7所述的方法，其特征在于，所述步骤S3中所述时域内发电机d、q轴积分方程离散为矩阵形式为AX=Y； 

将所述方程（8）或方程（9）写成矩阵的形式如方程（10）： 

AX=Y                              （10）； 

其中对于d轴状态方程（8）： 

X＝[μλ,(μ+λ+k₂λ),μλX″_d,μX′_d+(1+k₂)λX″_d,X_d,(1+k₂)λu₀+μv₀,μλu₀]； 

对于q轴状态方程（9）： 

X＝[μλ,(μ+λ+k₂λ),μλX″_q,μX′_q+(1+k₂)λX″_q,X_q,(1+k₂)λu₀+μv₀,μλu₀]。 

9.如权利要求8所述的方法，其特征在于，所述步骤S4中： 

所述X与所述发电机d轴实用参数之间的关系为： 

所述X与所述发电机q轴实用参数之间的关系为： 

。

10.如权利要求7所述的方法，其特征在于，所述当试验前为稳态工况时，将d轴状态变量初值u₀＝E″_q(0)＝U_q0+X″_di_d0，v₀＝E′_q(0)＝U_q0+X′_di_d0或q轴状态变量初值u₀＝E″_d(0)＝U_d0-X″_qi_q0v₀＝E′_d(0)＝U_d0-X′_qi_q0代入方程（8）或（9）中从而消去两个初值u₀和v₀参数从而减少辨识向量X的个数。