基于改進CVA算法的系統辨識方法

臧春華 高偉

作者簡介：臧春華（1963-），教授，從事系統建模與優化控制的研究。

通訊作者：高偉（1998-），碩士研究生，從事系統辨識的研究，243679720@qq.com。

引用本文：臧春華，高偉.基于改進CVA算法的系統辨識方法的研究[J].化工自動化及儀表，2023，50（6）：000-000.

DOI：10.20030/j.cnki.1000-3932.202306000

摘? 要? 為了進一步提高規范變量分析（CVA）辨識方法的準確性，降低在干擾因素下辨識精度的損失，提出改進CVA辨識方法。控制策略上比較兩種經典的估計CVA型矩陣A和C方法，并分析這兩種方法所獲得的估計結果之間的差異。在此基礎上，為消除誤差，提出統一的估計系統矩陣A和C的方法，通過修正表達式的形式，將兩種經典方法聯系起來。仿真結果表明：辨識結果在保證精度的同時，得到的傳遞函數可以很好地擬合跟蹤工業數據，而且辨識方法具有較好的魯棒性。

關鍵詞? 改進CVA算法? 系統辨識? 狀態空間模型? 子空間辨識方法? 擬合跟蹤? 工業數據? 一階誤差

中圖分類號? TP18? ? ? 文獻標志碼? A? ? ? 文章編號? 1000-3932（2023）06-0000-00

子空間狀態空間基于加權矩陣的奇異值分解，在對多輸入多輸出系統進行簡單參數化的基礎上，能直接利用輸入輸出數據估計狀態序列，不受przori參數化和非線性優化引起的問題的影響，能廣泛應用于大規模、多工業過程領域，擁有易辨識、普遍適應工況的特點[1～5]。隨著現代控制技術的突飛猛進，控制對象的復雜性持續提高，多變量工業過程愈發常見，子空間辨識方法就被逐漸運用到越來越多的領域中[6～8]。規范變量分析（Canonical Variate Analysis，CVA）辨識方法是一類典型的子空間辨識方法，在子空間辨識方法的所有種類中，CVA辨識方法是在最佳統計推理原則的基礎上開發的，所以該方法具有高精度、易辨識的特點，一直是控制領域一種重要的辨識方法。實際生產中，CVA與N4SID（Numerical Algorithms for Subspace State Space System Identification）[9]和PO-MOESP[10]都作為子空間辨識方法的分支而被高度應用，國內外眾多學者也都提出了關于采用此類方法進行辨識的研究[11，12]。CVA構建目標函數時能夠在預測未來數據空間的基礎上，自高維系統數據中分離出低維狀態變量，從而估計最優模型階次，全過程基于Akaike信息判據的基礎實現[13]。在子空間辨識方法的所有種類中，由于CVA辨識方法是在最佳統計推理原則[14]的基礎上開發的，故相比之下它能始終提供更好的估計參數，成功率要優于其他子空間辨識方法。相較于傳統的使用回歸識別模型的方法，在實際工況中當延時環節已知或者能夠準確預測時，子空間辨識方法的準確性更好。

