金秀章，李 京

(華北電力大學(xué) 控制與計(jì)算機(jī)工程學(xué)院，河北 保定 071003)

1 引 言

發(fā)展經(jīng)濟(jì)的同時(shí)，環(huán)保達(dá)標(biāo)已經(jīng)成為了一項(xiàng)重要指標(biāo)[1]。近年來(lái)我國(guó)對(duì)使用火電機(jī)組的電廠SO2排放的要求[2]為小于35 mg/m3。所以，火電廠需要嚴(yán)格控制SO2排放，SO2濃度的檢測(cè)也更需要達(dá)到準(zhǔn)確、迅速以及對(duì)入口SO2濃度進(jìn)行提前控制。

目前，由于測(cè)量SO2的儀器測(cè)量時(shí)間較長(zhǎng)且火電廠的SO2測(cè)量?jī)x的安裝位置距煙氣管道較遠(yuǎn)等原因，無(wú)法及時(shí)測(cè)得SO2值，使得存在煤燃燒后要經(jīng)過(guò)一定延遲才會(huì)發(fā)現(xiàn)SO2濃度的變化，對(duì)脫硫作業(yè)以及提前噴漿控制都十分不利。因此，需要對(duì)SO2入口濃度進(jìn)行預(yù)測(cè)，進(jìn)行噴漿提前控制。

傳統(tǒng)的SO2排放預(yù)測(cè)是通過(guò)CEMS儀表進(jìn)行檢測(cè)。這些儀表價(jià)格昂貴且需要人工安裝以及硬件的維護(hù)，建設(shè)與維護(hù)費(fèi)用較高。通過(guò)本文提出的軟測(cè)量預(yù)測(cè)模型，預(yù)測(cè)精度貼近實(shí)際儀表測(cè)量值，可以節(jié)約建設(shè)成本與維護(hù)費(fèi)用。

軟測(cè)量技術(shù)是通過(guò)計(jì)算機(jī)，利用實(shí)際中易測(cè)量的輔助變量，對(duì)難測(cè)量或無(wú)法測(cè)量的主導(dǎo)變量進(jìn)行預(yù)測(cè)[3]。軟測(cè)量建模有基于機(jī)理建模法[4]、數(shù)模建模法[5]、回歸分析法[6]、灰色預(yù)測(cè)法[7]、人工神經(jīng)網(wǎng)絡(luò)建模法(artificial neural network,ANN)[8]、支持向量機(jī)建模法[9]等。其中，在支持向量機(jī)分類法基礎(chǔ)上發(fā)展出了最小二乘支持向量機(jī)(least squares support vector machines，LSSVM)，預(yù)測(cè)結(jié)果的好壞，由LSSVM模型的參數(shù)決定。本文采用粒子群算法(particle swarm optimization，PSO)進(jìn)行參數(shù)選取的優(yōu)化，提高了模型的預(yù)測(cè)精度。

2 互信息理論

在機(jī)器學(xué)習(xí)中，信息熵(information entropy)常用來(lái)描述一個(gè)事物所包含信息量的多少，因此信息熵可以成為函數(shù)的優(yōu)化目標(biāo)或者方程式中參數(shù)選取的最佳值。信息熵的定義公式為：

(1)

若兩個(gè)時(shí)間不相互獨(dú)立，滿足:

H(A,B)=H(A)+H(B)-I(A,B)

(2)

式中：I(A,B)為A與B之間的互信息量[10]，表示A變量中含有B變量信息的多少或者B變量中含有A變量信息的多少。

事件中一個(gè)隨機(jī)事件關(guān)于另一個(gè)隨機(jī)事件的信息量稱為互信息量，互信息量I(A,B)是定量地研究信息流通問(wèn)題的重要基礎(chǔ)。通常用事件A的后驗(yàn)概率p(A|B)與先驗(yàn)概率p(A)比值的對(duì)數(shù)表示互信息量，即：

(3)

同理，可以定義B對(duì)A的互信息量：

(4)

而p(AB)=p(A)p(B|A)=p(B)p(A|B)，可得：

=I(A)-I(A|B)

=I(B)-I(B|A)

=I(A)+I(B)-I(AB)

(5)

算法性能的表現(xiàn)由評(píng)價(jià)函數(shù)決定，因此在互信息選擇特征時(shí)該函數(shù)起著至關(guān)重要的作用。通常直接選取輸出變量B與輸入變量A互信息值作為評(píng)價(jià)指標(biāo)，引入BIF算法進(jìn)行評(píng)價(jià)，BIF算法第一步會(huì)將所有候選特征變量f進(jìn)行其評(píng)價(jià)函數(shù)(J(f))計(jì)算，然后依據(jù)函數(shù)值大小按降序排序，最終選擇前k個(gè)特征變量作為子集S。BIF算法的優(yōu)點(diǎn)是效率高，因此它適用于高維數(shù)據(jù)處理中。其中BIF算法的評(píng)價(jià)函數(shù)如式(6)所示。

