CN104158191A - 一种多机电力系统稳定运行的分散协调控制方法 - Google Patents

Abstract

本发明公开了一种多机电力系统稳定运行的分散协调控制方法，以IEEE三机九节点系统模型为例，对多机电力系统进行数学建模并对模型进行分析与研究，将多机电力系统转子运动方程偏差化得到标准的线性状态变量方程，在线性化的基础上，通过直接迭代法求解得到的Levine-Athans方程组得到最优分散协调控制器，设计出能适应静态稳定情况的最优分散协调控制规律，能够保证电力系统的稳定运行，提高电网稳定性，对多机电力系统控制规律的研究与验证对提高电力系统稳定运行具有理论价值和实际意义。

Description

一种多机电力系统稳定运行的分散协调控制方法

技术领域

本发明属于多机电力系统控制技术领域，涉及一种对多机电力系统的控制方法，具体是指一种对多机电力系统有效控制的最优分散协调控制方法。

背景技术

多机组的互联运行带来显著的经济与社会效益的同时，对电网的稳定性运行也带来了新的挑战。由于电网的互联，各机组间存在着一定的耦合关系，单一的故障扰动不仅会影响到本系统机组的正常运行，也会牵扯到其他机组的运行状态，从而可能使整个电网的潮流、电压以及频率发生变化，导致自动保护装置出现连锁反应引发大面积停电的事故。由此，多机组电力系统的互联的实现解决各地区能源负荷分布不均，减少系统装机容量与备用容量等问题的同时，也带来了一系列的电网稳定性的问题。

在电网互联形成的多机电力系统中，各台发电机之间存在着相互的联系与影响。对多机电力系统中的发电机进行完整建模后发现，每台发电机的数学模型都是由一组动力学方程与一组非动力学耦合方程组成，一台发电机运行参数的变化势必影响其它互联的发电机的运行状态，从而影响到整个电网的稳定运行。

发明内容

为提高电网稳定性，本发明提供了一种多机电力系统稳定运行的分散协调控制方法。

本发明所采用的技术方案是：一种多机电力系统稳定运行的分散协调控制方法，其特征在于，包括以下步骤：

步骤1：对多机电力系统被控对象进行数学建模；其具体实现包括以下子步骤：

步骤1.1：首先从力学的角度推导同步发电机转子运动方程，即经典的摇摆方程；

步骤1.2：接着对同步电机转子运动方程中涉及的发电机输出功率进行推导与求解；

步骤1.3：最后从电路的角度对同步发电机的励磁绕组电磁动态方程进行分析与求解，得到多机电力系统第i台同步发电机的数学模型表达式为：

$\left\{\begin{array}{c}\frac{d{\mathrm{\δ}}_i}{\mathrm{dt}}=({\mathrm{\ω}}_i-1){\mathrm{\ω}}_{0i}\\ \frac{d{\mathrm{\ω}}_i}{\mathrm{dt}}=\frac{1}{T_{\mathrm{Ji}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}})-\frac{D_i}{T_{\mathrm{Ji}}}({\mathrm{\ω}}_i-1)\\ \frac{d{E^{\′}}_{\mathrm{qi}}}{\mathrm{dt}}=\frac{1}{{T^{\′}}_{d0i}}[E_{\mathrm{fi}}-E_{\mathrm{qi}}^{\′}-(x_{\mathrm{di}}-x_{\mathrm{di}}^{\′})i_{\mathrm{di}}]\end{array}\right.$

多机电力系统： $\begin{array}{c}{P}_{\mathrm{ei}}=E_{\mathrm{qi}}^2{G}_{\mathrm{ii}}+E_{\mathrm{qi}}{{\mathrm{\Σ}}_{\underset{j\≠i}{j=1}}^n{E}_{\mathrm{qj}}(G_{\mathrm{ij}}\cos {\mathrm{\δ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\δ}}_{\mathrm{ij}})}_{\begin{array}{c}\\ \end{array}}i=\{\mathrm{1,2},\·\·\·n\};\end{array}$ 式中，δ_i为第i台发电机转子运行角，即q轴与同步参考轴S间的夹角，单位rad，ω_i表示转子运动角速度(标幺值)，ω_0i为转子运动角速度的初始值，单位rad/s，T_Ji为转子转动惯量时间常数，单位s，P_mi表示第i台发电机输入的机械功率(标幺值)，P_ei表示第i台发电机输出的电磁功率(标幺值)，D_i表示阻尼系数(标幺值)，E’_qi表示第i台发电机暂态电动势(标幺值)，E’_fi表示稳态磁链在d轴定子侧产生的电势(标幺值)，i_di表示第i台发电机输电流瞬态值(标幺值)，x_di、x’_di分别为第i台发电机直轴同步电抗与瞬态电抗(标幺值)；且I＝YU，Y＝G+jB，为电力系统简化导纳矩阵，G_ii＝Y_iisinα_ii和B_ii＝Y_iicosα_ii为第i节点的自电导，G_ij＝Y_ijsinα_ii和B_ij＝Y_ijcosα_ii为第i节点与第j节点间的互电导及互导纳；

步骤2：对多机电力系统数学模型进行非线性系统线性化，其具体实现包括以下子步骤：

步骤2.1：由于对多机电力系统的最优分散协调控制设计要求系统为线性系统，因此为建立多机系统线性化的状态空间方程，首先对其转子运动方程进行线性化，即：

$\left\{\begin{array}{c}\mathrm{\Δ}{\stackrel{\·}{\mathrm{\δ}}}_i=\mathrm{\Δ}{\mathrm{\ω}}_i\\ \mathrm{\Δ}{\stackrel{\·}{\mathrm{\ω}}}_i=-\frac{{\mathrm{\ω}}_{0i}}{T_{\mathrm{Ji}}}\mathrm{\Δ}{P}_{\mathrm{ei}}-\frac{D_i}{T_{\mathrm{Ji}}}\mathrm{\Δ}{\mathrm{\ω}}_i+\frac{{\mathrm{\ω}}_{0i}}{T_{\mathrm{Ji}}}\mathrm{\Δ}{P}_{\mathrm{mi}}\end{array}\right.---\left(1\right)$

其中，其中，δ_i为第i台发电机转子运行角，单位rad，ω_i为发电机转子运动角速度(标幺值)，ω_0i为第i台发电机转子运动角速度的初始值，单位rad/s，T_Ji为第i台发电机转子转动惯量时间常数，单位s，P_mi表示第i台发电机输入的机械功率(标幺值)，P_ei表示第i台发电机输出的电磁功率(标幺值)，D_i为第i台发电机阻尼系数(标幺值)；

步骤2.2：对发电机输出的电磁功率进行线性化；

发电机输出的电磁功率

$\left[\begin{array}{cc}{P}_{\mathrm{ei}}=E_{\mathrm{qi}}^2{G}_{\mathrm{ii}}+E_{\mathrm{qi}}\underset{j\≠i}{\overset{n}{\underset{j=1}{\mathrm{\Σ}}}}{E}_{\mathrm{qj}}(G_{\mathrm{ij}}\cos {\mathrm{\δ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\δ}}_{\mathrm{ij}})& i=\mathrm{1,2,3}\end{array}\right]$

其中G_ij＝Y_ijcosψ_ij表示连接线间的电导，且实际中阻抗角ψ_ij接近90°，因此式

中电导G_ij可忽略不计，即

$\begin{array}{cc}{P}_{\mathrm{ei}}=E_{\mathrm{qi}}^2{G}_{\mathrm{ii}}+E_{\mathrm{qi}}\underset{j\≠i}{\overset{n}{\underset{j=1}{\mathrm{\Σ}}}}{E}_{\mathrm{qj}}{B}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)& i=\mathrm{1,2},\·\·\·n\end{array}$

