第一部分

模型預(yù)測控制

01

最優(yōu)控制理論處理的問題通常是找到一個滿足容許控制的 u*，把它作用于系統(tǒng)（被控對象）?(t)=f(x(t),u(t),t) 從而可以得到系統(tǒng)的狀態(tài)軌跡 x(t)，使得目標(biāo)函數(shù)最優(yōu)。對于軌跡跟蹤問題，那目標(biāo)函數(shù)就是使得這個軌跡在一定的時間范圍[t0tf]內(nèi)與我們期望的軌跡（目標(biāo)）x*(t) 越近越好。最優(yōu)控制問題更一般的表達(dá)如下：在被控對象符合動力學(xué)原理（狀態(tài)方程）和狀態(tài)約束

的條件下，求解控制函數(shù) u(t) 以使得連續(xù)時間性能指標(biāo)

最小。其中 t0?是初始時刻，tf?是終端時刻，E 是終端時刻代價，g 是運行時刻代價。例如，更具體的場景，對于時間最短問題（例如控制電流使得最短時間充電到 SOC100%），性能指標(biāo)即：

對于最小燃料消耗問題（接下來本文中主要使用的例子，假設(shè)燃料消耗與控制量 u 成正比，并且 u∈[0 1]），性能指標(biāo)變?yōu)椋?/p>

最優(yōu)控制的性能指標(biāo)函數(shù)和約束函數(shù)都是泛函，因為他們的自變量 u(t) 本身也是個函數(shù)。于是，求解性能指標(biāo)函數(shù)最優(yōu)的問題是一個求泛函極值的問題，類比于函數(shù)極值通過求導(dǎo)來獲得極值條件，泛函極值則通過變分法來獲得極值條件。而基于變分的兩大方法: 龐特里亞金極小值原理和動態(tài)規(guī)劃（求解 Hamilton-Jacobi-Bellman 方程），則是最優(yōu)控制問題求解的主要方法。對于簡化的最優(yōu)控制問題，例如被控對象為線性系統(tǒng)

性能指標(biāo)是二次型的無限時域的連續(xù)時間泛函：

（其中第一項是對狀態(tài)（例如追蹤誤差）的要求，第二項是對控制能量的要求，Q 和 R 是權(quán)重矩陣），優(yōu)化區(qū)間考慮整個時間域 [t0=0 tf=∞]（也就是 Linear-Quadratic Regulator (LQR)?）, 可以使用極小值原理或動態(tài)規(guī)劃或 Ricatti 代數(shù)方程求得閉環(huán)形式的最優(yōu)控制。例如，對一個車載倒立擺系統(tǒng)[1]：

其中系統(tǒng)輸出是車輛位移 x 和倒立擺角度 θ，控制輸出 u 是作用在車上的水平力?？刂剖菢?biāo)量，狀態(tài)是 4 維

Q = [1,0,0,0;...?

0,0,0,0;...?

0,0,1,0;...

0,0,0,0];

R = 1;?

[K,S,P] = lqr(sys,Q,R)

利用 lqr 函數(shù)就可以求得最優(yōu)增益矩陣 K。但對于更一般的非線性系統(tǒng)，使用極小值原理或動態(tài)規(guī)劃求解解析解幾乎不可能，數(shù)值解計算復(fù)雜度較高，所以可以考慮一定的問題簡化，僅考慮未來的一小段時間最優(yōu)（而非全時間域“最優(yōu)”），將最優(yōu)控制問題轉(zhuǎn)化未來幾個時間步的在線數(shù)值優(yōu)化問題，例如 MPC，可以大大降低“最優(yōu)”控制計算的復(fù)雜度。

我們通過和人的行為進(jìn)行類比來解釋 MPC 的思想[2]：設(shè)想開車的場景，我們的目標(biāo)是讓車按一定的軌跡行駛，人的控制方式和模型預(yù)測控制器的工作方式很像，司機(jī)的大腦中可能有一些經(jīng)驗（ 類似虛擬的動力學(xué)模型），于是他會利用這些經(jīng)驗（虛擬模型）在大腦中去預(yù)先“仿真”這個過程，可以預(yù)測假如他采取了某些動作（油門，剎車，轉(zhuǎn)向等等）之后的一段時間內(nèi)汽車會有多快或有多少轉(zhuǎn)向，從而他會從這些可選的一系列動作中選擇一個最能使汽車在未來一段時間內(nèi)接近期望軌跡的控制動作作為當(dāng)前時刻的行為，例如踩油門，剎車等等，并在每個時刻不斷重復(fù)類似的思想，從而驅(qū)動車輛到期望的軌跡。

算法介紹

（圖表1）