在實際工況中，噪聲的影響是不能忽視的，噪聲是伴隨工業生產過程而存在的。CVA方法能夠有效辨識過程數據整體的特征，同時有利于削弱過程數據中噪聲的影響。文獻[15]提出一種基于正交分解的新型遞推子空間辨識方法，在處理數據變量時合理利用了投影理論，并且結合遺忘因子的引入，在優化傳統方法的計算問題的同時解決了閉環狀態下的噪聲影響，能一致估計系統矩陣的特征值。文獻[16]提出一種閉環二維可分性系統子空間辨識方法SIMPCA-AR，通過代數矩陣重排方法預測馬爾可夫參數，該方法能在噪聲影響下和缺乏控制器先驗信息的情況下，直接預測狀態空間模型，還能同時利用輔助變量修正偏差問題。文獻[17]針對在不同工況時域下子空間辨識過程復雜的問題，提出時域非線性子空間重構方法，基于測量到的非線性位置響應數據訓練神經網絡模型，最后通過載荷識別提升了辨識速度。上述辨識方法能夠有效提升子空間辨識速度，減弱實際工況噪聲的影響。然而，辨識準確度仍然不夠優秀。文獻[18]提出基于機器學習的隨機子空間辨識的方法，借助外嵌的數字特征模態辨識神經網絡的深層結構，并獨立地通過在線過程獲取模態參數，在求解矩陣的奇異值的過程中能避免人工干預，提高了辨識的穩定性。文獻[19]提出基于高光譜圖像子空間辨識的FPGA實現方法，以計算噪聲估計矩陣，通過數值穩定的正交變換，避免了求解最小二乘方法的數值困難，同樣能夠保證辨識過程穩定運行。值得注意的是，上述方法基本能夠保證辨識過程的快速性或者穩定性，然而，由于階次過高，有時會使辨識結果仍然在精準度上有欠缺。筆者就此問題對比CVA方法中估計矩陣A和C的傳統方法，并分析傳統方法所獲得的兩種估計值之間的差異。為減小這種差異，在整合并簡化傳統方法中的表達式后，得到了近似于一階的誤差表達式，構建出統一的改進CVA算法。

1? CVA型矩陣A、C的經典估計方法

直接利用輸入輸出數據估計得到系統狀態空間模型的方法被稱為子空間辨識方法，對于多輸入多輸出（Multiple-Input Multiple-Output）系統，子空間辨識方法的簡單參數化是它能夠作為輸入輸出方法的有效替代方法的主要原因。在描述子空間辨識方法前，首先利用投影定理[20]，使用矩陣進行以下定義。

將一個矩陣的行空間投射到矩陣B的行空間上的算符為：

（1）

其中，

A的行空間在B的行空間上的投影為：

（2）

將A的行空間投影到B的行空間的正交補上的算符是：

（3）

定義（為Hilbert 空間），A的行空間沿B的行空間在C的行空間上的斜投影定義為：

（4）

將矩陣A和矩陣B之間的協方差定義為：

其中，為期望值運算符。

由于測量值有限，將協方差近似為：

其中，是測量的數量。

子空間辨識方法能夠從測量的輸入/輸出數據序列（和）中估計出下述系統[21]：

（5）

其中，假定輸入數據序列是一個持續激勵的準平穩確定性序列；；序列是一個具有的白噪聲隨機過程。

考慮以下輸入輸出數據：

其中，為采樣次數。

要想從收集到的數據中確定系統矩陣（A、B、C、D），就要定義：

輸入數據可以被寫在以下block Hankel矩陣[21，22]中：

（6）

其中，前兩個block Hankel矩陣的下標表示歷史數據，后兩個block Hankel矩陣的下標表示未來數據。

同樣，輸出數據可以按照上述方法構建成block Hankel矩陣，在中表示為，。

同時，定義輸入和輸出矩陣、。定義正交投影矩陣、，其中，表示矩陣P的行空間到M的行空間的正交投影，矩陣M的上標的加號符號表示Moore-Penrose偽逆[23]。

當輸入輸出數據被組織成具有這個特殊的block Hankel結構的矩陣時，式（5）可以寫為：

（7）

（8）

注意，是擴展的可觀矩陣。

同時，有：

（9）

其中，和為矩陣在運算的過程中產生的常數。

矩陣和是由一個卡爾曼濾波器組生成的狀態序列。

并行處理歷史輸入和輸出的block Hankel矩陣的每列：

（10）

（11）

是一個block Toeplitz[22，23]矩陣，如下所示：

（12）

在這個過程中，輸入持續地被激勵為f+p階，而且它與過程噪聲工作和測量噪聲不相關，測量次數也是無窮大的。過程噪聲和測量噪聲不會同時為0，加權矩陣W1和W2使得W1不是奇異的。定義作為斜投影：

