基于數據驅動的一類復雜系統動態特性分析

郭玲利,王目樹,王 燁

(1.北京科技大學 自動化學院,北京 100083; 2.長治學院 數學系, 山西 長治 046000)

基于數據驅動的一類復雜系統動態特性分析

郭玲利1,2,王目樹1,王 燁1

(1.北京科技大學 自動化學院,北京 100083; 2.長治學院 數學系, 山西 長治 046000)

針對一類復雜的無法對其機理建模的離散時間系統，根據采集的兩年工藝參數數據，結合復雜工藝特點，提出了基于數據驅動的系統動態特性建模方法，構建了時間序列受控回歸滑動平均(CARMA)胞映射模型;在模型結構確定的基礎上，采用改進的量子行為粒子群優化(IQPSO)算法對系統參數進行辨識;算法通過設計新的粒子更新式增加了粒子的多樣性，避免了算法的早熟收斂；算法通過在后期將搜索到的最優值傳遞給神經元作為初始權值，利用神經元增強算法的局部搜索能力，實現了算法探索與開發的平衡，達到對模型參數進行快速精確辨識的目的;在轉化為狀態空間模型基礎上，根據胞映射理論對系統進行了穩定性分析，通過對胞映射作圖快速獲得平衡胞，利用動態優化原理，找到所有的周期胞和吸引域，達到對系統穩定性分析的目的;利用現場工藝數據進行仿真，結果證明了所提方法的有效性。

數據驅動；建模；粒子群算法；胞映射；穩定性

0 引言

隨著科學技術、特別是信息科學技術的快速發展, 冶金、化工、機械、電子、電力、交通運輸和物流等企業發生了重大變化。企業的規模越來越大, 生產設備、生產工藝和生產過程越來越復雜[1-3]。 傳統方法, 即依據物理化學機理建立精確數學模型, 已變得越來越困難[4-9]。同時相當多的企業每天都在產生大量的生產和過程數據, 這些數據隱含著某些工藝變動或設備運行等信息。理論和工程界迫切需要解決的問題是，如何有效利用這些大量的數據和知識, 在難以建立受控系統較準確機理模型的條件下, 還可實現對生產過程的預報、評價和優化控制。因此, 發展基于數據驅動的控制理論與方法是新時期理論發展與應用的必然要求, 具有極為重要的理論與現實意義[10-12]。

本文針對一類復雜系統進行基于數據驅動方法的建模和控制正是在這一思想理念下產生的。由于這類復雜系統如燒結過程中發生的一系列的物理化學變化,使得其動態特性難以用精確的機理模型進行描述。但通過考察工藝特點，受工藝現場看火師傅根據火焰類型進行溫度調節從而使系統工況維持正常運行，進而保證良好的產品質量的思想啟發，我們嘗試通過已采集的工藝過程參數數據來構建這樣一個可以代替人工判斷實現復雜系統的自動控制模型。在模型結構確定的基礎上，采用改進的量子行為粒子群算法(IQPSO)進行參數辨識。在模型基礎上，我們針對相應的自治系統，利用胞映射理論對系統穩定性加以分析。最后，利用采集到的2年的現場實際數據，對改進的算法和胞映射理論加以驗證和應用，從而對這一類復雜系統有較為深入了解。

1 復雜系統建模

針對類似燒結機煅燒過程，各種燒結原料按一定比例混合布在臺車上，臺車一邊接受布料、同時向前行進，運行至點火器下，點火，抽風，空氣從料層表面吸入，遇到混合料中的燃料發生燃燒反應，在燃料燃燒產生高溫的同時發生一系列的物理化學反應，部分混合料顆粒表面發生軟化，熔化，產生一定量的液相，并潤濕其他未熔化的顆粒，當冷卻后，液相將礦粉顆粒粘結成塊，這個過程稱為燒結。

針對這一類機理難建模系統的復雜特性，通過分析工藝特點，得知燒結質量的好壞取決于燒結終點，因此我們利用歷時2年從工業現場采集到靠近燒結終點的燒結臺車兩邊的共8269組6維風箱溫度數據，而為了消除采樣過程中的干擾和測量誤差，因此對數據進行濾波，濾波前后的風箱溫度如圖2所示。