在操作点处进行线性化，用与分别替代P_ei、δ_i、δ_j，则

$\begin{array}{c}{P}_{\mathrm{ei}}^0+\mathrm{\Δ}{P}_{\mathrm{ei}}=E_{\mathrm{qi}}^2{G}_{\mathrm{ii}}+E_{\mathrm{qi}}\underset{j\≠i}{\overset{3}{\underset{j=1}{\mathrm{\Σ}}}}{E}_{\mathrm{qj}}{B}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i^0+\mathrm{\Δ}{\mathrm{\δ}}_i-{\mathrm{\δ}}_j^0-\mathrm{\Δ}{\mathrm{\δ}}_j)\\ =E_{\mathrm{qi}}^2{G}_{\mathrm{ii}}+E_{\mathrm{qi}}\underset{j\≠i}{\overset{3}{\underset{j=1}{\mathrm{\Σ}}}}{E}_{\mathrm{qj}}{B}_{\mathrm{ij}}[\sin ({\mathrm{\δ}}_i^0-{\mathrm{\δ}}_j^0)\cos (\mathrm{\Δ}{\mathrm{\δ}}_i-\mathrm{\Δ}{\mathrm{\δ}}_j)+\cos ({\mathrm{\δ}}_i^0-{\mathrm{\δ}}_j^0)\sin (\mathrm{\Δ}{\mathrm{\δ}}_i-\mathrm{\Δ}{\mathrm{\δ}}_j)]\end{array}$

假设增量很小，即偏差接近于零，有

cos(△δ_i-△δ)≈0

sin(△δ_i-△δ_j)≈△δ_i-△δ_j

则

$\mathrm{\Δ}{P}_{\mathrm{ei}}=E_{\mathrm{qi}}\underset{j\≠i}{\overset{3}{\underset{j=1}{\mathrm{\Σ}}}}{E}_{\mathrm{qj}}{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i^0-{\mathrm{\δ}}_j^0)(\mathrm{\Δ}{\mathrm{\δ}}_i-\mathrm{\Δ}{\mathrm{\δ}}_j)---\left(2\right)$

这里将发电机输入的机械功率的偏差△P_mi作为控制向量U，发电机转子转动时功角差△δ_i、转速差△ω_i作为状态向量X，将式(2)代入式(1)，整理得：

$\stackrel{\·}{X}=\mathrm{AX}+\mathrm{BU}---\left(3\right)$

其中状态向量X、控制向量U：

X＝[△δ₁,△ω₁,△δ₂,△ω₂,△δ₃,△ω₃]^T，U＝[△P_m1,△P_m2,△P_m3]^T；

且有

$\begin{array}{cc}A=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ -k_{12}-k_{13}& -\frac{D_1}{T_{j1}}& k_{12}& 0& k_{13}& 0\\ 0& 0& 0& 1& 0& 0\\ k_{21}& 0& -k_{21}-k_{23}& -\frac{D_2}{T_{j2}}& k_{23}& 0\\ 0& 0& 0& 0& 0& 1\\ k_{31}& 0& k_{32}& 0& -k_{31}-k_{32}& -\frac{D_3}{T_{j3}}\end{array}\right]& \end{array}B=\left[\begin{array}{ccc}0& 0& 0\\ \frac{{\mathrm{\ω}}_{0i}}{T_{j1}}& 0& 0\\ 0& 0& 0\\ 0& \frac{{\mathrm{\ω}}_{0i}}{T_{j2}}& 0\\ 0& 0& 0\\ 0& 0& \frac{{\mathrm{\ω}}_{0i}}{T_{j3}}\end{array}\right]$

$\begin{array}{ccc}{k}_{\mathrm{ij}}=\frac{{\mathrm{\ω}}_{0i}{E}_i{E}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i\left(0\right)-{\mathrm{\δ}}_j\left(0\right))}{T_{\mathrm{ji}}}& i=\mathrm{1,2,3}& i\≠j\end{array};$

步骤3：规定二次型性能指标，首先选用二次型性能指标作为系统的性能指标，即

$J={\∫}_0^{\∞}\frac{1}{2}(X^T\left(t\right)\mathrm{QX}\left(t\right)+U^T\left(t\right)\mathrm{RU}\left(t\right))\mathrm{dt}---\left(4\right)$

其中，X为系统的状态向量，这里X被选取为能够描述电力系统动态运行过程中变量的偏差，如△δ，△ω，△E’_q，△E_f，△P_m等；U为控制向量，如汽轮机调速汽门控制量偏差△P_m；Q为正定或半正定状态权矩阵，R为正定控制量权矩阵；

步骤4：对于多机电力系统的各子系统，按子系统状态向量反馈进行最优分散协调控制的设计；其具体实现包括以下子步骤：

步骤4.1：对于多机电力系统的各子系统，设计具有分散控制结构的形如U_i(t)＝K_iX_i(t)(i＝1,2,…,N)的控制器，即各子系统控制器仅反馈本子系统的状态向量，使全系统的二次型性能指标达到最小；

对于多机电力系统的各子系统：

