基于二次型規劃考慮網絡丟包的魯棒狀態估計

王中杰，易總根

（同濟大學電子與信息工程學院，上海201804）

隨著通信網絡的發展，基于網絡的反饋控制系統（NCS）受到越來越多的關注［1－2］.在NCS中，從傳感器節點到控制節點和從控制節點到執行器節點之間的數據傳輸與通信都必須通過網絡進行.而與此同時，由于微機電系統（MEMS）技術的發展，無線傳感網絡（WSN）的應用也越來越廣泛［3］.在WSN中，傳感器節點感知并且收集數據然后發送給中心節點進行處理從而實現對物理系統的感知.由于網絡的廣泛應用，對這種類型的系統進行物理狀態的估計與傳統方法有許多不同之處.

由于網絡環境的復雜性，比如網絡帶寬的有限性和網絡的共享性，使得網絡通信中總會存在不可靠性和不確定性.通常這些不確定性包括隨機的延時、數據包丟失和測量失效.本文只考慮網絡數據包丟包情形下的狀態估計問題.對于該問題，文獻［4］提出了一個用于間斷觀測的Kalman濾波器，對估計誤差的統計收斂屬性進行了研究，并且證明了在一定觀測到達率下的存在性，而超過一定的閾值后誤差將會發散.文獻［5］從部分觀測損失的方面提出了一種Kalman濾波，并且導出了適用于該情況下的濾波更新過程.文獻［6］對數據包丟失進行了較為詳細的描述與定義，并給出了一個考慮數據丟包的數學模型，基于這個模型設計了最優Η2濾波器.文獻［7］將隨機丟包率轉換為系統中的一個隨機參數，推導出了一種考慮隨機參數和隨機輸入的通用Η∞范數濾波方法.在噪聲是獨立的時不變情況下，文獻［8］設計了最優的全階和降階的狀態估計器對具有數據包丟失的離散時間系統進行狀態估計.文獻［9］考慮了多數據包丟失的情況，并且在離散時間隨機線性系統的框架下將這種情況轉換為一個測量延時和滑動平均測量噪聲，并且基于此設計了無偏的最優濾波器.文獻［10］從概率角度考慮具有網絡丟包的離散時間系統的狀態估計問題.文獻［11］從Markov鏈的角度提出了一個模型，該模型綜合考慮了隨機延時、數據包丟失和測量失效，并且基于此設計得到最優Kalman濾波器.文獻［12］通過概率權值的方法得到了一個適用于具有多數據包丟失特性的非線性系統的非線性濾波器.文獻［13］中將隨機多數據包丟失假定為服從Bernoulli分布，并且設計了一種基于零均值高斯白噪聲的Kalman濾波器.

盡管已經提出了很多關于網絡數據包丟失的處理算法，但是大多數算法都不具有魯棒性或者魯棒性還不夠，比如有些必須知道丟包率大?。?2］.當然也有一些具有魯棒特性較好的算法，比如文獻［6］和［7］中的算法，但是這兩種算法的求解過程必須計算相對比較復雜的線性矩陣不等式（LMI）問題，尤其是當狀態維數較大的時候，計算量會比較大，在實際應用中會受到一定的限制.即使是線性模型，當有數據包丟失時魯棒狀態估計會使得其變成一個非線性問題，也會造成算法的計算量上升.因此，為了在滿足一定魯棒性要求的前提下獲得更高的計算效率，可以將魯棒估計問題轉化為具有更高求解效率的二次型規劃問題（Quadratic Programming，QP）［14－16］.

本文給出了一種網絡數據包丟失情況下狀態估計的數學模型描述，目的是將該魯棒狀態估計問題轉化為向量優化問題，通過標量化方法將其轉化為一個二次型規劃問題.考慮到求解算法上的高效性，最后將該問題轉化為l1正則化最小平方問題（l1－regularized least squares problem）求解.進一步加入動態更新過程就形成本文提出的基于二次型規劃的魯棒狀態估計算法（Robust State Estimation Based on Quadratic Programming，RSE＿QP）.

1 問題描述

在狀態估計過程中，首先要假定物理過程數學模型是已知的.通過這個數學模型，能夠獲得狀態預測值.但是在實際的動態系統中，許多噪聲信號會使真實的狀態值與理論值不同.因此就需要一個對系統狀態進行動態估計的濾波器，并且必須保證對系統噪聲影響的魯棒性.而其中最基本的思想是通過設計一個濾波器使得狀態預測值與估計值之間誤差最小，這樣就能得到最優的狀態估計.圖1對狀態估計原理進行了簡單描述.