圖2 濾波前后時間序列

進一步我們對濾波后的工藝參數數據利用數據驅動和胞映射的建模方法，并進行自相關函數和偏自相關函數的分析，同時通過F檢驗法對模型階數和延時進行顯著性檢驗(表1)，構建了如下時間序列受控回歸滑動平均(CARMA)胞映射模型結構。

表1 殘差平方和檢驗

y(k+3)=-a2yc(k+2)-a1yc(k+1)-

a0yc(k)+b2u(k+2)+b1u(k+1)

(1)

在模型確定的基礎上，我們采用改進的量子行為粒子群優化(IQPSO)與神經元相結合的算法對模型參數進行辨識，該算法是針對粒子群(PSO)算法缺點而提出的。粒子群算法是由美國電氣工程師Eberhart博士和社會心理學家Kennedy博士[13]于1995年提出，算法的主要思想源于對鳥類群體行為的研究, 利用生物學家Heppner提出的模型構建粒子群模型和算法。在PSO算法中，搜索域空間中的鳥被作為優化問題的解看待，稱為”粒子”。每一個粒子對應著優化問題的一個適應值，優化問題的最優解對應著適應值最好的粒子，適應值根據具體問題而定，比如在控制問題中，往往把具體性能指標作為適應值。在該算法中，粒子的速度決定其飛行方向和距離，粒子通過追尋群體中的最優粒子最終完成解空間的搜索。PSO算法自提出以來，由于其具有計算簡單，控制參數少，易于實現等特點，引起了國內外很多學者的關注，開展了一系列研究，在算法的理論分析、性能改進以及應用等方面取得了豐碩成果。

但在PSO算法中，粒子的運動狀態由位置和速度確定，隨著時間的推移，粒子的運動軌跡是既定的；同時粒子的速度受到一定限制，使得粒子的搜索空間最終是一個有限的并逐漸減小的區域，因此不能覆蓋整個可行解空間，所以該算法不能保證全局收斂。

針對這一缺點，根據粒子群的基本收斂性質[14]，我國學者孫俊[15]利用量子力學中的相關理論為基礎，設計了一種基于全局水平的參數控制方法，提出了量子行為粒子群優化(QPSO)算法。算法基本原理是通過模擬量子力學中粒子在勢場中向勢能最低點的移動建立搜索機制。本文在該算法的基礎上加以改進，并將之用于我們構建的非線性系統的參數辨識。改進的基本思想是通過增加粒子的多樣性，使得算法在前期具有更好的全局探索能力，算法后期，在基本確定待辨識參數的最優值區域的基礎上，為提高辨識精度，將搜索到的最優值傳遞給神經元作為初始權值，利用神經元增強算法的局部開發能力。具體的改進算法如下：

1)確定辨識參數個數，數據長度以及參數初始化。

2)通過如下公式計算群體中每一個粒子的適應度。

3)判斷是否收斂或達到最大迭代次數。

4)根據如下改進公式進行等概率的進化粒子。

Xi,j(t+1)=Xi,j(t)±α·|Cj(t)-Xi,j(t)|·ln[1/ui,j(t)]

Xi,j(t+1)=pi,j(k)±α·|pi,j(t)-Xi,j(t)|·ln[1/ui,j(t)]

其中：

ui,j(t)～U(0,1)φj(t)～U(0,1)

Pi(t)=[Pi,1(t),Pi,2(t),...Pi,N(t)]表示個體最好位置。它由下式確定:

(2)

Xi(t)=[Xi,1(t),Xi,2(t),...Xi,N(t)]表示在t時刻第i個粒子的位置，其中i=1,2,..M.G(t)=[G1(t),G2(t),...GN(t)]表示群體的全局最好位置，且G(t)=Pg(t)，g為處于全局最好位置的粒子的下標，g∈{1,2,...M}。群體的全局最好位置由如式(3)、(4)確定：

(3)

G(t)=Pg(t)

(4)

5)重復步驟2)～步驟4)共M次。