J(fi)=I(fi;c)

(6)

式中：fi為輔助變量；c為預(yù)測(cè)模型的主導(dǎo)變量。

這種算法最終能從較多輔助變量中篩選出與主導(dǎo)變量相關(guān)性最大的輔助變量，卻未考慮到輔助變量之間的信息冗余。MIFS算法將懲罰因子引入了評(píng)價(jià)函數(shù)中，大大降低了輔助變量之間的信息冗余。MIFS算法的評(píng)價(jià)函數(shù)為：

(7)

式中：fi∈f為輔助變量；β為懲罰因子；Si∈S為已選變量。

利用式(7)的評(píng)價(jià)函數(shù)，既能篩選出與輸出變量相關(guān)性最大的輔助變量，還剔除了輔助變量之間的冗余性，達(dá)到了輔助變量之間相關(guān)性最大、信息冗余性最低的目的。β越大，則評(píng)價(jià)函數(shù)反映出的輔助變量之間的信息冗余越多。

參數(shù)β的選取很大程度上影響MIFS算法選取輔助變量，而最小冗余最大相關(guān)性(minimum redundancy and maximum relevance, mRMR)算法很好地解決了β難以確定的問(wèn)題。其評(píng)價(jià)函數(shù)為：

(8)

mRMR算法的具體步驟如下：

Step1： 初始化輔助變量中已選變量與未選變量的集合S和f。

Step2： 確定未選變量fi∈f，同時(shí)分別求得fi與主導(dǎo)變量c之間的互信息。

Step3： 確定第一輔助變量。用上述BIF算法，求得與主導(dǎo)變量之間的互信息，將輔助變量按互信息數(shù)值從大到小依次排列。

Step4： 貪心搜索。循環(huán)計(jì)算第一輔助變量fS與未選變量fi之間的互信息數(shù)值，依據(jù)評(píng)價(jià)函數(shù)J，將未選輔助變量中與第一輔助變量相關(guān)性最大的變量fi作為下一個(gè)變量；與此同時(shí)，f=f-{fi},S=S+{fi}，繼續(xù)循環(huán)，直到篩選的變量個(gè)數(shù)達(dá)到初始設(shè)置的變量個(gè)數(shù)n。

3 LSSVM理論

f(x)=WTφ(xk)+b+ek

(9)

式中：W∈Rdn；b∈R；ek∈R，k=1,2,…,n。

根據(jù)結(jié)果風(fēng)險(xiǎn)最小化原則，并引入最小二乘法，式(9)必須滿足式(10)：

(10)

式中：e為誤差；γ為正則化參數(shù)。

依據(jù)拉格朗日法，可寫(xiě)成下式：

(11)

式中αk(k=1,2,…,n)為拉格朗日因子。

消除參數(shù)W和e后，有：

(12)

式中：y=[y1,y2,…,yn];l=[1,1,…,1];α=[α1,α2,…,αn];Ωkl=φ(xk)Tφ(xl);k=1,2,…,n。

尋找滿足式(11)的α和b是LSSVM建模的目標(biāo)。最小二乘支持向量機(jī)核函數(shù)定義公式為：

k(xk,xl)=φ(xk)Tφ(xl)

k(xk,xl)是Mercer約束條件的任意對(duì)稱函數(shù)。引入核函數(shù)后的函數(shù)如式(13)。

(13)

核函數(shù)的選擇對(duì)最小二乘支持向量機(jī)的預(yù)測(cè)有較大影響，但核函數(shù)選取至今沒(méi)有好的方法及經(jīng)驗(yàn)。徑向基函數(shù)(RBF)、多項(xiàng)式函數(shù)和S函數(shù)等較為普遍使用,本文選用徑向基函數(shù)(RBF)。

4 粒子群(PSO)理論

粒子群算法源于對(duì)鳥(niǎo)群中飛行動(dòng)作的模仿[12]，較遺傳算法等尋優(yōu)算法簡(jiǎn)單，可用來(lái)優(yōu)化多種函數(shù)組合下的參數(shù)，從而尋找最佳工作參數(shù)。對(duì)于計(jì)算過(guò)程而言，使用粒子群算法首先需要初始化粒子群的各項(xiàng)參數(shù)，例如粒子的初始速度、位移及其決定的適應(yīng)值等，然后進(jìn)行迭代計(jì)算。迭代過(guò)程中，每個(gè)粒子都以2個(gè)最優(yōu)值為標(biāo)準(zhǔn)框定自己的行進(jìn)路線，1個(gè)是粒子本身所找到的目前最優(yōu)解，該最優(yōu)解將被粒子記錄；剩下的1個(gè)粒子在所有種群中尋找到最優(yōu)解。若粒子群間有相互通訊的方式，則各粒子都能尋找到自己所在的最優(yōu)位置。其算法如下：

