遺傳編程在低濃打漿過程建模中的應(yīng)用

張政江

（浙江大學(xué)，浙江杭州，310027）

打漿過程是制漿造紙生產(chǎn)過程中極為重要的一環(huán)。造紙所用的漿料需要通過打漿的機(jī)械作用，使?jié){料獲得一些特定的性質(zhì)，以滿足生產(chǎn)紙張的質(zhì)量要求。由于未經(jīng)打漿的原漿多種多樣，同時(shí)打漿過程中還存在許多擾動(dòng)，這些因素都會(huì)影響成漿質(zhì)量的穩(wěn)定性，因此有必要對(duì)打漿過程進(jìn)行自動(dòng)控制［1］。

目前國外廣泛用于打漿過程控制的機(jī)理假說是比邊緣負(fù)荷理論（Specific Edge Load Theory，簡稱SEL理論），這一假說是針對(duì)盤磨機(jī)打漿提出的。根據(jù)這一理論［2］，通過控制比能量、比負(fù)荷就可以得到質(zhì)量穩(wěn)定的成漿，因此為了實(shí)現(xiàn)打漿過程的控制就需要建立比能量、比負(fù)荷的動(dòng)態(tài)模型。但比能量、比負(fù)荷與相關(guān)變量之間的定量關(guān)系還有待于進(jìn)一步研究。

20世紀(jì)90年代發(fā)展起來的遺傳編程（Genetic Programming，簡寫為GP）［3-4］是一種自動(dòng)編程技術(shù)。遺傳編程是一種與遺傳算法（Genetic Algorithm）很相似的進(jìn)化算法，也是一種基于生物自然選擇與遺傳機(jī)理的隨機(jī)搜索算法，僅僅告訴計(jì)算機(jī)要做什么，而沒有明確告訴計(jì)算機(jī)如何逐步地去完成所要求的任務(wù)。自Koza提出遺傳編程以來，遺傳編程已經(jīng)在許多領(lǐng)域得到應(yīng)用，如用于一類化工過程靜態(tài)輸入輸出模型的建立［5］、動(dòng)態(tài)模型的建立［6］等。其結(jié)果表明，遺傳編程具有自動(dòng)選擇輸入變量、自動(dòng)確定模型結(jié)構(gòu)的特性，是一種新的建模方法。但遺傳編程存在進(jìn)化速度慢，搜索效率低，進(jìn)化到最優(yōu)解需要時(shí)間長等缺陷。

本研究包括兩部分內(nèi)容：①在分析了遺傳編程進(jìn)化速度慢的原因后，采用新的適應(yīng)值調(diào)整方法，由此改善了進(jìn)化速度；②進(jìn)行基于遺傳編程的低濃打漿過程動(dòng)態(tài)建模工作，利用遺傳編程具有自動(dòng)選擇輸入變量、自動(dòng)確定模型結(jié)構(gòu)的特性，自動(dòng)確定低濃打漿過程比能量、比負(fù)荷與各相關(guān)變量之間的定量關(guān)系。

1 遺傳編程算法及改進(jìn)

遺傳編程是對(duì)個(gè)體的集合——群體進(jìn)行遺傳操作。個(gè)體是一棵有根、有結(jié)點(diǎn)、帶標(biāo)記、有序分支的樹，即用樹表示一個(gè)問題的可能解或一個(gè)表達(dá)式。遺傳編程的搜索空間是問題所有可能的函數(shù)和端點(diǎn)組合的集合。函數(shù)集可以由數(shù)學(xué)函數(shù)、條件算子或任意定義的特殊函數(shù)等組成。端點(diǎn)集可以由常數(shù)、傳感器的輸出、系統(tǒng)的狀態(tài)變量、系統(tǒng)的控制輸入等組成。

遺傳編程首先從函數(shù)集和端點(diǎn)集中隨機(jī)選擇元素組成初始群體，然后確定每個(gè)個(gè)體的適應(yīng)值。適應(yīng)值是反應(yīng)個(gè)體與問題實(shí)際值之間的接近程度，通常是基于問題實(shí)際值與個(gè)體返回值之間的誤差計(jì)算的。確定了初始群體后，利用復(fù)制、交叉和變異操作從當(dāng)前群體中產(chǎn)生新的群體。復(fù)制操作首先按基于適應(yīng)值的某種選擇方法從當(dāng)前群體中選擇1個(gè)個(gè)體，然后把它復(fù)制到下一代群體中。再按交叉率從當(dāng)前群體中隨機(jī)選擇2個(gè)個(gè)體，通過交叉操作得到2個(gè)新的個(gè)體。交叉操作生成的子代個(gè)體不同于父代個(gè)體，但它們完全是根據(jù)父代個(gè)體的基因生成的。如果父代個(gè)體是問題的一個(gè)可能解，那么父代個(gè)體就會(huì)有與問題解有關(guān)的基因，由這些基因相結(jié)合而生成的子代個(gè)體就有可能更逼近于問題的最優(yōu)解。如此循環(huán)，直至滿足準(zhǔn)則或遺傳代數(shù)已達(dá)最大數(shù)，把當(dāng)前代中最好的個(gè)體作為遺傳編程設(shè)計(jì)的最好結(jié)果。算法詳細(xì)介紹見文獻(xiàn)［5-8］。

遺傳編程存在遺傳進(jìn)化速度慢，搜索效率低，進(jìn)化到最優(yōu)解需要時(shí)間長的缺陷。傳統(tǒng)的遺傳編程是通過復(fù)制、交叉等操作生成新的下一代。由復(fù)制操作選擇適應(yīng)值好的個(gè)體復(fù)制到交叉池中，再按交叉率隨即選擇個(gè)體進(jìn)行交叉，產(chǎn)生新的一代，循環(huán)執(zhí)行這個(gè)過程，直至得到最優(yōu)解。從這個(gè)過程可以看出：①遺傳編程只是將適應(yīng)值高的個(gè)體復(fù)制，也就是只著眼于當(dāng)前代中好的個(gè)體，并沒有考慮個(gè)體的遺傳進(jìn)化速度；②遺傳編程的目標(biāo)是平均適應(yīng)值的提高，隨著平均適應(yīng)值的提高，最好的個(gè)體也就越接近問題的解。不難想象，增加對(duì)適應(yīng)值較好的、進(jìn)化速度快的個(gè)體的復(fù)制數(shù)，會(huì)有助于改善遺傳進(jìn)化的收斂速度。文獻(xiàn)［9］根據(jù)群體的可進(jìn)化性（Evolvability）概念對(duì)交叉或變異操作后產(chǎn)生的群體中適應(yīng)值好于群體平均適應(yīng)值的個(gè)體進(jìn)行適應(yīng)值激勵(lì)，以改善進(jìn)化速度，提高遺傳算法的計(jì)算效率。而在遺傳編程中，由于交叉后的個(gè)體在計(jì)算過程動(dòng)態(tài)響應(yīng)時(shí)有可能產(chǎn)生溢出，為了維持群體的規(guī)模，在進(jìn)化過程中需不斷用新產(chǎn)生的個(gè)體代替溢出的個(gè)體，而新個(gè)體的適應(yīng)值往往是比較差的，這樣會(huì)使交叉后群體平均適應(yīng)值降低。如果仍采用文獻(xiàn)［9］的對(duì)好于交叉后適應(yīng)度平均值的個(gè)體進(jìn)行適應(yīng)值激勵(lì)，就有可能是對(duì)大多數(shù)個(gè)體進(jìn)行適應(yīng)值激勵(lì)，失去了進(jìn)行適應(yīng)值激勵(lì)的意義。因此，本研究采取如下方法：對(duì)好于復(fù)制后平均適應(yīng)值的個(gè)體，根據(jù)其交叉前后適應(yīng)值的變化量對(duì)適應(yīng)值進(jìn)行適當(dāng)調(diào)整的方法。在交叉產(chǎn)生的群體中，考慮好于復(fù)制后群體平均適應(yīng)值的個(gè)體O。假如個(gè)體O是由父代個(gè)體p1、p2交叉產(chǎn)生的，則個(gè)體O的調(diào)整后適應(yīng)值計(jì)算為：

其中，f′（·）為調(diào)整后的適應(yīng)值，f（·）為原適應(yīng)值；如果個(gè)體O是從復(fù)制后群體中直接復(fù)制的，則f′（O）=f（O）。對(duì)于低于復(fù)制后群體平均適應(yīng)值的個(gè)體O，f′（O）=f（O）。

2 基于遺傳編程的低濃打漿過程動(dòng)態(tài)建模

本研究建模所用的低濃打漿工藝流程如圖1所示，為國內(nèi)典型的φ450mm雙盤盤磨機(jī)打漿流程。流量為F0的清水和濃度為C1、流量為F1的原漿經(jīng)漿泵混合后成為流量為F2、濃度為C2的漿進(jìn)入盤磨機(jī)，在盤磨機(jī)中經(jīng)打漿處理后成為成漿送入成漿池。文獻(xiàn)［10］根據(jù)SEL理論建立了該流程原漿流量、清水流量、盤磨機(jī)功率與比能量、比負(fù)荷之間的機(jī)理動(dòng)態(tài)模型。

圖1 低濃打漿工藝流程圖

采用遺傳編程計(jì)算方法進(jìn)行建模工作。根據(jù)打漿過程機(jī)理模型的響應(yīng)，將所建模型的時(shí)域長度取為14。端點(diǎn)集選擇為4個(gè)輸入變量、1個(gè)輸出變量和一些任意取定的常數(shù)，即｛u1（k），u1（k-1），…，u1（k-14），u2（k），u2（k-1），…，u2（k-14），u3（k），u3（k-1），…u3（k-14），u4（k），u4（k-1），…，u4（k-14），yi（k-1），…，yi（k-14），4.1，4.5，5.6，7.8｝；其中，u1表示清水流量，u2表示原漿流量，u3表示盤磨機(jī)功率，y1為比能量、y2為比負(fù)荷。為了檢驗(yàn)遺傳編程是否具有選擇輸入變量的能力，加入了任意的變量u4。函數(shù)集取為｛+，-，×，pow2，sqrt｝。群體規(guī)模設(shè)為50個(gè)個(gè)體，初始個(gè)體的結(jié)點(diǎn)數(shù)為33，進(jìn)化過程中個(gè)體的最大結(jié)點(diǎn)數(shù)為35。適應(yīng)值函數(shù)取為：

式中，RMS（O）為某個(gè)體在所有試驗(yàn)點(diǎn)上返回值與模型參考值之間的均方根誤差。適應(yīng)值越大，則均方根誤差越小，即個(gè)體越接近問題的最優(yōu)解。復(fù)制操作時(shí)采取競(jìng)賽規(guī)模為2的聯(lián)賽選擇方法，將2個(gè)個(gè)體中適應(yīng)值好的個(gè)體復(fù)制到下一代群體。交叉率為0.8，在進(jìn)行交叉操作時(shí)除保留最佳個(gè)體，即最佳個(gè)體直接復(fù)制到下一代，最佳個(gè)體也參與交叉。遺傳計(jì)算停止條件為均方根誤差小于0.06，或進(jìn)化代數(shù)大于600。

個(gè)體中含有一些常數(shù)結(jié)點(diǎn)，如果對(duì)這些常數(shù)進(jìn)行優(yōu)化辨識(shí)，會(huì)有助于改善個(gè)體適應(yīng)值。遺傳編程雖然也可以對(duì)這些常數(shù)進(jìn)行辨識(shí)［3］，但遺傳編程更善于結(jié)構(gòu)辨識(shí)，因此有必要采用其他的系數(shù)辨識(shí)方法［7］。本研究采用單純形法對(duì)常數(shù)結(jié)點(diǎn)系數(shù)進(jìn)行優(yōu)化。將某個(gè)個(gè)體的系數(shù)（設(shè)總數(shù)為M）作為一個(gè)頂點(diǎn)（M維），再隨機(jī)產(chǎn)生M個(gè)頂點(diǎn)，利用單純形法進(jìn)行系數(shù)優(yōu)化計(jì)算，優(yōu)化目標(biāo)是頂點(diǎn)在各個(gè)試驗(yàn)點(diǎn)上的個(gè)體返回值與模型參考值之間的均方根誤差為最小。當(dāng)各個(gè)頂點(diǎn)的均方根誤差與均方根誤差平均值的差值小于設(shè)定值，或優(yōu)化迭代次數(shù)已達(dá)最大值，則結(jié)束優(yōu)化迭代，將均方根誤差最小的頂點(diǎn)作為該個(gè)體的優(yōu)化系數(shù)。

考慮了適應(yīng)值調(diào)整和系數(shù)優(yōu)化后，動(dòng)態(tài)建模的遺傳編程計(jì)算步驟如下：

（1）從函數(shù)集、端點(diǎn)集中隨機(jī)選擇元素生成初始個(gè)體，遍歷個(gè)體計(jì)算出個(gè)體的適應(yīng)值；

（2）依照聯(lián)賽選擇方法進(jìn)行復(fù)制操作，并計(jì)算復(fù)制后群體的平均適應(yīng)度；

（3）滿足遺傳停止條件，則轉(zhuǎn)至（8）；

（4）如果最好的適應(yīng)值連續(xù)6次不變，則加入隨機(jī)產(chǎn)生的新個(gè)體；

（5）根據(jù)交叉率進(jìn)行交叉操作；

（6）個(gè)體常數(shù)結(jié)點(diǎn)優(yōu)化計(jì)算，個(gè)體適應(yīng)值計(jì)算，如有個(gè)體產(chǎn)生溢出現(xiàn)象，則以隨機(jī)產(chǎn)生的新樹代替；

（7）根據(jù)第2部分的個(gè)體適應(yīng)值調(diào)整方法計(jì)算個(gè)體的適應(yīng)值，返回（2）；

（8）將當(dāng)前代中適應(yīng)值最好的個(gè)體作為問題的解，并結(jié)束計(jì)算。

將隨機(jī)產(chǎn)生的100點(diǎn)u1、u2、u3和u4的動(dòng)態(tài)序列組作為文獻(xiàn)［10］機(jī)理模型的輸入，將所得到的動(dòng)態(tài)響應(yīng)y1、y2和u1、u2、u3、u4序列作為本研究遺傳編程的訓(xùn)練數(shù)據(jù)；并分別產(chǎn)生y1、y2對(duì)于單位階躍輸入的響應(yīng)，作為遺傳編程所得模型的校正數(shù)據(jù)，以檢驗(yàn)?zāi)Ｐ偷木取?/p>

