一個模型預測控制（MPC）的簡單實現

1 生活中的啟示

情景如下：你們團隊每天早晨開一次例會，主要會議內容是你匯報工作進度，領導根據工作目標和工作進度，制定當天的工作計划，你領到工作計划后開始工作。每天都這樣周而復始，從領導的角度看，這件工作實現了“閉環”，工作進度“可控”，這就是閉環控制系統。

圖1 閉環控制系統

不同的領導（控制器）水平有高又低，按照介紹控制器先踩一腳PID的國際慣例，設:

則PID控制器的控制律為：

從該控制律中我們可以看到PID的兩個問題：

1， PID控制器不具有“前瞻性”：參與計算的各個量，有當前的  ，上個控制周期的  ，以及之前所有的  累計和，偏偏沒有未來的  。這個“領導”目光過於短淺，只求今天能完成任務，哪怕第二天公司就倒閉了，那也是第二天要解決的麻煩，今天該干啥干啥。

2， PID屬於無模型控制。作為一個“領導”，PID僅僅通過工作目標和工作進度上的差距，以及三個近乎魔法般的數字，就制定了工作計划，完全不考慮你的工作能力和這項工作的難度，這是非常失職的。

為了提高工作的的可行性，經過思考，我們還可以有另一種方案：

• 領導聽完你的匯報后，根據工作進度（系統狀態）、工作目標（參考值）、你的工作能力和任務難度（系統模型），制定了未來10天的工作計划，每一天要干什么都寫得清清楚楚明明白白。
• 領導將這十天計划里面第一天的計划交給你，剩下九天的計划銷毀。
• 第二天重復上述過程。

顯然第二種方案更具有可行性，它能根據任務的完成情況及時調整工作計划，同時兼顧了未來的工作計划、你的工作能力和任務的難度，更符合我們的認知。

當然這種方案還有一點小缺陷，比如“剩下九天計划銷毀”，領導看了這個方案肯定心里不舒服，辛辛苦苦做的計划憑什么銷毀呢？為了彌補這個缺陷，我們提出了“控制時域”和“預測時域”的概念。控制時域是指領導做多少天的計划，預測時域是指領導思考多長遠的進度。比如，領導考慮到計划跟不上變化，想太長遠沒意義，就先做了五天的計划（控制時域是五天），然后第6-10天的計划和第5天的計划一模一樣，領導預測了下按照這樣的十天計划，十天后進度完成的還算可以（預測時域是十天）。這樣依然給你第一天的計划，銷毀剩下九天的計划，領導的工作量卻少了許多。

模型預測控制的基本思想就蘊含在上述過程中，它利用一個已有的模型、系統當前的狀態和未來的控制量，來預測系統未來的輸出，然后與我們期望的系統輸出做比較，得到一個損失函數，即：

由於上式中模型、當前狀態、期望輸出都是已知的，因此只有未來控制量一個自變量。采用二次規划的方法求解出某個未來控制量，使得損失函數最小，這個未來控制量的第一個元素就是當前控制周期的控制量。

2 實際控制的例子

2.1 問題描述

在無限光滑的一維水平直線上有一個質量為  的滑塊，初始位置與初始速度都為0，現需要設計控制器，在傳感器測得滑塊位置  的基礎上，為滑塊提供外力  ，使其跟隨參考點 

圖2 例圖

2.2 預測模型

首先建立動力學方程：

選取狀態向量  （除非特殊說明，后文中x表示狀態向量，而不是滑塊位置)，構建系統狀態方程為：

其中  。

2.3 預測模型離散化

采用前向歐拉法將狀態方程離散化：

其中  ，T為控制周期。

2.4 預測

MPC方法的一個獨特之處就是需要對未來系統狀態進行預測，我們記未來  個控制周期內預測的系統狀態為：

 稱為預測時域，括號中  表示在當前  時刻預測  時刻的系統狀態，以此類推。另外，預測動態系統未來狀態時，還需要知道預測時域內的控制量  ：

這是我們接下來將要求解的優化問題的獨立變量。

現在，我們可以通過離散化狀態方程依次對未來  個控制周期的系統狀態進行預測：

整合成矩陣形式：

其中，psi ：     ，             theta ：  

上式中的下三角形式，直接反映了系統在時間上的因果關系，即  時刻的輸入對  時刻的輸出沒有影響，  時刻的輸入對  和  時刻沒有影響，等等。

2.5 優化

