葉偉琴，戚志東，田家欣，孫成碩

(南京理工大學自動化學院，江蘇南京 210094)

1 引言

質子交換膜燃料電池(proton exchange membrane fuel cell，PEMFC)是一種清潔、高效的氫燃料電池，具有啟動溫度低、功率密度高、響應速度快、效率高等突出優(yōu)點[1].PEMFC發(fā)電過程包括流動、傳熱、傳質和電化學反應等多種物理化學現(xiàn)象，具有典型的分數(shù)階特性[2-3].建立準確的電特性模型是對其進行電特性控制研究的基礎.

近年來，不少研究人員對PEMFC的電特性建模作出了相應探索，建模方法分為機理建模和辨識建模，機理建模過程復雜、模型參數(shù)眾多且獲取難度大，因此考慮辨識建模.然而，經(jīng)典的辨識方法在處理多變量系統(tǒng)辨識方面具有局限性和不足，考慮到辨識算法的實用性、計算量、收斂性等問題，子空間辨識方法是最合適的選擇[4].如Sami等人[5]采用子空間辨識方法建立了PEMFC/超級電容混合供電系統(tǒng)的狀態(tài)空間模型，用于研發(fā)一種動態(tài)響應速度快、可持續(xù)發(fā)電的高效電源管理單元.徐夏吟[6]采用多變量輸出誤差狀態(tài)方 程(multivariable output error state space，MOESP)子空間辨識算法建立了PEMFC溫度與冷卻水流量之間的狀態(tài)空間模型，該模型可以很好地描述PEMFC的溫度動態(tài)特性.

研究表明，PEMFC發(fā)電過程中氣體擴散[7]、對流擴散[8]、熱量傳導[9]以及電化學反應等動態(tài)過程均存在分數(shù)階特性.為此，一些研究人員在建模方法中加入分數(shù)階微積分理論，以求更加準確地描述PEMFC的動態(tài)特性.Taleb[8]建立了PEMFC分數(shù)階等效電路模型.Leng等人[10]和Qi等人[9]采用分數(shù)階微分方程描述PEMFC反應氣體擴散的動態(tài)過程，能更加精確的刻畫PEMFC內部的氣體擴散現(xiàn)象.胡聰[11]采用分數(shù)階子空間辨識方法建立了空冷型PEMFC的分數(shù)階狀態(tài)空間模型(fractional order state space，F(xiàn)OSS)，但沒有采用全局優(yōu)化算法對辨識算法中較多的未知參數(shù)進行尋優(yōu).戈衛(wèi)平[12]采用時域分數(shù)階子空間辨識方法建立了PEMFC分數(shù)階電特性模型，該辨識方法得到的模型能夠較好的描述PEMFC的電特性.但上述方法都是基于時域，計算時域信號的分數(shù)階導數(shù)時，存在較大的計算復雜度.目前用頻域分數(shù)階子空間辨識建立PEMFC電特性模型的研究較少.本文考慮將計算從時域轉換到頻域，時域中的微分將在頻域中轉化為乘積的形式.此外，只要覆蓋系統(tǒng)的帶寬，頻域辨識算法所需的數(shù)據(jù)量也較少，這些優(yōu)點使得頻域子空間辨識方法在分數(shù)階系統(tǒng)的辨識中具有更大的優(yōu)勢.

本文考慮到質子交換膜燃料電池PEMFC系統(tǒng)發(fā)電過程中的分數(shù)階特性，將分數(shù)階微積分理論與頻域子空間辨識結合，提出一種頻域分數(shù)階子空間辨識方法建立PEMFC的FOSS.為了得到輸入輸出的頻率響應數(shù)據(jù)，采用隨機多頻正弦激勵信號對時域采集的信號進行處理;構造實、虛部矩陣，通過RQ分解、SVD分解以及最小二乘法求解系統(tǒng)系數(shù)矩陣.考慮到頻域分數(shù)階子空間辨識過程中，存在多個未知參數(shù)，本文將粒子群算法(particle swarm optimizatio，PSO)作為主線，加入遺傳算法(genetic algorithm，GA)中的選擇、交叉和變異操作，提出了一種遺傳-粒子群GAPSO混合優(yōu)化算法對參數(shù)進行尋優(yōu)，以獲得更加準確的PEMFC的電特性FOSS模型.

2 PEMFC

