加權有向圖社區(qū)發(fā)現(xiàn)的子系統(tǒng)劃分

楊曉峰,謝 巍,張浪文?

(1.華南理工大學自動化科學與工程學院,廣東廣州 510640;2.廣州市標準化研究院,廣東廣州 510110)

1 引言

現(xiàn)代化學過程通常由相互連接的操作單元組成,它們通過能量和信息流緊密地集成在一起.先進控制算法在應用于這些過程時應滿足預定的目標,達到預期的安全,可持續(xù)性和經濟性水平[1].然而,傳統(tǒng)的集中或分散控制架構難以適用于復雜過程系統(tǒng)[2].

近十年來,為提高計算效率,維護靈活性和容錯性,分布式控制架構得到快速發(fā)展[3].分布式控制系統(tǒng)設計包括兩個關鍵步驟[4]:1)合理地將整個系統(tǒng)劃分成較小的子系統(tǒng),需要保證子系統(tǒng)的能控性、能觀性和靈活性等;2)針對每個子系統(tǒng)設計局部控制器及其協(xié)調算法.在過去的十年里,分布式控制算法(尤其是分布式模型預測控制)吸引了非常廣泛的研究關注[5].然而,與其同等重要的子系統(tǒng)劃分問題受到的關注相對較少[6].

近幾年來,子系統(tǒng)劃分方法取得了一些成果[7–9].文獻[7]提出了一種基于層次聚類的分散控制系統(tǒng)分解方法.文獻[8]提出了一種既考慮耦合強度又考慮結構緊密性的輸入–輸出配對方法.另外,文獻[9]研究了基于無權重有向圖的社區(qū)發(fā)現(xiàn)子系統(tǒng)劃分算法.在社區(qū)發(fā)現(xiàn)算法中,將一個大規(guī)模系統(tǒng)看作一個網(wǎng)絡,通過將網(wǎng)絡分解為較小的子網(wǎng)絡,建立面向分布式控制或估計的子系統(tǒng)模型.

考慮到不同狀態(tài)和輸出變量之間的結構緊密性,已有研究嘗試將非線性系統(tǒng)分解為較小的單元,然后設計分布式狀態(tài)估計[10–11].文獻[10]提出了一種面向分布式狀態(tài)估計的過程網(wǎng)絡子系統(tǒng)劃分方法.在文獻[11],分布式狀態(tài)估計和控制在一個架構內被共同考慮.上述面向分布式狀態(tài)估計的子系統(tǒng)劃分方法均只考慮過程系統(tǒng)的物理拓撲,而未考慮不同變量之間的連接強度.研究發(fā)現(xiàn),如果不考慮變量間連通性的強弱,可能導致子系統(tǒng)劃分的次優(yōu)或不適當,特別是對于連接強度差異較大的系統(tǒng)影響更大[12].

基于上述分析,本文提出一種面向分布式狀態(tài)估計的子系統(tǒng)劃分方法,構建了加權有向圖,既考慮結構的緊密性,又考慮系統(tǒng)狀態(tài)/輸出之間的連接強度.基于加權有向圖,利用社區(qū)發(fā)現(xiàn)算法將所有系統(tǒng)變量劃分為較小的群組,使得劃分得到的子系統(tǒng)內部關聯(lián)較強,而子系統(tǒng)之間的耦合強度較弱.

2 問題描述

考慮如下復雜非線性系統(tǒng)的子系統(tǒng)劃分問題:

其中:x ∈Rnx表示狀態(tài)向量,y ∈Rny表示輸出向量,f和h分別表示非線性系統(tǒng)狀態(tài)和輸出方程.

本研究的目標是將系統(tǒng)(1)劃分成多個子系統(tǒng),得到第i個子系統(tǒng)的模型,以設計分布式狀態(tài)估計器

其中:x(i)∈Rnx(i)表示第i個子系統(tǒng)的狀態(tài)向量,y(i)∈Rny(i)表示第i個子系統(tǒng)的測量輸出,X(i)表示所有與子系統(tǒng)i直接關聯(lián)的子系統(tǒng)狀態(tài),i=1,···,p,p是劃分得到子系統(tǒng)的個數(shù).在本設計中,假設p是已知的且數(shù)量少于測量輸出的個數(shù),即pny.

由于本文子系統(tǒng)的劃分是面向分布式狀態(tài)估計設計,對于劃分的子系統(tǒng),需要進行能觀性的判斷.為了檢驗非線性系統(tǒng)(1)的能觀性,可以判斷該系統(tǒng)的能觀矩陣是否滿秩[13],能觀矩陣通過如下方式構建:

