復(fù)雜非線性系統(tǒng)的子系統(tǒng)分解方法

謝苗苗，張浪文，2，謝巍，2

（1 華南理工大學自動化科學與工程學院，廣東廣州510641； 2 廣東道氏技術(shù)股份有限公司，廣東江門529400）

引 言

信息技術(shù)和工業(yè)技術(shù)的發(fā)展推動著整體工業(yè)體系發(fā)生巨大的變化，諸多系統(tǒng)朝著越來越復(fù)雜的方向發(fā)展，如化工系統(tǒng)、灌溉系統(tǒng)、智能交通系統(tǒng)等，其典型特征為結(jié)構(gòu)復(fù)雜、維數(shù)高、約束多、模型不確定、非線性。同時，對于這類系統(tǒng)控制的要求越來越高，要求獲得更好控制效果的同時滿足經(jīng)濟效益指標，通常傳統(tǒng)的PID 控制難以滿足系統(tǒng)對于工藝水平的要求，為滿足安全、可持續(xù)性，先進的控制系統(tǒng)對于這些過程來說是必要的[1]。

模型預(yù)測控制（model predictive control, MPC）憑借其在處理多變量、狀態(tài)/輸入約束、非線性等問題上的優(yōu)勢，成為系統(tǒng)先進控制方法的首選[2-3]。通常，由于結(jié)構(gòu)的復(fù)雜性、容錯性以及計算問題，集中式MPC 通常存在控制時效性差、靈活性差等問題，因而在復(fù)雜過程中難以適用[4-5]。對于大規(guī)模過程的整體控制，分散式MPC 是解決方案之一[6]。為了實現(xiàn)分散式MPC，首先將一個大規(guī)模的過程分解成多個子系統(tǒng)，然后對每個子系統(tǒng)獨立地設(shè)計控制，以達到減少控制器設(shè)計的復(fù)雜度[7-8]。分散式MPC的一個關(guān)鍵特性是子系統(tǒng)的控制器之間沒有或很少有信息交換，由于通?；み^程單元之間的耦合是不可忽略的，使得分散控制下的全局性能受到限制[9]。在此背景下，分布式MPC 被提出，在協(xié)作分布式MPC 中，子系統(tǒng)通常使用全局代價函數(shù)[10]，在每個采樣時刻，一個子系統(tǒng)MPC 都需要與其他所有子系統(tǒng)MPCs 通信，子系統(tǒng)MPCs 以迭代的方式計算它們的控制動作[11]，從而提高分布式MPC 的性能，諸多實際應(yīng)用說明分布式控制框架能夠提高計算效率、維護靈活性和容錯性的發(fā)展的趨勢。分布式控制系統(tǒng)的設(shè)計主要包括兩個關(guān)鍵步驟[12]：(1)將整個過程適當?shù)胤纸獬筛〉淖酉到y(tǒng)；(2)提出一種經(jīng)濟有效的控制算法，在此基礎(chǔ)上設(shè)計本地控制器并確定通信協(xié)議。

在過去的十五年里，分布式控制（特別是分布式MPC）算法的發(fā)展受到了極大的關(guān)注[13-14]。同時，同等重要的子系統(tǒng)分解得到的關(guān)注相對較少[15]。在文獻中，有一些關(guān)于分布式控制系統(tǒng)分解算法的重要結(jié)果[16]。文獻[17-18]針對過程網(wǎng)絡(luò)的控制分解，采用了層次聚類的方法；文獻[19]提出了一種面向系統(tǒng)分解的輸入-輸出配對方法，該方法能夠同時考慮耦合結(jié)構(gòu)和強度；文獻[20]則基于社區(qū)發(fā)現(xiàn)算法，提出了一種非線性模型預(yù)測控制的分布式優(yōu)化問題分解方法；文獻[21]提出了一種相對時間平均增益序列的網(wǎng)絡(luò)分解方法；文獻[22-23]提出了面向分布式滾動時域估計的子系統(tǒng)分解方法，并在污水處理系統(tǒng)中進行了應(yīng)用；文獻[24]則將現(xiàn)有的社區(qū)發(fā)現(xiàn)算法推廣到分布式狀態(tài)估計和控制的共同框架內(nèi)進行考慮。然而，上述基于社區(qū)發(fā)現(xiàn)算法的子系統(tǒng)劃分均只考慮了系統(tǒng)的關(guān)聯(lián)度，而忽略了不同變量之間的連接強度。先前的工作[25-26]提出了一種基于加權(quán)邊的群組檢測的子系統(tǒng)劃分方法，并研究了其在分布式狀態(tài)估計中的應(yīng)用，而對于面向分布式控制的加權(quán)社區(qū)發(fā)現(xiàn)算法尚未得到解決。

