CN102403731B - 微型燃气轮机发电系统的模拟方法 - Google Patents

Abstract

本发明公开了一种微型燃气轮机发电系统的模拟方法，在该微型燃气轮机发电系统中，微型燃气轮机发电系统通过整流逆变装置向电网输出电能，将微型燃气轮机发电系统用以电感电流、电容电压为状态变量的三阶动态微分方程来表征；并最终求出微型燃气轮机发电系统的有功响应和无功响应：本发明的微型燃气轮机发电系统的模拟方法，能准确地模拟微型燃气轮机发电系统在并网运行条件下的稳态和暂态特性，可作为研究含有微型燃气轮机发电系统的研究工具。

Description

微型燃气轮机发电系统的模拟方法

技术领域

本发明涉及一种微型燃气轮机发电系统的模拟方法。

背景技术

人们对电网安全、经济、交互等方面的需要和对电能质量高标准的要求，促使传统电网亟需向智能电网转变。分布式发电以其污染少、可靠性高、能源利用率高、安装地点灵活等优点而备受关注，将成为智能电网的重要组成部分。分布式电源大多接入电网10kV及以下的低压配电网中，改变了传统配电网的负荷组成和拓扑结构，使得传统配电网综合负荷模型已不再能很好的描述配电网综合负荷特性，考虑分布式电源影响的综合负荷建模已成为电力系统负荷建模发展的主要趋势之一。

微型燃气轮机作为目前最具有商业竞争力的分布式电源之一，具有可使用燃料丰富、体积小、重量轻、污染小、运行维护简单、易于实现冷热电联产等优点，因而受到人们的日益关注。因此，在鼓励分布式电源高渗透率的大环境下，研究含有微型燃气轮机发电的配电网综合负荷模型具有重要意义。

如何构建总体测辨法建模所需的微型燃气轮机仿真系统及如何对其进行等效以满足电网仿真计算对等效模型的要求，是研究含微型燃气轮机的配电网综合负荷建模面临两个必须解决的问题。国内外对微型燃气轮机发电系统仿真控制技术进行了相关研究。有文献采用斩波器来模拟微型燃气轮机发电系统直流侧的电特性，并通过直流侧模型进行近似线性化的和降阶处理，得到其传递函数；有文献通过将用永磁同步发电机及整流器部分作为一个整体来建模，建立了微型燃气轮机发电系统的整体模型；有文献构建了采用阻尼控制器控制的微型燃气轮机发电系统，分析了微型燃气轮机发电系统的动态性能，验证了所提出的阻尼控制器的有效性；有文献构建了简化的微型燃气轮机发电系统模型，分析了其在并网和孤网运行方式下的动态性能；有文献对并网的微型燃气轮机发电系统的各模块的数学模型进行了建模仿真。这些文献偏重于微型燃气轮机发电系统的动态性能及并网控制方法的研究上，缺乏既能反映自身动态特性又能满足负荷建模要求的等效模型，而针对问题二的研究未见报道。

发明内容

本发明所要解决的技术问题是提供一种微型燃气轮机发电系统的模拟方法，该方法提供了一种微型燃气轮机发电系统的等效模型，能准确地模拟微型燃气轮机发电系统在并网运行条件下的稳态和暂态特性，可作为研究含有微型燃气轮机发电系统的研究工具。

发明的技术解决方案如下：

一种微型燃气轮机发电系统的模拟方法，在该微型燃气轮机发电系统中，微型燃气轮机发电系统通过整流逆变装置向电网输出电能，将微型燃气轮机发电系统用以下以电感电流、电容电压为状态变量的三阶动态微分方程来表征：