其中:Lfh 表示方程h 對方程f 的Lie導數(shù),其定義為表示h對f的r階Lie導數(shù),定義為如果對于x ∈X,有rank(Q(x))=nx,那么,系統(tǒng)(1)在范圍X內是局部可觀的[13],否則系統(tǒng)(1)是不可觀的.

社區(qū)發(fā)現(xiàn)算法是一種將大規(guī)模復雜網(wǎng)絡劃分成多個子網(wǎng)絡的有效工具[14],通過由文獻[15]定義的模塊度來評價一個網(wǎng)絡劃分的好壞.考慮一個帶有N個節(jié)點的網(wǎng)絡,其不加權的鄰接矩陣A是一個N ×N矩陣,定義元素Aij為

通常,不考慮節(jié)點自身的回路,也就是說鄰接矩陣的主對角線上的元素為0.用模塊度來評價劃分好壞,模塊度Q的范圍為0到1,社區(qū)發(fā)現(xiàn)算法通過最大化模塊度來找到最優(yōu)群組(子系統(tǒng))結構,使得劃分的子系統(tǒng)內部的連接數(shù)最多,而群組之間的連接數(shù)最少.一般模塊度Q在0.3到0.7之間時被認為是較好的劃分方法,否則表示該系統(tǒng)不適合再進行劃分.

然而,實際應用的子系統(tǒng)劃分中,僅考慮劃分得到的子系統(tǒng)內部連接邊數(shù)最多是不夠的,還需要考慮連接邊的連接權重.本文將對系統(tǒng)構造帶權重的有向圖,并利用模塊度最大化方法對系統(tǒng)進行劃分,以保證各個劃分子系統(tǒng)之間的連接權重最小.

3 基于加權有向圖的社區(qū)發(fā)現(xiàn)子系統(tǒng)劃分

本節(jié)將提出一種系統(tǒng)化的子系統(tǒng)劃分方法,考慮將非系統(tǒng)(1)劃分成多個子系統(tǒng),并設計滾動時域估計算法對子系統(tǒng)狀態(tài)進行分布式狀態(tài)估計.

本文方法的主要結構如圖1所示.首先,將非線性系統(tǒng)(1)用加權有向圖進行描述,其中系統(tǒng)狀態(tài)和測量輸出變量看成網(wǎng)絡中的節(jié)點,這些節(jié)點通過加權連接邊進行連接,連接邊上的權重反映了節(jié)點之間的連接強度.然后,基于加權有向圖,利用社區(qū)發(fā)現(xiàn)算法找到最大模塊度,將系統(tǒng)劃分成多個子系統(tǒng)群組.最后,通過檢測劃分子系統(tǒng)的能觀性,獲得子系統(tǒng)劃分結果.

3.1 加權連接邊的定義

本文將利用加權有向圖對系統(tǒng)(1)進行描述,具體來說,將所有的狀態(tài)和測量輸出變量當做有向圖的節(jié)點,這些節(jié)點通過有向邊進行連接.令fi表示向量方程f的第i(i=1,···,nx)行,hj表示向量方程h的第j(j=1,···,ny)行,xi表示狀態(tài)向量x的第i(i=1,···,nx)行,yj(j=1,···,ny)表示輸出變量y的第j行.

在已有的無權重有向圖中,不同變量(節(jié)點)之間的連接強度未被考慮.為獲得權重,對非線性系統(tǒng)在工作點xs求一階偏導得到有向圖中連接邊的敏感度

圖1 本文子系統(tǒng)劃分結構Fig.1 Flowchart of the proposed subsystem decomposition

定義如下連接邊的權重:

狀態(tài)變量xi到另一個狀態(tài)變量xl連接邊的權重定義為(l,i=1,···,nx)

狀態(tài)變量xl到測量輸出變量yj連接邊的權重定義為(l=1,···,nx,j=1,···,ny)

其中α是一個介于[0 1]的參數(shù),用于調節(jié)連接邊的權重在子系統(tǒng)劃分過程中的重要性,當α=0時,所有連接邊的權重為1,即等效于已有的無權重有向圖;當α增加時,連接邊的權重將在子系統(tǒng)劃分中更多地被考慮.當S(xl,yj)=0或S(xi,xl)=0時,連接邊的權重為無窮,即兩個節(jié)點之間沒有直接連接.

由于不同的α值將影響子系統(tǒng)的劃分結果,如何選擇合適的參數(shù)α對子系統(tǒng)劃分相當重要.總體而言,當有向圖中所有連接邊的敏感度差異很小時,在子系統(tǒng)劃分時應重點考慮連接邊的數(shù)量,此時選擇較小的α即可.相反,當所有連接邊的敏感度差異很大時,應更多地考慮連接邊的權重對子系統(tǒng)劃分的影響[16].

3.2 最短路徑