模型預(yù)測控制 (MPC) 是一種最優(yōu)控制技術(shù)，在每個控制周期 tk , MPC 控制器獲得系統(tǒng)（被控對象）當(dāng)前的狀態(tài)。接下來它利用基于系統(tǒng)當(dāng)前狀態(tài)和系統(tǒng)動態(tài)模型通過在有限時間上預(yù)測系統(tǒng)未來 p 步狀態(tài)軌跡，并與目標(biāo)軌跡構(gòu)建代價函數(shù)和約束，進(jìn)行一次開環(huán)優(yōu)化問題（要優(yōu)化的變量是作用在被控對象上的控制輸入序列）求解，得到未來一段時間 [tk,tk+1,···,tk+p] 的控制輸入序列?[uk,uk+1,···,uk+p] 。當(dāng)然對于求解得到的控制輸入序列，控制器只將最近時刻 tk 的控制?uk 作用于系統(tǒng)（被控對象）而忽略掉控制序列的其他值[uk+1,···,uk+p], 在下一個時刻tk+1，MPC 控制器會重復(fù)上述優(yōu)化求解過程重新計算控制序輸入列并只將第一個控制值作為當(dāng)前時刻的控制量作用于系統(tǒng)，依次類推，重復(fù)上述過程進(jìn)行下一時刻的求解。所以 MPC 在每個時刻都在線進(jìn)行一個含約束的優(yōu)化問題的求解（滾動優(yōu)化，特殊情況除外）。接下來我們看一下這個優(yōu)化問題是如何定義以及求解的。MPC 的優(yōu)化問題 (更具體的二次規(guī)劃，QP) 包含這幾項主要內(nèi)容[3]：

代價函數(shù): 用于度量控制器性能，目標(biāo)是求它的最小值。

約束：目標(biāo)解必須滿足的條件，例如控制量（manipulated variables, MV）和被控對象狀態(tài)輸出的物理邊界。

控制與決策：得到優(yōu)化的控制量（MV）

代價函數(shù)

一個典型的代價函數(shù)由四部分組成，每一部分關(guān)注控制器性能的一個方面。

這里面的 zk?就是 QP 問題的決策變量（控制變量）。其中每一項量化一個性能指標(biāo)。

Jy(zk) 量化系統(tǒng)輸出與目標(biāo)軌跡的跟蹤效果

Ju(zk) 用于量化控制輸入與目標(biāo)控制變量跟蹤效果（在很多應(yīng)用中，控制器還需要保證控制變量（MV）保持在某個目標(biāo)附近，尤其在控制量數(shù)目遠(yuǎn)多于系統(tǒng)輸出量數(shù)目的情況下）

JΔu(zk) 控制變量波動約束（多數(shù)應(yīng)用場景下一般都希望控制變量有較小的變動或調(diào)整，而不希望較大的突變）。

例如對于常用的 Jy(zk)?，即量化圖表1中黑色雙箭頭代表的區(qū)域?？刂破鞯哪繕?biāo)是將被控對象保持在指定的參考軌跡附近。MPC 控制器通過如下計算得到一個標(biāo)量作為性能度量來實現(xiàn)目標(biāo)軌跡跟蹤:

其中，底標(biāo) i 是對 p 個預(yù)測步的循環(huán)，j 是對 ny 個輸出的循環(huán)。

k:當(dāng)前控制周期

p:預(yù)測區(qū)間（Prediction Horizon）

ny:被控對象輸出變量的數(shù)量

zk:QP 問題決策變量，也就是對應(yīng)的時間步上的控制器輸出序列，其中

yj(k+i|k):從 k 時刻開始的未來第 i 個（共 p 個）時間步長被控對象第 j 個（共 ny?個）輸出的預(yù)測值

rj(k+i|k):從時刻開始的未來第個（共個）時間步長被控對象第個（共個）輸出的目標(biāo)參考值

：未來第個時間步，被控對象第個（共個）輸出的權(quán)重因子（無量綱）

約束

MPC 中最常見的約束就是邊界約束，例如針對被控對象、控制變量（MV）、控制變量變化量的邊界約束，如下：

yj,min(i),yj,max(i):從 k 時刻開始的未來第 i 個時間步長被控對象第 j 個（共 ny?個）輸出的下界和上界

uj,min(i),uj,max(i):從 k 時刻開始的未來第 i 個時間步長第 j 個（共 nu?個）控制變量的下界和上界

Δuj,min(i),Δuj,max(i):從 k 時刻開始的未來第 i 個時間步長第 j 個（共 nu?個）控制變量的下界和上界

預(yù)測與決策

