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

摘要

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

主权项

1.一种多机电力系统稳定运行的分散协调控制方法，其特征在于，包括以下步骤：步骤1：对多机电力系统被控对象进行数学建模；其具体实现包括以下子步骤：步骤1.1：首先从力学的角度推导同步发电机转子运动方程，即经典的摇摆方程；步骤1.2：接着对同步电机转子运动方程中涉及的发电机输出功率进行推导与求解；步骤1.3：最后从电路的角度对同步发电机的励磁绕组电磁动态方程进行分析与求解，得到多机电力系统第i台同步发电机的数学模型表达式为：$\left\{\begin{array}{c}\frac{d{\delta }_i}{\mathrm{dt}}=({\omega }_i-1){\omega }_{0i}\\ \frac{d{\omega }_i}{\mathrm{dt}}=\frac{1}{T_{\mathrm{Ji}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}})-\frac{D_i}{T_{\mathrm{Ji}}}({\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}}{{\Sigma }_{\underset{j\ne i}{j=1}}^n{E}_{\mathrm{qj}}(G_{\mathrm{ij}}\cos {\delta }_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\delta }_{\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}\Delta {\stackrel{·}{\delta }}_i=\Delta {\omega }_i\\ \Delta {\stackrel{·}{\omega }}_i=-\frac{{\omega }_{0i}}{T_{\mathrm{Ji}}}\Delta P_{\mathrm{ei}}-\frac{D_i}{T_{\mathrm{Ji}}}\Delta {\omega }_i+\frac{{\omega }_{0i}}{T_{\mathrm{Ji}}}\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：对发电机输出的电磁功率进行线性化，即：$\Delta P_{\mathrm{ei}}=E_{\mathrm{qi}}\underset{j\ne i}{\overset{3}{\underset{j=1}{\Sigma }}}{E}_{\mathrm{qj}}{B}_{\mathrm{ij}}\cos ({\delta }_i^0-{\delta }_j^0)(\Delta {\delta }_i-\Delta {\delta }_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{{\omega }_{0i}}{T_{j1}}& 0& 0\\ 0& 0& 0\\ 0& \frac{{\omega }_{0i}}{T_{j2}}& 0\\ 0& 0& 0\\ 0& 0& \frac{{\omega }_{0i}}{T_{j3}}\end{array}\right]$$\begin{array}{ccc}{k}_{\mathrm{ij}}=\frac{{\omega }_{0i}{E}_i{E}_j{B}_{\mathrm{ij}}\cos ({\delta }_i\left(0\right)-{\delta }_j\left(0\right))}{T_{\mathrm{ji}}}& i=\mathrm{1,2,3}& i\ne j\end{array};$步骤3：规定二次型性能指标，首先选用二次型性能指标作为系统的性能指标，即$J={\int }_0^{\infty }\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{·}{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{·}{X}\left(t\right)=(A+B{K}_d)X\left(t\right)---\left(7\right)$其解为：$\stackrel{·}{X}\left(t\right)=\phi \left(t\right)X_0---\left(8\right)$其中$\phi \left(t\right)=e^{(A+B{K}_d)t}.;$将式(8)代入式(4)，则原系统性能指标为：$J=X_0^T\left[{\int }_0^{\infty }{\phi }^T\left(t\right)(Q+K_d^{T}R{K}_d)\phi \left(t\right)\mathrm{dt}\right]X_0=X_0^{T}P{X}_0$其中$P={\int }_0^{\infty }{\phi }^T\left(t\right)(Q+K_d^{T}R{K}_d)\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}{\Sigma }}\underset{j=1}{\overset{n}{\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-\theta )K_d^{\left(i\right)}+\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。