古志，曾云，李敏，吳一凡，侯睿，錢晶

(昆明理工大學(xué)冶金與能源工程學(xué)院， 云南 昆明 650093)

調(diào)速器PID參數(shù)是影響水力機(jī)組調(diào)節(jié)性能的關(guān)鍵因素，其合理的選擇和優(yōu)化不僅可以確保機(jī)組安全穩(wěn)定運(yùn)行[1]，還能抑制電力系統(tǒng)的低頻振蕩現(xiàn)象[2].因此，調(diào)速器PID參數(shù)優(yōu)化的研究具有重要意義.

隨著計(jì)算機(jī)技術(shù)和控制理論的不斷發(fā)展，智能優(yōu)化算法被廣泛應(yīng)用于水輪機(jī)調(diào)速器的參數(shù)優(yōu)化，例如粒子群算法(particle swarm optimization, PSO)[3]、遺傳算法[4]、引力搜索算法[5]、人工魚群算法[6]和狼群算法[7]等.然而，上述方法也存在一些不足，如：遺傳算法的性能對(duì)參數(shù)依賴性強(qiáng)；人工魚群算法的收斂速度慢；狼群算法對(duì)一種工況的理想控制方式很難滿足其他工況的控制要求.同時(shí)，這些算法還具有計(jì)算過程復(fù)雜、耗時(shí)長、程序編寫煩瑣困難等問題.

人群搜索算法(seeker optimization algorithm, SOA)是一種新型啟發(fā)式智能搜索優(yōu)化算法[8].SOA及其變體已經(jīng)被證明能解決各種優(yōu)化問題，例如優(yōu)化傳輸切換問題[9]、優(yōu)化電力系統(tǒng)無功功率的分配[10]、調(diào)整神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)和參數(shù)[11].文中正是基于SOA具有能解決系統(tǒng)控制問題和優(yōu)化PID控制參數(shù)的特點(diǎn)，將其應(yīng)用到水輪機(jī)調(diào)速系統(tǒng)中，以獲得最優(yōu)調(diào)節(jié)參數(shù).

1 水輪機(jī)調(diào)速系統(tǒng)模型

1.1 水輪機(jī)非線性模型

文中采用IEEE Working Group推薦的非線性水輪機(jī)模型.在該模型中，非線性主要通過水輪機(jī)水頭和流量的變化來體現(xiàn)，而出力采用近似線性的公式進(jìn)行計(jì)算.水輪機(jī)出力、流量、水頭的表達(dá)式為

(1)

式中：Pm為水輪機(jī)出力相對(duì)值；At為增益因子，At=1/hr(qr-qnl).其中，hr為水輪機(jī)額定水頭相對(duì)值；qr為水輪機(jī)額定流量相對(duì)值；qnl為空載流量相對(duì)值；Dt為水力阻尼因子，并網(wǎng)運(yùn)行條件下該項(xiàng)數(shù)值很小，文中忽略其影響；Δω為機(jī)組角速度偏差相對(duì)值；q為水輪機(jī)流量相對(duì)值；h為水輪機(jī)水頭相對(duì)值；h0為靜水頭相對(duì)值；y為水輪機(jī)接力器位移的相對(duì)值；fp為水頭損失系數(shù)；Δh為流量變化引起的管道末端水頭變化相對(duì)值，也稱暫態(tài)水頭相對(duì)值.

1.2 水力系統(tǒng)水擊模型

水輪機(jī)引水系統(tǒng)流量變化到水頭變化的傳遞函數(shù)為[12]

(2)

式中：Zn為管道的水力浪涌阻抗系數(shù)；Te為水力彈性時(shí)間，s；s為拉普拉斯算子.

文中采用的是單機(jī)單管簡單水力系統(tǒng)，傳遞函數(shù)為剛性水擊模型，即式(2)中n=0的情形.

1.3 調(diào)速器與電液隨動(dòng)系統(tǒng)的傳遞函數(shù)

水電站微機(jī)調(diào)速器最廣泛使用的是典型的并聯(lián)PID控制單元，電液隨動(dòng)系統(tǒng)采用的是一階慣性環(huán)節(jié)，兩者的傳遞函數(shù)為

(3)

式中：Kp，Ki，Kd分別為PID調(diào)節(jié)器的比例、積分、微分的增益系數(shù)；Td為常數(shù)；Ty為接力器的時(shí)間常數(shù).

1.4 水輪機(jī)調(diào)速系統(tǒng)仿真框架