$\left\{\begin{array}{c}\stackrel{\·}{X}\left(t\right)=\mathrm{AX}\left(t\right)+\mathrm{BU}\left(t\right)\\ X\left(0\right)=X_0\end{array}\right.---\left(5\right)$

U_i(t)＝K_iX_i(t)

写成集中形式即为：

或

U(t)＝K_dX(t)  (6)

其中K_d＝blockdiag{K₁K₂...K_N}；

将式(5)代入式(6)中，得：

$\stackrel{\·}{X}\left(t\right)=(A+B{K}_d)X\left(t\right)---\left(7\right)$

其解为：

$\stackrel{\·}{X}\left(t\right)=\mathrm{\φ}\left(t\right)X_0---\left(8\right)$

其中 $\mathrm{\φ}\left(t\right)=e^{(A+B{K}_d)t}.;$

将式(8)代入式(4)，则原系统性能指标为：

$J=X_0^T\left[{\∫}_0^{\∞}{\mathrm{\φ}}^T\left(t\right)(Q+K_d^{T}R{K}_d)\mathrm{\φ}\left(t\right)\mathrm{dt}\right]X_0=X_0^{T}P{X}_0$

其中

$P={\∫}_0^{\∞}{\mathrm{\φ}}^T\left(t\right)(Q+K_d^{T}R{K}_d)\mathrm{\φ}\left(t\right)\mathrm{dt}$

为矩阵方程

$P(A+B{K}_d)+{(A+B{K}_d)}^{T}P+Q+K_d^{T}R{K}_d=0---\left(9\right)$

的解，且由式(9)矩阵方程可知P^T也为方程的解，因此，P为对称矩阵，有P^T＝P；

步骤4.2：求解最优控制规律，即使性能指标函数J达到极小值对系统进行线性最优控制，即使二次型性能指标函数J达到极小值；

当矩阵A、B维数不同时，有下式成立

B^TA＝tr(AB^T)

其中tr()表示矩阵的迹，即矩阵对角线元素之和；

则性能指标函数J可改写成：

$J=\mathrm{tr}\left(P{X}_0{X}_0^T\right)$

假设初始状态变量X₀为均匀分布与n维单位球面上的随机向量，则上式又可写成：

J＝tr(P)

令

$G(P,K_d)=P(A+B{K}_d)+{(A+B{K}_d)}^{T}P+Q+K_d^{T}R{K}_d---\left(10\right)$

则约束条件式(10)为

G(P,K_d)＝0

做拉格朗日函数

$L=\mathrm{trP}+\underset{i=j}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{v}_{\mathrm{ij}}{g}_{\mathrm{ij}}---\left(11\right)$

式中g_ij为G(P,K)的第i行，第j列元素；g_ij为对应元素的拉格朗日乘子；

根据矩阵的迹定义，式(11)可写成

L＝trP+tr[V^TG(P,K)]    (12)

步骤4.3：最后利用矩阵迹的基本运算法则推导出如下方程组：

$\left\{\begin{array}{c}R{K}_d{V}_d+{\left(B^T\mathrm{PV}\right)}_d=0\\ P(A+B{K}_d)+{(A+B{K}_d)}^{T}P+Q+K_d^T{\mathrm{RK}}_d=0\\ (A+B{K}_d)V+V{(A+B{K}_d)}^T+I=0\end{array}\right.---\left(13\right)$

(13)式为Levine-Athans方程组，对该方程组求解得到反馈增益K_d即为最优分散协调控制规律；

步骤5：用直接迭代法求解最优分散协调控制规律，其具体实现包括以下子步骤：

步骤5.1：对于式(13)第一个方程式求解可以得到：

步骤5.2：选择初始K_d ⁰值，使初始稳定；用Matlab求解矩阵A的特征值，发现所有特征根具有负的实部，因此矩阵A初始稳定，此时取初始K_d ⁰为零矩阵；

步骤5.3：给定系统允许误差ε，取迭代初始步长θ，令i＝0，

步骤5.4：由式 $P^{\left(i\right)}(A+B{K}_d^{\left(i\right)})+{(A+B{K}_d^{\left(i\right)})}^T{P}^{\left(i\right)}+Q+{K_d^{\left(i\right)}}^{T}R{K}_d^{\left(i\right)}=0$ 求得P⁽ⁱ⁾；由式 $(A+B{K}_d^{\left(i\right)})V^{\left(i\right)}+V^{\left(i\right)}{(A+B{K}_d^{\left(i\right)})}^T+I=0$ 求得V⁽ⁱ⁾；

步骤5.5：计算 $K_{d0}^{(i+1)}=-R^{-1}{\left(B^T{P}^{\left(i\right)}{V}^{\left(i\right)}\right)}_d{V}_d^{\left(i\right)-1};$

步骤5.6：判断迭代精度，若则停止迭代，为所求解，否则进入下一步；

步骤5.7：计算 $K_d^{(i+1)}=(1-\mathrm{\θ})K_d^{\left(i\right)}+\mathrm{\θ}{K}_{d0}^{(i+1)};$

步骤5.8：由式 $P^{(i+1)}(A+B{K}_d^{(i+1)})+{(A+B{K}_d^{(i+1)})}^T{P}^{(i+1)}+Q+{K_d^{(i+1)}}^{T}R{K}_d^{(i+1)}=0$ 求得P⁽ⁱ⁺¹)；

步骤5.9：判断trP⁽ⁱ⁺¹⁾<trP⁽ⁱ⁾，若成立，则令θ变为1.2θ，i变为i+1，返回到步骤5.4；否则，θ变为0.5θ，返回到步骤5.8。

作为优选，步骤2.2中所述的对发电机输出的电磁功率进行线性化，其具体实现过程是对发电机输出的电磁功率在操作点处进行线性化，用与分别替代P_ei、δ_i、δ_j，则

$\begin{array}{c}{P}_{\mathrm{ei}}^0+\mathrm{\&Delta;}{P}_{\mathrm{ei}}=E_{\mathrm{qi}}^2{G}_{\mathrm{ii}}+E_{\mathrm{qi}}\underset{j\&NotEqual;i}{\overset{3}{\underset{j=1}{\mathrm{\&Sigma;}}}}{E}_{\mathrm{qj}}{B}_{\mathrm{ij}}\sin ({\mathrm{\&delta;}}_i^0+\mathrm{\&Delta;}{\mathrm{\&delta;}}_i-{\mathrm{\&delta;}}_j^0-\mathrm{\&Delta;}{\mathrm{\&delta;}}_j)\\ =E_{\mathrm{qi}}^2{G}_{\mathrm{ii}}+E_{\mathrm{qi}}\underset{j\&NotEqual;i}{\overset{3}{\underset{j=1}{\mathrm{\&Sigma;}}}}{E}_{\mathrm{qj}}{B}_{\mathrm{ij}}[\sin ({\mathrm{\&delta;}}_i^0-{\mathrm{\&delta;}}_j^0)\cos (\mathrm{\&Delta;}{\mathrm{\&delta;}}_i-\mathrm{\&Delta;}{\mathrm{\&delta;}}_j)+\cos ({\mathrm{\&delta;}}_i^0-{\mathrm{\&delta;}}_j^0)\sin (\mathrm{\&Delta;}{\mathrm{\&delta;}}_i-\mathrm{\&Delta;}{\mathrm{\&delta;}}_j)]\end{array}$

假设增量很小，即偏差接近于零，有

cos(△δ_i-△δ)≈0

sin(△δ_i-△δ_j)≈△δ_i-△δ_j。

作为优选，步骤4.2中所述的使性能指标函数J达到极小值，即使L取得极值，L取得极值必要条件为：

$\frac{\&PartialD;L}{\&PartialD;V}=0;\frac{\&PartialD;L}{\&PartialD;P}=0;\frac{\&PartialD;L}{\&PartialD;K_d}=0.$

作为优选，步骤5.2中所述的选择初始K_d ⁰值，有如下两类：

(1)当系统A矩阵初始稳定时，可将K_d ⁰取为零矩阵；

(2)当系统A矩阵初始值不稳定时，可适当降低系统运行点，或者设定一个机械阻尼(例如D_i＞8.2279)，使系统由不稳定变为稳定，先用零矩阵作初始反馈矩阵经上述步骤解得该运行点处系统反馈矩阵K′_d，再将该反馈矩阵K′_d作为初始反馈矩阵求解原系统的反馈矩阵K_d，从而得到原系统最优分散协调控制器的解。

作为优选，步骤5.3中所述的给定系统允许误差ε，取迭代初始步长θ，其中ε＝0.01，θ＝0.5。

本发明以多机电力系统中同步发电机为研究对象，对多机电力系统进行数学建模，对模型进行分析与研究，针对电力系统不同的稳定情况，利用线性控制理论，设计出能适应静态稳定情况的最优分散协调控制规律控制方法，再对控制规律进行离线与实时仿真研究进行验证，得出解决电力系统稳定运行直接而有效的方法，对多机电力系统控制规律的研究与验证对提高电力系统稳定运行具有理论价值和实际意义。

附图说明

图1：是本发明实施例所采用的IEEE三机九节点模型图。

图2：是本发明实施例所采用的三机九节点模型的Matlab仿真模型。

图3：是本发明实施例的Matlab仿真得出的最优分散协调控制与全局线性最优控制功角差比较图。

图4：是本发明实施例的Matlab仿真得出的以2号机组为例的最优分散协调控制与全局线性最优控制转速比较图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

下面以三机九节点模型为例，对多机电力系统分散协调最优控制进行设计。图1为IEEE三机九节点模型，其中所有阻抗都是以100MVA为基准的标幺值。

将多机电力系统转子运动方程列出并写成如下形式：

$\left\{\begin{array}{c}\frac{{\mathrm{d\&delta;}}_i}{\mathrm{dt}}={\mathrm{\&omega;}}_i-{\mathrm{\&omega;}}_{0i}\\ \frac{{\mathrm{d\&omega;}}_i}{\mathrm{dt}}=\frac{{\mathrm{\&omega;}}_{0i}}{T_{\mathrm{Ji}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}})-\frac{D_i}{T_{\mathrm{Ji}}}({\mathrm{\&omega;}}_i-{\mathrm{\&omega;}}_{0i})\end{array}\right.---\left(1\right)$

其中ω_i为有名值，单位rad/s。

三机九节点实验系统参数采用IEEE标准的三机九节点模型参数：

ω₀＝2πf₀＝314.159rad/s；

发电机阻尼系数D₁＝D₂＝D₃＝5；

发电机惯量时间常数T_jiT_j1＝47.28s T_j2＝12.8s T_j3＝6.02s；

发电机电动势E_iE₁＝1.0558 E₂＝0.9885 E₃＝0.9902；

机械功率P_miP_m1＝0.71641 P_m2＝1.63 P_m3＝0.85；

稳定运行时系统导纳矩阵：

$Y=\left[\begin{array}{ccc}0.846-2.988i& 0.287+1.513i& 0.210+1.226i\\ 0.287+1.513i& 0.420-2.724i& 0.213+1.088i\\ 0.210+1.226i& 0.213+1.088i& 0.277-2.368i\end{array}\right]$

接下来对上述系统进行最优分散协调控制的设计，具体步骤如下：

步骤1:非线性对电力系统模型系统线性化

这里采用IEEE的如图1所示的三机九节点数学建模。然后对该模型进行分散协调控制的设计，由于最优分散协调控制设计要求系统为线性系统，因此首先需对转子运动方程进行线性化与偏差化，即有：

$\left\{\begin{array}{c}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{\mathrm{\&delta;}}}_i=\mathrm{\&Delta;}{\mathrm{\&omega;}}_i\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_i=-\frac{{\mathrm{\&omega;}}_{0i}}{T_{\mathrm{Ji}}}\mathrm{\&Delta;}{P}_{\mathrm{ei}}-\frac{D_i}{T_{\mathrm{Ji}}}\mathrm{\&Delta;}{\mathrm{\&omega;}}_i+\frac{{\mathrm{\&omega;}}_{0i}}{T_{\mathrm{Ji}}}\mathrm{\&Delta;}{P}_{\mathrm{mi}}\end{array}\right.---\left(2\right)$

$\mathrm{\&Delta;}{P}_{\mathrm{ei}}=E_{\mathrm{qi}}\underset{j\&NotEqual;i}{\overset{3}{\underset{j=1}{\mathrm{\&Sigma;}}}}{E}_{\mathrm{qj}}{B}_{\mathrm{ij}}\cos ({\mathrm{\&delta;}}_i^0-{\mathrm{\&delta;}}_j^0)(\mathrm{\&Delta;}{\mathrm{\&delta;}}_i-\mathrm{\&Delta;}{\mathrm{\&delta;}}_j)$

写成线性状态方程为：

$\stackrel{\&CenterDot;}{X}=\mathrm{AX}+\mathrm{BU}---\left(3\right)$

其中状态向量X、控制向量U、输出向量Y分别为：

$\begin{array}{cc}X=\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&delta;}}_1\\ \mathrm{\&Delta;}{\mathrm{\&omega;}}_1\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_2\\ \mathrm{\&Delta;}{\mathrm{\&omega;}}_2\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_3\\ \mathrm{\&Delta;}{\mathrm{\&omega;}}_3\end{array}\right]U=\left[\begin{array}{c}\mathrm{\&Delta;}{P}_{m1}\\ \mathrm{\&Delta;}{P}_{m2}\\ \mathrm{\&Delta;}{P}_{m3}\end{array}\right]& Y=\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&delta;}}_1\\ \mathrm{\&Delta;}{\mathrm{\&omega;}}_1\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_2\\ \mathrm{\&Delta;}{\mathrm{\&omega;}}_2\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_3\\ \mathrm{\&Delta;}{\mathrm{\&omega;}}_3\end{array}\right]\end{array}$