PEMFC是一個復雜系統(tǒng)，其發(fā)電過程包含有氣體擴散、質子傳輸、電化學以及熱量耗散等物理、化學過程.空冷型PEMFC測控系統(tǒng)的結構如圖1所示.氫氣進入PEMFC電堆與氧氣發(fā)生還原反應生成水、熱量和電能;風扇負責帶來電化學反應所需的氧氣并對電堆進行散熱;排氣閥負責排出水和多余的氫氣;測控系統(tǒng)通過采集負載電流、氫氣流量、氫氣壓力、輸出電壓、輸出功率、電堆溫度、環(huán)境溫度等變量的電信號傳送給控制器，進而調節(jié)風扇電壓、氫氣流量對PEMFC的發(fā)電性能進行調控.

圖1 空冷型PEMFC測控系統(tǒng)結構圖Fig.1 Structure diagram of air cooled PEMFC measurement and control system

PEMFC是一個典型的多輸入多輸出系統(tǒng)，影響PEMFC發(fā)電過程的因素有氣體流量、氣體壓力、堆溫、濕度和負載電流等，為了減小建模復雜度，本文采用文獻[13]中的典型相關分析法[13]選取與PEMFC電特性相關性最大的輸入變量，并采用相關性分析去除冗余輸入變量.最終選取氫氣流量和負載電流作為輸入變量，輸出電壓和功率作為輸出變量.

3 頻域分數(shù)階子空間辨識方法

子空間辨識方法可以分為時域和頻域，但是在計算時域信號的分數(shù)階導數(shù)時，存在較大的計算復雜度.而頻域辨識方法能夠將時域中的微分將在頻域中轉化為乘積的形式，此外，只要覆蓋系統(tǒng)的帶寬，頻域辨識算法所需的數(shù)據(jù)量也較少.

3.1 問題描述

線性時不變分數(shù)階系統(tǒng)的分數(shù)階狀態(tài)空間方程為

其 中:A ∈Rn×n，B ∈Rn×m，C ∈Rl×n，D ∈Rl×m為常數(shù)未知矩陣，需要通過辨識得到，α是微積分階次，Dα為微積分的操作算子，n為系統(tǒng)矩陣A維數(shù)，u(t)∈Rm為輸入變量，y(t)∈Rl為輸出變量，m為輸入變量個數(shù)，l為輸出變量個數(shù).假定系統(tǒng)是穩(wěn)定且可觀測的，系統(tǒng)系數(shù)矩陣A的維數(shù)已知.

本文選取Caputo 定義作為理論推導的基礎[14]，Caputo定義的分數(shù)階微積分表達式如下:

其中:a和t表示微積分的上下限，Γ(·)為Gamma函數(shù)，m-1 ≤α ＜m，m ∈N，其Laplace變換為

借助分數(shù)階微積分的Laplace變換，可以得到系數(shù)矩陣與傳遞函數(shù)的關系如下:

其中n為系統(tǒng)矩陣A維數(shù).

3.2 頻域分數(shù)階子空間辨識

3.2.1 激勵信號

在進行頻域子空間辨識之前，首先需要從時域采樣的數(shù)據(jù)中獲得頻率響應函數(shù)(frequency response function，F(xiàn)RF).一般情況下，用于系統(tǒng)辨識的激勵信號需要滿足一些激勵條件，即用于辨識的信號需要足夠的帶寬相對于被辨識的系統(tǒng)，可用的輸入信號有白噪聲、偽隨機二進制序列、多頻正弦等[15].對于FRF[16]的測量，需要選擇合適的激勵信號，除了具有足夠的頻譜寬度外，還應該具有較小的峰值因子/時間因子[16].本文采用的是多頻正弦激勵信號.多頻正弦是一組正弦諧波信號的和，其表達式如下:

其中:Ts是采樣周期，f0是基頻，根據(jù)相位可以選擇多頻正弦的形式，不同的多頻正弦形式φk不同，本文選擇隨機多頻正弦:φk～U(0，2π)，U(0，2π)表示在區(qū)間(0，2π)內均勻分布.

當獲得輸入輸出數(shù)據(jù)后，使用以下公式直接從時域信號中獲得FRF:

其中Y(jωk)可以通過離散傅里葉變換獲得，公式如下:

其中N表示采樣點個數(shù)，U(jωk)通過同樣的方式定義.由于這些直接使用的時域信號不是穩(wěn)態(tài)響應，因此這種方法不能產生可靠的頻域數(shù)據(jù).

