楊大州，張尋，姜昊，劉暢，李益國

(1.大唐環(huán)境產(chǎn)業(yè)集團(tuán)股份有限公司特許經(jīng)營分公司，南京 211106；2.東南大學(xué)能源與環(huán)境學(xué)院，南京 210096)

近年來，大氣污染問題日益嚴(yán)峻，電站煤燃燒排放的大量SO2則是大氣污染的主要污染之一[1]。煙氣脫硫技術(shù)是世界上唯一大規(guī)模商業(yè)化應(yīng)用的脫硫方法，其中，石灰石-石膏濕法煙氣脫硫技術(shù)是應(yīng)用最廣泛和最成熟的脫硫工藝技術(shù)。但是，脫硫系統(tǒng)是一個(gè)典型的大慣性和非線性的多變量控制系統(tǒng)，傳統(tǒng)的PID控制難以滿足脫硫系統(tǒng)的控制要求。

為了研究電廠脫硫系統(tǒng)控制策略，建立高精度數(shù)學(xué)模型是首要工作。目前最常見的系統(tǒng)建模方法主要有機(jī)理建模方法與數(shù)據(jù)辨識建模方法。機(jī)理建模需要熟知控制系統(tǒng)的內(nèi)在機(jī)理，且計(jì)算復(fù)雜，而電廠脫硫系統(tǒng)的被控對象往往具有非線性以及大慣性的特點(diǎn)，且脫硫工藝流程繁多，系統(tǒng)特性復(fù)雜，不適合應(yīng)用機(jī)理建模方法；其次，當(dāng)脫硫系統(tǒng)自動控制投入時(shí)獲得的運(yùn)行數(shù)據(jù)，其輸入變量和輸出變量之間存在相關(guān)性，采用最小二乘法等一般的開環(huán)辨識方法難以實(shí)現(xiàn)對模型參數(shù)的無偏估計(jì)；再者，從自適應(yīng)控制的角度看，為了應(yīng)對脫硫系統(tǒng)的慢時(shí)變特性，也需要研究系統(tǒng)的閉環(huán)辨識方法。

子空間辨識方法可以直接獲得被控對象的狀態(tài)空間模型，便于進(jìn)行多變量預(yù)測控制器的設(shè)計(jì)。本文采用兩種多變量閉環(huán)子空間辨識方法，建立脫硫系統(tǒng)的兩輸入兩輸出的狀態(tài)空間模型，并與開環(huán)子空間辨識方法進(jìn)行仿真比較與動態(tài)特性分析。

1 濕法石灰石脫硫系統(tǒng)

濕法石灰石煙氣脫硫技術(shù)是目前世界上最廣泛使用和最成熟的脫硫工藝技術(shù)，采用石灰石作為脫硫吸收劑，除去煙氣中的SO2[2]。該技術(shù)脫硫效率較高，運(yùn)行較為可靠穩(wěn)定，并且能夠適應(yīng)大容量機(jī)組和高濃度SO2煙氣條件。但是，濕法石灰石脫硫系統(tǒng)的控制問題是其發(fā)展的最大阻礙之一。

典型的濕法石灰石煙氣脫硫系統(tǒng)可劃分為煙氣系統(tǒng)、吸收塔系統(tǒng)、石灰石漿液制備系統(tǒng)與石膏脫水系統(tǒng)，其中最核心的子系統(tǒng)是吸收塔系統(tǒng)，吸收塔內(nèi)脫硫過程示意圖如圖1所示。幾乎整個(gè)脫硫反應(yīng)都是在吸收塔內(nèi)完成的。

圖1 吸收塔內(nèi)脫硫過程示意圖