在水輪機(jī)調(diào)速系統(tǒng)仿真研究中，采用一階發(fā)電機(jī)簡化模型來近似處理.采用Matlab/Simulink仿真平臺(tái)以兩路并聯(lián)模塊的形式建立水輪機(jī)調(diào)速系統(tǒng)的非線性仿真模型，如圖1所示.

對(duì)該仿真模型有幾點(diǎn)說明：

1) 圖1中，傳遞函數(shù)均采用線性增量化的形式，通過設(shè)置初值，使其可用于非線性情形，例如y0,q0,h0,mg0,mt0等也通過預(yù)設(shè)初值的形式來表達(dá)，其中mg0為發(fā)電機(jī)側(cè)阻力矩相對(duì)值，mt0為水輪機(jī)側(cè)動(dòng)力矩相對(duì)值.

2) 流量變化到管道末端水頭的變化傳遞函數(shù)為假分式，需分解為真分式形式才能使用Simulink模塊構(gòu)建.圖1中是按剛性水擊傳遞函數(shù)分解的，若采用彈性水力傳遞函數(shù)，只需分解為真分式形式，修改該模塊即可.

圖1 非線性水輪機(jī)調(diào)速系統(tǒng)仿真模塊

3) 時(shí)間乘以誤差值絕對(duì)積分(ITAE)是衡量優(yōu)化算法性能好壞的標(biāo)準(zhǔn)，以控制量、誤差和時(shí)間為約束條件，直接在Simulink模塊里與調(diào)速器輸入連接.

4) 頻率擾動(dòng)是通過預(yù)先設(shè)置階躍信號(hào)形式實(shí)現(xiàn)的；功率擾動(dòng)采用圖中mg0輸入階躍擾動(dòng)來模擬.

2 基于SOA的PID設(shè)計(jì)

2.1 參數(shù)的編碼

當(dāng)人群中的搜索個(gè)體數(shù)為n時(shí)，每個(gè)搜索個(gè)體的位置是由PID控制器的參數(shù)構(gòu)成，整個(gè)人群就可以表示成矩陣M[13]

(4)

2.2 搜索步長的確定

SOA模擬人類搜索行為，其搜索個(gè)體的位置更新依賴于搜索步長.搜索步長計(jì)算公式為

(5)

式中：αij為個(gè)體i在j維空間的搜索步長，其值為非負(fù)數(shù)；δij為個(gè)體i在j維空間的高斯隸屬度函數(shù)的參數(shù)變量；uij為個(gè)體i在j維空間中隨機(jī)、均勻分布在[ui,1]上的隸屬度值.

高斯隸屬度的參數(shù)變量δij可確定為

δij=abs(xmax-xmin)·ω,

(6)

ω=(ηmax-η)/ηmax,

(7)

式中：abs(x)表示為x的絕對(duì)值；xmax,xmin分別為同一種群中個(gè)體的最大函數(shù)值和最小函數(shù)值的位置；ω為慣性加權(quán)值；ηmax,η分別為算法的迭代最大次數(shù)和算法當(dāng)前迭代次數(shù).

SOA的隸屬度值u的分布由高斯隸屬函數(shù)來描述，即

(8)

式中：uA為高斯隸屬度；x為系統(tǒng)輸入變量；u，δ分別為高斯隸屬函數(shù)的參數(shù).

將式(8)采用線性化函數(shù)表示，則個(gè)體位置的隸屬度值[14]為

(9)

uij=rand(ui,1),j=1,2,…,D,

(10)

式中：i為種群中的個(gè)體序列號(hào)；Fi為種群中個(gè)體i的函數(shù)值以降序方式排列得到的序列號(hào)；ui為個(gè)體i的隸屬度值；D為種群搜索的空間維度.

2.3 搜索方向的確定

運(yùn)用隨機(jī)加權(quán)幾何平均法確定搜索方向，即[15]

dij(t)=sign(ωdij,pro+φ1dij,ego+φ2dij,alt),

(11)

式中：dij(t)為個(gè)體i在j維空間的搜索方向；sign(·)為符號(hào)函數(shù)；dij,ego(t),dij,alt(t),dij,pro(t)分別為個(gè)體i在j維空間的利己方向、利他方向和預(yù)動(dòng)方向；φ1，φ2為區(qū)間[0,1]內(nèi)的常數(shù).

2.4 個(gè)體位置的更新

搜索的步長與方向確定后，SOA還需進(jìn)行相應(yīng)的位置更新，具體的更新方式為

