四、滾動控制
盡管優化向量 Δ U \Delta U ΔU包含了 N c N_c Nc個控制變量序列，但是根據滾動控制原理，通常隻使用這一序列中的第一個變量，即 Δ u ( k i ) \Delta u(k_i) Δu(ki)，忽略這一序列中的其餘變量。當下一可采樣周期到來時，使用新觀測的狀态變量 x ( k i + 1 ) x(k_i+1) x(ki+1)，重複以上優化過程，計算出新的控制變量序列，進而使優化視窗不斷向前推進，每一采樣時刻都進行實時預測與優化。
例：繼續上例，一階系統模型狀态空間為：
考慮r_w=0的情況。初始條件為： x ( 10 ) = [ 0.1 0.2 ] T x(10)=\begin{bmatrix}0.1 & 0.2 \end{bmatrix}^T x(10)=[0.10.2]T， u ( 9 ) = 0 u(9)=0 u(9)=0。
解：在采樣時間 k i = 10 k_i=10 ki=10時，上例已經計算出來的控制量為 Δ u ( 10 ) = 7.2 \Delta u(10)=7.2 Δu(10)=7.2，因為 u ( 9 ) = 0 u(9)=0 u(9)=0，那麼這一時刻的控制量為 u ( 10 ) = u ( 0 ) + Δ u ( 10 ) = 7.2 u(10)=u(0)+\Delta u(10)=7.2 u(10)=u(0)+Δu(10)=7.2， y ( 10 ) = x m ( 10 ) = 0.2 y(10)=x_m(10)=0.2 y(10)=xm(10)=0.2，将這一時刻的控制量和狀态量代入狀态方程，可以計算下一采樣時刻的狀态變量為：
在采樣時刻 k i = 11 k_i=11 ki=11，新的狀态資訊為 Δ x m ( 11 ) = 0.88 − 0.2 = 0.68 \Delta x_m(11)=0.88-0.2=0.68 Δxm(11)=0.88−0.2=0.68， y ( 11 ) = 0.88 y(11)=0.88 y(11)=0.88，可以寫成 x ( 11 ) = [ 0.68 0.88 ] T x(11)=\begin{bmatrix}0.68 & 0.88 \end{bmatrix}^T x(11)=[0.680.88]T，可以得到新一輪的優化結果：
這一時刻的優化控制量為 u ( 11 ) = u ( 10 ) + Δ u ( 11 ) = 2.96 u(11)=u(10)+\Delta u(11)=2.96 u(11)=u(10)+Δu(11)=2.96。使用這個控制量代入狀态方程得到：
在采樣時刻 k i = 12 k_i=12 ki=12，新的狀态資訊為 Δ x m ( 12 ) = 1 − 0.88 = 0.12 \Delta x_m(12)=1-0.88=0.12 Δxm(12)=1−0.88=0.12， y ( 12 ) = 1 y(12)=1 y(12)=1，可以寫成 x ( 11 ) = [ 0.12 1 ] T x(11)=\begin{bmatrix}0.12 & 1 \end{bmatrix}^T x(11)=[0.121]T，可以得到新一輪的優化結果：
這一時刻的優化控制量為 u ( 12 ) = u ( 11 ) + Δ u ( 12 ) = 2 u(12)=u(11)+\Delta u(12)=2 u(12)=u(11)+Δu(12)=2。使用這個控制量代入狀态方程得到：
在 k i = 13 k_i=13 ki=13采樣時刻，新的狀态資訊為 Δ x m ( 13 ) = 1 − 1 = 0 \Delta x_m(13)=1-1=0 Δxm(13)=1−1=0， y ( 13 ) = 1 y(13)=1 y(13)=1，可以寫成 x ( 11 ) = [ 0 1 ] T x(11)=\begin{bmatrix}0 & 1 \end{bmatrix}^T x(11)=[01]T，可以得到新一輪的優化結果：
上述計算過程中，控制信号及狀态變量如下圖所示。
圖2 控制信号、狀态變化量和輸出量曲線。這個例子還說明 Δ U \Delta U ΔU向量在不同時刻的的差異執行個體。注意，輸出響應達到的設定點信号時， Δ U \Delta U ΔU向量為零。
4.1 閉環控制
在 k i k_i ki時刻，優化的向量參數 Δ U \Delta U ΔU通過下式計算：
其中
對應于期望輸入的變化，
對應于狀态回報控制。全都取決于系統參數。是以在時不變系統中它們是常數矩陣。由于滾動控制中，隻取 Δ U \Delta U ΔU中的第一個變量作為控制增量，是以：
其中， K y K_y Ky是下式的第一個元素：
K m p c K_{mpc} Kmpc是下式的第一行：
這是标準的時不變系統狀态回報控制。 K m p c K_{mpc} Kmpc是回報控制增益向量。代入增廣矩陣當中得到閉環方程為：
閉環特征值可以通過閉環特征方程得到：
例1.5：求上例中産生的閉環回報增益矩陣和權值 r w = 0 r_w=0 rw=0和 r w = 10 r_w=10 rw=10的閉環系統的特征值。
解：當 r w = 0 r_w=0 rw=0時，
閉環系統的特征值由矩陣 A − B K m p c A-BK_{mpc} A−BKmpc決定。從上例可知：
是以特征值為：
幾乎在複平面的原點上。
當 r w = 10 r_w=10 rw=10時，
其特征值為：
表明閉環系統的動力學響應比 r w = 0 r_w=0 rw=0時的響應要慢得多，這個結論也印證了之前例子中的分析。
例1.6：假設一個連續時間系統的傳遞函數如下：
其中， ω = 10 \omega = 10 ω=10。用采樣周期 Δ t = 0.01 \Delta t=0.01 Δt=0.01将其離散化。當 N c = 3 N_c=3 Nc=3、 R ˉ = 0.5 I \bar{R}=0.5I Rˉ=0.5I時，當參數 N p = 20 N_p=20 Np=20和 N p = 200 N_p=200 Np=200時的參數敏感程度。
解：使用如下MATLAB函數将系統轉換成連續時間狀态空間模型：
omega=10;
numc=omega^2;
denc=[1 0.1*omega omega^2];
[Ac,Bc,Cc,Dc]=tf2ss(numc,denc);
得到連續狀态空間方程為：
根據例1.1中的方法，得到增廣狀态空間方程為：
假設目前采樣時刻為 k = 10 k=10 k=10，初始條件為
當 N p = 20 N_p=20 Np=20時的 Δ U \Delta U ΔU為：
狀态回報控制增益為：
閉環特征值為：
當時的為： N p = 200 N_p=200 Np=200時的為:
狀态回報控制增益為：
狀态回報控制增益為：
閉環特征值為：
對比上述的求解結果可以看出，不同 N p N_p Np對優化的控制變量 Δ U \Delta U ΔU影響巨大。這一比較研究說明了設計中控制變量 Δ U \Delta U ΔU關于的選擇存在的敏感性。
仔細考察這一現象背後的原因，從研究海塞矩陣開始：
當 N p = 20 N_p=20 Np=20時有：
其條件數為：
然而當 N p = 200 N_p=200 Np=200時有：
其條件數為：
可以看出，當 N p N_p Np由20增大到200時，海塞矩陣的條件數會劇烈增加，這就是導緻控制量優化結果對參數 N p N_p Np敏感度高的深層原因。
如果預測和控制周期較短，閉環預測控制系統不一定是穩定的。通常，這些參數是需要根據閉環控制穩定性和控制效果進行整定。後面會提出如何使用較長的預測和控制周期來確定閉環控制的穩定性。
4.2 Matlab實作
如何通過時域滾動控制方法在Matlab中實作模型預測控制？這是本教程想要達到的目的。假設系統的狀态空間模型為：
首先需要将系統模型和控制域長度 N c N_c Nc、預測域長度 N p N_p Np輸入m檔案。将上述二階雙積分系統，以及 N c = 4 N_c=4 Nc=4、 N p = 4 N_p=4 Np=4輸入到m檔案，如下：
Ap=[1 1;0 1];
Bp=[0.5;1];
Cp=[1 0];
Dp=0;
Np=20;
Nc=4;
然後調用mpcgain.m檔案中的函數得到增廣模型矩陣，以及預測控制需要的各種矩陣。初始的狀态變量為0，初始的狀态回報變量也為0，給定的期望輸入為1，仿真采樣個數為100，初始的控制量和輸出量都是0：
[Phi_Phi,Phi_F,Phi_R,A_e, B_e,C_e] = mpcgain(Ap,Bp,Cp,Nc,Np);
[n,n_in]=size(B_e);
xm=[0;0];
Xf=zeros(n,1);
N_sim=100;
r=ones(N_sim,1);
u=0; % u(k-1) =0
y=0;
當采樣時間為kk時，向量 Δ U \Delta U ΔU根據期望輸入r(kk)和狀态向量xf計算得到。控制量的增量隻取 Δ U \Delta U ΔU的第一個元素，控制量為：
預設的權重系數為0.1。利用所産生的控制信号對系統的狀态和輸出進行仿真，得到的狀态變量更新到變量xf中：
for kk=1:N_sim;
DeltaU=inv(Phi_Phi+0.1*eye(Nc,Nc))*(Phi_R*r(kk)-Phi_F*Xf);
deltau=DeltaU(1,1);
u=u+deltau;
u1(kk)=u;
y1(kk)=y;
xm_old=xm;
xm=Ap*xm+Bp*u;
y=Cp*xm;
Xf=[xm-xm_old;y];
end
将系統輸出、控制量曲線顯示出來：
k=0:(N_sim-1);
figure
subplot(211)
plot(k,y1)
xlabel('Sampling Instant')
legend('Output')
subplot(212)
plot(k,u1)
xlabel('Sampling Instant')
legend('Control')
運作程之後的結果如下：
可以将權重系數修改到5做一下對比，可以發現，閉環響應速度會明顯變慢：
通過修改 A p A_p Ap, B p B_p Bp, C p C_p Cp可以用來仿真其他系統模型，修改不同參數 N p N_p Np, N c N_c Nc, r w r_w rw來對比不同參數下的控制效果。
參考文獻：
[1] Liuping Wang. Model Predictive Control System Design and Implementation Using MATLAB.