$\left\{\begin{array}{c}\frac{{\mathrm{dI}}_{L.d}}{\mathrm{dt}}=\frac{1}{L}(U_{i.d}-U_{g.d}-I_{L.d}R)-\mathrm{\ω}{I}_{L.q}\\ \frac{{\mathrm{dI}}_{L.q}}{\mathrm{dt}}=\frac{1}{L}(U_{i.q}-U_{g.q}-I_{L.q}R)+\mathrm{\ω}{I}_{L.d}\\ \frac{{\mathrm{dU}}_{\mathrm{dc}}}{\mathrm{dt}}=\frac{1}{C}(\frac{P_s}{U_{\mathrm{dc}}}-S_d{I}_{L.d}-S_q{I}_{L.q})\end{array}\right.;$

其中，I_L.d、I_L.q是逆变器输出电流的d轴分量和q轴分量，PCC为公共连接点，逆变器的输出电流也即微型燃气轮机发电系统由PCC点注入电网的电流；U_i是逆变器的出口电压；U_g.d、U_g.q是微型燃气轮机发电系统与电网的公共连接点电压的d轴分量和q轴分量；R、L、C分别为等值电阻、等值电感和等值电容器；P_S为微型燃气轮机原动机输入功率；U_dc为直流母线电压；S_d、S_q为逆变器同步坐标系下的平均开关函数，通过以下修正公式来不断修正：

$\left\{\begin{array}{c}\mathrm{\Δ}{U}_{\mathrm{dc}}=U_{\mathrm{dc}}\left(t\right)-U_{\mathrm{dc}}^{*}\\ \mathrm{\Δ}{I}_d=I_q^{*}-I_{L.d}\left(t\right)\\ \mathrm{\Δ}{I}_q=I_q^{*}-I_{L.q}\left(t\right)\\ I_d^{*}\left(t\right)=I_d^{*}+K_{p1}\mathrm{\Δ}{U}_{\mathrm{dc}}+K_{i1}{\∫}_{t_0}^t\mathrm{\Δ}{U}_{\mathrm{dc}}\mathrm{dt}\\ S_d\left(t\right)=S_d+K_{p2}\mathrm{\Δ}{I}_d+K_{i2}{\∫}_{t_0}^t\mathrm{\Δ}{I}_d\mathrm{dt}\\ S_q\left(t\right)=S_q+K_{p2}\mathrm{\Δ}{I}_q+K_{i2}{\∫}_{t_0}^t\mathrm{\Δ}{I}_q\mathrm{dt}\end{array}\right.;$

其中，k_p1、k_i1、k_p2、k_i2为控制环节参数；ΔU_dc、ΔI_d、ΔI_q分别为相应变量在单位时刻的变化量；I_L.d(t)、I_L.q(t)、S_d(t)、S_q(t)、

分别为I_L.d、I_L.q、S_d、S_q、

在t时刻的值；分别为U_dc、I_L.d、I_L.q、I_d(t)的参考值，解出电流I_L.d、I_L.q的初值【I_L.d、I_L.q的初值就是I_L.d(0)、I_L.q(0)】之后，通过给定S_d＝0.87、S_q＝0，通过下式求得U_dc的初值，将U_dc的初值赋予

即可获取

$\left\{\begin{array}{c}{P}_s=U_{i.d}{I}_{L.d}+U_{i.q}{I}_{L.q}\\ \frac{P_s}{U_{\mathrm{dc}}}=S_d{I}_{L.d}+S_q{I}_{L.q}\end{array}\right.;$

在后续辨识中保持不变，

取值为0以使得逆变器传输无功为0，通过初值求解求得I_L.d的初值【参见后文的上述状态变量值I_L.x(0)、I_L.y(0)、U_g.x(0)、U_g.y(0)、U_i.x(0)、U_i.y(0)通过派克变化到同步坐标系d、q轴，即求得I_L.d(0)，I_L.q(0)，】，将I_L.d的初值作为

的初始值，在后续过程中动态修正；

求解前述的三阶动态微分方程求得任意时刻的电流的d、q轴分量I_L.d和I_L.q，再对I_L.d和I_L.q进行派克反变换，即可求出同步坐标下的电流I_L.x、I_L.y，进而求出微型燃气轮机发电系统的有功响应和无功响应：

$\left\{\begin{array}{c}P=U_{g.x}{I}_{L.x}+U_{g.y}{I}_{L.y}\\ Q=U_{g.y}{I}_{L.x}-U_{g.x}{I}_{L.y}\end{array}\right.,$ 其中，U_g.x和U_g.y分别为同步坐标系下的公共连接点处的x电压分量和y电压分量，【U_g.x和U_g.y的取值与给定的电压激励有关，对于负荷建模研究来言，一般取给定的电压激励的相角为参考相角，由此建立的同步坐标下，U_g.x等于给定的电压激励值，而U_g.y等于0；】【向量包括幅值和相角，给定的电压机理的相角作为参考值，即其他向量的相角与参考值之间的差值作为该向量的相角。】

获取初始值并求解微分方程的步骤为：

所述的初始值包括已知的初值u₀＝[U_g.x(0)U_g.y(0)]^T、y₀＝[P(0)Q(0)]^T和求得的初值U_i.x(0)U_i.y(0)、I_L.x(0)、I_L.y(0)；

其中： $U_{i.x}\left(0\right)=U_{g.x}\left(0\right)+\frac{P\left(0\right)R+Q\left(0\right)\mathrm{\ωL}}{U_g\left(0\right)};$

$U_{i.y}\left(0\right)=U_{g.y}\left(0\right)+\frac{P\left(0\right)\mathrm{\ωL}-Q\left(0\right)R}{U_g\left(0\right)};$

$I_{L.x}\left(0\right)=\frac{(P\left(0\right)R+Q\left(0\right)\mathrm{\ωL})R}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}+\frac{(P\left(0\right)\mathrm{\ωL}-Q\left(0\right)R)\mathrm{\ωL}}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)};$

$I_{L.y}\left(0\right)=\frac{(P\left(0\right)\mathrm{\ωL}-Q\left(0\right)R)R}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}-\frac{(P\left(0\right)R+Q\left(0\right)\mathrm{\ωL})\mathrm{\ωL}}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)};$

其中，P(0)、Q(0)即微型燃气轮机发电系统上网的有功和无功的初始值；求解微分方程的过程为：将上述初始值进行派克变换，得到所有初值的d、q轴分量，作为求解微分方程的初始值，最后用四阶龙格库塔法解微分方程求解所述的三阶动态微分方程。

参数辨识方法如下：

通过仿真或实测在微型燃气轮机发电系统并网处设置不同程度的短路故障，获取建模样本数据，包括微型燃气轮机发电系统上网点的电压U、有功功率P、无功功率Q数据；

采用遗传算法进行等效模型参数辨识，得到最终的等效模型-即由三阶动态微分方程和修正公式表征的微型燃气轮机发电系统的等效模型；

等效模型参数包括k_p1、k_i1、k_p2、k_i2、R、L、C。

具体的辨识步骤如下：

步骤1：给定实测样本u(k)，P(k)，Q(k)(k＝0，2，…，L-1)，对遗传寻优迭代次数i进行初始化，令i＝0；

步骤2：随机产生独立待辨识参数样本R、L、C、k_p1、k_i1、k_p2、k_i2；

步骤3：设定遗传寻优迭代次数i＝i+1；

步骤4：当k＝0时，即为初值求解，有：将实测样本代入下列公式求变量初值I_L.x(k)、I_L.y(k)、U_dc(k)：【这里k＝0是初始时刻的求解过程，其与一般的点的求解过程是不一样的。对于初值点来说，U、P、Q都是已知的；而对于非初始时刻，只有电压激励U是已知的，P、Q是未知的。】