經(jīng)遺傳計(jì)算所得的低濃打漿過程動(dòng)態(tài)模型為：

式中，y1的RMS誤差為0.0615，y2的RMS誤差為0.0631。

原漿流量階躍作用下的新建模型比能量、比負(fù)荷響應(yīng)曲線如圖2、圖3所示；盤磨機(jī)功率階躍作用下的新建比能量、比負(fù)荷響應(yīng)曲線如圖4、圖5所示。結(jié)果說明，基于遺傳編程所建立的模型較好地反映了打漿過程中比能量、比負(fù)荷的變化。

在上述模型的輸入項(xiàng)中沒有出現(xiàn)清水流量u1，這跟工藝過程分析的情況不同。從過程的動(dòng)態(tài)數(shù)據(jù)中發(fā)現(xiàn)，當(dāng)清水流量發(fā)生變化時(shí)，對(duì)比能量、比負(fù)荷的影響比較小。而遺傳編程在搜索解空間時(shí)是依據(jù)適應(yīng)值函數(shù)進(jìn)行的，當(dāng)某個(gè)輸入變量對(duì)輸出變量的作用比較小時(shí)，在個(gè)體中有這個(gè)輸入變量與沒有這個(gè)輸入變量的適應(yīng)值相差很小，則在進(jìn)化過程中這個(gè)輸入變量就有可能被當(dāng)作無關(guān)的變量而舍去，由此使上述模型中不含清水流量。另外，上述模型中也沒有出現(xiàn)u4，因?yàn)閡4是任意加入的變量，與y1、y2并沒有關(guān)聯(lián)，故遺傳編程也沒有選擇u4。