（13）

的奇異值分解：

（14）

其中，矩陣等價于；系統的階數n等于方程（14）中奇異值的個數；擴展的可觀矩陣[1，22]等價于；狀態序列為。

辨識時，在確定完模型的階數后，通過加權投影的奇異值分解估計可觀矩陣。其中，，再聯系到，可知與二者可以表示零空間的正交投影。假定是可觀的，那么由以及投影定理可知，并且能夠得到以下關系：

（15）

其中，是奇異值為的對角矩陣；是正交矩陣。

系統的階數n由S的非零奇異值的個數決定。然而，在許多實際工況下，噪聲是無法避免的，而噪聲又會影響測量值，進而無法直接區分非零奇異值和零奇異值。此時，比較這些值或假設不同的順序以及比較模擬誤差都是必要的。和的列空間是相同的，所以，從中刪除最后1行后，可計算出=。然后可以得到系統矩陣，以求解以下最小二乘問題：

（16）

其中，表示矩陣的Frobenius范數[24]。

值得注意的是，卡爾曼濾波器狀態序列可以在不知道系統矩陣的情況下，能直接從擴展可觀矩陣的估計和輸入輸出數據確定，如式（1）～（4）。

綜上，提出了經典的估計A、C矩陣的方法。

2.1? 移位不變性法

在移位不變性法中。可以通過求解得到。其中，結合到的移位不變性，有，，因此，可得：

（17）

2.2? 最小二乘法

最小二乘法采用卡爾曼濾波狀態序列和求解最小二乘問題。投影定理指出，和與正交投影和有關，由投影定理公式可以將卡爾曼濾波器狀態序列表示為

，將（9）重寫為最小二乘問題（18），并得到矩陣C、A和的估計值：

+? ? ? ? ? ? ? ? ? ? ?（18）

3? 改進的CVA型矩陣A、C的統一估計方法

如果給定奇異值分解（式（15）），那么可得：

（19）

（20）

（21）

如果不是一個標準正交矩陣，那么：

（22）

再考慮方程：

（23）

方程（9）的解和式（23）的解是相同的。考慮到矩陣（）和（）的行空間是正交的，那么給定，可得：

）? ? ? ? ? （24）

給定，可知的行空間和））正交。在此基礎上，給定噪聲矩陣則Moore-Penrose偽逆可以近似為一階：

（25）

接下來，通過保持二階項，有：

（26）

對于CVA算法來說，如果對給定式（15）中描述的SVD進行分解，由最小二乘問題給出的解和由移位不變性方法得到的解會略有不同，情況如下：

（27）

（28）

（29）

（30）

證明

首先，方程（9）的解可以寫成：

（24）

然后，如果是標準正交的，那么，根據式（20）～（22），有，同時，由的行空間和（））正交可知：

（25）

（26）

（27）

（28）

在CVA算法中，有以及。再聯系式（25）、（26），可得：

（29）

因此有：

（30）

（31）

式（31）能夠在兩種經典方法存在估計差值的基礎上，更準確地估計出A、C的近似值，即式（31）為統一兩種經典方法提供了公式基礎。事實上，改進后的CVA方法就是要取代計算的偽逆，轉變為計算，這個矩陣的行比要少。因此，新的算法在數值上比前一代更高效。

4? 仿真實例

為了驗證筆者所述方法的性能，本研究展示的仿真案例有著動態過程，并且質量相關，系統階數為4。輸出的開關概率為0.2，范圍為[-1，1]，白噪聲以單位矩陣為協方差。在實驗中，收集了700個數據點，500個用于辨識步驟、200個用于驗證。在辨識出模型之后，使用相同的輸入序列計算所辨識系統的輸出，并與原系統以及以經典CVA方法辨識出的結果進行比較。歷史和未來的上限是f= 20，p=20。

使用經典位移不變性方法得到的估計值如下：