原煙氣降溫后從吸收塔底部進(jìn)入。新鮮石灰石漿液從漿液槽送入，與反應(yīng)后下落至漿液槽的石灰石漿液相互混合?；旌虾蟮臐{液通過循環(huán)泵實(shí)現(xiàn)循環(huán)使用，而經(jīng)過結(jié)晶沉淀后的石膏等副產(chǎn)品被不斷吹掃防止聚集，并被運(yùn)送到儲蓄罐中。吸收區(qū)內(nèi)，循環(huán)漿液從吸收塔頂部由噴嘴噴出，自上而下的漿液與自下上升的煙氣逆向接觸，發(fā)生一系列化學(xué)反應(yīng)，除去煙氣中的SO2。凈化后的煙氣經(jīng)升溫后從煙囪排出。

2 閉環(huán)子空間辨識方法

閉環(huán)子空間辨識方法越來越受到關(guān)注，相較于開環(huán)辨識方法，其用于脫硫系統(tǒng)辨識有以下幾點(diǎn)優(yōu)勢：

1)由于工業(yè)生產(chǎn)過程常常具有很大的干擾，開環(huán)辨識會使得系統(tǒng)在大范圍內(nèi)變得非線性；

2)閉環(huán)辨識有利于控制系統(tǒng)的設(shè)計(jì)，閉環(huán)辨識所得到的模型在控制算法設(shè)計(jì)相關(guān)的頻率范圍內(nèi)是相當(dāng)準(zhǔn)確的，而且對于一些特殊的控制系統(tǒng)設(shè)計(jì)方法，閉環(huán)辨識的方差也要小于開環(huán)辨識；

3)閉環(huán)子空間辨識大多數(shù)都有在未來輸入的補(bǔ)空間進(jìn)行正交投影計(jì)算的步驟，從而可將噪聲項(xiàng)消除[3]。

目前國內(nèi)外閉環(huán)子空間辨識方法應(yīng)用較多的方法主要有基于正交投影的閉環(huán)子空間辨識方法(Closed-loop Subspace Orthogonal Projection Identification Method，CSOPIM)[4]與基于新息估計(jì)的閉環(huán)子空間辨識方法(Parallel parsimonious SIM with Estimation，PARSIME)[5]。

2.1 子空間辨識的基礎(chǔ)知識

考慮如下的離散線性時(shí)不變狀態(tài)空間模型：

(1)

式(1)中，xk∈n為系統(tǒng)的狀態(tài)變量；uk∈nu為輸入觀測向量；yk∈ny為輸出觀測向量；wk∈n與vk∈ny分別為系統(tǒng)的測量噪聲與過程噪聲。

假設(shè)系統(tǒng)是能觀測的，可設(shè)計(jì)如下Kalman濾波器：

(2)

式(2)中，K為卡爾曼濾波增益。

(3)

式(3)中，新息ek為零均值白噪聲。

定義輸入數(shù)據(jù)Hankel矩陣Up與Uf：

定義狀態(tài)序列Xi=[xixi+1…xi+j-2xi+j-1]∈n×j，則過去狀態(tài)序列為Xp=X0=[x0x1…xj-2xj-1]，未來狀態(tài)序列為Xf=Xi=[xixi+1…xi+j-2xi+j-1]。

首先考慮以下狀態(tài)預(yù)測方程：

x1=Ax0+Bu0+Ke0

x2=Ax1+Bu1+Ke1

(4)

=A(Ax0+Bu0+Ke0)+Bu1+Ke1

類推可得到i時(shí)刻系統(tǒng)狀態(tài)變量與i，i+1，…，i+j-1時(shí)刻的狀態(tài)預(yù)測，并代入新息形式的狀態(tài)空間模型式(3)，可得到系統(tǒng)i時(shí)刻與0，1，…，i-1時(shí)刻的輸出。

Γi=[CCA…CAi-1]T

由此得到構(gòu)造的Hankel矩陣Yp，Yf以及Ep，Ef，下標(biāo)p和f分別表示過去和將來相對時(shí)間的概念：

Yp=ΓiXp+HiUp+GiEp

(5)

(6)

Yf=ΓiXf+HiUf+GiEf

(7)

2.2 基于正交投影的閉環(huán)子空間辨識方法

CSOPIM辨識方法是基于正交投影的辨識方法，從控制對象的輸入輸出著手辨識對象模型，通過正交投影將輸入輸出數(shù)據(jù)轉(zhuǎn)化為卡爾曼狀態(tài)序列，再求出系統(tǒng)矩陣[6]。關(guān)鍵是引入了包含設(shè)定值Hankel矩陣的輔助變量。

2.2.1 基于正交投影的子空間辨識方法