6)經過M次學習后將當前最優值傳給單神經元，按下列公式對最優解進行局部精確搜索。

(5)

(6)

其中:i=1,...,an+bn,η為學習率。wk={ai(i=1,...an),bi(i=1,...bn)}

[X1,...Xn,...Xan+bn]=[Y(k-1),...Y(k-an),...

U(k-1),...U(k-bn)]

an=3,bn=2

7)轉至步驟3)。

8)輸出辨識結果，即參數a2,a1,a0,b2,b1的值給系統。

2 模型的穩定性分析

系統的穩定性是系統的固有屬性，線性系統由系統的結構和參數決定，與系統的初始條件和外界擾動的大小無關。而非線性系統的穩定性則與初始條件與外界擾動的大小有關。一個系統能夠正常工作的前提條件是系統必須是穩定的。在此，針對我們所構建的模型(1)特點，我們選用胞映射為工具討論系統的穩定性。胞映射的理論基礎是馬爾科夫鏈理論、龐加萊點映射原理和相空間重構原理[16-17],它是目前處理和研究時間序列問題，尤其是非線性動力系統全局性能問題的成功方法之一。胞映射方法最先在20世紀80年代初由Hsu[16]提出，方法以其快速，準確和適用范圍廣等特點很快成為研究的熱點，而被很多研究者所關注。它的基本思想是將動力系統狀態空間離散化形成許多小的幾何體(胞)，全部胞組成的集合形成胞空間。通過劃分將動力系統狀態空間轉化為胞空間后，動力系統中狀態的轉移就對應為胞之間的轉移。進而通過對胞之間轉移關系的研究完成對原動力系統的相應研究。這就是胞映射方法的基本思想。胞映射的實現方式是通過構造局部的轉移矩陣，以迭代算法進行長時間演化，實際上是用多周期運動來逼近所要描述的非線性運動，最終得到系統狀態變量在研究區域的極限概率分布。胞映射的突出特點是只要計算內存容量和速度允許，對非線性動力學系統能夠以任意尺度來劃分胞元并進行胞映射全局分析[18-21]。通過胞映射分析系統的穩定性，為進一步設計控制器作準備。

針對本文所提的模型特點，用已知輸出胞和對應時刻的輸入，預測未來系統輸出所在胞，胞劃分的依據是通過聚類算法得到的，而不同于以往的對狀態空間等分劃分。由于涉及到依賴未來系統輸出尋找其對應的胞，進而通過其胞標號作為反饋信息找到所對應的胞中心從而代入模型(1)得到下一時刻的系統輸出胞，因此在此構建的模型不同于以往的線性CARMA模型。針對這種非線性的胞模型特點，我們利用胞映射具有處理非線性系統穩定性的優勢特點，優勢在于平衡胞的確定和其相應的吸引域的確定比較容易。通過轉化為其對應的狀態空間模型進而構建胞映射，通過尋找平衡胞以及搜索到達以平衡胞為目標胞的吸引域，對上述特點的系統開展穩定性分析。首先對上述得到的離散時間模型，通過[12]，轉化為如下的狀態空間胞映射模型。

yc(k)=x1c(k)

通過考察胞空間中每一個胞，在系統對應的自治模型映射下的像胞，得到所有的胞對。對于這些胞對，找到所有的平衡胞，搜索以平衡胞為目標胞，到達目標胞的所有路徑。路徑上的所有胞構成了吸引域，對應于系統的穩定域。利用動態規劃的方法，通過尋找到達以平衡胞為目標胞的所有路徑，從而將狀態空間中的所有胞進行分類，分類后的胞與穩定性的對應關系如下：

1)平衡胞。對一簡單胞映射C，如果胞元素Ze滿足：

Ze=C(Ze)

則稱Ze為系統的平衡胞。由于系統的穩定性是指：處于平衡點的系統經擾動后，偏離平衡位置，在擾動消失后，系統最終回到平衡點?？疾煜到y的平衡點，對應胞空間中的平衡胞，因此需要對所有的胞對。找到平衡胞，并將之作為在系統穩定性分析中的目標胞，所有目標胞構成目標胞集。從此胞集出發搜索到達此胞集的所有路徑，從而確定系統相應的的穩定域。