為了獲得準確的無頻譜泄漏誤差的FRF，需要對系統(tǒng)進行特殊的激勵信號處理.本文采用的方法是利用整數(shù)周期的u(t)，如圖2所示，這表示整個激勵信號u(tall)包括幾個完整的信號，即u(t)，u[1](t)，u[2](t)，···，u[i](t).然后，每一個輸入周期都對應一個輸出周期，最后一個輸出周期y[i](t)對應的u[i](t)可以看作是相對準確的穩(wěn)態(tài)響應.最后，利用時域中的這部分信號，通過計算式(5)-(6)得到FRF.

圖2 周期性多頻正弦激勵信號示意圖Fig.2 Schematic diagram of periodic multi-frequency sinusoidal excitation signal

利用設計的周期性激勵信號得到輸入輸出數(shù)據(jù)的FRF后，開始進行頻域分數(shù)階子空間辨識過程，最后期望得到一個FOSS模型.

3.2.2 系統(tǒng)系數(shù)矩陣辨識

1) 將時域模型轉換到頻域模型.

首先，通過Laplace變換，將分數(shù)階系統(tǒng)模型轉換到s域，令sjω，可得

其中q是一個輔助階次，表示擴展可觀測矩陣的行數(shù)，且滿足條件q≥n，其矩陣構成為

且Oq為下三角Toeplitz矩陣，其矩陣構成為

對于式(8)中所有的ω，都是使用不同采樣頻率ωk下的采樣數(shù)據(jù)，其中k1，2，···，M，這M個向量可以被合并為

其中G表示為

Gk表示為

引入術語Gre，定義如下:

3) 對R22進行SVD分解

其中Σ+的對角線元素是R22的非零奇異值，且Σ+的選擇取決于系統(tǒng)矩陣A的維數(shù).

5) 計算矩陣A和C

其中符號?表示矩陣的廣義逆.

6) 用最小二乘法計算矩陣B和D

其中:?是Kronecker積，vec(·)表示將矩陣的每一列合并排成一列，且有

3.2.3 未知參數(shù)優(yōu)化

在上述辨識過程中，同元分數(shù)階次α(0＜α ＜2)、輔助階次q(2＜q ＜10)以及頻域采樣點數(shù)M(2＜M＜20)未知的，因此需采用優(yōu)化算法對參數(shù)進行優(yōu)化.本文將GA和PSO算法結合起來，提出GA-PSO混合優(yōu)化算法，將PSO算法作為主線，其中加入GA中的選擇、交叉和變異操作，以進一步提高種群中個體的自適應調整搜索方向、增強全局尋優(yōu)的能力.GA-PSO混合優(yōu)化算法具體步驟如下:

1) 初始化GA-PSO混合優(yōu)化算法中的粒子種群和參數(shù);

2) 根據(jù)構造的適應度函數(shù)，計算每個粒子的適應度值，將每個粒子作為粒子組中當前種群的局部最優(yōu)粒子，記為Ppbest，將適應度最好的粒子作為全局種群中的全局最優(yōu)粒子，記為Pgbest;

3) 設置當前迭代次數(shù)為g，根據(jù)以下公式更新粒子的速度和位置;

4) 從均勻分布U(0，1)中選取隨機數(shù)，判斷其是否大于預設的閾值，若滿足條件，則對種群進行編碼，并執(zhí)行種群的選擇、交叉、變異操作，否則，將粒子的速度歸零，位置保持不變，跳過該步驟;

5) 重新計算粒子的適應度值，將當前粒子組中適應度最好的粒子作為局部最優(yōu)，記為Gpbest，將粒子群中適應度最好的粒子作為全局最優(yōu)，記為Ggbest;

6) 比較Ppbest和Gpbest選出局部最優(yōu)粒子記為Gpbest，比較Pgbest和Ggbest選出全局最優(yōu)粒子記為Ggbest;

7) 根據(jù)式(25)-(26)更新粒子的速度和位置，判斷是否滿足算法終止條件，即當前迭代次數(shù)達到最大進化代數(shù)或粒子群的適應度值滿足優(yōu)化目標，若未滿足終止條件則令迭代次數(shù)gg+1并返回步驟2)，直到滿足終止條件，輸出粒子群的全局最優(yōu)解Ggbest即為待優(yōu)化問題的最優(yōu)解.

為了對上述提出的GA-PSO混合優(yōu)化算法的尋優(yōu)性能進行評估，與GA、PSO算法進行尋優(yōu)性能的對比仿真測試.選取4個常用的標準測試函數(shù)作為優(yōu)化對象，其表達式如下:

1) Sphere函數(shù)，如式(27)所示:

2) Ackley函數(shù)，如式(28)所示:

3) Rastrigin函數(shù)，如式(29)所示:

4) Griewank函數(shù)，如式(30)所示:

其中Sphere函數(shù)是簡單的單峰值函數(shù)，而Ackley函數(shù)、Rastrigin函數(shù)以及Griewank函數(shù)是含有許多局部最優(yōu)值的多峰值函數(shù)，可以用來測試優(yōu)化算法的全局收斂性能和尋優(yōu)精度.

標準測試函數(shù)的優(yōu)化參數(shù)設置:維數(shù)為10，種群規(guī)模為100，適應度函數(shù)為測試函數(shù)本身，最優(yōu)值為0，最大進化代數(shù)為500，學習因子c1的變化范圍為[1，2.5]，學習因子c2的變化范圍為[1.5，2.75]，慣性權重的變化范圍為[0.4，0.9].Sphere函數(shù)的取值范圍為[-10，10]，Ackley函數(shù)的取值范圍為[-40，40]，Rastrigin函數(shù)的取值范圍為[-5.12，5.12]，Griewank函數(shù)的取值范圍為[-600，600].3種優(yōu)化算法對4種標準測試函數(shù)的尋優(yōu)曲線圖如圖3所示.

圖3 適應度函數(shù)的尋優(yōu)曲線Fig.3 Optimization curves of fitness functions

如圖3，對于單峰值的Sphere函數(shù)，GA的搜索速度較慢，PSO的收斂精度不足，而GA-PSO混合優(yōu)化算法的收斂速度和收斂精度都很高;對于存在多個局部極值的多峰函數(shù)Ackley 函數(shù)、Rastrigin 函數(shù)以及Griewank函數(shù)，PSO尋優(yōu)時很容易陷入局部最優(yōu)解，收斂精度和尋優(yōu)成功率都很低;而GA-PSO混合優(yōu)化算法的尋優(yōu)成功率很高，這是因為該算法含有交叉、變異操作，能夠幫助粒子跳出局部最優(yōu)解，并收斂于全局最優(yōu)解.

綜上所述，分數(shù)階系統(tǒng)的FOSS模型辨識的流程圖如圖4所示.

圖4 分數(shù)階系統(tǒng)的FOSS模型辨識流程圖Fig.4 FOSS model identification flow chart of fractional order system

設置同元分數(shù)階次、輔助階次、頻域采樣點數(shù)的初始值，得到系統(tǒng)的系數(shù)矩陣后，將預測輸出與實際輸出的差值作為目標函數(shù)，判斷其是否小于等于容許誤差，即≤ε;若滿足上述條件，則辨識結束;若不滿足則需用優(yōu)化算法對分數(shù)階次進行尋優(yōu).定義目標函數(shù)如下:

轉化為一個非線性優(yōu)化問題:

4 PEMFC電特性FOSS模型

4.1 采集時域輸入/輸出數(shù)據(jù)，設計輸入激勵信號

為PEMFC的輸入變量負載電流和氫氣流量設計周期性多頻正弦信號作為激勵信號.輸入信號為隨機相位的多頻正弦，系統(tǒng)由5個周期的輸入信號激勵，選擇采樣頻率fs5 Hz，頻譜分辨率f00.05 Hz，則采樣間隔Ts1/fs0.2 s，測量采集輸入/輸出信號N100個點對應一個激勵周期.輸入u(t)和輸出y(t)信號曲線如圖5-6所示.

圖5 負載電流和氫氣流量的周期性多頻正弦激勵信號Fig.5 Periodic multi-frequency sinusoidal excitation signal of load current and hydrogen flow

4.2 獲取頻域數(shù)據(jù)，計算FRF

對最后一個周期的輸入/輸出時域信號按照式(3)作離散傅里葉變換，得到頻域數(shù)據(jù)U(jωk)和Y(jωk).然后按照式(6)求得系統(tǒng)的FRF，即G(jωk)矩陣.

圖6 PEMFC系統(tǒng)輸出電壓和功率Fig.6 Output voltage and power curve of PEMFC system

4.3 頻域分數(shù)階子空間辨識

設置初始值，迭代次數(shù)i0，同元分數(shù)階次α0.9，輔助階次q5，頻域采樣點數(shù)M7.系統(tǒng)奇異值分布的柱形圖如圖7所示，系統(tǒng)階次取n2;得到系統(tǒng)的系數(shù)矩陣后，將預測輸出與實際輸出的差值作為目標函數(shù)，判斷其是否小于等于容許誤差，即≤ε;若滿足上述條件，則辨識結束，輸出模型參數(shù){α，q，M，A，B，C，D};否則令迭代次數(shù)ii+1，重復頻域分數(shù)階子空間辨識步驟，直到結束條件滿足為止.