$U_{i.x}\left(0\right)=U_{g.x}\left(0\right)+\frac{P\left(0\right)R+Q\left(0\right)\mathrm{\ωL}}{U_g\left(0\right)}---\left(15\right)$

$U_{i.y}\left(0\right)=U_{g.y}\left(0\right)+\frac{P\left(0\right)\mathrm{\ωL}-Q\left(0\right)R}{U_g\left(0\right)}---\left(16\right)$

$U_i\left(0\right)=\sqrt{U_{i.x}{\left(0\right)}^2+U_{i.y}{\left(0\right)}^2}---\left(17\right)$

$I_{L.x}\left(0\right)=\frac{(P\left(0\right)R+Q\left(0\right)\mathrm{\ωL})R}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}+\frac{(P\left(0\right)\mathrm{\ωL}-Q\left(0\right)R)\mathrm{\ωL}}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}---\left(18\right)$

$I_{L.y}\left(0\right)=\frac{(P\left(0\right)\mathrm{\ωL}-Q\left(0\right)R)R}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}-\frac{(P\left(0\right)R+Q\left(0\right)\mathrm{\ωL})\mathrm{\ωL}}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}---\left(19\right)$

将上述状态变量值I_L.x(0)、I_L.y(0)、U_g.x(0)、U_g.y(0)、U_i.x(0)、U_i.y(0)通过派克变化到同步坐标系d、q轴，忽略逆变器功率损耗，依据逆变器两侧功率守恒原则，求得逆变器输入功率，进而求取机侧输入功率和稳压电容电压初值，在后续过程中保持不变，有：

P_s＝P_i(0)＝U_i.d(0)I_L.d(0)+U_i.q(0)I_L.q(0)    (20)

I_dc(0)＝S_d(0)I_L.d(0)+S_q(0)I_L.q(0)           (21)

U_dc(0)＝P_s/I_dc(0)                           (22)

当k＞0，即非初值时刻，有：

