[發(fā)明專利]一種多機電力系統(tǒng)穩(wěn)定運行的分散協(xié)調(diào)控制方法

有效

• 申請?zhí)枺?/td>: 201410455919.X
• 申請日: 2014-09-09
• 公開（公告）號: CN104158191A
• 公開（公告）日: 2014-11-19
• 發(fā)明（設(shè)計）人: 來金鋼;周洪;胡文山;宋春德;路曉慶;高星冉;江朝強
• 申請（專利權(quán)）人: 武漢大學(xué)
• 主分類號: H02J3/00
• 分類號: H02J3/00
• 代理公司: 武漢科皓知識產(chǎn)權(quán)代理事務(wù)所(特殊普通合伙) 42222
• 代理人: 薛玲

1.一種多機電力系統(tǒng)穩(wěn)定運行的分散協(xié)調(diào)控制方法，其特征在于，包括以下步驟：

步驟1：對多機電力系統(tǒng)被控對象進行數(shù)學(xué)建模；其具體實現(xiàn)包括以下子步驟：

步驟1.1：首先從力學(xué)的角度推導(dǎo)同步發(fā)電機轉(zhuǎn)子運動方程，即經(jīng)典的搖擺方程；

步驟1.2：接著對同步電機轉(zhuǎn)子運動方程中涉及的發(fā)電機輸出功率進行推導(dǎo)與求解；

步驟1.3：最后從電路的角度對同步發(fā)電機的勵磁繞組電磁動態(tài)方程進行分析與求解，得到多機電力系統(tǒng)第i臺同步發(fā)電機的數(shù)學(xué)模型表達式為：

dδidt=(ωi-1)ω0idωidt=1TJi(Pmi-Pei)-DiTJi(ωi-1)dE′qidt=1T′d0i[Efi-Eqi′-(xdi-xdi′)idi]]]>

多機電力系統(tǒng)：Pei=Eqi2Gii+EqiΣj=1j≠inEqj(Gijcosδij+Bijsinδij)i={1,2,···n};]]>式中，δ_i為第i臺發(fā)電機轉(zhuǎn)子運行角，即q軸與同步參考軸S間的夾角，單位rad，ω_i表示轉(zhuǎn)子運動角速度(標幺值)，ω_0i為轉(zhuǎn)子運動角速度的初始值，單位rad/s，T_Ji為轉(zhuǎn)子轉(zhuǎn)動慣量時間常數(shù)，單位s，P_mi表示第i臺發(fā)電機輸入的機械功率(標幺值)，P_ei表示第i臺發(fā)電機輸出的電磁功率(標幺值)，D_i表示阻尼系數(shù)(標幺值)，E’_qi表示第i臺發(fā)電機暫態(tài)電動勢(標幺值)，E’_fi表示穩(wěn)態(tài)磁鏈在d軸定子側(cè)產(chǎn)生的電勢(標幺值)，i_di表示第i臺發(fā)電機輸電流瞬態(tài)值(標幺值)，x_di、x’_di分別為第i臺發(fā)電機直軸同步電抗與瞬態(tài)電抗(標幺值)；且I＝Y(jié)U，Y＝G+jB，為電力系統(tǒng)簡化導(dǎo)納矩陣，G_ii＝Y(jié)_iisinα_ii和B_ii＝Y(jié)_iicosα_ii為第i節(jié)點的自電導(dǎo)，G_ij＝Y(jié)_ijsinα_ii和B_ij＝Y(jié)_ijcosα_ii為第i節(jié)點與第j節(jié)點間的互電導(dǎo)及互導(dǎo)納；

步驟2：對多機電力系統(tǒng)數(shù)學(xué)模型進行非線性系統(tǒng)線性化，其具體實現(xiàn)包括以下子步驟：

步驟2.1：首先對其轉(zhuǎn)子運動方程進行線性化，即：

Δδ·i=ΔωiΔω·i=-ω0iTJiΔPei-DiTJiΔωi+ω0iTJiΔPmi---(1)]]>

其中，δ_i為第i臺發(fā)電機轉(zhuǎn)子運行角，單位rad，ω_i為發(fā)電機轉(zhuǎn)子運動角速度(標幺值)，ω_0i為第i臺發(fā)電機轉(zhuǎn)子運動角速度的初始值，單位rad/s，T_Ji為第i臺發(fā)電機轉(zhuǎn)子轉(zhuǎn)動慣量時間常數(shù)，單位s，P_mi表示第i臺發(fā)電機輸入的機械功率(標幺值)，P_ei表示第i臺發(fā)電機輸出的電磁功率(標幺值)，D_i為第i臺發(fā)電機阻尼系數(shù)(標幺值)；

步驟2.2：對發(fā)電機輸出的電磁功率進行線性化，即：

ΔPei=EqiΣj=1j≠i3EqjBijcos(δi0-δj0)(Δδi-Δδj)---(2)]]>

這里將發(fā)電機輸入的機械功率的偏差ΔP_mi作為控制向量U，發(fā)電機轉(zhuǎn)子轉(zhuǎn)動時功角差Δδ_i、轉(zhuǎn)速差Δω_i作為狀態(tài)向量X，將式(2)代入式(1)，整理得：

X·=AX+BU---(3)]]>

其中狀態(tài)向量X、控制向量U：

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

且有

A=010000-k12-k13-D1Tj1k120k130000100k210-k21-k23-D2Tj2k230000001k310k320-k31-k32-D3Tj3B=000ω0iTj1000000ω0iTj2000000ω0iTj3]]>