圖7 PEMFC系統(tǒng)奇異值分布圖Fig.7 Singular value distribution of PEMFC system

4.4 GA-PSO優(yōu)化過程

對GA-PSO混合優(yōu)化算法中的粒子種群和參數(shù)初始化，其中種群規(guī)模N20，慣性權重的上下限值ωmax0.9，ωmin0.4，最大迭代次數(shù)Gmax100，交叉概率pc0.5，變異概率pm0.05，粒子速度的上下限值Vmax1，Vmin-1，粒子位置的上下限值Xmax10，Xmin-10，學習因子c1和c2的初始值c1s1，c2s1.5，終止值c1e2.5，c2e2.75.另外，適應度函數(shù)為模型實際輸出與預測輸出的均方根誤差(root mean square error，RMSE)，遺傳算法的選擇操作采用輪盤賭法，學習因子的取值方法采用非對稱線性變化取值法，并選用線性慣性權重遞減策略.

5 仿真結果

為了驗證頻域分數(shù)階子空間辨識算法及改進優(yōu)化算法的有效性，用輸入/輸出數(shù)據(jù)進行辨識和驗證.以式(22)為目標函數(shù)，待優(yōu)化參數(shù)為同元分數(shù)階次α(0＜α ＜2)、輔助階次q(2＜q ＜10)以及頻域采樣點數(shù)M(2＜M ＜20).經(jīng)過GA-PSO的尋優(yōu)和頻域分數(shù)階子空間辨識過程，得到優(yōu)化參數(shù)和PEMFC電特性FOSS模型的值如下:

基于GA 優(yōu)化[17]、PSO 優(yōu)化[18]、GA-PSO 優(yōu)化的PEMFC頻域分數(shù)階子空間辨識模型的辨識結果如圖8-11所示.

從圖8-9中可以看出，帶GA優(yōu)化、帶PSO優(yōu)化與帶GA-PSO優(yōu)化的頻域分數(shù)階子空間辨識方法辨識得到的輸出均能夠跟隨實測數(shù)據(jù)的變化趨勢，但都存在一定的誤差.通過圖10-11可以看出，帶GA-PSO優(yōu)化的頻域分數(shù)階子空間辨識的輸出電壓誤差絕對值小于0.25 V，輸出功率誤差絕對值小于0.63 W.與帶GA優(yōu)化、帶PSO優(yōu)化的頻域分數(shù)階子空間辨識相比，誤差較小，精度更高.帶GA-PSO優(yōu)化的頻域分數(shù)階子空間辨識能夠更精確的描述PEMFC的電特性.表1為頻域分數(shù)階子空間辨識與時域分數(shù)階子空間辨識[12]的辨識時間比較.

表1 頻域與時域分數(shù)階子空間辨識的辨識時間比較Table 1 Time comparison of fractional subspace identification in time domain and frequency domain

圖8 輸出電壓Fig.8 Output voltage

圖9 輸出功率Fig.9 Output power

圖10 輸出電壓誤差絕對值Fig.10 Absolute value of output voltage error

圖11 輸出功率誤差絕對值Fig.11 Absolute value of output power error

從表1可以看出頻域分數(shù)階子空間辨識的辨識時間明顯比時域分數(shù)階子空間辨識短.當辨識數(shù)據(jù)組數(shù)取100，200，300時，頻域分數(shù)階子空間辨識的辨識時間大約是時域分數(shù)階子空間辨識三分之二;當辨識數(shù)據(jù)組數(shù)取400，500時，頻域分數(shù)階子空間辨識的辨識時間大約是時域分數(shù)階子空間辨識二分之一.由此可以看出，頻域分數(shù)階子空間辨識能夠大大減小計算復雜度.

6 結論

本文針對PEMFC發(fā)電過程中存在的分數(shù)階特性，建立了PEMFC電特性的FOSS模型.在頻域子空間辨識方法中引入分數(shù)階微積分理論，提出一種頻域分數(shù)階子空間辨識方法.對影響PEMFC系統(tǒng)電特性的因素進行分析，選取負載電流、氫氣流量為模型輸入，輸出電壓和功率為模型輸出，采用上述頻域分數(shù)階子空間辨識方法得到系統(tǒng)的系數(shù)矩陣，并在辨識過程中采用一種GA-PSO混合優(yōu)化算法對辨識算法中的未知參數(shù)進行尋優(yōu).仿真結果驗證了算法的有效性，頻域分數(shù)階子空間辨識得到的模型能夠準確地描述PEMFC的電特性變化過程.