Δxij(t+1)=αij(t)dij(t),

(12)

xij(t+1)=xij(t)+Δxij(t+1).

(13)

3 優(yōu)化過程

3.1 目標(biāo)函數(shù)

在迭代過程中，SOA只需評(píng)價(jià)目標(biāo)函數(shù)的適應(yīng)值W來判斷種群個(gè)體或解的優(yōu)劣.選擇性能指標(biāo)較好的積分時(shí)間絕對(duì)誤差(integral time absolute error, ITAE)作為該算法的目標(biāo)函數(shù)

(14)

式中：e為系統(tǒng)誤差；n為仿真計(jì)算中的采樣點(diǎn)數(shù)；T為仿真計(jì)算的步距.

3.2 參數(shù)優(yōu)化的流程

運(yùn)用SOA優(yōu)化調(diào)節(jié)參數(shù)的步驟具體如下：

步驟1：算法初始化.設(shè)置算法初始參數(shù)，包括種群規(guī)模數(shù)n、最大迭代次數(shù)ηmax、最大隸屬度umax、最小隸屬度umin、慣性權(quán)重最大值ωmax、慣性權(quán)重最小值ωmin、最小適應(yīng)值Wmin、空間維度D；確定優(yōu)化參數(shù)的控制范圍[xmin,xmax]，并在此范圍內(nèi)初始化種群個(gè)體的位置，初始化歷史最優(yōu)目標(biāo)函數(shù)值Wbest和最優(yōu)解xij_best；設(shè)置迭代次數(shù)η=0.

步驟2：計(jì)算個(gè)體的目標(biāo)函數(shù)值.調(diào)用水力機(jī)組調(diào)速系統(tǒng)控制仿真平臺(tái)，采用智能化的搜尋策略，計(jì)算每個(gè)個(gè)體的搜索方向dij(t)和步長αij(t)，將得到的系統(tǒng)響應(yīng)按式(14)進(jìn)行目標(biāo)函數(shù)的計(jì)算.

步驟3：求出種群中最優(yōu)目標(biāo)函數(shù)值Wbest_current和對(duì)應(yīng)的個(gè)體位置xij_current.當(dāng)Wbest_current

步驟4：按照式(12)，(13)計(jì)算更新每個(gè)個(gè)體的位置.

步驟5：η=η+1；當(dāng)η<ηmax時(shí)，轉(zhuǎn)至步驟2，否則結(jié)束算法循環(huán)，得出最優(yōu)目標(biāo)函數(shù)值Wbest和最優(yōu)解xij_best.

4 仿真計(jì)算

4.1 仿真基本參數(shù)

以文中所述方法對(duì)所建立的非線性水輪機(jī)調(diào)速系統(tǒng)模型進(jìn)行仿真模擬，機(jī)組的相關(guān)參數(shù)：Td=0.02，水流慣性時(shí)間常數(shù)Tw=1.2，fp=0.01，At=1.053，Ty=0.3，qnl=0.05，h0=1，發(fā)電機(jī)機(jī)械慣性時(shí)間常量Ta=7.7，發(fā)電機(jī)自調(diào)節(jié)系數(shù)eg=5.

目前已使用的SOA參數(shù)都是事先人為給定的，對(duì)于非線性水輪機(jī)調(diào)節(jié)系統(tǒng)不具備普適性.文中經(jīng)過大量仿真試算，有以下幾點(diǎn)發(fā)現(xiàn)：

1) 非線性水輪機(jī)調(diào)節(jié)系統(tǒng)仿真中，ITAE指標(biāo)是與最小適應(yīng)值Wmin聯(lián)系起來進(jìn)入算法循環(huán)的，Wmin的設(shè)置值過大或過小都會(huì)影響SOA算法的迭代尋優(yōu)過程.取Wmin=0.1能較好地滿足系統(tǒng)調(diào)節(jié)需要.

2) 搜索方向中適當(dāng)增加利己方向權(quán)重φ1和利他方向的權(quán)重φ2有利于加快SOA搜索速度.取φ1=0.9，φ2=0.8能使得SOA在非線性水輪機(jī)調(diào)節(jié)系統(tǒng)仿真中的迭代次數(shù)較少.

仿真其他SOA的參數(shù)設(shè)置：D=3,umax=0.950 0，umin=0.011 1，ωmax=0.9，ωmin=0.1，n=100，ηmax=100，Kp∈[0, 6.5]，Ki∈[0, 4.5]，Kd∈[0, 5].