為實施社區(qū)發(fā)現(xiàn)算法,需要構造加權鄰接矩陣,因此首先應找到有向圖中的一個節(jié)點到另一個節(jié)點的最短路徑.令P表示一個狀態(tài)變量到另一個狀態(tài)變量或輸出變量的路徑,在這條路徑P中的所有連接邊表示為e ∈P.路徑L(P)的長度可以計算得到

其中:l=1,···,nx, j=1,···,ny,Plj,Plj表示從狀態(tài)xl到輸出yj所有路徑和其中一條路徑.本文利用Dijkstra’s算法[17]找到任意兩個節(jié)點的最短路徑,其計算復雜度為O(E log V),其中V 和E分別是節(jié)點的數(shù)量和連接的邊數(shù).

3.3 構建鄰接矩陣

本節(jié)將構建非線性系統(tǒng)(1)的鄰接矩陣.對于系統(tǒng)(1)所形成的有向圖,總共有na(其中na=nx+ny)個頂點.令ca為狀態(tài)和測量輸出向量,定義為ca=[x1,···,xnx,y1,···,yny].

那么,所構建的加權鄰接矩陣Aw將被用于社區(qū)發(fā)現(xiàn)算法,以對系統(tǒng)(1)進行劃分.上述定義的加權鄰接矩陣(15)是無權重鄰接矩陣[14]的擴展,其主要區(qū)別是:a)加權鄰接矩陣考慮了兩個節(jié)點之間連接的權重,而文獻[14]只考慮了兩個節(jié)點之間的連接度;b)本文方法中,考慮了兩個節(jié)點之間的最短路徑,因此,鄰接矩陣中的0元素表征了兩個節(jié)點之間沒有連接.

3.4 尋找最大模塊度

模塊度Qw既考慮了靈敏度,又考慮了邊的數(shù)量.子系統(tǒng)劃分問題等價于通過尋找最大化的模塊化Qw進行社區(qū)結構發(fā)現(xiàn).首先,需要進行初步的群組配置,群組結構的初始化基于如下假設:i)每個群組至少有一個測量輸出;ii)應將直接影響yi的系統(tǒng)狀態(tài)分配給屬于yi同一個群組.初始化群組結構步驟如下:

步驟1定義關于復雜系統(tǒng)的群組結構sw(0)=[cx1(0) ··· cxnx(0) cy1(0) ··· cyny(0)],其中cη(0)表示每個節(jié)點xi(i=1,···,nx)或yj(j=1,···,ny)的群體標簽.應用快速二分法(如文獻[18])尋找最大模塊度Qw時,在sw中的第i個變量被分配到第i個群體,i=1,···,na,即(0)=[1 ··· na].

步驟2對于每個狀態(tài)變量xi,i=1,···,nx,找到所有和狀態(tài)變量xi關聯(lián)的測量輸出,在ˉsw(0)中將相應的元素標為相同編號.

步驟3對于每個測量輸出yj,j=1,···,ny,找到所有和輸出方程hj(·)關聯(lián)的狀態(tài)變量,并將這些狀態(tài)標量標為和yj相同的編號.

步驟4所得到的分組(0)即可作為實施社區(qū)發(fā)現(xiàn)算法的初始群組.

得到的初始群組結構將用于尋找子系統(tǒng)分解中模塊度的最大值.在所提出的方法中,子系統(tǒng)(即式(2)中的p)的數(shù)量應預先設定.在本文中,快速二分法用于最大化模塊度[18],其步驟如下:

步驟1對每個節(jié)點i實施如下算法:

步驟1.1對于節(jié)點i的每個鄰居節(jié)點j,通過將節(jié)點i從當前群組移動到節(jié)點j所在的群組,計算模塊度的變化值?Qw.

步驟1.2找到上一步中最大的?Qw>0,將頂點i放入到相應的群組中(即將節(jié)點i編入最大?Qw相應的分組),得到新的分組.

步驟2令nc(k+1)作為節(jié)點集結后的群組數(shù)量,如果nc(k+1)

3.5 分布式滾動時域估計設計

為對子系統(tǒng)的劃分結果進行分析,本文將對所得到的子系統(tǒng)應用于分布式滾動時域估計算法,并說明考慮加權有向圖的必要性.本文提出分布式滾動時域估計算法,局部估計器通過共享的通信網(wǎng)絡進行信息交互.子系統(tǒng)i的局部估計器為求解如下優(yōu)化問題:

4 實例應用

本節(jié)考慮將所提出的子系統(tǒng)劃分及分布式滾動時域估計算法應用于反應–分離過程系統(tǒng)[19],該過程包括兩個連續(xù)攪拌反應釜和一個快速分離釜(見圖2).第一個連續(xù)攪拌反應釜包含反應物J,其給定蒸汽為F10,然后產出產品K.然后,產品K在第2個反應釜中轉換成最終產品L.