圖1 狀態估計原理圖Fig.1 State estimation process

1.1 模型

式中：A∈Rn×n；C∈Rm×n；x（k）∈Rn是狀態向量；y（k）∈Rm是量測輸出；｛w（k）｝和｛v（k）｝是穩態、零均值離散高斯白噪聲過程，w（k）～N（0，W）是輸入噪聲，v（k）～N（0，V）是量測輸出噪聲，其中W和V為對應的方差值.

當物理過程與濾波器之間存在網絡環境時，數據包的丟失必然存在.上面的模型是無法描述這樣的情況的.如果要將網絡丟包考慮到數學模型中，必須對上面的模型進行一定的修正，使用文獻［14］中類似的方法，可以得到下面的模型：

式中增加項u（k）∈Rm是為了描述網絡丟包的.在這個模型中，u（k）的值是可以調整的，這樣可以模擬網絡丟包現象.在時刻k，當沒有丟包時，u（k）＝0；當有丟包時，u（k）的值可以使得y（k）的某些分量的值為0，即存在丟包.由于一般情況下網絡丟包的概率相對比較小，所以可以假定向量u（k）是稀疏的.

1.2 問題轉化

存在這樣一個事實，當一個系統可觀測時，必須通過量測輸出y（k）來對系統的狀態進行估計.為了使狀態估計更加準確，除了要考慮狀態誤差，y（k）理論值與實測值之間的誤差應該盡可能小.因此，當要對一個動態系統的狀態進行準確的估計與跟蹤時，必須考慮輸出誤差的衡量和狀態誤差的衡量.定義輸出誤差的目標函數為J1，狀態誤差的目標函數為J2.

根據文獻［14］，輸出誤差的目標函數J1可以表示為

式中：h為一個調整參數.

由于v（k）＝y（k）－C′x′（k）

因此優化問題（6）可以被用來衡量輸出誤差的最優值.

在標準的Kalman濾波中，存在一個如下形式的狀態更新過程：

其中xk｜k－1是時刻k時的狀態預測值，而xk－1｜k－1是時刻k－1時的狀態估計值.為了得到時刻k時的狀態估計值，可以使用下面形式的目標函數來衡量：

其中Pk｜k－1是時刻k時的狀態估計誤差的先驗協方差矩陣.通過上面的目標函數，狀態的估計值可以通過下面的二次型優化問題求解得到.

通過上面的討論，如果必須同時考慮輸出誤差與狀態誤差兩部分，即必須要獲得這兩個目標函數的最小值，這個問題就成為一個多目標優化問題，也是一種向量優化的問題［17］.因此，考慮網絡數據包丟失的魯棒估計問題可以描述成一個向量優化問題.下面用相應的向量優化問題來近似描述該魯棒狀態估計問題，首先定義一個向量函數J，有

然后可以得到該魯棒狀態估計問題的向量優化表達形式

為了得到時刻k的狀態估計值，必須在已經得到的時刻k－1的狀態估計值xk－1｜k－1基礎上對上述向量優化問題進行求解，得到向量優化問題的Pareto最優值［17］.必須指出的是，通過該形式獲得的估計值無法從理論上證明其為最優估計值，只是一個近似最優估計值，但可以使得這個魯棒估計問題在求解效率上有很大的提高.

2 RSE＿QP算法

對于一個完整的狀態估計與跟蹤算法，必須包括當前時刻的狀態估計與狀態更新過程.其中當前時刻狀態估計就是要對上述向量優化問題模型進行求解，首先要將向量優化問題轉化為標量優化問題，然后對該標量優化問題求解最優值得到當前狀態估計值.

2.1 標量化

有很多種方法可以對向量優化問題進行求解.而在這些方法中，標量化（scalarization）方法較簡單也較實用.通過標量化方法可以將向量優化問題轉化為一個普通標量形式的優化問題，并且可以得到該問題的Pareto最優值.下面是轉化的過程.

式中：θ1和θ2可以被認為是在優化過程中目標函數J中兩部分的權重值.定義λ＝θ1h，可得

上述的向量優化問題（10）可以等價為下面的標量二次型優化問題.

2.2 更新過程

進行上述標量化過程以后，得到的優化問題（13）是一個典型的二次型規劃問題.而對該二次型優化問題求解時，只能得到時刻k的狀態估計值，相當于一個靜態估計過程.如果要對一個動態系統的狀態進行狀態跟蹤與估計，必須要有每個時刻的更新過程.本文借助標準的Kalman濾波更新過程來對每個時刻的狀態誤差協方差Pk｜k和增益矩陣K進行更新，更新過程的方程如下.

狀態預測方程為