因此，本文同時考慮狀態(tài)、輸入和輸出之間的關(guān)聯(lián)性，研究面向化工過程分布式狀態(tài)估計和控制的子系統(tǒng)分解方法。具體地，首先基于化工過程形成有向網(wǎng)絡(luò)，在構(gòu)建的網(wǎng)絡(luò)中，狀態(tài)變量、輸入變量和測量輸出變量都被視為節(jié)點，提出基于社區(qū)結(jié)構(gòu)發(fā)現(xiàn)算法的子系統(tǒng)分解方法，并在化工過程實例的分解中進行驗證，設(shè)計分布式MPC 算法實現(xiàn)來說明方法優(yōu)勢。

1 預(yù)備工作

1.1 問題描述

針對復(fù)雜非線性系統(tǒng)的子系統(tǒng)問題，通?？紤]如下動態(tài)模型：

其中，x ∈Rnx、u ∈Rnu、y ∈Rny分別表示狀態(tài)、輸入和輸出向量，f 表示系統(tǒng)的狀態(tài)方程，h 表示非線系統(tǒng)的輸出方程。為實現(xiàn)分布式控制，需要將系統(tǒng)(1)劃分成如下描述的子系統(tǒng)模型：

其中，i=1,…,p，p 為子系統(tǒng)的個數(shù)，x(i)∈Rnx(i)表示第i 個子系統(tǒng)的狀態(tài)向量，u(i)∈Rnu(i)表示第i 個子系統(tǒng)的輸入向量，y(i)∈Rny(i)表示第i 個子系統(tǒng)的測量輸出，X(i)表示所有與子系統(tǒng)直接關(guān)聯(lián)的子系統(tǒng)狀態(tài)，U(i)表示所有與子系統(tǒng)直接關(guān)聯(lián)的子系統(tǒng)對應(yīng)的輸入。

1.2 有向圖

本文將構(gòu)建加權(quán)有向圖對系統(tǒng)(1)進行描述，有向圖的構(gòu)建方法可參考文獻[24]，將系統(tǒng)的所有狀態(tài)、輸入和測量輸出變量當作有向圖的節(jié)點，這些節(jié)點通過有向邊進行連接。令hj,j = 1,…,ny表示向量方程h 的第j 行；fi,i = 1,…,nx表示向量方程f 的第i 行；xi,i = 1,…,nx表示狀態(tài)向量x 的第i 行；uk,k = 1,…,nu表示輸出向量u 的第k 行；yj,j = 1,…,ny表示輸出向量y的第j行。

根據(jù)如下規(guī)則定義有向圖：

1.3 能觀性與能控性

對于分解得到的子系統(tǒng)(2)，需要進行能觀能控性檢測。對于非線性系統(tǒng)，能觀性的判斷可以檢驗如下能觀矩陣是否滿秩：

類似地，對于得到的子系統(tǒng)(2)，能控性的判斷可以檢驗如下能觀矩陣是否滿秩：

1.4 社區(qū)發(fā)現(xiàn)算法

社區(qū)發(fā)現(xiàn)算法用于將復(fù)雜網(wǎng)絡(luò)劃分成子網(wǎng)絡(luò)，其構(gòu)建方法可參考文獻[27-28]，假設(shè)一個網(wǎng)絡(luò)具有N個節(jié)點，其鄰接矩陣A 的各個元素Aij可以通過如下方式構(gòu)建：

同時，定義模塊度進行分解結(jié)果的評價：

1.5 分布式Lyapunov MPC算法

在對系統(tǒng)進行子系統(tǒng)劃分后，對其進行分布式模型預(yù)測控制來驗證子系統(tǒng)劃分的好壞，本文引入基于迭代子系統(tǒng)的分布式Lyapunov MPC[11]，即在時刻tk的實現(xiàn)如下步驟：

（1）控制器i接收和發(fā)送x(i)(tk)給它的下游鄰居子系統(tǒng)。同時，它通過上游鄰居的通信網(wǎng)絡(luò)獲得子系統(tǒng)狀態(tài)x(l)(tk)；

（2）迭代控制器評估。

①在迭代c(c ≥1)中，控制器i 基于其上游鄰居在迭代c-1 中評估的狀態(tài)軌跡評估其未來的輸入軌跡；