使用經典最小二乘方法得到的估計值如下：

使用改進統一估計方法得到的估計值如下：

同時，能夠得到：

6.3563

6.3531

6.3553

其中，trace為矩陣的跡，表明改進的方法更好地縮小了傳統的兩種方法預估C、A矩陣的差異。

值得注意的是，改進方法在預估矩陣A時，雖然已經比傳統方法效果更好，但橫向對比發現，它仍然比自己預估矩陣C時產生的誤差略大，這是因為在式（25）、（26）中，某些二階項被忽略了。辨識系統的輸入和輸出信號如圖1、2所示。

為了測試模型的質量，現進行驗證。驗證過程分別在原始系統中、用經典CVA方法辨識得到的系統中以及用改進后的CVA方法辨識得到的系統中輸入相同的序列，觀察并比較其各自的輸出情況。驗證輸入序列如圖3所示，輸出序列如圖4所示。

辨識性能通過使用解釋方差（EV）評估：

（27）

其中，代表輸出的平均值；代表第i個輸出的平均值；k代表第k個采樣時間；是估計殘差。

當一個模型返回一個輸出的時，解釋方差被近似為零。通常，若某方法的解釋方差越接近1，則代表該方法擬合系統的效果越好。考慮4種不同級別的白噪聲。計算了一些解釋方差并取平均值，辨識和驗證數據集的平均解釋方差EV的結果見表1。

可以看出，改進后的CVA算法在辨識階段具有較好的性能，而且抗噪聲能力更強。同時，從圖4上也能看出，新方法能夠更好地擬合原系統的曲線，展現出了很好的性能和更高的準確性。

在此基礎上，更換一組輸入序列，只使用改進的CVA算法進行辨識，以便更加直觀地觀察它的性能，原始信號和所辨識系統的輸出如圖5、6所示。

可見，改進后的CVA方法能夠很好地擬合系統，反映系統的性能。

5? 結束語

筆者提出一種統一的估計系統矩陣A和C的算法，將傳統的兩種估計A、C矩陣的方法最小二乘法與移位不變性法相聯系，得到了近似于一階的誤差表達式，減小了傳統的兩種估計A、C矩陣方法所估算結果的誤差。在此基礎上，進行仿真實驗。通過圖像以及平均解釋方差可以看出，改進后的CVA方法能夠更好地展現原系統的性能，并且在抗噪聲污染方面表現出了更好的性能，多次驗證過程也驗證了這個結論。證實筆者提出的方法，在保證精度的同時，得到的傳遞函數可以很好地擬合跟蹤工業數據，辨識方法有著較好的魯棒性。

參? 考? 文? 獻

[1] 尹宏鵬，陳波，柴毅，等.基于視覺的目標檢測與跟蹤綜述[J].自動化學報，2016，42（10）：1466-1489.

[2] YANG C，ZHANG T，ZHANG Z L，et al.MLD-MPC for Ultra-Supercritical Circulating Fluidized Bed Boiler Unit Using Subspace Identification[J].Energies，2022，15（15）：5476-5499.

[3] NI Z，LIU J，WU Z，et al.Identification of the state-space model and payload mass parameter of a flexible space manipulator using a recursive subspace tracking method[J].Chinese Journal of Aeronautics，2019，32（2）：513-530.

[4] 朱銳.基于子空間的非線性系統辨識研究[D].南京：東南大學.

[5] SZYMON G，RICCARDO R，S?LEYMAN CEM YENICELI，et al.Uncertainty quantification of modal parameter estimates obtained from subspace identification：An experimental validation on a laboratory test of a large-scale wind turbine blade[J].Engineering Structures，2022，256：114001.1-114001.13.

[6] 李明珠，占永寧，張守興.子空間辨識方法在污水處理中的應用[J].時代農機，2019，46（10）：139-140.