圖2 反應–分離過程系統(tǒng)Fig.2 Reactor-separator process

其中:xAi,xBi表示產品J,K在容器i(i=1,2,3)中的質量分數(shù);xC3表示產品L在反應釜3中的質量分數(shù);Ti容器i(i=1,2,3)中液體的溫度;T10,T20表示輸入到容器1和2的蒸汽溫度;F1,F2表示物料從容器1和2流出的流速;F10,F20表示物料從容器1和2流出的穩(wěn)態(tài)流速;Fr,Fp分別表示產品的循環(huán)和廢棄流速;Vi表示容器i(i=1,2,3)的體積;E1和E2表示反應過程1和2的催化劑;k1和k2表示反應過程1和2的指前因子值;?H1和?H2表示反應過程1和2的熱量;αA,αB,αC表示物質J,K,L的相對揮發(fā)率;Q1,Q2,Q3分別表示容器i(i=1,2,3)的熱量輸入;Cp熱容量;R表示氣體常數(shù);ρ表示溶液濃度.

相關模型參數(shù)參見表1.

表1 模型參數(shù)Table 1 Model parameters

該反應–分離過程的其中一個穩(wěn)定工作點為

為實現(xiàn)利用測量輸出Ti(i=1,2,3)對系統(tǒng)狀態(tài)xAi,xBi,Ti(i=1,2,3)的估計,首先對子系統(tǒng)進行劃分,并和已有的子系統(tǒng)劃分方法(文獻[14])進行比較,得到如表2所示的劃分結果.

表2 子系統(tǒng)劃分結果Table 2 Decompositions for the CSTR process

在應用分布式估計算法時,如果進行實時信息交互,分布式估計器的通信成本較高,因此希望子估計器之間可以在較少通信交互的同時獲得較好的估計效果.如圖3所示,應用分布式滾動時域估計算法時(通訊間隔n=1),基于子系統(tǒng)劃分2的估計效果比基于子系統(tǒng)劃分1獲得了更好的效果,本文提出的基于加權有向圖的子系統(tǒng)劃分方法比未考慮連接權重的分解算法具有更小的估計誤差.

圖3 當n=1時,實際狀態(tài)軌跡(藍色實線)、基于分解1的狀態(tài)估計(綠色虛線)和基于分解2的狀態(tài)估計(紅色點劃線)Fig.3 The trajectories for n=1 of the actual states(blue solid lines),estimated states based on Decomposition 1(green dashed-dot lines)and Decomposition 2(red dashed lines)

圖4 當n=3時,實際狀態(tài)軌跡(藍色實線)、基于分解1的狀態(tài)估計(綠色虛線)和基于分解2的狀態(tài)估計(紅色點劃線)Fig.4 The trajectories for n=3 of the actual states(blue solid lines),estimated states based on Decomposition 1(green dashed-dot lines)and Decomposition 2(red dashed lines).

考慮不同的通信間隔n=1,2,3,對分布式滾動時域估計算法進行測試.如表3和圖4(通信間隔n=3)所示,隨著通信間隔的增加,基于分解1的狀態(tài)估計效果受到較大的影響.而對于基于本文方法的分解2,由于子系統(tǒng)之間的關聯(lián)較弱,因此較大的通信間隔對估計性能的影響并不大.其現(xiàn)實意義是在保證估計效果的同時,可以減輕分布式狀態(tài)估計設計過程中子系統(tǒng)通信的次數(shù),降低通信成本.

表3 不同通信間隔對估計性能的影響Table 3 Summary of the estimation performance

5 總結

本文提出一種基于加權有向圖的社區(qū)發(fā)現(xiàn)子系統(tǒng)劃分方法,并應用于分布式狀態(tài)估計設計中.基于耦合關聯(lián)的強度研究了一種系統(tǒng)化的復雜大系統(tǒng)分解方法,根據(jù)系統(tǒng)的動態(tài)模型自動構建系統(tǒng)的加權有向圖,進而獲得系統(tǒng)的劃分,得到的子系統(tǒng)內部的關聯(lián)遠強于子系統(tǒng)之間的關聯(lián).這類分解方能夠為在復雜系統(tǒng)的分布式狀態(tài)估計和控制奠定基礎,同時也對研究復雜網(wǎng)絡的特性具有重要作用.未來可擴展到帶系統(tǒng)輸入的復雜系統(tǒng)劃分方法,及其對分布式控制策略下的性能提升研究,通過分布式狀態(tài)估計和控制,對該子系統(tǒng)劃分方法進行系統(tǒng)的評估.