(8)

對Z進(jìn)行SVD分解如下：

(9)

理論上Z的秩應(yīng)該為nui+n，所以根據(jù)上述SVD分解結(jié)果進(jìn)行階次判定。

根據(jù)式(9)可知Z的正交補(bǔ)空間為U2。取M為單位矩陣，矩陣U2M可分解為U2M=(P1P2)T，簡單變形得到：

(10)

但是上述方法用于系統(tǒng)閉環(huán)辨識效果并不理想，一種可行的解決方法是引入包含設(shè)定值Hankel矩陣的輔助變量[8]：

(11)

需要指出的是，輔助變量的構(gòu)造方法并不是唯一的。

2.2.2 系統(tǒng)矩陣的計(jì)算方法

卡爾曼濾波狀態(tài)可根據(jù)式(12)計(jì)算得到：

(12)

(13)

根據(jù)式(13)可求出Hi1，并構(gòu)造出Hi。

定義：Rf=Ri+1|2i-1Up=U0|iYp=Y0|i

Uf=Ui+1|2i-1Yf=Yi+1|2i-1，從而得到：

(14)

式(14)中，Γi為Γi去掉最后ny行得到，Hi由Hi去掉最后ny行與nu列得到。

系統(tǒng)的參數(shù)矩陣A，B，C，K可根據(jù)以下公式求解得到：

(15)

但是該方法在閉環(huán)系統(tǒng)中應(yīng)用效果可能不理想，由于反饋控制律的存在，Ui|i可能與Xi相關(guān)，從而導(dǎo)致有偏估計(jì)問題。為了解決這一問題，可以通過下述步驟求解：

2)假定系統(tǒng)的參數(shù)矩陣中D=0；

3)參數(shù)矩陣B與K可以通過下式求解：

(16)

2.3 基于新息估計(jì)的閉環(huán)子空間辨識方法

在閉環(huán)辨識環(huán)境下，未來輸入與過去噪聲是相關(guān)的，使得正交投影不能消除噪聲項(xiàng)，所以基于正交投影的子空間辨識在閉環(huán)下會存在有偏估計(jì)。而基于新息估計(jì)的閉環(huán)子空間辨識方法(PARSIME)的基本思想為由上一行的估計(jì)新息序列結(jié)果帶到下一行中，下一行估計(jì)得到的值繼續(xù)應(yīng)用到后面的一行，即PARSIME算法從第一行塊中估計(jì)新息序列，然后認(rèn)為是已知新息去估計(jì)其它模型參數(shù)。

2.3.1 基于新息估計(jì)的子空間辨識方法

閉環(huán)子空間辨識與開環(huán)子空間辨識最大的不同在于未來輸入數(shù)據(jù)與過去噪聲數(shù)據(jù)存在相關(guān)性，故需要先將相關(guān)性消除，然后再進(jìn)行無偏差辨識。在PARSIME方法中考慮將Hankel矩陣Yf按行分塊：

Yf=[Yf1Yf2…Yfi…Yff]T，i=1，…，f

(17)

采用2.1節(jié)中類似的方法，得到如下的矩陣等式：

(18)

式(18)中，

Γfi=CAi-1

Hfi=[CAi-2B…CBD]?[Hi-1…H1H0]

Gfi=[CAi-2K…CKI]?[Gi-1…G1G0]

(19)

每一行需要計(jì)算Efi的估計(jì)值，計(jì)算如下：

(20)

當(dāng)i=1時(shí)有：

(21)

(22)

因?yàn)樾孪⑿蛄幸呀?jīng)得到估計(jì)值替代原始數(shù)據(jù)，故未來輸入數(shù)據(jù)不再與新息序列相關(guān)，直接使用最小二乘求解。

借鑒其它子空間辨識方法，可確定權(quán)矩陣如下[9]：

(23)

(24)

2.3.2 系統(tǒng)矩陣的計(jì)算

關(guān)于系統(tǒng)矩陣的計(jì)算，根據(jù)上述計(jì)算結(jié)果重構(gòu)系統(tǒng)狀態(tài)，再進(jìn)行系統(tǒng)矩陣的計(jì)算[10]。其具體計(jì)算步驟如下：