vij(t+1)=wvij(t)+c1r1j(t)(pij(t)-Xij(t))+

c2r2j(t)(pgj(t)-Xij(t))

(14)

Xij(t+1)=Xij(t)+vij(t+1)

(15)

式中各參數(shù)的含義如下：

1)i為第i個(gè)粒子，i=1,2,…,M。M是該群體粒子的總數(shù)，本文選取M=300。

2)j為各粒子中尋優(yōu)的第j個(gè)參數(shù)。

3)w為慣性權(quán)值因子，其大小影響整體的尋優(yōu)能力，為了保證其收斂到全局最優(yōu)，避免出現(xiàn)早熟收斂，一般采用權(quán)值w在wmax與wmin之間線性遞減的方法。本文選取wmax=0.9，wmin=0.4。

4)t為此時(shí)優(yōu)化的代數(shù)。

5)vij(t)為第i個(gè)粒子在j維空間中的速度。

6)c1和c2為加速因子，本文中選取c1=1.7，c2=1.5。

7)r1j和r2j為相對(duì)獨(dú)立的兩個(gè)隨機(jī)參數(shù)，其值變化范圍在[0～1]之間。

8)pij(t)為粒子i的歷史最優(yōu)解的j維值。

9)pgj(t)=min{pij(t)}為各粒子在t時(shí)刻的歷史最優(yōu)解的j維值，即所有粒子在所優(yōu)化的第j個(gè)參數(shù)中的歷史最優(yōu)解。

10)Xij(t)為粒子i處于j維空間的位置。

根據(jù)這2個(gè)方程，進(jìn)行迭代計(jì)算動(dòng)態(tài)更新粒子的位置和速度，搜索迭代過(guò)程中粒子的各個(gè)適應(yīng)值的最優(yōu)值，根據(jù)目標(biāo)函數(shù)計(jì)算適應(yīng)度值。整個(gè)搜索過(guò)程中獲得的最佳方程解是全局最優(yōu)解。當(dāng)群體搜索到滿足最小適應(yīng)值的最優(yōu)位置或者達(dá)到迭代次數(shù)之后即可結(jié)束計(jì)算。

5 入口SO2濃度預(yù)測(cè)模型的建立及仿真

對(duì)陜西某電廠所采樣的數(shù)據(jù)進(jìn)行優(yōu)化處理后，利用互信息的方法篩選出輔助變量，并采用PSO算法對(duì)LSSVM的參數(shù)進(jìn)行尋優(yōu)，建立出SO2濃度的預(yù)測(cè)模型。

5.1 輔助變量的選擇

本文通過(guò)機(jī)理分析[13]，篩選出了36個(gè)輔助變量，其中主要有機(jī)組負(fù)荷、總煤量、總風(fēng)量、石灰石漿液濃度、石灰石漿液循環(huán)泵電流等。通過(guò)互信息(mutual information,MI)和最大相關(guān)-最小冗余(mRMR)結(jié)合的方法，以預(yù)測(cè)模型的輸出精度作為篩選標(biāo)準(zhǔn)來(lái)確定最終輔助變量的個(gè)數(shù)，建立了以入口SO2濃度為主導(dǎo)變量的預(yù)測(cè)模型。篩選出的變量有漿液循環(huán)泵電流、氧化風(fēng)機(jī)電流、空預(yù)器含氧量等18個(gè)輔助變量作為模型輸入。輔助變量及預(yù)測(cè)流程如圖1所示。

圖1 SO2入口濃度排放模型結(jié)構(gòu)圖Fig.1 Structure diagram of SO2 inlet concentration emission model

5.2 數(shù)據(jù)預(yù)處理

在分析數(shù)據(jù)的相關(guān)性方面，采用皮爾遜(Pearson)相關(guān)性[14]分析變量的相關(guān)性，得出數(shù)據(jù)的相關(guān)性非常低，不能較好地作為模型輔助變量的輸入。通過(guò)對(duì)10 000個(gè)樣本進(jìn)行分析，輔助變量與主導(dǎo)變量之間的相關(guān)性系數(shù)大都在0.5以下，很難確定出相關(guān)性較高的輔助變量。因此，本文采用互信息方法篩選輔助變量，從現(xiàn)場(chǎng)的36個(gè)輔助變量中最終篩選出18個(gè)輔助變量作為模型的輸入，再通過(guò)mRMR方法[15]排列出與主導(dǎo)變量信息量最大的輔助變量，經(jīng)過(guò)調(diào)試選出最合理的輔助變量。為了防止各輔助變量之間因數(shù)值差異過(guò)大而產(chǎn)生的樣本間的相互影響，對(duì)數(shù)據(jù)進(jìn)行了歸一化和反歸一化。對(duì)輔助變量進(jìn)行皮爾遜相關(guān)性分析，相關(guān)性較大的輔助變量如表1所示?；バ畔⒑Y選出的輔助變量如表2所示。