$\left\{\begin{array}{c}\mathrm{\Δ}{U}_{\mathrm{dc}}(k-1)=U_{\mathrm{dc}}(k-1)-U_{\mathrm{dc}}^{*}\\ \mathrm{\Δ}{I}_d(k-1)=I_d^{*}(k-1)-I_{L.d}(k-1)\\ \mathrm{\Δ}{I}_q(k-1)=I_q^{*}(k-1)-I_{L.q}(k-1)\\ I_d^{*}\left(k\right)=I_d^{*}(k-1)+K_{p1}\mathrm{\Δ}{U}_{\mathrm{dc}}(k-1)+K_{i1}{\∫}_{t_0}^t\mathrm{\Δ}{U}_{\mathrm{dc}}(k-1)\mathrm{dt}\\ S_d\left(k\right)=S_d(k-1)+K_{p2}\mathrm{\Δ}{I}_d(k-1)+K_{i2}{\∫}_{t_0}^t\mathrm{\Δ}{I}_d(k-1)\mathrm{dt}\\ S_q\left(k\right)=S_q(k-1)+K_{p2}\mathrm{\Δ}{I}_q(k-1)+K_{i2}{\∫}_{t_0}^t\mathrm{\Δ}{I}_q(k-1)\mathrm{dt}\end{array}\right.---\left(\text{23}\right)$

$\left\{\begin{array}{c}{U}_{i.d}\left(k\right)=S_d\left(k\right)U_{\mathrm{dc}}(k-1)\\ U_{i.q}\left(k\right)=S_q\left(k\right)U_{\mathrm{dc}}(k-1)\end{array}\right.---\left(24\right)$

通过采用四阶龙格库塔法求解下面的微分方程组，求得k时刻U_dc、I_L.d、I_L.q相应的值U_dc(k)、I_L.d(k)、I_L.q(k)；

$\frac{{\mathrm{dU}}_{\mathrm{dc}}}{\mathrm{dt}}=\frac{1}{C}(\frac{P_s}{U_{\mathrm{dc}}}-S_d{I}_{L.d}-S_q{I}_{L.q})---\left(25\right)$

$\left\{\begin{array}{c}\frac{{\mathrm{dI}}_{L.d}}{\mathrm{dt}}=\frac{1}{L}(U_{i.d}-U_{g.d}-I_{L.d}R)-\mathrm{\ω}{I}_{L.q}\\ \frac{{\mathrm{dI}}_{L.q}}{\mathrm{dt}}=\frac{1}{L}(U_{i.q}-U_{g.q}-I_{L.q}R)+\mathrm{\ω}{I}_{L.d}\end{array}\right.---\left(26\right)$

步骤5：将步骤4中求出d、q坐标下的电流I_L.d(k)、I_L.q(k)进行派克反变换，得到响应序列P_m(k)，Q_m(k)，k＝0，1，2，…，L-1，L为实测样本长度；

$\left\{\begin{array}{c}{P}_m{\left(k\right)=U}_{g.x}\left(k\right)I_{L.x}\left(k\right)+U_{g.y}{\left(k\right)I}_{L.y}\left(k\right)\\ Q_m\left(k\right)=U_{g.y}\left(k\right)I_{L.x}\left(k\right)-U_{g.x}\left(k\right)I_{L.y}\left(k\right)\end{array}\right.---\left(27\right)$

步骤6：计算单个序列误差jks。

jks＝(P-P_m)²+(Q-Q_m)²    (28)

步骤7：判断目标函数值即累计误差minJ(x)是否满足终止条件，终止条件为累计误差在设定的范围之内，满足终止条件则转步骤8，不满足终止条件则利用遗传算法进行i+1次遗传寻优迭代，求解最优化问题minJ→m(i)，转步骤2；

$\min J\left(x\right)=\min \underset{k=0}{\overset{L-1}{\mathrm{\Σ}}}[{(P\left(k\right)-P_m\left(k\right))}^T\·(P\left(k\right)-P_m\left(k\right))+{(Q\left(k\right)-Q_m\left(k\right))}^T\·(Q\left(k\right)-Q_m\left(k\right))]---\left(29\right)$

步骤8：输出最终的辨识参数即等效模型参数。

有益效果：

本发明提出了一种微型燃气轮机发电系统的模拟方法，能满足电网计算要求的MT发电系统(即微型燃气轮机发电系统)暂态仿真。

现实意义在于微型燃气轮机发电系统的并网容量不断增加，影响了配电网侧综合负荷特性，使得传统的负荷模型描述能力下降，将本发明的微型燃气轮机发电系统的模拟方法与电力系统仿真分析软件对接，可以研究MT发电系统对电力系统综合负荷特性的影响。

本发明提出的微型燃气轮机发电系统的模拟方法填补了满足电网计算要求的MT发电系统暂态仿真模拟方法的空白。

本发明不仅可以描述MT发电系统的稳态特性而且可以描述其暂态特性，对MT发电系统具有普遍的实用性。

仿真表明采用本方法构建的MT等效模型与Matlab/Simuink中MT仿真系统的曲线是吻合的，能准确模拟微型燃气轮机发电系统的外特性。

实验表明，采用本发明建模方法构造的模型具有较好的描述能力、泛化能力、模型参数稳定性。

附图说明

图1为本发明中使用的MT整流逆变系统的控制原理图。

图2为本发明中使用的4节点系统单线图。图2是在Matlab/Simulink工具箱里搭建的仿真系统，本发明的建模样本是在这个系统里的PCC处测量的。

图3为本发明中并网的微型燃气轮机发电系统等值电路图。

图4为本发明中等效描述模型拟合效果；其中，(a)电压跌落20％时有功响应；(b)电压跌落20％时无功响应；

图5为本发明中内插外推能力验证的效果图。其中，(a)电压跌落10％时有功响应；(b)电压跌落10％时无功响应；(c)电压跌落30％时有功响应；(d)电压跌落30％时无功响应。

图6为总体测辨负荷建模实现过程的流程图。

图7为仿真模型和本发明的构造的模型在电压跌落时响应结果示意图。图a为电压跌落20％时电压激励的曲线，图b、c分别为此情况下的有功响应和无功响应曲线。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明：

实施例1：

一种微型燃气轮机发电系统的模拟方法，步骤如下：

(1)研究MT系统的发电原理，本发明以单轴微型燃气轮机发电的数学模型为代表，在MATLAB/Simulink中搭建MT发电系统的数学模型，研究MT发电系统的特性。

对于MT发电系统的数学模型，国内外已有相关的研究，本发明只是引用IEEE已收录论文中的公式搭建的模块，参见以下参考文献：Saha A K，Chowdhury S，Chowdhury SP，et al.Modeling and performance analysis of amicroturbine as a distributed energy resource[J].IEEE Transactions on EnergyConversion，2009，24(2)：529-538.

基于上述MT仿真系统的数学模型，采用MATLAB/Simulink工具进行仿真，为本领域的普通技术人员的常用技术手段。

选择适合的整流逆变控制策略，在MATLAB/Simulink中搭建MT发电系统，建立MT发电系统的详细数字仿真模型。

(2)通过对MT发电系统工作机理分析以及发电系统运行特性仿真分析，提出能准确模拟MT发电系统外特性并满足电网计算要求的MT发电系统机电暂态等效仿真模型，该模型的数学模型以PCC点电压和电网频率为激励，以MT发电系统向电网注入的d、q轴电流、及直流母线电压为状态变量的三阶微分-代数方程组；其物理模型如图3所示。

(3)在不同的扰动下，对MT发电系统进行一系列的数字仿真，通过在MT发电系统并网处设置不同程度的短路故障，获取大量的建模样本数据，包括MT发电系统上网点的电压U、有功P、无功Q数据，采样时间间隔为2ms，采用综合改进的遗传算法对模型进行辨识参数建模，获得等效模型参数，包括k_p1、k_i1、k_p2、k_i2、R、L、C。其等效电路如图3所示，图2所示电路是模拟微型燃气轮机发电系统接入的配电网络，事实上，其他配电网络都可以简化成相似的等效网络，只是网络参数(线路阻抗、负荷等)和复杂程度不同而已，且采集的建模数据是源自微型燃气轮机发电系统的并网母线的数据，因此，在其他网络结构及参数下，依然能按照本发明模拟微型燃气轮机的运行特性。图7为电压跌落20％条件下的建模样本及辨识结果实例。

图3为图1的等效电路，图2作为数据获取的外部条件，只是将我所构建的模型放在一个配网系统中，所得模型与网络的具体结构并无特殊关系，也即即使不是这个网络结构，只要是图1所示的微型燃气轮机发电系统结构是一致的，就能按此方法得到所需的模型。

(4)对模型的自描述能力(内插、外推)进行检验，检验结果如图5所示，并研究了MT发电系统等效模型在不同负荷水平条件下适用性，表明模型具有较好的描述能力和泛化能力。

图1中，I_L是逆变器的输出电流，即MT发电系统由公共连接点(Point ofCommon Coupling，PCC)注入电网的电流，U_i是逆变器的出口电压，U_g是PCC点处电压，R、L分别为等值电阻、等值电感。

根据图1构造的微型燃气轮机发电的等值电路图，令状态向量x＝[I_L.x I_L.y U_dc]^T；模型参数相量θ＝[R L C k_p1 K_i1 k_p2 k_i2]^T。输入相量u＝[U_x U_y]^T；输出相量y＝[P Q]^T，则微型燃气轮机发电系统等效模型写成状态方程的一般形式如式(1)所示。

$\left\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=f(x,\mathrm{\θ},u)\\ y=h(x,\mathrm{\θ},u)\end{array}\right.---\left(1\right)$

式(1)中第一式是MT发电系统的动态微分方程，第二式是MT发电系统的输出方程。动态微分方程的具体形式如式(2)所示。

$\frac{d{\stackrel{\·}{I}}_L}{\mathrm{dt}}=f(x,\mathrm{\θu})=\frac{1}{L}({\stackrel{\·}{U}}_i-{\stackrel{\·}{U}}_g-{\stackrel{\·}{I}}_{L}R)---\left(2\right)$

其中：

${\stackrel{\·}{U}}_i=\stackrel{\·}{S}{U}_{\mathrm{dc}}---\left(3\right)$

U_dc是MT发电系统整流器的输出电压，

为逆变器开关函数。

MT发电系统的入网电流和PCC电压用电网同步坐标下的相量表示，分别如式(4)、(5)所示。

${\stackrel{\·}{I}}_L=I_{L.x}+j{I}_{L.y}---\left(4\right)$

${\stackrel{\·}{U}}_g=U_{g.x}+j{U}_{g.y}---\left(5\right)$

取PCC电压为参考相量，即PCC点的电压初值为

同时已知稳态时的激励u₀＝[U_g.x(0)U_g.y(0)]^T和响应y₀＝[P(0)Q(0)]^T，逆变器的出口电压表达如式(6)。

${\stackrel{\·}{U}}_i\left(0\right)=U_{i.x}\left(0\right)+j{U}_{i.y}\left(0\right)---\left(6\right)$

利用电压降落公式，可求得

如式(7)、(8)所示。

$U_{i.x}\left(0\right)=U_{g.x}\left(0\right)+\frac{P\left(0\right)R+Q\left(0\right)\mathrm{\ωL}}{U_g\left(0\right)}---\left(7\right)$

$U_{i.y}\left(0\right)=U_{g.y}\left(0\right)+\frac{P\left(0\right)R-Q\left(0\right)\mathrm{\ωL}}{U_g\left(0\right)}---\left(8\right)$

进而可以进一步求出电路中的电流，如式(9)、(10)、(11)所示。

${\stackrel{\·}{I}}_L\left(0\right)=I_{L.x}\left(0\right)+j{I}_{L.y}\left(0\right)=\frac{1}{R+\mathrm{j\ωL}}({\stackrel{\·}{U}}_i\left(0\right)-{\stackrel{\·}{U}}_g\left(0\right))---\left(9\right)$

$I_{L.x}\left(0\right)=\frac{(P\left(0\right)R+Q\left(0\right)\mathrm{\ωL})R}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}+\frac{(P\left(0\right)\mathrm{\ωL}-Q\left(0\right)R)\mathrm{\ωL}}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}---\left(10\right)$

$I_{L.y}\left(0\right)=\frac{(P\left(0\right)\mathrm{\ωL}-Q\left(0\right)R)R}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}-\frac{(P\left(0\right)R+Q\left(0\right)\mathrm{\ωL})\mathrm{\ωL}}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}---\left(11\right)$