且有

$\begin{array}{cc}A=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ -k_{12}-k_{13}& -\frac{D_1}{T_{j1}}& k_{12}& 0& k_{13}& 0\\ 0& 0& 0& 1& 0& 0\\ k_{21}& 0& -k_{21}-k_{23}& -\frac{D_2}{T_{j2}}& k_{23}& 0\\ 0& 0& 0& 0& 0& 1\\ k_{31}& 0& k_{32}& 0& -k_{31}-k_{32}& -\frac{D_3}{T_{j3}}\end{array}\right]& \end{array}B=\left[\begin{array}{ccc}0& 0& 0\\ \frac{{\mathrm{\&omega;}}_{0i}}{T_{j1}}& 0& 0\\ 0& 0& 0\\ 0& \frac{{\mathrm{\&omega;}}_{0i}}{T_{j2}}& 0\\ 0& 0& 0\\ 0& 0& \frac{{\mathrm{\&omega;}}_{0i}}{T_{j3}}\end{array}\right]$

$\begin{array}{ccc}{k}_{\mathrm{ij}}=\frac{{\mathrm{\&omega;}}_{0i}{E}_i{E}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\&delta;}}_i\left(0\right)-{\mathrm{\&delta;}}_j\left(0\right))}{T_{\mathrm{ji}}}& i=\mathrm{1,2,3}& i\&NotEqual;j\end{array}$