所采用的Lyapunov MPC i 的設(shè)計，在采樣時間tk、迭代c中，被表述為以下優(yōu)化問題：

2 子系統(tǒng)分解方法

本文復(fù)雜非線性系統(tǒng)的子系統(tǒng)劃分主要流程如圖1 所示。首先，將非線性系統(tǒng)(1)用加權(quán)有向圖進行描述；然后，構(gòu)建系統(tǒng)的加權(quán)有向圖，實施社區(qū)發(fā)現(xiàn)算法找到最大模塊Q，得到推薦的子系統(tǒng)劃分結(jié)果；最后，對所得到的子系統(tǒng)進行能觀能控性判斷。

圖1 復(fù)雜非線性系統(tǒng)的子系統(tǒng)的劃分流程Fig.1 Flowchart of the proposed subsystem decomposition method of nonlinear systems.

2.1 加權(quán)連接邊的定義

為考慮節(jié)點之間的連接強度，構(gòu)建如下有向圖連接邊的敏感度，可以通過對非線性系統(tǒng)在工作點(xs,us)求一階偏導(dǎo)得到：

其中，S(uk,xi)表示從輸入uk到狀態(tài)xi的敏感度,S(xi,xl)表示從狀態(tài)xi到另一個狀態(tài)xl的敏感度，S(xl,yj)表示從狀態(tài)xl到輸出yj的敏感度。

進而，構(gòu)建敏感度矩陣：

本節(jié)中，定義如下連接邊的權(quán)重。

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

輸入變量uk到狀態(tài)變量xi連接邊的權(quán)重定義為(k = 1,…,nu,i = 1,…,nx)：

狀態(tài)變量xl到輸出變量yj連接邊的權(quán)重定義為：

當S(xi,xl)= 0，S(uk,xi)= 0 或S(xl,yj)= 0 時，連接邊的權(quán)重為無窮，即兩個節(jié)點之間沒有直接連接。

2.2 最短路徑

為構(gòu)造加權(quán)鄰接矩陣，需要找到兩個節(jié)點之間的最短路徑。假設(shè)P 為兩個節(jié)點之間的路徑，那么所有連接邊表示為e ∈P，路徑L(P)的長度為：

假設(shè)xi到xl的所有路徑集合表示為，其最短路徑為：

類似地，從輸入變量uk到狀態(tài)變量xi的路徑長度Lki(Pki)為：

相應(yīng)從輸入變量uk到狀態(tài)變量xi的最短路徑(uk,xi)為：

其中，k = 1,…,nu,i = 1,…,nx，P′ki、Pki分別為從uk到xi的所有路徑和其中一條路徑。

從狀態(tài)變量xl到輸出變量yj的路徑長度Llj(Plj)為：

相應(yīng)從狀態(tài)變量xl到輸出變量yj的最短路徑(xl,yj)為：

其中，l = 1,…,nx,j = 1,…,ny，、Plj表示從狀態(tài)變量xl到輸出變量yj的所有路徑和其中一條路徑。

從輸入變量uk到輸出變量yj的路徑長度Lkj(Pkj)為：

相應(yīng)從輸入變量uk到輸出變量yj的最短路徑(uk,yj)為：

其 中，k = 1,…,nu,j = 1,…,ny，P′kj、Pkj表 示 從 輸入變量uk到輸出變量yj的所有路徑和其中一條路徑。

從輸出變量yj到另一個輸出變量ym的路徑長度Ljm(Pjm)為：

相應(yīng)從輸出變量yj到另一個輸出變量ym的最短路徑(yj,ym)為：

其中，j,m = 1,…,ny，P′jm、Pjm表示從輸出變量yj到另一個輸出變量ym的所有路徑和其中一條路徑，節(jié)點到自身的最短路徑為(yj,yj)= ∞。

2.3 構(gòu)建鄰接矩陣

那么，所構(gòu)建的加權(quán)鄰接矩陣Aw將被用于社區(qū)發(fā)現(xiàn)算法，以對系統(tǒng)(1)進行劃分。

那么，可以通過如下步驟構(gòu)建加權(quán)鄰接矩陣Aw。

2.4 尋找最大模塊度

（1）對每個頂點i重復(fù)以下步驟：

①對于節(jié)點i 的每個鄰居節(jié)點j，通過將節(jié)點i從當前群組移動到節(jié)點j 所在的群組，計算模塊度的變化值ΔQw；

②找到上一步中最大的ΔQw＞0，將頂點i放入到相應(yīng)的群組中，得到新的分組。