圖6為采用文獻(xiàn)［9］的激勵(lì)方法和本文第2部分的計(jì)算方法對(duì)上述動(dòng)態(tài)過程中y1建模的遺傳編程計(jì)算結(jié)果對(duì)比。圖中縱座標(biāo)為個(gè)體的RMS（O），曲線所示為一代中最小的RMS（O）。因此，該均方根誤差越小，所得的個(gè)體越接近問題的最優(yōu)解。結(jié)果說明，本研究的計(jì)算方法有助于改善群體的收斂速度。

3 結(jié)語

本研究的結(jié)果表明，所采用的改進(jìn)算法有助于改善遺傳編程的收斂速度，所建立的低濃打漿過程動(dòng)態(tài)模型較好地反映了比能量、比負(fù)荷的變化。因此，遺傳編程適合于建立變量相互關(guān)系不明的打漿過程的模型。同時(shí)也要注意到，如果某一個(gè)輸入變量對(duì)輸出變量的作用不大，則遺傳編程有可能無法對(duì)該變量做出正確的選擇。

圖6 遺傳編程計(jì)算結(jié)果比較

［1］陳蘇平，皮道映，孫優(yōu)賢.打漿過程控制的現(xiàn)狀和趨勢(shì)［J］.中國造紙，1994，13（3）：51.