步骤2:规定二次型性能指标

根据线性最优控制设计原理，首先选用二次型性能指标作为系统的性能指标，即

$J={\&Integral;}_0^{\&infin;}\frac{1}{2}(X^T\left(t\right)\mathrm{QX}\left(t\right)+U^T\left(t\right)\mathrm{RU}\left(t\right))\mathrm{dt}---\left(4\right)$

其中，X为系统的状态向量，例如对于电力系统来说，X常取能够描述电力系统动态运行过程中变量的偏差，如△δ，△ω，△E’_q，△E_f，△P_m等；U为控制向量，如汽轮机调速汽门控制量偏差△P_m；Q为正定或半正定状态权矩阵；R为正定控制量权矩阵。对系统进行线性最优控制，即使二次型性能指标J达到极小值。

步骤3:对于多机电力系统的各子系统，进行子系统状态量反馈的最优分散协调控制的设计；

对于多机电力系统：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{X}}_i\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{\mathrm{ij}}{X}_j\left(t\right)+B_i{U}_i\left(t\right)\\ X_i\left(0\right)=X_{i0},i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,N\end{array}\right.$

其中X_i为状态向量，且X_i∈Rⁿ _i；U_i为控制向量，且U_i∈R^r _i；

上式各子系统可写成集合的形式：

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{X}\left(t\right)=\mathrm{AX}\left(t\right)+\mathrm{BU}\left(t\right)\\ X\left(0\right)=X_0\end{array}\right.---\left(5\right)$

状态量反馈分散控制所要解决的问题可描述如下：设计具有分散控制结构的，形如U_i(t)＝K_iX_i(t)(i＝1,2,…,N)的控制器，即各子系统控制器仅反馈本子系统的状态量，使全系统的二次型性能指标达到最小。U_i(t)＝K_iX_i(t)写成集中形式即为：

或

U(t)＝K_dX(t)    (6)

其中K_d＝blockdiag{K₁ K₂…K_N}

将式(6)代入式(5)中，得：

$\stackrel{\&CenterDot;}{X}\left(t\right)=(A+B{K}_d)X\left(t\right)---\left(7\right)$

其解为：

$\stackrel{\&CenterDot;}{X}\left(t\right)=\mathrm{\&phi;}\left(t\right)X_0---\left(8\right)$

其中

$\mathrm{\&phi;}\left(t\right)=e^{(A+B{K}_d)t}$

将式(8)代入式(4)，则原系统性能指标为：

$J=X_0^T\left[{\&Integral;}_0^{\&infin;}{\mathrm{\&phi;}}^T\left(t\right)(Q+K_d^{T}R{K}_d)\mathrm{\&phi;}\left(t\right)\mathrm{dt}\right]X_0=X_0^{T}P{X}_0$

其中

$P={\&Integral;}_0^{\&infin;}{\mathrm{\&phi;}}^T\left(t\right)(Q+K_d^{T}R{K}_d)\mathrm{\&phi;}\left(t\right)\mathrm{dt}$

为矩阵方程

$P(A+B{K}_d)+{(A+B{K}_d)}^{T}P+Q+K_d^{T}R{K}_d=0---\left(9\right)$

的解，且有(9)式矩阵方程可知P^T也为方程的解，因此，P为对称矩阵，有P＝P^T。

求解最优控制规律，即使性能指标函数J达到极小值对系统进行线性最优控制，即使二次型性能指标函数J达到极小值。

当矩阵A、B不同维，有下式成立

B^TA＝tr(AB^T)

其中tr()表示矩阵的迹，即矩阵对角线元素之和。

则性能指标函数J可改写成：

$J=\mathrm{tr}\left(P{X}_0{X}_0^T\right)$

假设初始状态变量X₀为均匀分布与n维单位球面上的随机向量，则上式又可写成：

J＝tr(P)

令

$G(P,K_d)=P(A+B{K}_d)+{(A+B{K}_d)}^{T}P+Q+K_d^{T}R{K}_d---\left(10\right)$

则约束条件式(10)为

G(P,K_d)＝0

做拉格朗日函数