（2）令nc(k + 1)作為節(jié)點集結(jié)后的群組數(shù)量，如果nc(k + 1) ＜nc(k)，回到步驟（1）。否則，找到最大化的模塊度Q。

那么，最大化模塊度Q 對應(yīng)的分組即為本文的子系統(tǒng)分解結(jié)果，然后，對所得到的子系統(tǒng)進行能觀能控性判斷，如果滿足，則當前分解為最優(yōu)分解，如果不滿足，則找到第二大的模塊度Q 對應(yīng)的分解方案，以此類推。

圖2 兩個CSTR和一個分離釜的復(fù)雜反應(yīng)分離器系統(tǒng)Fig.2 Two CSTRs and a flash tank with recycle stream

3 仿真驗證

本研究考慮的過程是由兩個連續(xù)攪拌反應(yīng)釜(continuous stirred tank reactor, CSTR)和一個快速分離釜組成的三個容器的復(fù)雜反應(yīng)分離器系統(tǒng)（圖2）[30]。第一個CSTR 包含反應(yīng)物A，轉(zhuǎn)化為所需的產(chǎn)物B。物種A也可以反應(yīng)成副產(chǎn)品C。溶劑不反應(yīng)，標記為D。第一個CSTR 的出水和額外的新鮮飼料構(gòu)成了第二個CSTR 的入口。反應(yīng)A→B 和A→C（分別稱為1 和2）在兩個CSTR 中串聯(lián)進行。假設(shè)三個容器都有靜態(tài)滯留。

在標準建模假設(shè)下，通過材料和能量平衡獲得的描述系統(tǒng)行為的動力學方程如下：

式(29)中各參數(shù)值見表1。每個儲罐都有一個外部熱量輸入/排出。所需的工作點xs是過程的不穩(wěn)定穩(wěn)態(tài)點，為xs=[369.5,3.318,0.172,0.042,435.3,2.751,0.446,0.111,435.3,2.882,0.497,0.120]T。

閃蒸罐分離器模型的運行假設(shè)是：在閃蒸罐的工作溫度范圍內(nèi)，每種物質(zhì)的相對揮發(fā)度保持恒定。該假設(shè)允許根據(jù)容器液體部分的質(zhì)量分數(shù)計算塔頂?shù)馁|(zhì)量分數(shù)。還假設(shè)分離器中發(fā)生的反應(yīng)量可以忽略不計。對塔頂流的組成與閃蒸罐中持液率的組成進行建模：

式中，xD是閃蒸罐液含率中溶劑的質(zhì)量分數(shù)。

對該化工過程實施本文提出的分解方法，并和文獻[11,30]中的方法進行比較，得到了三個分解：文獻[30]中根據(jù)物質(zhì)種類劃分方法的分解1；根據(jù)文獻[11]的結(jié)構(gòu)劃分建議得到的分解2；應(yīng)用本方法的分解3。這三個分解如表2 所示。實現(xiàn)本文子系統(tǒng)分解的計算時間0.068 s。

表1 式（29）、式（30）中的參數(shù)值Table 1 Parameter values in Eq.（29）and Eq.（30）

表2 系統(tǒng)（29）的三種分解Table 2 Three decompositions for system（29）

分 別 對 表2 中 的 分 解1、2、3 進 行 分 布 式Lyapunov MPC 控制實現(xiàn)。在不同的Lyapunov MPC實現(xiàn)下，每次運行將使用相同的初始條件和干擾序列。初始系統(tǒng)狀態(tài)設(shè)定為：x0=[356.9, 3.227, 0.030,0.015, 446.5, 2.696, 0.465, 0.125, 446.5, 2.821,0.530,0.139]T。

定義性能指標如下：

首先，比較了無干擾的情況，兩個分解模型在一次評估中，分布式Lyapunov MPC 迭代次數(shù)為2 次實現(xiàn)時，標稱閉環(huán)系統(tǒng)的溫度和輸入軌跡如圖3 所示，性能指標和平均CPU 時間如表3 所示。圖3 和表3 的結(jié)果表明，考慮耦合強度的分解3 在一定程度上優(yōu)化了性能，且沒有過多消耗時間，在同樣的控制策略下，本文所提出的分解方法子系統(tǒng)內(nèi)部耦合更為緊密，能夠提升控制效果。

表3 不同分解模型的迭代分布式MPC和集中式MPC的性能指標（J）和CPU時間Table 3 Performance index（J）and CPU time of iterative distributed MPC with different decomposition models and centralized MPC