4)計(jì)算系統(tǒng)參數(shù)矩陣A，B與K，用MATLAB矩陣形式表示如下：

(25)

PARSIME方法將整體辨識模型分解為若干子模型，通過對子模型進(jìn)行多步計(jì)算，避免了回歸框架中對冗余非因果參數(shù)的重復(fù)估計(jì)，同時(shí)對新息模型計(jì)算得到新息估計(jì)值來代替系統(tǒng)噪聲估計(jì)值，消除噪聲的影響，解決了基于正交投影的閉環(huán)子空間辨識方法的有偏估計(jì)問題。

3 脫硫系統(tǒng)的閉環(huán)子空間辨識

3.1 閉環(huán)數(shù)據(jù)生成

針對某脫硫系統(tǒng)的傳遞函數(shù)模型[11]式(26)，采用基于擴(kuò)增狀態(tài)空間模型的預(yù)測控制器對該模型進(jìn)行控制，生成閉環(huán)數(shù)據(jù)的預(yù)測控制器參數(shù)見表1。

(26)

式(26)中，pH表示漿液的pH值；SO2表示出口的SO2濃度，10-6；Q表示石灰石供給泵轉(zhuǎn)速，r/min；L表示漿液循環(huán)泵電機(jī)頻率，Hz。

以脫硫系統(tǒng)某運(yùn)行平衡點(diǎn)為基準(zhǔn)，令pH和SO2設(shè)定值按方波形式變化，同時(shí)在輸入和輸出變量上疊加白噪聲信號，以模擬測量噪聲，閉環(huán)數(shù)據(jù)生成相關(guān)參數(shù)見表2。仿真時(shí)間為15 000 s，最終生成的閉環(huán)輸入和輸出數(shù)據(jù)如圖2與圖3所示。

表1 生成閉環(huán)數(shù)據(jù)的預(yù)測控制器參數(shù)

圖2 生成的閉環(huán)輸入數(shù)據(jù)

圖3 生成的閉環(huán)輸出數(shù)據(jù)

3.2 閉環(huán)子空間辨識

3.2.1 辨識結(jié)果

分別采用CSOPIM與PARSIME辨識方法，計(jì)算得到的狀態(tài)空間模型的系統(tǒng)矩陣如下：

CSOPIM方法的辨識結(jié)果：

PARSIME方法的辨識結(jié)果：

3.2.2 數(shù)據(jù)擬合效果對比分析

將兩種閉環(huán)子空間辨識方法CSOPIM和PARSIME以及開環(huán)子空間辨識方法N4SID，與原始閉環(huán)數(shù)據(jù)的擬合效果進(jìn)行對比，基于閉環(huán)數(shù)據(jù)的三種辨識方法辨識效果比較如圖4所示。將每種辨識方法得到的兩個(gè)輸出量分別與原始閉環(huán)數(shù)據(jù)作差，得到輸出量的誤差，再計(jì)算輸出量的絕對誤差平均值與誤差標(biāo)準(zhǔn)差，基于閉環(huán)數(shù)據(jù)的三種辨識方法的辨識效果比較指標(biāo)見表3。由圖表可知，閉環(huán)子空間辨識方法的效果總體比開環(huán)子空間辨識方法好，PARSIME方法的辨識效果優(yōu)于CSOPIM方法。

表3 基于閉環(huán)數(shù)據(jù)的三種辨識方法的辨識效果比較指標(biāo)

從理論角度來說，N4SID方法適用于開環(huán)數(shù)據(jù)的辨識，且易受干擾影響，因此N4SID方法用于閉環(huán)辨識效果不好。CSOPIM方法雖然屬于閉環(huán)辨識方法，但是由于正交投影無法消除噪聲項(xiàng)的影響，存在有偏估計(jì)的問題，而PARSIME方法通過將未知新息序列視為已知參數(shù)，從而能夠得到系統(tǒng)參數(shù)的一致估計(jì)值，通過估計(jì)值來代替系統(tǒng)噪聲估計(jì)值，以消除噪聲的影響。因此PARSIME方法的辨識效果會更好。

圖4 基于閉環(huán)數(shù)據(jù)的三種辨識方法辨識效果比較