$L=\mathrm{reP}+\underset{i=j}{\overset{n}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{v}_{\mathrm{ij}}{g}_{\mathrm{ij}}---\left(11\right)$

式中g_ij为G(P，K)的第i行，第j列元素；g_ij为对应元素的拉格朗日乘子。

根据矩阵的迹定义，上式(11)可写成

L＝trP+tr[V^TG(P,K)]    (12)

为使性能指标J取得极值即使L取得极值，L取得极值必要条件为：

$\frac{\&PartialD;L}{\&PartialD;V}=0;\frac{\&PartialD;L}{\&PartialD;P}=0;\frac{\&PartialD;L}{\&PartialD;K_d}=0;$

最后利用矩阵迹的基本运算法则推导方程组：

$\left\{\begin{array}{c}R{K}_d{V}_d+{\left(B^T\mathrm{PV}\right)}_d=0\\ P(A+B{K}_d)+{(A+B{K}_d)}^{T}P+Q+K_d^T{\mathrm{RK}}_d=0\\ (A+B{K}_d)V+V{(A+B{K}_d)}^T+I=0\end{array}\right.---\left(13\right)$

(13)式为著名的Levine-Athans方程组，对该方程组求解得到反馈增益Kd即为最优分散协调控制规律。

步骤4:最优分散协调控制规律求取；

根据最优分散协调控制的设计方法，K_d可通过求解以下Levine-Athans方程组得到：

$\left\{\begin{array}{c}R{K}_d{V}_d+{\left(B^T\mathrm{PV}\right)}_d=0\\ P(A+B{K}_d)+{(A+B{K}_d)}^{T}P+Q+K_d^T{\mathrm{RK}}_d=0\\ (A+B{K}_d)V+V{(A+B{K}_d)}^T+I=0\end{array}\right.$

由最优分散协调控制的设计原理可知，需找到U＝K_d·X，即求解分块对角反馈增益矩阵K_d使全系统性能指标

$J=\underset{0}{\overset{T}{\&Integral;}}[X^T\left(t\right)\mathrm{QX}\left(t\right)+U^T\left(t\right)\mathrm{RU}\left(t\right)]\mathrm{dt}$

达极小值。

步骤5:三机九节点电力系统最优分散协调控制Matlab仿真

利用IEEE标准的三机九节点模型参数，采用直接迭代法求解下式Levine-Athans方程组

$\left\{\begin{array}{c}R{K}_d{V}_d+{\left(B^T\mathrm{PV}\right)}_d=0\\ P(A+B{K}_d)+{(A+B{K}_d)}^{T}P+Q+K_d^T{\mathrm{RK}}_d=0\\ (A+B{K}_d)V+V{(A+B{K}_d)}^T+I=0\end{array}\right.$

其中取初始控制器K_d ⁰＝diag(100 100,100 100,100 100)；取状态权矩阵Q＝diag(1,1,1,1,1,1)，取控制量权矩阵R＝diag(1,1,1)。

利用Matlab编写直接迭代法程序求解控制器K_d。

求得最终控制器K_d：

$K_{-1.2405d}=\left[\begin{array}{cccccc}-0.8780& -1.1799& 0& 0& 0& 0\\ 0& 0& -1.0508& -1.1972& 0& 0\\ 0& 0& 0& 0& -1.1567& -1.2405\end{array}\right]$

即有：

△P_m1＝-0.8780△δ₁-1.1799△ω₁

△P_m2＝-1.0508△δ₂-1.1972△ω₂

△P_m3＝-1.1567△δ₃-1.2405△ω₃

图2为Matlab环境下搭建的三机九节点仿真模型，下面利用Matlab对三机九节点最优分散协调控制策略进行仿真，并与全局线性最优控制策略比较分析，以二号机组为例，结果如图3和图4所示。由图3和图4可知，分散协调控制对系统扰动同样有很好的抑制作用，控制作用与全状态量反馈线性最优控制类似，但综合考虑，分散协调控制效果不及全状态量反馈线性最优控制效果好，主要体现在调节时间比全状态量反馈线性最优控制稍长。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (5)

1.一种多机电力系统稳定运行的分散协调控制方法，其特征在于，包括以下步骤：

步骤1：对多机电力系统被控对象进行数学建模；其具体实现包括以下子步骤：

步骤1.1：首先从力学的角度推导同步发电机转子运动方程，即经典的摇摆方程；

步骤1.2：接着对同步电机转子运动方程中涉及的发电机输出功率进行推导与求解；

步骤1.3：最后从电路的角度对同步发电机的励磁绕组电磁动态方程进行分析与求解，得到多机电力系统第i台同步发电机的数学模型表达式为：

$\left\{\begin{array}{c}\frac{d{\mathrm{\&delta;}}_i}{\mathrm{dt}}=({\mathrm{\&omega;}}_i-1){\mathrm{\&omega;}}_{0i}\\ \frac{d{\mathrm{\&omega;}}_i}{\mathrm{dt}}=\frac{1}{T_{\mathrm{Ji}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}})-\frac{D_i}{T_{\mathrm{Ji}}}({\mathrm{\&omega;}}_i-1)\\ \frac{d{E^{\&prime;}}_{\mathrm{qi}}}{\mathrm{dt}}=\frac{1}{{T^{\&prime;}}_{d0i}}[E_{\mathrm{fi}}-E_{\mathrm{qi}}^{\&prime;}-(x_{\mathrm{di}}-x_{\mathrm{di}}^{\&prime;})i_{\mathrm{di}}]\end{array}\right.$

多机电力系统： $\begin{array}{c}{P}_{\mathrm{ei}}=E_{\mathrm{qi}}^2{G}_{\mathrm{ii}}+E_{\mathrm{qi}}{{\mathrm{\&Sigma;}}_{\underset{j\&NotEqual;i}{j=1}}^n{E}_{\mathrm{qj}}(G_{\mathrm{ij}}\cos {\mathrm{\&delta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\&delta;}}_{\mathrm{ij}})}_{\begin{array}{c}\\ \end{array}}i=\{\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;n\};\end{array}$ 式中，δ_i为第i台发电机转子运行角，即q轴与同步参考轴S间的夹角，单位rad，ω_i表示转子运动角速度(标幺值)，ω_0i为转子运动角速度的初始值，单位rad/s，T_Ji为转子转动惯量时间常数，单位s，P_mi表示第i台发电机输入的机械功率(标幺值)，P_ei表示第i台发电机输出的电磁功率(标幺值)，D_i表示阻尼系数(标幺值)，E’_qi表示第i台发电机暂态电动势(标幺值)，E’_fi表示稳态磁链在d轴定子侧产生的电势(标幺值)，i_di表示第i台发电机输电流瞬态值(标幺值)，x_di、x’_di分别为第i台发电机直轴同步电抗与瞬态电抗(标幺值)；且I＝YU，Y＝G+jB，为电力系统简化导纳矩阵，G_ii＝Y_iisinα_ii和B_ii＝Y_iicosα_ii为第i节点的自电导，G_ij＝Y_ijsinα_ii和B_ij＝Y_ijcosα_ii为第i节点与第j节点间的互电导及互导纳；

步骤2：对多机电力系统数学模型进行非线性系统线性化，其具体实现包括以下子步骤：

步骤2.1：首先对其转子运动方程进行线性化，即：

$\left\{\begin{array}{c}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{\mathrm{\&delta;}}}_i=\mathrm{\&Delta;}{\mathrm{\&omega;}}_i\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_i=-\frac{{\mathrm{\&omega;}}_{0i}}{T_{\mathrm{Ji}}}\mathrm{\&Delta;}{P}_{\mathrm{ei}}-\frac{D_i}{T_{\mathrm{Ji}}}\mathrm{\&Delta;}{\mathrm{\&omega;}}_i+\frac{{\mathrm{\&omega;}}_{0i}}{T_{\mathrm{Ji}}}\mathrm{\&Delta;}{P}_{\mathrm{mi}}\end{array}\right.---\left(1\right)$

其中，δ_i为第i台发电机转子运行角，单位rad，ω_i为发电机转子运动角速度(标幺值)，ω_0i为第i台发电机转子运动角速度的初始值，单位rad/s，T_Ji为第i台发电机转子转动惯量时间常数，单位s，P_mi表示第i台发电机输入的机械功率(标幺值)，P_ei表示第i台发电机输出的电磁功率(标幺值)，D_i为第i台发电机阻尼系数(标幺值)；

步骤2.2：对发电机输出的电磁功率进行线性化，即：

$\mathrm{\&Delta;}{P}_{\mathrm{ei}}=E_{\mathrm{qi}}\underset{j\&NotEqual;i}{\overset{3}{\underset{j=1}{\mathrm{\&Sigma;}}}}{E}_{\mathrm{qj}}{B}_{\mathrm{ij}}\cos ({\mathrm{\&delta;}}_i^0-{\mathrm{\&delta;}}_j^0)(\mathrm{\&Delta;}{\mathrm{\&delta;}}_i-\mathrm{\&Delta;}{\mathrm{\&delta;}}_j)---\left(2\right)$

这里将发电机输入的机械功率的偏差ΔP_mi作为控制向量U，发电机转子转动时功角差Δδ_i、转速差Δω_i作为状态向量X，将式(2)代入式(1)，整理得：

$\stackrel{\&CenterDot;}{X}=\mathrm{AX}+\mathrm{BU}---\left(3\right)$

其中状态向量X、控制向量U：

X＝[Δδ₁,Δω₁,Δδ₂,Δω₂,Δδ₃,Δω₃]^T，U＝[ΔP_m1,ΔP_m2,ΔP_m3]^T；

且有

$\begin{array}{cc}A=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ -k_{12}-k_{13}& -\frac{D_1}{T_{j1}}& k_{12}& 0& k_{13}& 0\\ 0& 0& 0& 1& 0& 0\\ k_{21}& 0& -k_{21}-k_{23}& -\frac{D_2}{T_{j2}}& k_{23}& 0\\ 0& 0& 0& 0& 0& 1\\ k_{31}& 0& k_{32}& 0& -k_{31}-k_{32}{}_{}& -\frac{D_3}{T_{j3}}\end{array}\right]& \end{array}B=\left[\begin{array}{ccc}0& 0& 0\\ \frac{{\mathrm{\&omega;}}_{0i}}{T_{j1}}& 0& 0\\ 0& 0& 0\\ 0& \frac{{\mathrm{\&omega;}}_{0i}}{T_{j2}}& 0\\ 0& 0& 0\\ 0& 0& \frac{{\mathrm{\&omega;}}_{0i}}{T_{j3}}\end{array}\right]$

$\begin{array}{ccc}{k}_{\mathrm{ij}}=\frac{{\mathrm{\&omega;}}_{0i}{E}_i{E}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\&delta;}}_i\left(0\right)-{\mathrm{\&delta;}}_j\left(0\right))}{T_{\mathrm{ji}}}& i=\mathrm{1,2,3}& i\&NotEqual;j\end{array};$