這一節我們將求解預測時域內的控制輸出  ，在求解優化問題之前，我們首先明確優化問題的數學描述。

我們的控制目標是使系統的狀態跟蹤期望的一條軌跡，通常稱為參考值，定義預測時域內的參考值序列：

注意，在  時刻進行控制的時候，控制器就必須已經得到了  時刻到  時刻的參考值，而PID就不需要這么多信息，這是MPC的一個缺點。

我們希望尋找最佳的控制量  ，使得預測時域內的狀態向量與參考值越接近越好，這是一個開環最優控制問題。為此，我們用預測狀態向量與參考值之間的累計誤差定義一個簡單的優化目標函數：

經常地，我們不希望控制動作太大，優化目標函數再添加一項對控制量的約束：

因此，該優化問題可以描述如下：

我們將優化函數  展開后合並同類項：

注意：：：：：：：：

　　在這里如果你仔細去化簡上面的公式，會發現一個問題 2 ET* Q * thita * U是怎么來的，化簡的結果應該是ET* Q * thita * U和UT*thitaT*Q*E，為什么上面的結果是 2 ET* Q * thita * U，難道說ET* Q * thita * U和UT*thitaT*Q*E相等不成，是的，這兩個式子確實相等，但這兩個是各自的轉置后的形式啊（Q是對角矩陣，轉置和原來一樣這個應該沒啥疑問），怎么相等呢，這和我們認知的矩陣轉置一般不等於原矩陣不符啊，等等因為這個特殊唄，特殊形式中對角矩陣算一個，那這個是什么特殊情況呢？別急，我推導一遍你就清楚了。

上式中  是常數項，對“  為何值時  取得最小值”這一問題沒有影響，因此直接舍去。

如圖3，matlab輸入 “help quadprog”查看二次型優化函數quadprog的說明文檔，令：

可得最終優化目標函數，至此可直接調用matlab quadprog函數求解  ，將  的第一個元素提取出來，作為本控制周期的控制量。

圖3 matlab quadprog函數

2.6 仿真

對於2.1中例子的動力學方程：

兩邊同時拉普拉斯變換：

可得傳遞函數：

在simulink中搭建仿真環境如圖4，並編寫MPC控制器（2.2—2.5）：

圖4 simulink中仿真

function u = MPCcontroller(pos_ref, pos, vel)  2 %參數設置  3 m    = 1.05;          %滑塊質量,增加了5%作為建模誤差  4 T    = 0.01;          %控制周期10ms //采樣時間  5 p    = 45;            %控制時域（預測時域）  6 Q    = 10*eye(2*p);   %狀態誤差權重 //因為Xk中有兩個狀態變量所以P乘以2以下道理相同  7 W    = 0.0001*eye(p); %控制輸出權重  8 umax = 100;           %控制量限制，即最大的力  9 Rk   = zeros(2*p,1);  %參考值序列 10 Rk(1:2:end) = pos_ref;%參考位置由函數參數指定 11 Rk(2:2:end) = vel;    %參考速度跟隨實際速度 12 %構建中間變量 13 xk    = [pos;vel];    %xk 14 A_    = [1 T;0 1];    %離散化預測模型參數A 15 B_    = [0;T/m];      %離散化預測模型參數B 16 psi   = zeros(2*p,2); %psi
　　%構建psi 17 for i=1:1:p 18     psi(i*2-1:i*2,1:2)=A_^i; 19 end
　　%構建theta 20 theta = zeros(2*p,p);     %theta 21 for i=1:1:p 22     for j=1:1:i 23         theta(i*2-1:i*2,j)=A_^(i-j)*B_; 24  end 25 end 26 E = psi*xk-Rk;            %E 27 H = 2*(theta'*Q*theta+W); %H
f = (2*E'*Q*theta)';      %f 29 %優化求解 30 coder.extrinsic('quadprog'); 31 Uk=quadprog(H,f,[],[],[],[],-umax,umax); 32 %返回控制量序列第一個值 33 u = 0.0;                %指定u的類型 34 u = Uk(1);              %提取控制序列第一項

p=40，仿真結果

p=60，仿真結果

 

參考文獻

陳虹. 模型預測控制[M]. 科學出版社, 2013.

龔建偉, 姜岩, 徐威. 無人駕駛車輛模型預測控制[M]. 北京理工大學出版社, 2014.

轉自：https://zhuanlan.zhihu.com/p/141871796?utm_source=wechat_session