每次滾動優(yōu)化計算過程中，模型預(yù)測控制器可以獲得整個預(yù)測區(qū)間 p 時間步上的參考軌跡 rj(k+i|k)。模型預(yù)測控制器中的“模型”包括被控對象模型、擾動模型和噪聲模型，如圖表2，控制器每次滾動優(yōu)化計算過程中，使用這些模型加上可調(diào)（可優(yōu)化）的控制變量 zk?來預(yù)測被控對象的輸出?yj(k+i|k)。

（圖表2）

上述預(yù)測過程可以通過簡化的狀態(tài)空間離散預(yù)測模型表示：

其中 v(k) 是可測干擾輸入，d(k) 是不可測干擾輸入白噪聲。

我們詳細(xì)描述一下在 k=0 時刻預(yù)測模型在預(yù)測區(qū)間 p 上的軌跡預(yù)測過程：對所有預(yù)測瞬時 i, 將 d(i) 設(shè)為0，可得，問題變成一個遞歸計算問題(從 k=1到 k=0)即可得到 y(i|0) 序列：

整理成矩陣形式：

其中（此處不詳細(xì)列出這些矩陣表達(dá)式，可以根據(jù)上式遞歸自行推導(dǎo)）

組合各個元素得到優(yōu)化問題定義

將上述的代價函數(shù)，約束，預(yù)測模型矩陣結(jié)合起來構(gòu)建 MPC 的開環(huán)優(yōu)化問題如下：

狀態(tài)變量 x(k)∈Rnx，控制變量 u(k)∈Rnu，滿足上述系統(tǒng)動力學(xué)方程以及時域約束，通過預(yù)測計算，求解如下最小化性能指標(biāo)對應(yīng)的優(yōu)化變量zk

其中

即在預(yù)測區(qū)間 p 上的待優(yōu)化的控制輸入序列。

示例：

非線性 MPC 的飛行機(jī)器人軌跡開環(huán)規(guī)劃和閉環(huán)反饋控制?

我們通過一個非線性 MPC 的飛行機(jī)器人軌跡優(yōu)化與控制控制器設(shè)計示例來主要介紹介紹 MPC 設(shè)計方法[4]。MPC 即可以用于軌跡追蹤控制（在線實時），也可以用于規(guī)劃（進(jìn)行一次開環(huán)的優(yōu)化）。本示例先用 MPC 進(jìn)行了一次開環(huán)優(yōu)化求解完成了規(guī)劃，找到了從一個位置到另一個位置的最佳軌跡，然后再結(jié)合狀態(tài)估計器：擴(kuò)展卡爾曼濾波器，來控制機(jī)器人沿規(guī)劃的最優(yōu)軌跡進(jìn)行閉環(huán)控制。

這個簡化的例子中的飛行機(jī)器人有四個推進(jìn)器在平面空間中移動。

該模型動力學(xué)系統(tǒng)如下:

function dxdt = FlyingRobotStateFcn(x, u)

% 狀態(tài)方程：x 是六個狀態(tài)，u是推力控制。

%x(1)???– 質(zhì)心x坐標(biāo)

%x(2)??– 質(zhì)心y坐標(biāo)

%x(3)? – theta, 機(jī)器人(推力)方向

%x(4)??– x方向的速度vx

%x(5)??– y方向的速度vy

%x(5)??– omega, 角速度

…

dxdt = zeros(6,1);

dxdt(1) = x(4);

dxdt(2) = x(5);

dxdt(3) = x(6);

dxdt(4) = (T1 + T2)*cos(theta);

dxdt(5) = (T1 + T2)*sin(theta);

dxdt(6) = alpha*T1 - beta*T2;

end

示例假設(shè)機(jī)器人有四個物理推力 [u(1)u(2)u(3)u(4)]，范圍從0到1。

軌跡規(guī)劃

前面提到了，MPC 用于軌跡規(guī)劃其實是求解了一次開環(huán)的帶約束的優(yōu)化問題，區(qū)別于后面的針對規(guī)劃好的軌跡的路徑跟蹤（利用反饋狀態(tài)的滾動優(yōu)化）。

先定義本示例的軌跡規(guī)劃問題：機(jī)器人最初停留在 [-10, -10] 位置，方向角 theta 為 pi/2。本例中的飛行任務(wù)是在 12s 內(nèi)移動機(jī)器人并停在最終位置 [0, 0]，方向角 theta 為 0。目標(biāo)是找到最優(yōu)路徑，使推進(jìn)器在控制過程中消耗的燃料總量最小。

在此例中，對于規(guī)劃軌跡，我們開環(huán)優(yōu)化需要使用全部時間長度：12s。設(shè)定一個采樣時間 Ts=0.4s，所以對應(yīng)的預(yù)測區(qū)間為 30 步（12/0.4）。規(guī)劃的運行周期通?？梢员确答伩刂浦芷诟?，或更低頻，所以有更多計算資源允許計算一個相對較大的優(yōu)化問題。對于軌跡規(guī)劃我們可以創(chuàng)建一個多級非線性MPC對象，它允許對每個預(yù)測步都定義不同的代價函數(shù)和約束。