［2］皮道映，孫優(yōu)賢.打漿過程機(jī)理模型的研究現(xiàn)狀與展望［J］.中國造紙學(xué)報(bào)，1995，10：86.

［3］Koza J R.Genetic programming：on the programming of computers by means of natural selection［C］.USA：TheM IT Press，1992.

［4］Koza J R.Genetic programmingⅡ：automatic discovery of reusable programs［C］.USA：TheM IT Press，1994.

［5］McKay B，Willis M，Barton G.Steady-state modeling of chemical process systems using genetic programming［J］.Computers chem.Engng.，1997，21（9）：981.

［6］W illisM，Hiden H，HinchliffeM.et al.Systemsmodeling using genetic programming［J］.Computers chem.Engng.，1997，21（Suppl）：1161.

［7］Gray G J，Murray-Smith D J，Li Y，et al.Nonlinearmodel structural identification using genetic programming［J］.Control Engineering Practice，1998，6：1341.

［8］劉 勇，康立山，陳毓屏.非數(shù)值并行算法——遺傳算法［M］.北京：科學(xué)出版社，1998.

［9］林 丹，李敏強(qiáng)，寇紀(jì)淞.用適應(yīng)值激勵(lì)機(jī)制提高遺傳算法的效率［J］.控制與決策，2000，15（6）：759.

［10］皮道映.多模型控制及應(yīng)用研究［D］.杭州：浙江大學(xué)，1995.