狀態誤差協方差預測方程為

Kalman增益方程為

狀態更新方程為

狀態估計誤差后驗協方差矩陣更新方程為

2.3 優化問題求解

對于上述二次型優化問題，求解的方法可以有很多種.為了尋求一種效率更高，通用性更好的方法.本文通過一定的轉化，使得求解變得相對直接并且高效.在本文中，優化問題（13）將轉化為標準的l1正則化最小平方問題.

首先，定義

從公式（17）和（19）可以得到

根據式（3）

將公式（20）和（21）代入公式（12）中，目標函數可以做如下變換：

定義

可以證明Q是一個正定矩陣，下面是證明過程.

設存在任意的一個非零向量α，有

由于V和Pk｜k－1都是正定的，因此

這樣就得到Q是一個正定矩陣.對其進行Cholesky分解，可以表示成如下形式：

式中：M為上三角矩陣.通過公式（22）和（27），可以得到目標函數的描述.

通過上述變換，優化問題（13）已轉化為一個標準的l1正則化最小平方問題.

使用一種截短牛頓內點法（truncated Newton interior－point method，TNIPM）［18］進行求解.

求解過程中，參數λ的值是用來設定向量u（k）的稀疏程度的，其值由數據丟包的多少決定.當沒有數據丟包時，可以設定λ的值非常大；而當存在數據包丟失時，可以通過下面的方法進行自適應設定：

其中0＜k＜1，依據一定的經驗，k的值應該要遠小于1，在本文中設定k＝0.01.

2.4 算法流程

通過上述討論，當前時刻的狀態估計可以通過求解相應的l1正則化最小平方問題得到，然后通過Kalman動態更新過程動態進入下一步的狀態估計，以下是整個RSE＿QP算法運行流程.

Step 1 初始化x0，P0｜0＝P0.

Step 2 當k≥1，根據公式（14），（15）和（16），計算xk｜k－1，Pk｜k－1，P－1k｜k－1和K（k）的值.

Step 3 獲得時刻k的測量值y（k），然后利用公式（19）計算e（k）.

Step 4 利用l1正則化最小平方問題求解器求解優化問題（29）得到u（k）的值.

Step 5 根據公式（17）和（18）更新xk｜k和Pk｜k.

Step 6 讓k＝k＋1，跳轉到step 2.

3 仿真實驗

為了驗證本文所提出的對于網絡數據包丟失的算法的有效性，本文分別給出了利用二階和高階（以四階為例）離散線性時不變系統的模型對算法進行仿真實驗.

3.1 二階系統仿真

利用二階離散線性時不變系統的模型對算法進行仿真實驗，其數學模型形式如下［6］：

圖2與圖3分別是仿真過程中兩個狀態x1和x2的仿真結果.從圖中可以看到，本文算法（RSE＿QP）能夠很好地滿足估計和跟蹤實際狀態的任務，并且跟蹤效果比標準的Kalman濾波（KF）好，表明本文算法體現出一定的魯棒性，能夠較好地處理具有網絡數據包丟失的情況.

在圖4中，使用狀態誤差的均方根（RMS）來評價估計結果，在每一時刻k，對狀態所有分量進行計算得到其均方根.圖5是在整個仿真過程中網路丟包的分布圖.結合這兩個圖，可以看到當存在數據丟包時，本文算法估計過程中的誤差均方根要遠小于標準的Kalman濾波算法，反映了在具有丟包情況下本文算法的優越性.

通過改變仿真過程中網絡丟包率的大小進一步驗證本文算法的魯棒性.圖6顯示了仿真實驗結果，整個仿真過程中，設置了不同的丟包率，對于每種丟包率下都仿真100個時間間隔并且求取100個均方根的平均值得到仿真實驗結果圖.從圖6中可以看出，隨著丟包率的增加，平均均方根總體趨勢是增加的，但從增長的速度來看，本文算法明顯慢于標準的Kalman濾波算法，也反映本文算法是具有魯棒性的.

圖6 不同數據丟包率下的平均誤差均方根（二階）Fig.6 Average RMS under different network packet dropout（2nd order）

3.2 高階系統仿真

利用四階離散線性時不變系統的模型對算法進行仿真實驗，其數學模型形式如下：

與二階系統仿真類似，設定W＝diag（1，1，1，1），V＝diag（0.01，0.01，0.01，0.01），P0｜0＝P0＝diag（1，1，1，1），x0＝［1 1 1 1 ］T，參數θ1＝θ2＝1.其他仿真參數與二階系統一致，并且丟包分布圖也一致.同樣在整個仿真過程中，將本文提出的算法與標準的Kalman濾波算法進行比較分析.