其中，P(0)、Q(0)即MT发电系统上网的有功和无功的初始值。

把已知的初值u₀＝[U_g.x(0)U_g.y(0)]^T、y₀＝[P(0)Q(0)]^T和求得的初值U_i.x(0)U_i.y(0)、I_L.x(0)、I_L.y(0)进行派克变换，得到U_i、I_L的d、q轴初始分量，忽略逆变器功率损耗，依据逆变器两侧功率守恒原则，通过求取逆变器出口侧功率，可以求得逆变器输入功率，进而可以求取机侧输入功率和稳压电容电压，由于整流逆变器的隔离作用，微型燃气轮机原动机能在一定程度上保持原有运行状态，故在后续计算过程中保持机侧输入功率恒定。

P_s＝P_i(0)＝U_i.d(0)I_L.d(0)+U_i.q(0)I_L.q(0)    (12)

I_dc(0)＝S_d(0)I_L.d(0)+S_q(0)I_L.q(0)           (13)

U_dc(0)＝P_s/I_dc(0)                           (14)

式中S_d(0)、S_q(0)值是给定的，在后续求解过程中依据式(15)不断动态修正。

$\left\{\begin{array}{c}\mathrm{\Δ}{U}_{\mathrm{dc}}=U_{\mathrm{dc}}\left(t\right)-U_{\mathrm{dc}}^{*}\\ \mathrm{\Δ}{I}_d=I_q^{*}-I_{L.d}\left(t\right)\\ \mathrm{\Δ}{I}_q=I_q^{*}-I_{L.q}\left(t\right)\\ I_d^{*}\left(t\right)=I_d^{*}+K_{p1}\mathrm{\Δ}{U}_{\mathrm{dc}}+K_{i1}{\∫}_{t_0}^t\mathrm{\Δ}{U}_{\mathrm{dc}}\mathrm{dt}\\ S_d\left(t\right)=S_d+K_{p2}\mathrm{\Δ}{I}_d+K_{i2}{\∫}_{t_0}^t\mathrm{\Δ}{I}_d\mathrm{dt}\\ S_q\left(t\right)=S_q+K_{p2}\mathrm{\Δ}{I}_q+K_{i2}{\∫}_{t_0}^t\mathrm{\Δ}{I}_q\mathrm{dt}\end{array}\right.---\left(15\right)$

至此，把已知的初值u₀＝[U_g.x(0)U_g.y(0)]^T、y₀＝[P(0)Q(0)]^T和求得的初值U_i.x(0)U_i.y(0)、I_L.x(0)、I_L.y(0)进行派克变换，得到所有初值的d、q轴分量，作为求解微分方程的初始值，最后用四阶龙格库塔法解微分方程(16)：

$\left\{\begin{array}{c}\frac{{\mathrm{dI}}_{\mathrm{dc}}}{\mathrm{dt}}=\frac{1}{C}\left(\frac{P_s}{U_{\mathrm{dc}}}-S_d{I}_{L.d}-S_q{I}_{L.q}\right)\\ \frac{{\mathrm{dI}}_{L.d}}{\mathrm{dt}}=\frac{1}{L}(U_{i.d}-U_{g.d}-I_{L.d}R)-\mathrm{\ω}{I}_{L.q}\\ \frac{{\mathrm{dU}}_{L.q}}{\mathrm{dt}}=\frac{1}{L}(U_{i.q}-U_{g.q}-I_{L.q}R)+\mathrm{\ω}{I}_{L.d}\end{array}\right.---\left(16\right)$

从而可以求得任意时刻的电流的d、q轴分量I_L.d、I_L.q及U_dc，再对I_L.d、I_L.q进行派克反变换，即可求出同步坐标下的电流I_L.x、I_L.y；进而求出本发明所构造的模型的有功响应和无功响应：

$\left\{\begin{array}{c}P=U_{g.x}{I}_{L.x}+U_{g.y}{I}_{L.y}\\ Q=U_{g.y}{I}_{L.x}-U_{g.x}{I}_{L.y}\end{array}\right.---\left(17\right)$

辨识步骤如下，框图见附图6：

步骤1：给定实测样本u(k)，P(k)，Q(k)(k＝0，1，2，…，L-1)，对遗传寻优迭代次数i进行初始化，令i＝0；

步骤2：随机产生独立待辨识参数样本R、L、C、k_p1、k_i1、k_p2、k_i2；

步骤3：设定遗传寻优迭代次数i＝i+1；

步骤4：当k＝0时，即为初值求解，有：将实测样本代入下列公式求变量初值I_L.x(k)、I_L.y(k)、U_dc(k)：

$U_{i.x}\left(0\right)=U_{g.x}\left(0\right)+\frac{P\left(0\right)R+Q\left(0\right)\mathrm{\ωL}}{U_g\left(0\right)}---\left(15\right)$

$U_{i.y}\left(0\right)=U_{g.y}\left(0\right)+\frac{P\left(0\right)\mathrm{\ωL}-Q\left(0\right)R}{U_g\left(0\right)}---\left(16\right)$

$U_i\left(0\right)=\sqrt{U_{i.x}{\left(0\right)}^2+U_{i.y}{\left(0\right)}^2}---\left(17\right)$

$I_{L.x}\left(0\right)=\frac{(P\left(0\right)R+Q\left(0\right)\mathrm{\ωL})R}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}+\frac{(P\left(0\right)\mathrm{\ωL}-Q\left(0\right)R)\mathrm{\ωL}}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}---\left(18\right)$

$I_{L.y}\left(0\right)=\frac{(P\left(0\right)\mathrm{\ωL}-Q\left(0\right)R)R}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}-\frac{(P\left(0\right)R+Q\left(0\right)\mathrm{\ωL})\mathrm{\ωL}}{[R^2+{\left(\mathrm{\ωL}\right)}^2]U_g\left(0\right)}---\left(19\right)$