% 具有 nx=6 個狀態(tài)和 nu=4 個控制輸入（mv）

% 預(yù)測區(qū)間 p = 30

nlobj=nlmpcMultistage(p,nx,nu);

nlobj.Ts=Ts;

for?ct=1:p

%并且每個 stage/ 預(yù)測步有自己的代價函數(shù)

其中代價函數(shù)

包含 stage 這個參數(shù)，

見下面 FlyingRobotCostFcn 函數(shù)

nlobj.Stages(ct).CostFcn='FlyingRobotCostFcn';

end

%?代價函數(shù)示例

function J=FlyingRobotCostFcn(stage, x, u)

%本例中因為推力與燃油正相關(guān)，所以這個示例每個預(yù)測步（1到p）簡化燃油消耗的表達(dá)為上四個推力之

%和。本示例沒有直接使用 stage 參數(shù)。

J = u(1) + u(2) + u(3) + u(4);

end

指定預(yù)測模型的狀態(tài)轉(zhuǎn)移方程和對應(yīng)的雅可比函數(shù)解析形式，可以大幅提升計算效率

nlobj.Model.StateFcn="FlyingRobotStateFcn";

nlobj.Model.StateJacFcn=@

FlyingRobotStateJacobianFcn;

控制目標(biāo)是在第12秒將機(jī)器人停在 [0,0] 處，角度為0弧度。將這個指定為終端狀態(tài)約束，其中最后一個預(yù)測步驟(第p+1步)的每個位置和速度狀態(tài)都應(yīng)該為零。

nlobj.Model.TerminalState = zeros(6,1);

如果可以為每個 stage 都提供代價函數(shù)和約束函數(shù)的雅可比函數(shù)解析式會大大提升優(yōu)化效率。本例沒有提供，所以它們的雅可比矩陣由 MPC 控制器使用內(nèi)置的數(shù)值差分方法計算。

每個推力的工作范圍在0到1之間，也就是 MV 的下限和上限。

for?ct = 1:nu

nlobj.MV(ct).Min = 0;

nlobj.MV(ct).Max = 1;

end

指定飛行器初始條件：

x0=[-10;-10;pi/2;0;0;0];% 機(jī)器人停在 [-10, -10], 方向角為 pi/2

u0 = zeros(nu,1);????%初始推力為0

通過調(diào)用 nlmpcmove 命令，就可以完成一次優(yōu)化，找到最佳狀態(tài)和控制（MV）軌跡。

[~,~,info] = nlmpcmove(nlobj,x0,u0);

圖表 3規(guī)劃（優(yōu)化求解）的最優(yōu)軌跡對應(yīng)六個狀態(tài)量的預(yù)測值

圖表4對應(yīng)的4個推力在最優(yōu)軌跡規(guī)劃求解時的最優(yōu)解序列（MV）

圖表5 機(jī)器人最優(yōu)軌跡的坐標(biāo)和方位信息，從?[-10 -10 pi/2]?到?[0 0 0].

軌跡跟蹤的反饋控制

在通過開環(huán)規(guī)劃得到最優(yōu)軌跡后，需要一個反饋控制器來使機(jī)器人沿著路徑移動。理論上，可以將開環(huán)規(guī)劃的到的最優(yōu)控制輸入 MV（圖表4）直接應(yīng)用于推進(jìn)器實現(xiàn)前饋控制。但在實際應(yīng)用中，通常會需要一個反饋控制器修正規(guī)劃時的建模誤差和抑制干擾，如下圖，所以接下來的任務(wù)我們就是設(shè)計 MPC Controller 和 State Estimator

圖表 6

設(shè)計MPC控制器?

本示例使用典型的非線性 MPC 控制器將機(jī)器人沿著參考軌跡移動到最終目標(biāo)位置。在這個軌跡跟蹤問題中，參考軌跡包含所有六個狀態(tài)的軌跡(輸出的數(shù)量等于狀態(tài)的數(shù)量，ny=nx)。

ny = 6;

nlobj_tracking = nlmpc(nx,ny,nu);

% 對于狀態(tài)方程和雅可比函數(shù)，我們和上面用于軌跡規(guī)劃的MPC中使用同樣的動力學(xué)模型

nlobj_tracking.Model.StateFcn=nlobj.、Model.StateFcn;

nlobj_tracking.Jacobian.StateFcn = nlobj.Model.StateJacFcn;

% 指定滾動優(yōu)化的預(yù)測區(qū)間 Prediction Horizon (不需要考慮太遠(yuǎn)的未來時間)和控制區(qū)間 Control