表1 皮爾遜相關(guān)性分析得到的相關(guān)性較大的輔助變量Tab.1 The auxiliary variables with high correlation were obtained by Pearson correlation analysis

火電廠鍋爐的燃燒工作過(guò)程是一個(gè)非常復(fù)雜的過(guò)程，存在多變量、多耦合以及遲延滯后的問(wèn)題，通過(guò)分析給煤量、總風(fēng)量以及機(jī)組負(fù)荷等變量的變化與入口二氧化硫濃度變化之間的變化趨勢(shì)做了時(shí)間序列的調(diào)整，通過(guò)調(diào)整入口二氧化硫濃度與輔助變量的時(shí)序進(jìn)一步提升了模型的預(yù)測(cè)精度，調(diào)整時(shí)間為600 s。

5.3 預(yù)測(cè)結(jié)果分析

利用上述方法選取1 100組變負(fù)荷時(shí)的數(shù)據(jù)進(jìn)行仿真預(yù)測(cè)實(shí)驗(yàn)，對(duì)800組采樣間隔為1 s的數(shù)據(jù)進(jìn)行訓(xùn)練，采用300組采樣間隔為1 s的具有代表性的升負(fù)荷時(shí)數(shù)據(jù)進(jìn)行預(yù)測(cè)。LSSVM的參數(shù)采用PSO算法進(jìn)行尋優(yōu)，LSSVM采用徑向基函數(shù)RBF作為核函數(shù)；通過(guò)PSO算法得到模型的正規(guī)化參數(shù)δ2和核函數(shù)參數(shù)C，分別為37.140 7和128。利用模型LSSVM與PSO-LSSVM得出的預(yù)測(cè)結(jié)果如圖2和圖3所示。

圖2 Pearson LSSVM模型預(yù)測(cè)結(jié)果Fig.2 Pearson LSSVM model prediction results

圖3 Pearson PSO-LSSVM模型預(yù)測(cè)結(jié)果Fig.3 Pearson PSO-LSSVM model prediction results

利用互信息理論分析得出的輔助變量如表2所示，利用模型LSSVM與PSO-LSSVM得出的預(yù)測(cè)結(jié)果如圖4和圖5所示。

圖4 互信息 LSSVM模型預(yù)測(cè)結(jié)果Fig.4 Mutual Information (MI) LSSVM model prediction results

圖5 互信息 PSO-LSSVM模型預(yù)測(cè)結(jié)果Fig.5 Mutual Information (MI) PSO-LSSVM model prediction results

經(jīng)對(duì)比可以發(fā)現(xiàn)，基于互信息方法篩選出的輔助變量利用PSO算法進(jìn)行參數(shù)尋優(yōu)后的LSSVM模型相對(duì)于李龍等[16]提出的未尋優(yōu)得到的LSSVM模型對(duì)SO2濃度的預(yù)測(cè)值更接近實(shí)際值，且均方根誤差與最大誤差均優(yōu)于普通LSSVM模型。綜合對(duì)比可以得出，基于互信息的PSO-LSSVM模型預(yù)測(cè)結(jié)果最好。各算法的均方根誤差(RMSE)、平均相對(duì)誤差(MRE)和最大誤差(ME)如表3所示。

表3 不同模型誤差對(duì)比Tab.3 Comparison of different model errors

6 結(jié) 論

本文提出了基于互信息的PSO-LSSVM模型預(yù)測(cè)二氧化硫濃度的方法，通過(guò)互信息與mRMR結(jié)合的方法篩選出的18個(gè)與主導(dǎo)變量相關(guān)性較大的輔助變量作為模型的輸入，減少了輔助變量，降低了算法的復(fù)雜性。通過(guò)粒子群(PSO)算法對(duì)LSSVM模型參數(shù)尋優(yōu)提升了模型的預(yù)測(cè)精度。對(duì)比LSSVM與PSO-LSSVM兩種模型的仿真結(jié)果，可以發(fā)現(xiàn)本文所提出的模型預(yù)測(cè)精度更貼近實(shí)際值，誤差小于未經(jīng)尋優(yōu)的LSSVM模型。