步骤3：规定二次型性能指标，首先选用二次型性能指标作为系统的性能指标，即

$J={\&Integral;}_0^{\&infin;}\frac{1}{2}(X^T\left(t\right)\mathrm{QX}\left(t\right)+U^T\left(t\right)\mathrm{RU}\left(t\right))\mathrm{dt}---\left(4\right)$

其中，X为系统的状态向量，这里X被选取为能够描述电力系统动态运行过程中变量的偏差，U为控制向量，Q为正定或半正定状态权矩阵，R为正定控制量权矩阵；

步骤4：对于多机电力系统的各子系统，按子系统状态向量反馈进行最优分散协调控制的设计；其具体实现包括以下子步骤：

步骤4.1：对于多机电力系统的各子系统，设计具有分散控制结构的控制器，即各子系统控制器仅反馈本子系统的状态向量，使全系统的二次型性能指标达到最小；

对于多机电力系统的各子系统：

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{X}\left(t\right)=\mathrm{AX}\left(t\right)+\mathrm{BU}\left(t\right)\\ X\left(0\right)=X_0\end{array}\right.---\left(5\right)$

U(t)＝K_dX(t)   (6)

其中K_d＝blockdiag{K₁K₂...K_N}；

将式(5)代入式(6)中，得：

$\stackrel{\&CenterDot;}{X}\left(t\right)=(A+B{K}_d)X\left(t\right)---\left(7\right)$

其解为：

$\stackrel{\&CenterDot;}{X}\left(t\right)=\mathrm{\&phi;}\left(t\right)X_0---\left(8\right)$

其中 $\mathrm{\&phi;}\left(t\right)=e^{(A+B{K}_d)t}.;$

将式(8)代入式(4)，则原系统性能指标为：

$J=X_0^T\left[{\&Integral;}_0^{\&infin;}{\mathrm{\&phi;}}^T\left(t\right)(Q+K_d^{T}R{K}_d)\mathrm{\&phi;}\left(t\right)\mathrm{dt}\right]X_0=X_0^{T}P{X}_0$

其中

$P={\&Integral;}_0^{\&infin;}{\mathrm{\&phi;}}^T\left(t\right)(Q+K_d^{T}R{K}_d)\mathrm{\&phi;}\left(t\right)\mathrm{dt}$

为矩阵方程

$P(A+B{K}_d)+{(A+B{K}_d)}^{T}P+Q+K_d^{T}R{K}_d=0---\left(9\right)$

的解，且由式(9)矩阵方程可知P^T也为方程的解，因此，P为对称矩阵，有P^T＝P；

步骤4.2：求解最优控制规律，即使性能指标函数J达到极小值对系统进行线性最优控制，即使二次型性能指标函数J达到极小值；

当矩阵A、B维数不同时，有下式成立

B^TA＝tr(AB^T)

其中tr()表示矩阵的迹，即矩阵对角线元素之和；

则性能指标函数J可改写成：

$J=\mathrm{tr}\left(P{X}_0{X}_0^T\right)$

假设初始状态变量X₀为均匀分布与n维单位球面上的随机向量，则上式又可写成：

J＝tr(P)

令

$G(P,K_d)=P(A+B{K}_d)+{(A+B{K}_d)}^{T}P+Q+K_d^{T}R{K}_d---\left(10\right)$