% Horizon (例如，僅在前幾個預(yù)測區(qū)間允許變量可優(yōu)化 來減少計算量)。

nlobj_tracking.PredictionHorizon = 10;

nlobj_tracking.ControlHorizon = 4;

非線性 MPC 的默認(rèn)代價函數(shù)是一個適合軌跡跟蹤的標(biāo)準(zhǔn)二次函數(shù)。對于軌跡跟蹤問題，代價函數(shù)中就 J(zk)=Jy(zk)+Ju(zk)+JΔu(zk), ?本示例對于跟蹤誤差 Jy(zk) 更看重，所以代價函數(shù)中的懲罰權(quán)重設(shè)置的大一些，控制抖動項 JΔu(zk)（MV 率）上的懲罰權(quán)重設(shè)置的小一些，如下：

nlobj_tracking.Weights.ManipulatedVariablesRate = 0.2*ones(1,nu);

nlobj_tracking.Weights.OutputVariables=5*ones(1,nx);

此外，在控制過程中，u1 和 u2 上下限的約束設(shè)定在 [0,1] 范圍內(nèi)，因為兩個是同側(cè)相反的方向，所以不需要同時工作，因此加一個等式約束，使 u(1)*u(2) 對所有預(yù)測區(qū)間都為 0，u3 和 u4 也是同樣。

nlobj_tracking.Optimization.CustomEqConFcn = ...

@(X,U,data) [U(1:end-1,1).*U(1:end-1,2); U(1:end-1,3).*U(1:end-1,4)];

▼

設(shè)計狀態(tài)估計器?

在這個例子中，只測量了三個位置狀態(tài) (x、y 和角度)。速度狀態(tài)是沒有測量，需要估計。此處使用擴(kuò)展卡爾曼濾波器 (EKF) 進(jìn)行非線性狀態(tài)估計。因為 EKF 需要離散時間模型，但我們之前的狀態(tài)空間模型是連續(xù)模型，所以使用梯形規(guī)則

通過解 nx 個非線性代數(shù)方程，對連續(xù)模型進(jìn)行轉(zhuǎn)換，得到離散狀態(tài)方程模型FlyingRobotStateFcnDiscreteTime.m。

其中m中的梯形規(guī)則如下實現(xiàn)：

% 利用上述梯形法則計算待定系數(shù)xk1

FUN=@(xk1) TrapezoidalRule(xk,xk1,uk,Ts,fk,ffun);

Options=optimoptions('fsolve','Display','none');

xk1 = fsolve(FUN,xk1,Options);

% 梯形法則函數(shù)

function f=TrapezoidalRule(xk,xk1,uk,Ts,fk,ffun)

% 計算梯形法則中的 f(ti+1,Yi+1)

fk1=ffun(xk1,uk);

% 通過使 f=0,求得待定系數(shù) xk1.

f=xk1-(xk + (Ts/2)*(fk1+fk));

接下來將上述離散后的模型給到 EKF 估計器用于狀態(tài)估計。

DStateFcn=@(xk,uk,Ts) FlyingRobotStateFcnDiscreteTime(xk,uk,Ts);

DMeasFcn = @(xk) xk(1:3);

% 創(chuàng)建EKF,并設(shè)定測量噪聲

EKF=extendedKalmanFilter(DStateFcn,DMeasFcn,x0);

EKF.MeasurementNoise = 0.01;

▼

閉環(huán)仿真?

現(xiàn)在圖表6中 MPC Controller, Plant, State Estimator 都有了，可以進(jìn)行追蹤控制的閉環(huán)仿真

設(shè)置好初始條件，仿真 32 個時間步，參考軌跡就是在規(guī)劃階段計算出來的最優(yōu)狀態(tài)軌跡。使用 nlmpcmove 和 nlmpcmoveopt 命令進(jìn)行閉環(huán)仿真。

Tsteps = 32;

for?k = 1:Tsteps

%獲取被控對象當(dāng)前時刻測量輸出（3維輸出數(shù)據(jù)：x、y和角度，含噪聲）.

yk = xHistory(k,1:3)' + randn*0.01;

%基于k時刻的測量值估計 k 時刻的狀態(tài)（得到6維狀態(tài)）.

xk = correct(EKF, yk);

%nlmpcmove 使用參考軌跡 Xref 進(jìn)行一次 p 個未來步長的優(yōu)化求解，得到當(dāng)前時刻的優(yōu)化控制uk.

[uk,options]=nlmpcmove(nlobj_tracking,xk,

lastMV,Xref(k:min(k+9,Tsteps),:),[],options);

%預(yù)測下一時刻的狀態(tài)和協(xié)方差矩陣，用于下步correct.

predict(EKF,uk,Ts);

%保存控制量并更新控制量用于作用在被控對象上