4.2 仿真結(jié)果與分析

仿真分為在剛性水擊模型下的頻率擾動(dòng)與負(fù)荷擾動(dòng).2組仿真中的仿真時(shí)間均為20 s.在第1組仿真中，施加5%的頻率擾動(dòng)，得到的目標(biāo)函數(shù)適應(yīng)度曲線如圖2所示，圖中η為程序迭代次數(shù).

從圖2可以看出，經(jīng)SOA優(yōu)化后目標(biāo)函數(shù)的適應(yīng)曲線在第28次迭代后開始趨于穩(wěn)定，這種現(xiàn)象表明，自第29次迭代起，SOA已收斂并達(dá)到了尋找PID控制最佳參數(shù)的目標(biāo).

圖2 在5%頻率擾動(dòng)下的適應(yīng)度曲線

在第2組仿真中，施加10%的孤網(wǎng)負(fù)荷擾動(dòng)，得到的目標(biāo)函數(shù)適應(yīng)度曲線如圖3所示.

從圖3可以看出，經(jīng)SOA優(yōu)化后目標(biāo)函數(shù)的適應(yīng)曲線在第24次迭代后開始趨于穩(wěn)定，這表明自第25次迭代起，SOA已收斂并達(dá)到了尋找PID控制最佳參數(shù)的目標(biāo).

圖3 在10%孤網(wǎng)負(fù)荷擾動(dòng)下的適應(yīng)度曲線

為驗(yàn)證SOA的優(yōu)化效果，在相同仿真條件下比較SOA-PID優(yōu)化方案和PSO-PID優(yōu)化方案，不同仿真方案得到的控制效果對(duì)比見圖4，5，所比較的PID控制參數(shù)和對(duì)應(yīng)的控制性能指標(biāo)見表1.

圖4 在5%頻率擾動(dòng)下的對(duì)比圖

通過圖4可以看出，同等仿真條件下，在5%頻率擾動(dòng)中，經(jīng)SOA優(yōu)化的系統(tǒng)能在8 s內(nèi)趨于穩(wěn)定，而經(jīng)PSO優(yōu)化的系統(tǒng)卻要11 s才能穩(wěn)定.同時(shí)，PSO優(yōu)化系統(tǒng)的超調(diào)量ψ比SAO優(yōu)化系統(tǒng)的超調(diào)量ψ大21.8%.通過圖5可以看出，同等仿真條件下，SOA-PID優(yōu)化方案與PSO-PID優(yōu)化方案在10%負(fù)荷擾動(dòng)下優(yōu)化效果基本一樣，均在10 s內(nèi)使系統(tǒng)穩(wěn)定，這也從側(cè)面驗(yàn)證了SOA優(yōu)化結(jié)果的正確性.從表1可以看出，在5%頻率擾動(dòng)工況下，SOA-PID優(yōu)化方案中的Kp，Ki，Kd數(shù)值均比PSO-PID優(yōu)化方案中的Kp，Ki，Kd數(shù)值小，而在10%孤網(wǎng)負(fù)荷擾動(dòng)工況下卻出現(xiàn)了相反的情況，但是SOA的ITAE指標(biāo)始終小于PSO的ITAE指標(biāo)，表明SOA具有更好的搜索功能.

圖5 在10%孤網(wǎng)負(fù)荷擾動(dòng)下的對(duì)比圖

表1 PID參數(shù)及控制性能指標(biāo)對(duì)比表

5 結(jié) 論

1) 在5%頻率擾動(dòng)下，經(jīng)SOA優(yōu)化的系統(tǒng)能在8 s內(nèi)趨于穩(wěn)定，此時(shí)系統(tǒng)的超調(diào)量ψ只有1.6%，滿足了系統(tǒng)穩(wěn)定要求，實(shí)現(xiàn)了PID參數(shù)的良好優(yōu)化.

2) 在10%負(fù)荷擾動(dòng)下，SOA優(yōu)化效果與PSO優(yōu)化效果基本相同，兩者均在10 s內(nèi)讓系統(tǒng)趨于穩(wěn)定，唯一的差別就是SOA優(yōu)化的ITAE指標(biāo)比PSO優(yōu)化的要小.

3) 驗(yàn)證了基于SOA的水輪機(jī)調(diào)速系統(tǒng)PID參數(shù)優(yōu)化方法的可行性和有效性，為以后SOA針對(duì)水電機(jī)組的改進(jìn)提供了基礎(chǔ)條件和方向.