2)簡單胞映射。對每個胞，經過胞映射后，對應的像胞只有一個，該映射的含義是，從原胞到像胞的映射關系是確定的。這種胞映射的算法簡單，也使得利用此方法對分析系統的全局性能很有效。但對于復雜系統，利用此方法所需要的胞數量很大。以此簡單胞映射為基礎，在以后的研究中，要考慮利用廣義胞映射、圖胞和馬爾科夫鏈理論為基礎進一步分析系統的穩定性。

2)周期運動和周期胞。設Cm為映射C進行m次映射，若k個不同胞:Z*(j)(j=1,2,…k)的一個序列滿足：

Z*(m+1)=Cm(Z*(1))

m=1,2,…k-1

Z*(1)=Ck(Z*(1))

則說該序列構成一個周期為k的周期運動，簡稱P-K運動。序列中的每一胞Z*(j)(j=1,2,…k)稱為周期為k的周期胞，簡稱P-K胞。周期胞對應于系統的周期運動。在非線性系統的分析中，實質是用周期胞的周期運動來逼近系統的動態軌跡。

3)暫態胞和吸引域(漸近穩定域)。經過有限次映射并最終到達周期胞的那種非周期胞稱之為暫態胞。如果r是一個使胞Z的r次映射成為一個P-K運動的P-K胞Z*(j)，即Cr(Z)=Z*(j)，則稱Z距離周期運動有r步遠。所有距離周期運動有r步遠或者小于r步遠的暫態胞的集合叫做周期運動的r步吸引域。而P-K運動的吸引域則是當r→∞時所有這種吸引域的總和。在本文中，我們要做的就是搜索所有能到達以周期為1的平衡胞為目標胞的暫態胞。所有的這些暫態胞集構成了系統的漸進穩定域。

因此，可通過胞映射的搜索對如圖3狀態胞進行分類，通過考察平衡胞，周期胞以及吸引域達到分析系統穩定性的目的。

圖3 胞映射圖

3 數值仿真

利用第1節改進的量子行為粒子群優化與神經元相結合的算法對構建的模型(1)，通過采集的歷時兩年的安陽燒結機的燒結臺車上靠近燒結終點的東西風箱溫度以及點火溫度進行仿真實驗，并利用第2節以胞映射為主要工具對模型相應的自治系統進行穩定性分析。IQPSO算法的各參數設定如下：辨識參數個數為N=5個，首先從最小二乘法中將參數的范圍確定，a2,a1,a0,b2,b1∈[-2,2]按數據長度為M=525規模隨機初始化5維粒子，參數范圍均為[-2，2]，Q=1 000，α=0.8。利用上述第二節的改進算法，得到如下(圖4，圖5)的辨識結果。圖4表示原始時間序列，QPSO辨識得到的參數對應的時間序列以及改進的的辨識參數對應的時間序列。圖5為改進前后得到的辨識參數對應的誤差時間序列。無論是從圖中定性觀察，還是從得到的數據定量計算，都可以得出改進的算法是有效的這樣的結論。

最后，經辨識得到的模型為:

y(k+3)-0.4231yc(k+2)-0.1430yc(k+1)-
0.1445yc(k)=-0.0591u(k+2)+0.0188u(k+1)

圖4 改進前后辨識參數對應的以及實際時間序列

圖5 改進前后辨識參數對應的誤差時間序列

通過差分方程向離散時間狀態空間模型的轉換，轉換后的狀態空間模型為:

y(k)=x(k)

要討論其系統的穩定性，則需要令u(k)=0考察其對應的自治系統。通過對胞空間中所有胞經過胞映射得到像胞，與其組成的胞映射對進行搜索算法，利用動態規劃最優原理確定系統的穩定域。首先對狀態空間進行胞劃分，由差分方程的輸出和點火溫度的輸入確定最終的狀態空間每一維的范圍，即對[-5 5]*[-5 5]*[-5 5]均勻劃分為10*10*10=1 000個胞。胞的劃分依據對輸出工況空間通過聚類反映工況特點的性質得到的。相應地，每一維對應的胞中心如下圖6。圖7表明對劃分后的每一個胞號作為橫坐標通過狀態空間模型得到的像胞號作為縱坐標得到的胞映射圖形。