则约束条件式(10)为

G(P,K_d)＝0

做拉格朗日函数

$L=\mathrm{trP}+\underset{i=j}{\overset{n}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{v}_{\mathrm{ij}}{g}_{\mathrm{ij}}---\left(11\right)$

式中g_ij为G(P,K)的第i行，第j列元素；g_ij为对应元素的拉格朗日乘子；

根据矩阵的迹定义，式(11)可写成

L＝trP+tr[V^TG(P,K)]   (12)

步骤4.3：最后利用矩阵迹的基本运算法则推导出如下方程组：

$\left\{\begin{array}{c}R{K}_d{V}_d+{\left(B^T\mathrm{PV}\right)}_d=0\\ P(A+B{K}_d)+{(A+B{K}_d)}^{T}P+Q+K_d^T{\mathrm{RK}}_d=0\\ (A+B{K}_d)V+V{(A+B{K}_d)}^T+I=0\end{array}\right.---\left(13\right)$

(13)式为Levine-Athans方程组，对该方程组求解得到反馈增益K_d即为最优分散协调控制规律；

步骤5：用直接迭代法求解最优分散协调控制规律，其具体实现包括以下子步骤：

步骤5.1：对于式(13)第一个方程式求解可以得到：

步骤5.2：选择初始K_d ⁰值，使初始稳定；用Matlab求解矩阵A的特征值，发现所有特征根具有负的实部，因此矩阵A初始稳定，此时取初始K_d ⁰为零矩阵；

步骤5.3：给定系统允许误差ε，取迭代初始步长θ，令i＝0，

步骤5.4：由式 $P^{\left(i\right)}(A+B{K}_d^{\left(i\right)})+{(A+B{K}_d^{\left(i\right)})}^T{P}^{\left(i\right)}+Q+{K_d^{\left(i\right)}}^{T}R{K}_d^{\left(i\right)}=0$ 求得P⁽ⁱ⁾；由式 $(A+B{K}_d^{\left(i\right)})V^{\left(i\right)}+V^{\left(i\right)}{(A+B{K}_d^{\left(i\right)})}^T+I=0$ 求得V⁽ⁱ⁾；

步骤5.5：计算 $K_{d0}^{(i+1)}=-R^{-1}{\left(B^T{P}^{\left(i\right)}{V}^{\left(i\right)}\right)}_d{V}_d^{\left(i\right)-1};$

步骤5.6：判断迭代精度，若则停止迭代，为所求解，否则进入下一步；

步骤5.7：计算 $K_d^{(i+1)}=(1-\mathrm{\&theta;})K_d^{\left(i\right)}+\mathrm{\&theta;}{K}_{d0}^{(i+1)};$

步骤5.8：由式 $P^{(i+1)}(A+B{K}_d^{(i+1)})+{(A+B{K}_d^{(i+1)})}^T{P}^{(i+1)}+Q+{K_d^{(i+1)}}^{T}R{K}_d^{(i+1)}=0$ 求得P⁽ⁱ⁺¹⁾；

步骤5.9：判断trP⁽ⁱ⁺¹⁾＜trP⁽ⁱ⁾，若成立，则令θ变为1.2θ，i变为i+1，返回到步骤5.4；否则，θ变为0.5θ，返回到步骤5.8。

2.根据权利要求1所述的多机电力系统稳定运行的分散协调控制方法，其特征在于：步骤2.2中所述的对发电机输出的电磁功率进行线性化，其具体实现过程是对发电机输出的电磁功率在操作点处进行线性化，用与 分别替代P_ei、δ_i、δ_j，则

$\begin{array}{c}{P}_{\mathrm{ei}}^0+\mathrm{\&Delta;}{P}_{\mathrm{ei}}=E_{\mathrm{qi}}^2{G}_{\mathrm{ii}}+E_{\mathrm{qi}}\underset{j\&NotEqual;i}{\overset{3}{\underset{j=1}{\mathrm{\&Sigma;}}}}{E}_{\mathrm{qj}}{B}_{\mathrm{ij}}\sin ({\mathrm{\&delta;}}_i^0+\mathrm{\&Delta;}{\mathrm{\&delta;}}_i-{\mathrm{\&delta;}}_j^0-\mathrm{\&Delta;}{\mathrm{\&delta;}}_j)\\ =E_{\mathrm{qi}}^2{G}_{\mathrm{ii}}+E_{\mathrm{qi}}\underset{j\&NotEqual;i}{\overset{3}{\underset{j=1}{\mathrm{\&Sigma;}}}}{E}_{\mathrm{qj}}{B}_{\mathrm{ij}}[\sin ({\mathrm{\&delta;}}_i^0-{\mathrm{\&delta;}}_j^0)\cos (\mathrm{\&Delta;}{\mathrm{\&delta;}}_i-\mathrm{\&Delta;}{\mathrm{\&delta;}}_j)+\cos ({\mathrm{\&delta;}}_i^0-{\mathrm{\&delta;}}_j^0)\sin (\mathrm{\&Delta;}{\mathrm{\&delta;}}_i-\mathrm{\&Delta;}{\mathrm{\&delta;}}_j)]\end{array}$

假设增量很小，即偏差接近于零，有

cos(Δδ_i-Δδ)≈0

sin(Δδ_i-Δδ_j)≈Δδ_i-Δδ_j。

3.根据权利要求1所述的多机电力系统稳定运行的分散协调控制方法，其特征在于：步骤4.2中所述的使指标函数J达到极小值，即使L取得极值，L取得极值必要条件为：

$\frac{\&PartialD;L}{\&PartialD;V}=0;\frac{\&PartialD;L}{\&PartialD;P}=0;\frac{\&PartialD;L}{\&PartialD;K_d}=0.$

4.根据权利要求1所述的多机电力系统稳定运行的分散协调控制方法，其特征在于：步骤5.2中所述的选择初始K_d ⁰值，有如下两类：

(1)当系统A矩阵初始稳定时，可将K_d ⁰取为零矩阵；

(2)当系统A矩阵初始值不稳定时，可通过降低系统运行点，或者设定一个机械阻尼，使系统由不稳定变为稳定，先用零矩阵作初始反馈矩阵经上述步骤解得该运行点处系统反馈矩阵K′_d，再将该反馈矩阵K′_d作为初始反馈矩阵求解原系统的反馈矩阵K_d，从而得到原系统最优分散协调控制器的解。

5.根据权利要求1所述的多机电力系统稳定运行的分散协调控制方法，其特征在于：步骤5.3中所述的给定系统允许误差ε，取迭代初始步长θ，其中ε＝0.01，θ＝0.5。