3.2.3 階躍響應(yīng)對比分析

以表2中平衡點(diǎn)為基準(zhǔn)，對石灰石供給泵轉(zhuǎn)速做+1 r/min的階躍，保持漿液循環(huán)泵頻率不變，石灰石供給泵階躍響應(yīng)曲線比較圖如圖5所示；另外對漿液循環(huán)泵頻率做+1Hz的階躍，保持石灰石供給泵轉(zhuǎn)速不變，漿液循環(huán)泵階躍響應(yīng)曲線比較如圖6所示。

綜合分析圖4、圖5和圖6可知，不管是從數(shù)據(jù)擬合效果，還是從階躍響應(yīng)曲線，對于閉環(huán)數(shù)據(jù)而言，閉環(huán)子空間辨識方法的效果都比開環(huán)子空間辨識方法好，同時(shí)PARSIME方法的辨識效果優(yōu)于CSOPIM方法。

圖5 石灰石供給泵階躍響應(yīng)曲線比較

圖6 漿液循環(huán)泵階躍響應(yīng)曲線比較

3.3 模型動態(tài)特性分析

分析圖5可知，當(dāng)石灰石供漿量增大時(shí)，由于新增的石灰石的溶解，漿液pH開始增大，同時(shí)吸收塔出口SO2濃度開始降低。當(dāng)石灰石溶解速率與SO2吸收速率相等時(shí)，系統(tǒng)達(dá)到新的穩(wěn)定狀態(tài)。兩者的變化方向相反，都屬于大慣性和有自平衡能力的被控對象，其上升時(shí)間達(dá)3 000 s左右。

分析圖6可知，當(dāng)漿液循環(huán)流量發(fā)生階躍變化時(shí)，由于液氣比增加，煙氣中更多的SO2被捕獲，因此漿液pH和出口SO2濃度同時(shí)減小。pH的減小由離子液體相的緩沖能力阻止，并且隨著漿液中石灰石的逐漸溶解pH降低幅度慢慢減小。因?yàn)镾O2的吸收速率大于石灰石溶解速率，所以離子液體相的緩沖能力耗盡，pH繼續(xù)減小。SO2吸收速率隨之降低，石灰石溶解速率提高。當(dāng)兩者速率相等時(shí)，系統(tǒng)達(dá)到新的平衡。兩者都屬于大慣性和有自平衡能力的被控對象，但漿液循環(huán)流量變化對pH的影響明顯快于對SO2濃度的影響。

通過分析，對濕法脫硫系統(tǒng)的主要動態(tài)特性總結(jié)如下：

1)由于脫硫系統(tǒng)的石灰石溶解速率低，出口SO2濃度和漿液pH的動態(tài)特性較慢。因此，需要在控制器設(shè)計(jì)中考慮輸出預(yù)測和動態(tài)超調(diào)等措施，以改善系統(tǒng)的調(diào)節(jié)品質(zhì)。

2)漿液循環(huán)流量變化對出口SO2濃度的影響在初期比較明顯，但是在低Ca/S(pH)比運(yùn)行時(shí)，漿液循環(huán)流量變化對出口SO2濃度幾乎沒有影響。這就表明，脫硫系統(tǒng)SO2吸收能力主要取決于外部石灰石的供給。

4 結(jié)語

本文采用基于正交分解和基于新息估計(jì)兩種典型的閉環(huán)子空間辨識方法，建立了脫硫系統(tǒng)的狀態(tài)空間模型，并對其有效性進(jìn)行了對比分析。通過將兩種閉環(huán)子空間辨識方法的辨識效果與開環(huán)子空間辨識方法N4SID進(jìn)行對比，驗(yàn)證了閉環(huán)子空間辨識方法的有效性。同時(shí)通過理論分析與辨識效果比較，發(fā)現(xiàn)以PARSIME方法建立的脫硫系統(tǒng)的動態(tài)數(shù)學(xué)模型精度優(yōu)于CSOPIM方法。此外，本文對濕法脫硫系統(tǒng)主要動態(tài)特性做出總結(jié)，為脫硫系統(tǒng)先進(jìn)控制策略的研究奠定基礎(chǔ)。