圖6 各維劃分10個胞時所對應的胞中心

圖7 各維劃分10個胞時對應胞映射

將胞映射的對應關系繪制成如圖7，更能直觀地看出系統的動態特性，并且能從對應的胞映射表中得到平衡胞，即原胞和像胞為相同胞。同時從圖8中也可以得到平衡胞，即只要將圖7中胞映射圖與對角線引出的直線相交，相交位置的胞即為平衡胞。例如，將原點所在胞[5,5,5]和[6,6,6]胞標號為445和556的平衡胞作為目標胞，找到目標胞的所有吸引域即為系統對應此平衡胞的穩定域。

圖8 平衡胞

從圖8可以看出相交的點對應的胞均為平衡胞，但結合胞映射表進行搜索，得知并不是所有平衡胞都存在吸引域，不存在吸引域的平衡胞對應于著系統的不穩定平衡點，這和系統存在不穩定的平衡點理論是相符的。除穩定的平衡胞，這類胞也稱為可達的平衡胞，和其對應的吸引域外，胞空間中其余平衡胞和非平衡胞對應于狀態空間中的不穩定區域。利用matlab進行編程，進一步得到穩定平衡胞和其吸引域對應的胞的個數為616，占整個胞數的74/125，穩定平衡胞和其吸引域對應于狀態空間中的穩定域。因此，可以看出，以胞映射為工具可以實現對這類復雜系統的穩定性分析，并且分析過程得到簡化。

4 結束語

針對如燒結機等一類復雜系統難建模的特點，本文利用數據驅動的建模方法，通過采集2年的工藝參數數據，建立受控回歸滑動平均(CARMA)胞映射模型。在確定模型結構的基礎上，利用改進的量子行為粒子群與神經元結構相結合的算法，對模型進行參數辨識。進一步，將離散時間的CARMA模型轉為狀態空間模型。在狀態空間中進行胞劃分，所有胞的集合構成胞空間。利用簡單胞映射，對胞空間中的每一個胞以胞中心作代表經過狀態空間模型得到像胞，構建胞映射。利用平衡胞定義，找到所有的平衡胞。以平衡胞為目標胞搜索所有到達平衡胞的路徑，對應系統的穩定域。實現以胞映射理論討論此類復雜系統的穩定性目的。最后，給出了仿真實例，實驗證明此方法是可行的，因此為這類復雜系統的動態特性分析提供了一條思路。在此基礎上，可進一步開展針對控制器的設計工作。

[1]HuangLG.Influenceofsinteringtemperatureoncrystallinephasessynthesizedfromfinesiliconpower[J].Chinaceramics,2007,43(9):23-25.

[2]KangGZ,JiaBG,WangLMetal.Theuseofrotarykilninferroalloysproduction[J].Ferro-alloys,2011,218(3).

[3] 蔡 健,王百紅,梁慧強,等.電解鋁生產電解電容器負極箔的關鍵技術[J].輕合金加工技術.2004,32(9):17-21.

[4]AdamyJ.A.flemming.softvariable-structurecontrols:asurvey[J].Automatica, 2004,40: 1821-1844.

[5]KonstantinosK,JohnL.ModelingandanalysisofDNAreplication[J].Automatica, 2011,47, 1156-1164.

[6]HenryY,GranthamK,NelsonH.Performancevaluationformotif-basedpatternedtexturedefectdetection[J].IEEETransactionsonAutomationScienceandEngineering, 2010, 7(1):58-72.

[7]StachW,WitoldP,LukaszA.Learningoffuzzycognitivemapsusingdensityestimate[J].IEEETransactionsonsystems,man,Cybernetics-partB:cybernetics, 2012, 42(3): 900-912.

[8]HouZS,JinST.Anoveldata-drivencontrolapproachforaclassofdiscrete-timenonlinearsystems[J].IEEEtransactionsoncontrolsystemstechnology, 2011, 19(6):1549-1558.

