復雜非線性系統的子系統分解方法

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

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

引 言

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

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

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

因此，本文同時考慮狀態、輸入和輸出之間的關聯性，研究面向化工過程分布式狀態估計和控制的子系統分解方法。具體地，首先基于化工過程形成有向網絡，在構建的網絡中，狀態變量、輸入變量和測量輸出變量都被視為節點，提出基于社區結構發現算法的子系統分解方法，并在化工過程實例的分解中進行驗證，設計分布式MPC 算法實現來說明方法優勢。

1 預備工作

1.1 問題描述

針對復雜非線性系統的子系統問題，通常考慮如下動態模型：

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

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

1.2 有向圖

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

根據如下規則定義有向圖：

1.3 能觀性與能控性

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

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

1.4 社區發現算法

社區發現算法用于將復雜網絡劃分成子網絡，其構建方法可參考文獻[27-28]，假設一個網絡具有N個節點，其鄰接矩陣A 的各個元素Aij可以通過如下方式構建：

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

1.5 分布式Lyapunov MPC算法

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

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

（2）迭代控制器評估。

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

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

2 子系統分解方法

本文復雜非線性系統的子系統劃分主要流程如圖1 所示。首先，將非線性系統(1)用加權有向圖進行描述；然后，構建系統的加權有向圖，實施社區發現算法找到最大模塊Q，得到推薦的子系統劃分結果；最后，對所得到的子系統進行能觀能控性判斷。

圖1 復雜非線性系統的子系統的劃分流程Fig.1 Flowchart of the proposed subsystem decomposition method of nonlinear systems.

2.1 加權連接邊的定義

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

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

進而，構建敏感度矩陣：

本節中，定義如下連接邊的權重。

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

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

狀態變量xl到輸出變量yj連接邊的權重定義為：

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

2.2 最短路徑

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

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

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

相應從輸入變量uk到狀態變量xi的最短路徑(uk,xi)為：

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

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

相應從狀態變量xl到輸出變量yj的最短路徑(xl,yj)為：

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

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

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

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

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

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

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

2.3 構建鄰接矩陣

那么，所構建的加權鄰接矩陣Aw將被用于社區發現算法，以對系統(1)進行劃分。

那么，可以通過如下步驟構建加權鄰接矩陣Aw。

2.4 尋找最大模塊度

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

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

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

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

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

圖2 兩個CSTR和一個分離釜的復雜反應分離器系統Fig.2 Two CSTRs and a flash tank with recycle stream

3 仿真驗證

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

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

式(29)中各參數值見表1。每個儲罐都有一個外部熱量輸入/排出。所需的工作點xs是過程的不穩定穩態點，為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。

閃蒸罐分離器模型的運行假設是：在閃蒸罐的工作溫度范圍內，每種物質的相對揮發度保持恒定。該假設允許根據容器液體部分的質量分數計算塔頂的質量分數。還假設分離器中發生的反應量可以忽略不計。對塔頂流的組成與閃蒸罐中持液率的組成進行建模：

式中，xD是閃蒸罐液含率中溶劑的質量分數。

對該化工過程實施本文提出的分解方法，并和文獻[11,30]中的方法進行比較，得到了三個分解：文獻[30]中根據物質種類劃分方法的分解1；根據文獻[11]的結構劃分建議得到的分解2；應用本方法的分解3。這三個分解如表2 所示。實現本文子系統分解的計算時間0.068 s。

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

表2 系統（29）的三種分解Table 2 Three decompositions for system（29）

分 別 對 表2 中 的 分 解1、2、3 進 行 分 布 式Lyapunov MPC 控制實現。在不同的Lyapunov MPC實現下，每次運行將使用相同的初始條件和干擾序列。初始系統狀態設定為：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 迭代次數為2 次實現時，標稱閉環系統的溫度和輸入軌跡如圖3 所示，性能指標和平均CPU 時間如表3 所示。圖3 和表3 的結果表明，考慮耦合強度的分解3 在一定程度上優化了性能，且沒有過多消耗時間，在同樣的控制策略下，本文所提出的分解方法子系統內部耦合更為緊密，能夠提升控制效果。

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

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

圖3 分解1（藍色點劃線）、分解2（紅色實線）和分解3（黑色虛線）的閉環系統的溫度和輸入軌跡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 進行不同迭代次數的分布式Lyapunov MPC 實現的比較，其參數值與第一組仿真相同，分別考慮了迭代次數為1 次、2 次和3 次的情況。閉環系統的溫度和輸入軌跡如圖4 所示，消耗的性能指標和平均CPU 時間如表5所示，結果表明，隨迭代次數增加，雖然計算時間增長，但控制性能提升，在迭代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 不同迭代次數的分布式LMPC的性能指標和CPU時間Table 5 Performance index and CPU time of distributed LMPC with different iteration times

4 結 論

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

符 號 說 明

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在循環中的濃度，kmol/m3

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

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

E1,E2——分別為反應1、2的活化能，kJ/kmol Fr——循環流量，m3/h

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

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

Hvap——汽化熱，kJ/kmol

ΔH1,ΔH2——分別為反應1、2的反應熱，kJ/kmol

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

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

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

R——氣體常數，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的相對揮發度