将上述状态变量值I_L.x(0)、I_L.y(0)、U_g.x(0)、U_g.y(0)、U_i.x(0)、U_i.y(0)通过派克变化到同步坐标系d、q轴，忽略逆变器功率损耗，依据逆变器两侧功率守恒原则，求得逆变器输入功率，进而求取机侧输入功率和稳压电容电压初值，在后续过程中保持不变，有：

P_s＝P_i(0)＝U_i.d(0)I_L.d(0)+U_i.q(0)I_L.q(0)    (20)

I_dc(0)＝S_d(0)I_L.d(0)+S_q(0)I_L.q(0)           (21)

U_dc(0)＝P_s/I_dc(0)                           (22)

当k＞0(即非初值时刻)，有：

$\left\{\begin{array}{c}\mathrm{\Δ}{U}_{\mathrm{dc}}(k-1)=U_{\mathrm{dc}}(k-1)-U_{\mathrm{dc}}^{*}\\ \mathrm{\Δ}{I}_d(k-1)=I_d^{*}(k-1)-I_{L.d}(k-1)\\ \mathrm{\Δ}{I}_q(k-1)=I_q^{*}(k-1)-I_{L.q}(k-1)\\ I_d^{*}\left(k\right)=I_d^{*}(k-1)+K_{p1}\mathrm{\Δ}{U}_{\mathrm{dc}}(k-1)+K_{i1}{\∫}_{t_0}^t\mathrm{\Δ}{U}_{\mathrm{dc}}(k-1)\mathrm{dt}\\ S_d\left(k\right)=S_d(k-1)+K_{p2}\mathrm{\Δ}{I}_d(k-1)+K_{i2}{\∫}_{t_0}^t\mathrm{\Δ}{I}_d(k-1)\mathrm{dt}\\ S_q\left(k\right)=S_q(k-1)+K_{p2}\mathrm{\Δ}{I}_q(k-1)+K_{i2}{\∫}_{t_0}^t\mathrm{\Δ}{I}_q(k-1)\mathrm{dt}\end{array}\right.---\left(\text{23}\right)$

$\left\{\begin{array}{c}{U}_{i.d}\left(k\right)=S_d\left(k\right)U_{\mathrm{dc}}(k-1)\\ U_{i.q}\left(k\right)=S_q\left(k\right)U_{\mathrm{dc}}(k-1)\end{array}\right.---\left(24\right)$

通过采用四阶龙格库塔法求解式(25)(26)所示的微分方程组，可以求得k时刻U_dc、I_L.d、I_L.q相应的值U_dc(k)、I_L.d(k)、I_L.q(k)；

$\frac{{\mathrm{dU}}_{\mathrm{dc}}}{\mathrm{dt}}=\frac{1}{C}(\frac{P_s}{U_{\mathrm{dc}}}-S_d{I}_{L.d}-S_q{I}_{L.q})---\left(25\right)$

$\left\{\begin{array}{c}\frac{{\mathrm{dI}}_{L.d}}{\mathrm{dt}}=\frac{1}{L}(U_{i.d}-U_{g.d}-I_{L.d}R)-\mathrm{\ω}{I}_{L.q}\\ \frac{{\mathrm{dI}}_{L.q}}{\mathrm{dt}}=\frac{1}{L}(U_{i.q}-U_{g.q}-I_{L.q}R)+\mathrm{\ω}{I}_{L.d}\end{array}\right.---\left(26\right)$

步骤5：将步骤4中求出d、q坐标下的电流I_L.d(k)、I_L.q(k)进行派克反变换，得到响应序列P_m(k)，Q_m(k)，k＝0，1，2，…，L-1，L为实测样本长度；

$\left\{\begin{array}{c}{P}_m{\left(k\right)=U}_{g.x}\left(k\right)I_{L.x}\left(k\right)+U_{g.y}{\left(k\right)I}_{L.y}\left(k\right)\\ Q_m\left(k\right)=U_{g.y}\left(k\right)I_{L.x}\left(k\right)-U_{g.x}\left(k\right)I_{L.y}\left(k\right)\end{array}\right.---\left(27\right)$

步骤6：计算单个序列误差jks。

jks＝(P-P_m)²+(Q-Q_m)²    (28)

步骤7：判断目标函数值即累计误差minJ(x)是否满足终止条件，终止条件为累计误差在设定的范围之内，满足终止条件则转步骤8，不满足终止条件则进行i+1次遗传寻优迭代，求解最优化问题minJ→m(i)，转步骤2；

$\min J\left(x\right)=\min \underset{k=0}{\overset{L-1}{\mathrm{\Σ}}}[{(P\left(k\right)-P_m\left(k\right))}^T\·(P\left(k\right)-P_m\left(k\right))+{(Q\left(k\right)-Q_m\left(k\right))}^T\·(Q\left(k\right)-Q_m\left(k\right))]---\left(29\right)$

求解模型响应序列P_m(k)，Q_m(k)，k＝0，1，2，…，L-1的平均误差Er(数据长度不一致时，每个序列的误差不一致，就应使用平均误差的概念)，用于后面的模型检验。

$E_r=\underset{k=0}{\overset{L-1}{\mathrm{\Σ}}}\sqrt{{(P\left(k\right)-P_m\left(k\right))}^T\·(P\left(k\right)-P_m\left(k\right))+{(Q\left(k\right)-Q_m\left(k\right))}^T\·(Q\left(k\right)-Q_m\left(k\right))}/L---\left(30\right)$

步骤8：输出电压激励和模型响应P_m(k)，Q_m(k)，完成。为验证上面建立的动态模型对MT发电系统等效描述模型的有效性，以图2所示系统为仿真实验对象，通过设置三相短路故障，使得B3节点电压跌落在10％～35％之间，测得PCC处母线电压和MT发电系统注入PCC母线的有功功率和无功功率功率共6组数据样本，分别作为建模激励和实测响应，对MT发电系统进行辨识建模。模型结构如式(26)、(27)所示，参数辨识采用综合改进的遗传算法。限于篇幅，给出母线B3电压跌落20％时的仿真响应及其对应的模型响应曲线如图4所示。

图4所示结果表明，暂态过程中，模型响应均能较好地拟合仿真实验数据，说明本发明提出的FC模型对数据样本的逼近效果较好，具有较强的自描述能力。因此，通过上面的分析充分说明，采用所建立的三阶微分方程描述的等效模型能较好地描述FC发电系统的特性。

将B3母线电压跌落分别为10％～30％时的电压激励依次施加于20％电压跌落时辨识所得模型，比较相应的模型响应对仿真实测响应的拟合程度。图5为10％内插和30％外推的响应曲线。检验结果表明，虽然拟合样本与建模样本的电压激励幅度相差较大，但模型具有良好的内插和外推特性，因而具有良好的泛化能力。

以电压跌落20％时的数据为例，结果见图7，图中的实测数据指从Simuink中得到的数据，通过前述的辨识方法，辨识参数值如下R＝0.02135(pu)L＝0.09310(pu)C＝0.00515(pu)k_p1＝0.0205 k_i1＝2.8899 k_p2＝0.2015 k_i2＝8.546，残差为0.00084。图7直观表明采用本方法构建的模型与Simuink中仿真的曲线能够较好的拟合，能准确模拟微型燃气轮机发电系统的响应过程。

Claims (4)

1.一种微型燃气轮机发电系统的模拟方法，其特征在于，在该微型燃气轮机发电系统中，微型燃气轮机发电系统通过整流逆变装置向电网输出电能，将微型燃气轮机发电系统用以下以电感电流、电容电压为状态变量的三阶动态微分方程来表征： 

其中，I_L.d、I_L.q是逆变器输出电流的d轴分量和q轴分量，PCC为公共连接点，逆变器的输出电流也即微型燃气轮机发电系统由PCC点注入电网的电流；U_i.d、U_i.q是逆变器的出口电压的d轴和q轴分量；U_g.d、U_g.q是微型燃气轮机发电系统与电网的公共连接点电压的d轴分量和q轴分量；R、L、C分别为等值电阻、等值电感和等值电容器；P_S为微型燃气轮机原动机输入功率；U_dc为直流母线电压；S_d、S_q为逆变器同步坐标系下的平均开关函数，通过以下修正公式来不断修正： 

其中，k_p1、k_i1、k_p2、k_i2为控制环节参数；ΔU_dc、ΔI_d、ΔI_q分别为相应变量在单位时刻的变化量；I_L.d(t)、I_L.q(t)、S_d(t)、S_q(t)、

分别为I_L.d、I_L.q、S_d、S_q、 

在t时刻的值；

分别为U_dc、I_L.d、I_L.q、I_L.d(t)的参考值，解出电流I_L.d、I_L.q的初值之后，通过给定S_d=0.87、S_q=0，通过下式求得U_dc的初值，将U_dc的初值赋予

即可获取

U_dc(t)为U_dc在t时刻的值； 

在后续辨识中保持不变，

取值为0以使得逆变器传输无功为0，通过初值求解求得I_L.d的初值，将I_L.d的初值作为

的初始值，在后续过程中动态修正； 

求解前述的三阶动态微分方程求得任意时刻的电流的d、q轴分量I_L.d和I_L.q，再对I_L.d和I_L.q进行派克反变换，即可求出同步坐标下的电流I_L.x、I_L.y，进而求出微型燃气轮机发电系统的有功响应和无功响应： 

其中，U_g.x和U_g.y分别为同步坐标系下的公共连接点处的x电压分量和y电压分量。 

2.根据权利要求1所述的微型燃气轮机发电系统的模拟方法，其特征在于，获取初始值并求解微分方程的步骤为： 

所述的初始值包括已知的初值u₀=[U_g.x(0)U_g.y(0)]^T、y₀=[P(0)Q(0)]^T和求得的初值U_i.x(0)、U_i.y(0)、I_L.x(0)、I_L.y(0)； 

其中：

其中，P(0)、Q(0)即微型燃气轮机发电系统上网的有功和无功的初始值；U_g(0)是微型燃气轮机发电系统与电网的公共连接点电压的初始值； 

求解微分方程的过程为：将上述初始值进行派克变换，得到所有初值的d、q轴分量，作为求解微分方程的初始值，最后用四阶龙格库塔法解微分方程求解所述的三阶动态微分方程。 

3.根据权利要求1或2所述的微型燃气轮机发电系统的模拟方法，其特征在于， 参数辨识方法如下： 

通过仿真或实测在微型燃气轮机发电系统并网处设置不同程度的短路故障，获取建模样本数据，包括微型燃气轮机发电系统上网点的电压U、有功功率P、无功功率Q数据； 

采用遗传算法进行等效模型参数辨识，得到最终的等效模型—即由三阶动态微分方程和修正公式表征的微型燃气轮机发电系统的等效模型； 

等效模型参数包括k_p1、k_i1、k_p2、k_i2、R、L、C。 

4.根据权利要求3所述的微型燃气轮机发电系统的模拟方法，其特征在于，具体的辨识步骤如下： 

步骤1：给定实测样本u(k),P(k),Q(k)(k=0,2,…,L-1)，对遗传寻优迭代次数i进行初始化，令i=0； 

步骤2：随机产生独立待辨识参数样本R、L、C、k_p1、k_i1、k_p2、k_i2； 

步骤3：设定遗传寻优迭代次数i=i+1； 

步骤4：当k=0时，即为初值求解，有：将实测样本代入下列公式求变量初值I_L.x(k)、I_L.y(k)、U_dc(k)： 

将上述状态变量值I_L.x(0)、I_L.y(0)、U_g.x(0)、U_g.y(0)、U_i.x(0)、U_i.y(0)通过派克变化到同步坐标系d、q轴，忽略逆变器功率损耗，依据逆变器两侧功率守恒原则，求得逆变器输入功率，进而求取机侧输入功率和稳压电容电压初值，在后续过程中保持不变，有： 

             P_s=P_i(0)=U_i.d(0)I_L.d(0)+U_i.q(0)I_L.q(0) 

             I_dc(0)=S_d(0)I_L.d(0)+S_q(0)I_L.q(0) 

               U_dc(0)=P_s/I_dc(0) 

当k>0，即非初值时刻，有： 

通过采用四阶龙格库塔法求解下面的微分方程组，求得k时刻U_dc、I_L.d、I_L.q相应的值U_dc(k)、I_L.d(k)、I_L.q(k)； 

步骤5：将步骤4中求出d、q坐标下的电流I_L.d(k)、I_L.q(k)进行派克反变换，得到响应序列P_m(k),Q_m(k),k=0,1,2,…,L-1，L为实测样本长度； 

步骤6：计算单个序列误差jks： 

          jks=(P-P_m)²+(Q-Q_m)²

步骤7：判断目标函数值即累计误差minJ(x)是否满足终止条件，终止条件为累计误差在设定的范围之内，满足终止条件则转步骤8，不满足终止条件则利用遗传算法进行i+1次遗传寻优迭代，求解最优化问题minJ→m(i)，转步骤2； 

步骤8：输出最终的辨识参数即等效模型参数。 

Images