圖7～10分別是仿真過程中4個狀態x1，x2，x3，x4的仿真結果.從圖中可以看到，本文算法針對高階系統同樣能夠較好地進行狀態的跟蹤，在丟包的情況下可以看出本文算法能很好地平滑跟蹤過程，較好地完成狀態估計任務.

同樣在圖11中，使用狀態誤差的均方根來評價估計結果，在每一時刻k，對狀態所有分量進行計算得到其均方根.同時，在圖12中，改變丟包率，驗證本文算法的魯棒性.從圖中也可以看出，對高階系統魯棒性也同樣得到了提高，丟包率變化的干擾很小.

4 結論

本文考慮了網絡數據包丟失情況下狀態估計問題，給出了一種新的包含數據包丟失描述的模型，并且將此魯棒狀態估計問題轉化為一個二次型規劃問題，最終經過一定的轉換利用標準的l1正則化最小平方問題求解器得到所要估計的狀態.通過仿真可以看到，本文算法具有較好的魯棒特性，能很好地進行動態系統的魯棒狀態估計，由于其屬于二次型規劃問題，容易實現，有較高的運算效率，因此具有較好的工程應用性.

［1］ Hespanha J P，Naghshtabrizi P，Xu Y.A survey of recent results in networked control systems［J］.IEEE Transactions on Automatic Control，2007，95（1）：138.

［2］ Zhang W，Branicky M S，Phillips S M.Stability of networked control systems［J］.IEEE Control System Magazine，2001，21：84.

［3］ Akyildiz I F，Su W，Sankarasubramaniam Y，et al.Wireless sensor networks：a survey［J］.Computer Networks，2002，38（1）：393.

［4］ Sinopoli B，Schenato L，Franceschetti M，et al.Kalman filtering with intermittent observations［J］.IEEE Transactions on Automatic Control，2004，49（9）：1453.

［5］ Liu X，Goldsmith A.Kalman filtering with partial observation losses［C］／／43rd IEEE Conference on Decision and Control.Atlantis：IEEE，2004：4180－4186.

［6］ Sahebsara M，Chen T，Shah S L.Optimal H2filtering with random sensor delay，multiple packet dropout and uncertain observations［J］.Int J Control，2007，80（2）：292.

［7］ Sahebsara M，Chen T，Shah S L.Optimal H∞filtering in networked control systems with multiple packet dropouts［J］.Systems ＆Control letters，2008，57（9）：696.

［8］ Sun S，Xie L，Xiao W.Optimal full－order and reduced－order estimators for discrete－time systems with multiple packet dropouts［J］.IEEE Transactions on Signal Processing，2008，56（8）：4031.

［9］ Sun S，Xie L，Xiao W.Optimal filtering for systems with multiple packet dropouts［J］.IEEE Transactions on Circuits and Systems—II：express briefs，2008，55（7）：695.

［10］ Shi L，Epstein M，Murray R M.Kalman filtering over a packet－dropping network：aprobabilistic perspective［J］.IEEE Transactions on Automatic Control，2010，55（3）：594.

［11］ Moayedi M，Soh Y C，Foo Y K.Optimal kalman filtering with random sensor delays，packet dropouts and missing measurements［C］／／American Control Conference.St.Louis：IEEE，2009：3405－3410.

［12］ Chen J，Li J，Ma L.Nonlinear filtering with multiple packet dropouts［C］／／International Conference on Image Analysis and Signal Processing.Hangzhou：IEEE，2010：83－88.

［13］ Xu Y，Wang W.Kalman filtering for systems with multiple packet dropouts［C］／／The 8th World Congress on Intelligent Control and Automation.Ji’nan：IEEE，2010：4996－5001.

［14］ Fuchs J J.An inverse problem approach to robust regression［C］／／IEEE International Conference on Acoustics，Speech，and Signal Processing.Phoenix：IEEE，1999：1809－1812.

［15］ Borguet S，Leonard O.A sensor－fault－tolerant diagnosis tool based on a quadratic programming approach［J］.Journal of Engineering for Gas Turbines and Power，2008，130（2）：021605.

［16］ Mattingley J，Boyd S.Real－time convex optimization in signal process［J］.IEEE Transactions on Signal Processing，2010，27（3）：50.

［17］ Boyd S，Vandenberghe L.Convex optimization［M］.Cambridge：Cambridge University Press，2004.

［18］ Kim S，Koh K，Lustig M，et al.An interior－point method for large－scale l1－regularized least squares problems［J］.IEEE Journal of Selected Topics in Signal Processing，2007，1（4）：606.