[9]BeghiA,BrignoliR,CecchinatoL,etal.Data-drivenfaultdetectionanddiagnosisforHVACwaterchillers[J].ControlEngineeringPractice,2016,53(8):79-91.

[10]WangJS,YangGH.Data-drivenoutput-feedbackfault-tolerantcontrolforunknowndynamicsystemswithfaultschangingsystemdynamics[J].JournalofProcessControl,2016,43(7):10-23.

[11]WangYJ,YangD,ZhangX,etal.Probabilitybasedremainingcapacityestimationusingdata-drivenandneuralnetworkmodel[J].Journalofpowersources,2016,315(5):199-208.

[12]LiuBao.ModernControlTheory[M].Beijing:ChinaMachinePress,2000.

[13]KenneyJ,EberhartRC.Particleswarmoptimization[C].ProceedingofIEEEInternationalConferenceonNeuralNetworks,1995,1942-1948.

[14]ClercM,KennedyJ.Theparticleswarm:explosion,stabilityandconvergencyamulti-dimensionalcomplexspace[J].IEEETransactionsonEvolutionaryComputation, 2002，6:58-73.

[15]SunJ,XuWB,FengB,Aglobalsearchstrategyofquantum-behavedparticleswarmoptimization[c].In:ProceedingsofIEEEConferenceonCyberneticsandIntelligentSystems,Singapore,2004, 111-116.

[16]HsuCS.Ageneralizedtheoryofcell-to-cellmappingfornonlineardynamicalsystems[J].ASMEJ.Appl.Mech, 1981,48:634-842.

[17]HengL,LieY,XieY.ThecontinuousPoincare-likecellmappingmethodanditsapplicationtononlineardynamicsanalysisofabearing-rotorsystem[J].TribologyInternational,1998,31(7):369-375.

[18]SunJQ,HsuCS.Globalanalysisofnonlineardynamicalsystemswithfuzzyuncertaintiesbythecellmappingmethod[J].Comput.Meth.Appl.Mech.Eng.,1990,83(2):109-120.

[19]SunJQ,HsuCS.Thegeneralizedcellmappingmethodinnonlinearrandomvibrationbaseduponshort-timeGaussianapproximation[J].J.Appl.Mech.,1990,57:1018-1025.

[20]HongL,XuJX.Crisesandchaotictransientsstudiedbythegeneralizedcellmappingdigraphmethod[J].Phys.Lett.A,1999,262:361-375.

[21]HongL,XuJX.Discontinuousbifurcationsofchaoticattractorsinforcedoscillatorsbygeneralizedcellmappingdigraph(GCMD)method[J].Int.J.BifurcatChaos,2001,11:723-736.

Dynamical Performance Analysis for a Class of Complex Systems Based on Data Driven

Guo Lingli1,2,Wang Mushu1，Wang Ye1

(1.School of Automation and Electrical Engineering, University of Science &Technology Beijing, Beijing 100083,China; 2.Department of Mathematics Changzhi College, Changzhi 046011,China)

A model and the stability analysis are proposed for a class of dynamical systems which is difficult in modeling due to the complex physicochemical change by data-driven method. A CARMA cell mapping model is established, whose parameters are identificated by quantum-behaved particle swarm optimization. To avoid the premature convergence of the IQPSO algorithm, the particles position are updated and accurate search capacity is enhanced by applying the neural networks that train the optimal value as the initial weighted value. Based on the model, the stability performance is analyzed by cell mapping theory and dynamic optimization principle, which is helpful to find the whole periodic cell and attractive domain. Simulation studies are included to demonstrate the effectiveness of the proposed approach.

data driven; modeling; particle swarm algorithm; parameter identification; cell mapping; stability

2016-12-12；

2017-01-19。

國家自然科學基金(61175122)。

郭玲利(1978-)，女，山西長治人，博士研究生，講師，主要從事控制科學與工程方向的研究。

1671-4598(2017)04-0220-05

10.16526/j.cnki.11-4762/tp.2017.04.060

TP273

A