一、MPC簡介
1.1 預測控制的日常應用
模型預測控制的設計目标是計算未來控制變量u的軌迹，以優化未來的系統輸出y。優化過程在一個有限的時間視窗進行，并且利用優化時間視窗開始時的系統資訊進行優化。為了了解預測控制的基本思想，以一個日常工作為例來進行說明。
工作從上午9點開始，團隊的工作目标是完成一個液體容器模型預測控制系統的設計和實作任務。我們會計劃好接下來8個小時的工作内容，但是隻按照接下來一個小時的計劃内容執行，一個小時之後又開始計劃接下來8個小時的工作内容，重複這一過程，知道工作完成為止。
在9點這一時刻，計劃工作内容的時候，會将9點之前已經完成的工作内容考慮在内。假設工作任務可以分為模組化、設計、仿真和實作，那麼完成這些任務可以看成是一些列變量的函數，這些變量包括：團隊付出的努力程度、團隊的工作默契程度、是否有外援幫助等。這些都是規劃問題中的可操縱變量。同時團隊能力也有局限性，比如了解設計問題的能力、計算機軟硬體工程的能力等。這些都是規劃問題中的硬限制和軟限制。我們已經擷取的關于這項工作的背景資訊（已完成的工作情況）至關重要。
在考慮了所有事情之後，将接下來8小時的任務定義成可操縱變量的函數。然後，計算出為了完成任務，接下來8小時中每小時所需要的做的事情。在計算中，我們将基于背景資訊，考慮到我們自身的局限性作為限制條件，尋找實作目标的最佳途徑。得到從上午9點到下午5點的工作計劃後，開始執行第一個小時的計劃，即9點到10點的工作計劃。
在上午10點，需要檢查在這一個小時内完成了多少工作。這些資訊用于計劃我們下一階段的工作。也許會因為各種原因導緻這一小時内沒有完成預計的工作。然而，在10點鐘，我們對我們所取得的成績進行了評估，并利用這一最新資訊來規劃我們接下來8個小時的工作。最終工作目标可能保持不變，也可能會改變。做計劃的時間一直是8小時不會改變。重複9點的計劃過程，然後給出接下來8小時的新的計劃。同樣隻執行計劃中的第一個小時的内容。11點，我們再次評估我們所取得的成績，并将最新的資訊用于接下來8小時的工作計劃。計劃和實施過程每小時重複一次，直到達到最終的目标。
在做計劃的過程中，有三點非常關鍵。第一是要能預測可能發生的事情（也就是要有模型）；第二個是能夠評估目前的狀态和進展（也就是要有測量方法）；第三個是要能夠實施計劃的内容（也就是要實作控制）。規劃工作主要包括：
- 規劃的時間固定在8個小時；
- 在規劃之前需要知道目前的狀态；
- 通過考慮限制條件，對8小時的工作采取最佳規劃，利用最新的可用資訊，在移動的時間視窗進行實時優化。
以上描述的規劃活動與MPC的原理相同。在這個例子中，提出了一些機率會經常用到：
- 移動時間視窗：從任意時刻ti到ti+Tp的時間視窗。視窗的長度Tp是不變的。在本例中每個小時會對接下來8小時的工作進行規劃，是以時間視窗長度為Tp=8，ti是時間視窗的開始點，每小時進行遞增，規劃從上午9點開始的，是以開始時ti=9。
- 預測時域：規定了希望未來能被預測到多久。這個參數與時間視窗長度Tp相同。
- 滾動時域控制：雖然未來控制信号的最優軌迹完全在移動時間視窗内描述，但實際的控制輸入隻取控制信号的第一個樣本，而忽略了軌迹的其餘部分。
- 目前狀态：在規劃過程中，需要時刻的資訊來預測未來。這個資訊描述為x(ti)，這是包含了相關因素的向量，通常是測量或預測得到的。
- 模型：在預測控制中，描述系統動态的模型至關重要。一個好的動态模型能對未來做出準确一緻的預測。
- 代價函數：為了做出最優決策，需要一個标準來衡量達到期望目标的程度。用來反映期望目标與實際響應之間差别的函數叫做目标函數，通常也叫做代價函數J，在優化視窗内使代價函數最小的控制軌迹，就是最優規劃軌迹。
1.2 預測控制中的模型
模型的表示方式有幾種：FIR（有限沖擊響應）模型、階躍響應模型、傳遞函數模型、狀态空間模型等。不同的模型表示方式，對應的預測控制算法不一樣。基于FIR或階躍響應模型的控制算法包括動态矩陣控制（DMC）以及二次型DMC。但是這種方法僅限于穩定的控制對象，而且需要非常大的模型階數和非常多的模型系數。傳遞函數模型能更簡潔的描述過程動力學，适用于穩定和不穩定對象。基于傳遞函數的控制算法代表包括Peterka的預測控制算法和Clarke等人的廣義預測控制算法。但是傳遞函數模型在處理多變量對象時效率較低。今年來，狀态空間模型用于預測控制設計越來越多。本書也将使用狀态空間模型。
二、從連續時間模型到增廣模型
2.1 連續模型到離散模型
連續時間狀态空間模型可以表示為：
對應的離散時間狀态空間模型可以表示為：
在MATLAB中可以通過c2dm函數，并指定一個采樣時間間隔Δt，将連續模型轉換成離散模型。
例：給定連續時間狀态空間模型如下，求出對應的離散時間狀态空間模型，采樣時間間隔為1。
MATLAB程式如下：
Ac = [0 1 0; 3 0 1; 0 1 0 ];
Bc= [1; 1; 3];
Cc=[0 1 0];
Dc=zeros(1,1);
Delta_t=1;
[Ad,Bd,Cd,Dd]=c2dm(Ac,Bc,Cc,Dc,Delta_t);
得到的離散模型矩陣如下：
2.2 離散模型到增廣模型
增廣模型如下：
其中，
Am, Bm, Cm分别對應離散模型中的Ad, Bd, Cd。(A,B,C)就是表示的增廣模型。
例：離散時間狀态空間模型的系統矩陣如下，求對應的增廣模型矩陣(A,B,C)。
首先要确定om的列數，由增廣模型可知，om的列數與Am的行數相同，那麼n1=2，om=[0 0]，将Am, Bm, Cm代入增廣模型矩陣方程中，就可以得到增廣模型矩陣：
用MATLAB實作上述過程如下：
Ad = [1 1;0 1];
Bd = [0.5;1];
Cd = [1 0];
[m1,n1]=size(Cd);
[n1,n_in]=size(Bd);
A_e=eye(n1+m1,n1+m1);
A_e(1:n1,1:n1)=Ad;
A_e(n1+1:n1+m1,1:n1)=Cd*Ad;
B_e=zeros(n1+m1,n_in);
B_e(1:n1,:)=Bd;
B_e(n1+1:n1+m1,:)=Cd*Bd;
C_e=zeros(m1,n1+m1);
C_e(:,n1+1:n1+m1)=eye(m1,m1);
得到的結果如下：
與直接代入計算的結果完全相同，驗證了MATLAB程式的正确性。
三、狀态預測與優化
3.1 狀态和輸出變量的預測
假設目前時刻為ki，x(ki)為通過檢測得到的系統狀态變量，這裡假設狀态變量可以檢測到。那麼Nc個未來的控制變量序清單示為：
（ki不是表示目前的時刻嗎，為什麼△u(ki)又表示的未來的控制量？因為目前時刻的控制量也需要通過優化才能得出，是以相對來說△u(ki)也是未來的控制量）
Np個預測的狀态變量表示為：
代入增廣模型中，得到Np個預測的狀态變量為：
那麼，Np個預測的輸出變量為：
可以看出，所有預測的變量都隻與目前觀測到的系統狀态x(ki)以及未來的控制變量序列
有關。
将未來的輸出變量序列和未來的控制變量序列寫成：
在單輸入，單輸出系統中，Y的維數是Np×1，△U的維數是Nc×1。将1.10到1.11合起來寫得到：
其中，
3.2 優化
在ki時刻的期望輸入值為r(ki)，需要通過MPC控制，将系統輸出盡可能的接近系統輸入。假設在優化視窗内，輸入是恒定值。問題就轉化成了需要在優化視窗内，找到一個最佳的控制變量序列，使得系統輸出與輸入的誤差最小。
将輸入r(ki)改寫成輸入序列：
定義代價函數為：
代價函數的第一項反映了輸入信号與預測的輸出信号之間的誤差，第二項則是為了限制控制量△U，避免控制量太大。
其中，
是對角矩陣：
其中，rw是影響閉環性能的可調參數。如果rw=0，則說明不用考慮控制變量△U有多大，隻考慮控制誤差越小越好。如果rw取值比較大，則表示需要仔細考慮△U的大小。
将J寫成：
求一階導：
當J的一階導為零時，可以求得控制變量的最優解：
其中假設了矩陣
存在。
因為：
将優化的△U寫成與r(ki)和x(ki)有關的形式為：
例：一階系統狀态空間方程如下：
其中a=0.8, b=0.1 ，求 (1) 增廣狀态空間模型。(2) 計算矩陣F、Φ，以及
(3) 假設在ki=10時刻，r(ki)=1，
當rw分别為0和10的時候，求最優解△U，并比較結果。
解：增廣狀态空間方程為：
矩陣F和Φ如下：
矩陣F中的系數按如下方式計算：
其中，
矩陣Φ中的系數按如下方式計算：
将a=0.8, b=0.1, 代入得：
向量
與矩陣
的最後一列相同，是因為矩陣F的最後一列和
相同都是1。
當rw=0時，即控制系統隻考慮減小輸入輸出的誤差，不考慮控制量大小的情況下，計算得到的Nc個控制量序列為：
可以看出控制量變化比較大。
當rw=10時，即限制了控制量變化大小的情況下，計算得到的Nc個控制量序列為：
與之前一種情況相比，控制量平穩很多。
将這兩種情況下的輸出量和狀态量曲線畫出來如下圖所示。可以看到圖a中當狀态變化量減少到0時，預測輸出可以到達期望輸入。而在圖b中，預測輸出緩慢向期望輸入靠近，在預測視窗内誤差比較大，收斂速度比a慢，但是狀态變化曲線更加平穩。
(a) rw=0 (b) rw=10
圖1 不同控制量權重下，輸出與狀态量變化曲線：1-△xm；2-y
由上述對比可以看出，如果想要控制量更平穩的變化，那麼控制信号達到穩定狀态需要更長的時間。如果将控制域從Nc=4增加到Nc=9，在rw=10的時候得到的最優控制量序列為：
可以看出，前四個數值都不一樣，而且都有所減少。在不同情況下，要根據實際需要選擇合适的Nc, Np, rw來進行MPC控制。
3.3 MATLAB實作
用MATLAB函數實作上述計算各矩陣的過程，最關鍵的是求矩陣F和Φ。其中Φ是托普利茲矩陣，通過定義其第一列，後面的列都可以通過移動前一列獲得。
定義一個函數，函數的輸入是離散時間狀态空間模型(Ap, Bp, Cp)和Np, Nc。MATLAB程式如下：
function [Phi_Phi,Phi_F,Phi_R,A_e, B_e,C_e]=mpcgain(Ap,Bp,Cp,Nc,Np)
[m1,n1]=size(Cp);
[n1,n_in]=size(Bp);
A_e=eye(n1+m1,n1+m1);
A_e(1:n1,1:n1)=Ap;
A_e(n1+1:n1+m1,1:n1)=Cp*Ap;
B_e=zeros(n1+m1,n_in);
B_e(1:n1,:)=Bp;
B_e(n1+1:n1+m1,:)=Cp*Bp;
C_e=zeros(m1,n1+m1);
C_e(:,n1+1:n1+m1)=eye(m1,m1);
n=n1+m1;
h(1,:)=C_e;
F(1,:)=C_e*A_e;
for kk=2:Np
h(kk,:)=h(kk-1,:)*A_e;
F(kk,:)= F(kk-1,:)*A_e;
end
v=h*B_e;
Phi=zeros(Np,Nc); %declare the dimension of Phi
Phi(:,1)=v; % first column of Phi
for i=2:Nc
Phi(:,i)=[zeros(i-1,1);v(1:Np-i+1,1)]; %Toeplitz matrix
end
BarRs=ones(Np,1);
Phi_Phi= Phi'*Phi;
Phi_F= Phi'*F;
Phi_R=Phi'*BarRs;
為了驗證程式的正确性，在MATLAB的work space輸入：
Ap = 0.8;
Bp = 0.1;
Cp = 1;
Nc = 4;
Np = 10;
[Phi_Phi, Phi_F, Phi_R, A_e, B_e, C_e] = mpcgain(Ap, Bp, Cp, Nc, Np);
運作後得到：
可以看出計算結果與上例中的相同，表明程式正确。
參考文獻：
[1] Liuping Wang. Model Predictive Control System Design and Implementation Using MATLAB.