[7] CHENGSHUO S，ZHIDONG Q，HAO Q，et al.PEMFC Fractional-order Subspace Identification Model[J].China Petroleum Processing & Petrochemical Technology，2022，24（3）：151-160.

[8] FAHIM A H，JAKIR H S，SUKUMAR K.A Novel Hybrid Deterministic-Stochastic Recursive Subspace Identification for Electromechanical Mode Estimation，Classification，and Control[J].IEEE Transactions on Industry Applications，2021，57（5）：5476-5787.

[9] 李正明，郭世偉.基于N4SID子空間辨識的線性自抗擾控制器[J].信息與控制，2017，46（2）：159-164.

[10] CABRERA R N，VM ALVARADO V M，GL LOPEZ G L，et al.An Alternative Approach to the Inference of the Extended Observability Matrix，and Its Relation With the PO-MOESP Algorithm[J].IEEE Transactions on Control Systems Technology，2017，25（3）：888-898.

[11] JEMMALI M A，OTIS M J-D，ELLOUZE M.Robust stabilization for discrete-time Takagi-Sugeno fuzzy system based on N4SID models[J].Engineering Computations，2019，36（4）：1400-1427.

[12] KENJI IKEDA.Frequency Domain Uncertainty in PO-MOESP Method[C]//Proceedings of the ISCIE International Symposium on Stochastic Systems Theory and Its Applications，2016.DOI：10.5687/sss.2016.104.

[13] KHALID A，SARWAT A I.Unified Univariate-Neural Network Models for Lithium-Ion Battery State-of-Charge Forecasting Using Minimized Akaike Information Criterion Algorithm[J].IEEE Access，2021（9）：39154-39170.

[14] HASHEMPOUR M，DOOSTPARAST M，PAKDAMAN Z.Statistical Inference on the Basis of Sequential Order Statistics under a Linear Trend for Conditional Proportional Hazard Rates[J].Statistics Optimization &；Information Computing，2020，8（2）：462-470.

[15] 周楚程.基于子空間辨識的工業過程建模及故障檢測[D].沈陽：東北大學，2017.

[16] 索寒生，程金旭.二維可分性系統遞推子空間辨識方法研究[J].控制工程，2018，25（9）：1713-1720.

[17] 龐志雅，馬志賽，丁千.基于神經網絡和子空間的非線性系統載荷識別[J].振動.測試與診斷，2022，42（5）：877-885；1031.

[18] LI Y C，LIN F K，LUO X G，et al.Application of an ensemble learning model based on random subspace and a J48 decision tree for landslide susceptibility mapping：a case study for Qingchuan，Sichuan，China[J].Environmental Earth Sciences，2022，81（9）：267-299.

[19] LE?N G，GONZ?LEZ C，MAYO R，et al.Noise estimation for hyperspectral subspace identification on FPGAs[J].The Journal of Supercomputing，2019，75（3）：1323-1335.

[20] 黃金杰，宮煜晴，郝現志.平穩切換LPV系統的H∞控制器設計[J].控制與決策，2022，37（5）：1167-1173.

[21] GUPTA S，NARANG R，HARISH M，et al.MDS block hankel-like rhotrices using conjugate elements and self-dual bases of finite fields[J]. Bulletin of Pure & Applied Sciences-Mathematics and Statistics，2022，41e（2）：184-198.

[22] ZHANG Y Y，CAO G F，LI H.Products of Toeplitz and Hankel Operators on Fock-Sobolev Spaces[J].Chinese Annals of Mathematics，Series B，2022，43（3）：401-416.

[23] HON S.Optimal block circulant preconditioners for block Toeplitz systems with application to evolutionary PDEs[J].Journal of Computational and Applied Mathematics，2022，407.

[24] Guo P C.A Frobenius Norm Regularization Method for Convolutional Kernel Tensors in Neural Networks[J].

Computational Intelligence and Neuroscience，2022，2022（48）：3277730.

（收稿日期：2023-01-11，修回日期：2023-04-04）