kij=ω0iEiEjBijcos(δi(0)-δj(0))Tjii=1,2,3i≠j;]]>

步驟3：規(guī)定二次型性能指標，首先選用二次型性能指標作為系統(tǒng)的性能指標，即

J=∫0∞12(XT(t)QX(t)+UT(t)RU(t))dt---(4)]]>

其中，X為系統(tǒng)的狀態(tài)向量，這里X被選取為能夠描述電力系統(tǒng)動態(tài)運行過程中變量的偏差，U為控制向量，Q為正定或半正定狀態(tài)權(quán)矩陣，R為正定控制量權(quán)矩陣；

步驟4：對于多機電力系統(tǒng)的各子系統(tǒng)，按子系統(tǒng)狀態(tài)向量反饋進行最優(yōu)分散協(xié)調(diào)控制的設(shè)計；其具體實現(xiàn)包括以下子步驟：

步驟4.1：對于多機電力系統(tǒng)的各子系統(tǒng)，設(shè)計具有分散控制結(jié)構(gòu)的控制器，即各子系統(tǒng)控制器僅反饋本子系統(tǒng)的狀態(tài)向量，使全系統(tǒng)的二次型性能指標達到最小；

對于多機電力系統(tǒng)的各子系統(tǒng)：

X·(t)=AX(t)+BU(t)X(0)=X0---(5)]]>

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

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

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

X·(t)=(A+BKd)X(t)---(7)]]>

其解為：

X·(t)=φ(t)X0---(8)]]>

其中φ(t)=e(A+BKd)t.;]]>

將式(8)代入式(4)，則原系統(tǒng)性能指標為：

J=X0T[∫0∞φT(t)(Q+KdTRKd)φ(t)dt]X0=X0TPX0]]>

其中

P=∫0∞φT(t)(Q+KdTRKd)φ(t)dt]]>

為矩陣方程

P(A+BKd)+(A+BKd)TP+Q+KdTRKd=0---(9)]]>

的解，且由式(9)矩陣方程可知P^T也為方程的解，因此，P為對稱矩陣，有P^T＝P；

步驟4.2：求解最優(yōu)控制規(guī)律，即使性能指標函數(shù)J達到極小值對系統(tǒng)進行線性最優(yōu)控制，即使二次型性能指標函數(shù)J達到極小值；

當矩陣A、B維數(shù)不同時，有下式成立

B^TA＝tr(AB^T)

其中tr()表示矩陣的跡，即矩陣對角線元素之和；

則性能指標函數(shù)J可改寫成：

J=tr(PX0X0T)]]>

假設(shè)初始狀態(tài)變量X₀為均勻分布與n維單位球面上的隨機向量，則上式又可寫成：

J＝tr(P)

令

G(P,Kd)=P(A+BKd)+(A+BKd)TP+Q+KdTRKd---(10)]]>

則約束條件式(10)為

G(P,K_d)＝0

做拉格朗日函數(shù)

L=trP+Σi=jnΣj=1nvijgij---(11)]]>

式中g(shù)_ij為G(P,K)的第i行，第j列元素；g_ij為對應(yīng)元素的拉格朗日乘子；

根據(jù)矩陣的跡定義，式(11)可寫成

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

步驟4.3：最后利用矩陣跡的基本運算法則推導(dǎo)出如下方程組：

RKdVd+(BTPV)d=0P(A+BKd)+(A+BKd)TP+Q+KdTRKd=0(A+BKd)V+V(A+BKd)T+I=0---(13)]]>

(13)式為Levine-Athans方程組，對該方程組求解得到反饋增益K_d即為最優(yōu)分散協(xié)調(diào)控制規(guī)律；

步驟5：用直接迭代法求解最優(yōu)分散協(xié)調(diào)控制規(guī)律，其具體實現(xiàn)包括以下子步驟：

步驟5.1：對于式(13)第一個方程式求解可以得到：

步驟5.2：選擇初始K_d⁰值，使初始穩(wěn)定；用Matlab求解矩陣A的特征值，發(fā)現(xiàn)所有特征根具有負的實部，因此矩陣A初始穩(wěn)定，此時取初始K_d⁰為零矩陣；

步驟5.3：給定系統(tǒng)允許誤差ε，取迭代初始步長θ，令i＝0，

步驟5.4：由式P(i)(A+BKd(i))+(A+BKd(i))TP(i)+Q+Kd(i)TRKd(i)=0]]>求得P⁽ⁱ⁾；由式(A+BKd(i))V(i)+V(i)(A+BKd(i))T+I=0]]>求得V⁽ⁱ⁾；

步驟5.5：計算Kd0(i+1)=-R-1(BTP(i)V(i))dVd(i)-1;]]>

步驟5.6：判斷迭代精度，若則停止迭代，為所求解，否則進入下一步；

步驟5.7：計算Kd(i+1)=(1-θ)Kd(i)+θKd0(i+1);]]>

步驟5.8：由式P(i+1)(A+BKd(i+1))+(A+BKd(i+1))TP(i+1)+Q+Kd(i+1)TRKd(i+1)=0]]>求得P⁽ⁱ⁺¹⁾；

步驟5.9：判斷trP⁽ⁱ⁺¹⁾＜trP⁽ⁱ⁾，若成立，則令θ變?yōu)?.2θ，i變?yōu)閕+1，返回到步驟5.4；否則，θ變?yōu)?.5θ，返回到步驟5.8。