uHistory(k,:) = uk'; %#ok<*SAGROW>

lastMV = uk;

%基于當(dāng)前的狀態(tài) xk, 將求解的最優(yōu)的uk作用在被控對象上（連續(xù)時間 ODE）并求解 ODE 得到下一時刻的狀態(tài).

ODEFUN=@(t,xk) FlyingRobotStateFcn(xk,uk);

[TOUT,YOUT] = ode45(ODEFUN,[0 Ts], xHistory(k,:)');

%將狀態(tài)結(jié)果保留下來用于下個循環(huán)進(jìn)行測量.

xHistory(k+1,:) = YOUT(end,:);

end

將規(guī)劃的參考軌跡和實際的閉環(huán)控制軌跡進(jìn)行比較。

?

非線性 MPC 反饋控制器成功地移動機(jī)器人(藍(lán)色塊)，沿著最優(yōu)軌跡(黃色塊)，并將其停在最后一個圖中的最終位置(紅色塊)。當(dāng)然計算出的實際燃油損失高于規(guī)劃的油量。產(chǎn)生此結(jié)果的主要原因是，由于我們在反饋控制器中使用了更小時間域的預(yù)測和控制范圍最優(yōu)求解，因此與規(guī)劃階段使用全時間域優(yōu)化問題相比，每個區(qū)間的控制決策都是次優(yōu)的。

第二部分

為被控對象辨識神經(jīng)網(wǎng)絡(luò)狀態(tài)

空間模型并用于 MPC

02

區(qū)別于第一性原理建模的另一種方法是基于實驗數(shù)據(jù)的黑盒建模。對于動態(tài)系統(tǒng)（被控對象）進(jìn)行數(shù)據(jù)建模的方式有很多，包括系統(tǒng)辨識（線性系統(tǒng)，Grey-Box, Nonlinear 等等，具體見前面文章[5] 數(shù)據(jù)驅(qū)動的動態(tài)系統(tǒng)建模-系統(tǒng)辨識），以及深度學(xué)習(xí)（LSTM, Neural ODE 等等，具體見前面的文章[6] 數(shù)據(jù)驅(qū)動的動態(tài)系統(tǒng)建模-深度學(xué)習(xí)）。

我們前面看到 MPC 控制器中進(jìn)行預(yù)測時需要被控對象的狀態(tài)空間方程（Plant Model），也就是飛行器的狀態(tài)空間方程腳本。在工業(yè)應(yīng)用中，有時很難利用第一性原理手動推導(dǎo)得到形如?的非線性狀態(tài)空間動態(tài)模型。在這種情況下我們可以訓(xùn)練一個神經(jīng)網(wǎng)絡(luò)狀態(tài)空間模型，，其中 F 和 G 是兩個神經(jīng)網(wǎng)絡(luò)模型，然后作為對象模型放到 MPC 中用于預(yù)測過程。示意如下

神經(jīng)網(wǎng)絡(luò) ODE 算法與訓(xùn)練

Neural ODE 思想介紹

神經(jīng)網(wǎng)絡(luò)狀態(tài)空間模型是基于神經(jīng)網(wǎng)絡(luò) ODE 實現(xiàn)的，這里可以簡單介紹一下實現(xiàn)思路(在前面的文章 數(shù)據(jù)驅(qū)動的動態(tài)系統(tǒng)建模-深度學(xué)習(xí)中簡單介紹了訓(xùn)練的思想。在視頻[7]

https://www.bilibili.com/video/BV1sU4y1z7Qb/中介紹了實現(xiàn)過程)。

對于一個被控對象初始狀態(tài) x(t0)，作用一個控制輸入序列 [u(t0)u(t1)···u(tend)]，通過實驗測量或者仿真計算（使用 ode 求解器，例如 ode45），可以得到系統(tǒng)的輸出 [y(t1)y(t2)···y(tend)]，作為真值。我們的目的是訓(xùn)練神經(jīng)網(wǎng)絡(luò)狀態(tài)空間模型的參數(shù)?[θf，θg] 使得在相同的初始狀態(tài)和作用相同的控制輸入序列的情況下，通過對神經(jīng)網(wǎng)絡(luò)狀態(tài)空間模型

進(jìn)行計算得到的輸出序列 [?(t1)?(t2)···?(tend)] 與系統(tǒng)的真值 [y(t1)y(t2)···y(tend)] 盡量接近，也就是loss=L{?(t),y(t)} 越小越好，其中 L 可以是任何自定義的損失函數(shù)。計算神經(jīng)網(wǎng)絡(luò)狀態(tài)空間模型系統(tǒng)輸出的過程可以使用 ode 求解器，例如 ode45, 先計算狀態(tài)轉(zhuǎn)移方程

然后計算測量輸出方程

這樣就可以計算損失函數(shù)，接下來可以按神經(jīng)網(wǎng)絡(luò)的訓(xùn)練步驟對參數(shù)?[θf，θg]?進(jìn)行訓(xùn)練。

Neural ODE 訓(xùn)練過程實現(xiàn)

我們通過示意腳本來介紹訓(xùn)練的過程，主要包含六步。

將整個數(shù)據(jù)集?(一輪 /epoch) 分成小的 chrunk(minibatches).

在每次迭代 (each minibatch), 進(jìn)行預(yù)處理.

將輸入數(shù)據(jù)傳給網(wǎng)絡(luò)進(jìn)行前饋計算?(inference).

和真值進(jìn)行比較?(i.e. compute loss function).

自動微分?(compute gradients of each layer or operation
in the network w.r.t. loss)

更新模型權(quán)重和偏置?(parameters)

%遍歷輪數(shù)

for?epoch= 1:numEpochs

%每一輪的迭代數(shù)，一次迭代一個minibatch用于更新一次模型參數(shù)

for i= 1:numIterperEpoch ? ?

iter=iter + 1;

????????%創(chuàng)建每次迭代用于訓(xùn)練模型參數(shù)的 minibatch?

????????數(shù)據(jù)集，包含控制序列和輸出序列（真值）

[dlU,dlY,timeIds]=

createMiniBatch(neuralOdeTimesteps,?

UtrainCell,YtrainCell);

…

%@modelGradients 函數(shù)為了計算損失對模型

參數(shù)的梯度，總共完成了三個任務(wù)：

· 將輸入數(shù)據(jù)傳給網(wǎng)絡(luò)進(jìn)行前饋計算(inference).

· 和真值進(jìn)行比較計算loss

(i.e. compute loss function).

· 計算梯度

[grads,loss]=dlfeval(@modelGradients,

dlX0,dlYTarge,dlU(或者插值體),

parameters,tspan);

%利用計算得到的梯度使用優(yōu)化求解器 adam?

更新模型權(quán)重和偏置

[parameters,averageGrad2,averageSqGrad2]=

adamupdate(parameters,grads,averageGrad,

averageSqGrad,iter,learnRate,gradDecay,

sqGradDecay);

end

end

%modelGradient 函數(shù)完成三個任務(wù)

function[grad,loss]=modelGradients(dlTarget,

dlX0, dlU, params, tspan)

%neuralStateMode l 函數(shù)（接下來列出來的函數(shù)）

就是利用初始狀態(tài)和控制序列進(jìn)行網(wǎng)絡(luò)推斷計算。

dlY= neuralStateModel(dlX0, dlU, params, tspan);

%利用模型推斷結(jié)果和真值進(jìn)行損失計算

loss= huber(dlY, dlTarget);

%利用損失進(jìn)行梯度計算

grad= dlgradient(loss, params);

end

function Y= neuralStateModel(stateX0, controlX, Params, tspan)

%odeModel 函數(shù)對應(yīng)著狀態(tài)方程中的F(x,u|θf)，基于初始狀態(tài)，和當(dāng)前神經(jīng)網(wǎng)絡(luò)參數(shù)下的狀態(tài)方程，利用 dlode45 進(jìn)行積分求解,得到對應(yīng)的模型推斷的狀態(tài)序列 stateX

%odeModel 對應(yīng) F

stateX= dlode45(@(t,y,p)odeModel

(t,y,p,controlX),[0tspan],

stateX0,Params.neuralode);

% 下面的代碼利用計算的狀態(tài)序列和輸入的控制序列，通過構(gòu)造G(x,u|θg)來計算得到測量輸出Y

stateControl = cat(1,stateX,controlX);

…(此處省去多層連接)

Y = fullyconnect(stateControl,Param.fc1.Weights,);

end

function dxdt = odeModel(t,y,thetaODE,controlu)

% 這個函數(shù)就是構(gòu)造狀態(tài)方程 ?=F(x,u|θf) 中的右半部分 F(x,u|θf)

dxdt = cat(1,u,y);

dxdt = thetaODE.fc1.Weights*dxdt + thetaODE.fc1.Bias;

… （可以自定義層邏輯）

end

其中利用損失進(jìn)行梯度計算時，可以直接用反向傳播方法進(jìn)行自動微分，根據(jù) BP 算法原理，要保留在前向計算過程中模型中的激活值，尤其進(jìn)行 ode 積分求解時，積分時長 tspan 較長,或者積分步長較小導(dǎo)致積分步數(shù)過多,都會使網(wǎng)絡(luò)前向計算保留較多的激活值（用于反向傳播進(jìn)行自動微分），這就對內(nèi)存有了較高的要求。為了解決這個問題，計算梯度的時候，原文[8]使用了伴隨梯度算法（adjoint method）, 把梯度的計算變成了求解帶初值的另外一個常微分方程問題，只需要當(dāng)前激活作為初值，直接可以使用 ODE 求解器計算這個新的常微分方程得到各個時刻損失對網(wǎng)絡(luò)參數(shù)的梯度，而不需要保留前向過程中的激活。具體的細(xì)節(jié)可以看原論文[8]證明。具體使用哪種微分計算可以通過設(shè)置 dlode45 這個針對神經(jīng)網(wǎng)絡(luò)微分方程的 ODE 求解器的屬性來實現(xiàn)[9]。?

工程師更友好的神經(jīng)網(wǎng)絡(luò)狀態(tài)空間模型辨識

通過上述腳本實現(xiàn)模型的訓(xùn)練，可以得到一個訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)狀態(tài)空間模型。對于領(lǐng)域工程師，MATLAB 提供了一種更集成的方式實現(xiàn)上述過程，將上述的六個步驟集成在一個 nlssest 函數(shù)中，我們通過這種簡單易用的方式對飛行機(jī)器人被控對象模型進(jìn)行神經(jīng)網(wǎng)絡(luò)狀態(tài)空間模型的辨識。?

%創(chuàng)建一個 idNeuralStateSpace 對象，包含6個狀態(tài)量和4個控制輸入.

obj=idNeuralStateSpace(6,NumInputs=4);

%系統(tǒng)因為狀態(tài)量就是測量輸出，所以只需要創(chuàng)建一個狀態(tài)方程 ?=F(x,u|θf)，測量方程是個恒等式 (G(x,u|θg)=I,即 Y=I*X）,自定義狀態(tài)方程的網(wǎng)絡(luò)節(jié)點數(shù)和層數(shù).

obj.StateNetwork = createMLPNetwork(obj,'state',...

LayerSizes=[32 32 32],...

Activations='tanh',...

WeightsInitializer="glorot",...

BiasInitializer="zeros");

% 設(shè)置訓(xùn)練選項

options = nssTrainingOptions('adam');

% 使用nlssest利用實驗數(shù)據(jù)/仿真數(shù)據(jù)對構(gòu)建的進(jìn)行參數(shù)的訓(xùn)練

% 完成了上述六個步驟

nss=nlssest(U,X,obj,options,'UseLastExperimentForValidation',true);

訓(xùn)練完成得到狀態(tài)空間模型 nss 后，可以使用 generateMATLABFunction 自動為這個神經(jīng)網(wǎng)絡(luò)狀態(tài)空間模型生成狀態(tài)方程和雅可比函數(shù)解析式。

generateMATLABFunction(nss,'nssStateFcn');

function [A, B] = nssStateFcnJacobian(x,u)

% 這是自動生成的神經(jīng)網(wǎng)絡(luò)狀態(tài)模型的雅可比函數(shù)

? ? %# codegen

persistent StateNetwork

MATname = 'nssStateFcnData';

if isempty(StateNetwork)

StateNetwork = coder.load(MATname);

end

out = [x;u];

J = eye(length(out));

% 第一層隱含層

Jfc=StateNetwork.fc1.Weights;

out=StateNetwork.fc1.Weights*out+StateNetwork.fc1.Bias;

Jac = deep.internal.coder.jacobian.tanh(out);

…

% output layer

J = StateNetwork.output.Weights*J;

end

將訓(xùn)練的神經(jīng)網(wǎng)絡(luò)狀態(tài)空間模型作為MPC中的預(yù)測模型

我們把自動生成生成的狀態(tài)方程和雅可比函數(shù)賦給 MPC 的模型。

nlobj_tracking.Model.StateFcn =

'nssStateFcn';

nlobj_tracking.Jacobian.StateFcn=

'nssStateFcnJacobian';

使用含有神經(jīng)網(wǎng)絡(luò)狀態(tài)空間模型的 MPC 再次進(jìn)行閉環(huán)仿真，代碼跟上述使用第一原理的 MPC 一樣，查看追蹤效果。?

將計劃的閉環(huán)軌跡與實際的閉環(huán)軌跡進(jìn)行比較。響應(yīng)與基于第一性原理的預(yù)測模型結(jié)果相對接近。

總結(jié)

本文旨在提供一種將 AI 和模型預(yù)測控制結(jié)合的方法，可以將訓(xùn)練的神經(jīng)網(wǎng)絡(luò)狀態(tài)空間方程的模型給到模型預(yù)測控制器中使用，從而在第一原理模型不存在的情況下也可以實現(xiàn)預(yù)測控制。對于模型預(yù)測控制本身比較復(fù)雜，包括類型選擇和加速等等。

編輯：黃飛

?