然后，為測試方法的抗干擾性，對三個子系統(tǒng)模型加入相同干擾。對每個子系統(tǒng)模型進行五次具有不同干擾序列的運行。性能指標和總CPU 時間如表4 所示，結(jié)果說明在相同的分布式控制器設(shè)計條件下，本文分解方法比分解1 和分解2 方法在性能上取得了更好的效果，說明考慮子系統(tǒng)之間的權(quán)重考慮，能夠獲得內(nèi)部連接更強的子系統(tǒng)劃分，提升分布式模型預(yù)測控制的性能。

圖3 分解1（藍色點劃線）、分解2（紅色實線）和分解3（黑色虛線）的閉環(huán)系統(tǒng)的溫度和輸入軌跡Fig.3 Temperature and input trajectories of the nominal closed-loop system under Decomposition 1(blue dash dotted line),Decomposition 2(red solid line),and Decomposition 3(black dotted line)

最后，針對分解3 進行不同迭代次數(shù)的分布式Lyapunov MPC 實現(xiàn)的比較，其參數(shù)值與第一組仿真相同，分別考慮了迭代次數(shù)為1 次、2 次和3 次的情況。閉環(huán)系統(tǒng)的溫度和輸入軌跡如圖4 所示，消耗的性能指標和平均CPU 時間如表5所示，結(jié)果表明，隨迭代次數(shù)增加，雖然計算時間增長，但控制性能提升，在迭代2次時，性能收斂。

表4 不同分解模型的分布式MPC和集中式MPC在5次運行中的性能指標（J）和CPU時間Table 4 Performance index（J）and CPU time of distributed MPC with different decomposition models and centralized MPC in five runs

圖4 分布式LMPC在迭代1次（紅色實線），2次（黑色虛線），3次（藍色點劃線）時的溫度和輸入軌跡Fig.4 The temperature and input trajectories of distributed LMPC when iteration times is 1(red solid line),2(black dotted line)and 3(blue dash dotted line)

表5 不同迭代次數(shù)的分布式LMPC的性能指標和CPU時間Table 5 Performance index and CPU time of distributed LMPC with different iteration times

4 結(jié) 論

提出了一種復(fù)雜化工過程的子系統(tǒng)分解，用于考慮耦合強度的分布式狀態(tài)估計和控制，在連續(xù)攪拌反應(yīng)釜過程的實例仿真表明，在不改變控制策略的前提下，本文提出的分解方法提升了過程控制性能。相對比只考慮系統(tǒng)連接度的子系統(tǒng)劃分方法，本文分解方法得到的子系統(tǒng)內(nèi)部之間連接度更強，而子系統(tǒng)之間的連接弱，因而更適合用于分布式控制設(shè)計。未來可以考慮改進連接邊權(quán)重的定義，使之能夠更好地反映系統(tǒng)的連接強度，應(yīng)用于更為復(fù)雜的系統(tǒng)分解，并通過設(shè)計分布式狀態(tài)估計進行綜合驗證。

符 號 說 明

CA1,CA2,CA3——分別為A在容器1、2、3中的濃度，kmol/m3

CB1,CB2,CB3——分別為B在容器1、2、3中的濃度，kmol/m3

CC1,CC2,CC3——分別為C在容器1、2、3中的濃度，kmol/m3

CAr,CBr,CCr——分別為A、B、C在循環(huán)中的濃度，kmol/m3

CA10,CA20——分別為容器1、2進料流中A的濃度，kmol/m3

CP——比熱容，kJ/（kg·K）

E1,E2——分別為反應(yīng)1、2的活化能，kJ/kmol Fr——循環(huán)流量，m3/h

F1,F2,F3——分別為容器1、2、3的出水流量，m3/h

F10,F20——分別為進料流至容器1、2的流量，m3/h

Hvap——汽化熱，kJ/kmol

ΔH1,ΔH2——分別為反應(yīng)1、2的反應(yīng)熱，kJ/kmol

k1,k2——分別為反應(yīng)1、2的指前因子，h-1

MWA,MWB,MWC——分別為A、B、C的摩爾質(zhì)量，kg/kmol

Q1,Q2,Q3——分別為容器1、2、3的熱量輸入/排出，kJ/h

R——氣體常數(shù)，kJ/(kmol·K)

T1,T2,T3——分別為容器1、2、3中的溫度，K

T10,T20——分別為到容器1、2的進料流溫度，K

V1,V2,V3——分別為容器1、2、3的體積，m3

αA,αB,αC,αD——分別為A、B、C、D的相對揮發(fā)度