第一部分

模型預測控制

01

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

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

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

對于最小燃料消耗問題（接下來本文中主要使用的例子，假設燃料消耗與控制量 u 成正比，并且 u∈[0 1]），性能指標變為：

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

性能指標是二次型的無限時域的連續時間泛函：

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

其中系統輸出是車輛位移 x 和倒立擺角度 θ，控制輸出 u 是作用在車上的水平力。控制是標量，狀態是 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 函數就可以求得最優增益矩陣 K。但對于更一般的非線性系統，使用極小值原理或動態規劃求解解析解幾乎不可能，數值解計算復雜度較高，所以可以考慮一定的問題簡化，僅考慮未來的一小段時間最優（而非全時間域“最優”），將最優控制問題轉化未來幾個時間步的在線數值優化問題，例如 MPC，可以大大降低“最優”控制計算的復雜度。

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

算法介紹

（圖表1）

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

代價函數: 用于度量控制器性能，目標是求它的最小值。

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

控制與決策：得到優化的控制量（MV）

代價函數

一個典型的代價函數由四部分組成，每一部分關注控制器性能的一個方面。

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

Jy(zk) 量化系統輸出與目標軌跡的跟蹤效果

Ju(zk) 用于量化控制輸入與目標控制變量跟蹤效果（在很多應用中，控制器還需要保證控制變量（MV）保持在某個目標附近，尤其在控制量數目遠多于系統輸出量數目的情況下）

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

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

其中，底標 i 是對 p 個預測步的循環，j 是對 ny 個輸出的循環。

k:當前控制周期

p:預測區間（Prediction Horizon）

ny:被控對象輸出變量的數量

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

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

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

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

約束

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?個）控制變量的下界和上界

預測與決策

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

（圖表2）

上述預測過程可以通過簡化的狀態空間離散預測模型表示：

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

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

整理成矩陣形式：

其中（此處不詳細列出這些矩陣表達式，可以根據上式遞歸自行推導）

組合各個元素得到優化問題定義

將上述的代價函數，約束，預測模型矩陣結合起來構建 MPC 的開環優化問題如下：

狀態變量 x(k)∈Rnx，控制變量 u(k)∈Rnu，滿足上述系統動力學方程以及時域約束，通過預測計算，求解如下最小化性能指標對應的優化變量zk

其中

即在預測區間 p 上的待優化的控制輸入序列。

示例：

非線性 MPC 的飛行機器人軌跡開環規劃和閉環反饋控制?

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

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

該模型動力學系統如下:

function dxdt = FlyingRobotStateFcn(x, u)

% 狀態方程：x 是六個狀態，u是推力控制。

%x(1)???– 質心x坐標

%x(2)??– 質心y坐標

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

%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

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

軌跡規劃

前面提到了，MPC 用于軌跡規劃其實是求解了一次開環的帶約束的優化問題，區別于后面的針對規劃好的軌跡的路徑跟蹤（利用反饋狀態的滾動優化）。

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

在此例中，對于規劃軌跡，我們開環優化需要使用全部時間長度：12s。設定一個采樣時間 Ts=0.4s，所以對應的預測區間為 30 步（12/0.4）。規劃的運行周期通常可以比反饋控制周期更大，或更低頻，所以有更多計算資源允許計算一個相對較大的優化問題。對于軌跡規劃我們可以創建一個多級非線性MPC對象，它允許對每個預測步都定義不同的代價函數和約束。

% 具有 nx=6 個狀態和 nu=4 個控制輸入（mv）

% 預測區間 p = 30

nlobj=nlmpcMultistage(p,nx,nu);

nlobj.Ts=Ts;

for?ct=1:p

%并且每個 stage/ 預測步有自己的代價函數

其中代價函數

包含 stage 這個參數，

見下面 FlyingRobotCostFcn 函數

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

end

%?代價函數示例

function J=FlyingRobotCostFcn(stage, x, u)

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

%和。本示例沒有直接使用 stage 參數。

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

end

指定預測模型的狀態轉移方程和對應的雅可比函數解析形式，可以大幅提升計算效率

nlobj.Model.StateFcn="FlyingRobotStateFcn";

nlobj.Model.StateJacFcn=@

FlyingRobotStateJacobianFcn;

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

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

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

每個推力的工作范圍在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];% 機器人停在 [-10, -10], 方向角為 pi/2

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

通過調用 nlmpcmove 命令，就可以完成一次優化，找到最佳狀態和控制（MV）軌跡。

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

圖表 3規劃（優化求解）的最優軌跡對應六個狀態量的預測值

圖表4對應的4個推力在最優軌跡規劃求解時的最優解序列（MV）

圖表5 機器人最優軌跡的坐標和方位信息，從?[-10 -10 pi/2]?到?[0 0 0].

軌跡跟蹤的反饋控制

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

圖表 6

設計MPC控制器?

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

ny = 6;

nlobj_tracking = nlmpc(nx,ny,nu);

% 對于狀態方程和雅可比函數，我們和上面用于軌跡規劃的MPC中使用同樣的動力學模型

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

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

% 指定滾動優化的預測區間 Prediction Horizon (不需要考慮太遠的未來時間)和控制區間 Control

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

nlobj_tracking.PredictionHorizon = 10;

nlobj_tracking.ControlHorizon = 4;

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

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

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

此外，在控制過程中，u1 和 u2 上下限的約束設定在 [0,1] 范圍內，因為兩個是同側相反的方向，所以不需要同時工作，因此加一個等式約束，使 u(1)*u(2) 對所有預測區間都為 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)];

▼

設計狀態估計器?

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

通過解 nx 個非線性代數方程，對連續模型進行轉換，得到離散狀態方程模型FlyingRobotStateFcnDiscreteTime.m。

其中m中的梯形規則如下實現：

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

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

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

xk1 = fsolve(FUN,xk1,Options);

% 梯形法則函數

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

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

fk1=ffun(xk1,uk);

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

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

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

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

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

% 創建EKF,并設定測量噪聲

EKF=extendedKalmanFilter(DStateFcn,DMeasFcn,x0);

EKF.MeasurementNoise = 0.01;

▼

閉環仿真?

現在圖表6中 MPC Controller, Plant, State Estimator 都有了，可以進行追蹤控制的閉環仿真

設置好初始條件，仿真 32 個時間步，參考軌跡就是在規劃階段計算出來的最優狀態軌跡。使用 nlmpcmove 和 nlmpcmoveopt 命令進行閉環仿真。

Tsteps = 32;

for?k = 1:Tsteps

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

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

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

xk = correct(EKF, yk);

%nlmpcmove 使用參考軌跡 Xref 進行一次 p 個未來步長的優化求解，得到當前時刻的優化控制uk.

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

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

%預測下一時刻的狀態和協方差矩陣，用于下步correct.

predict(EKF,uk,Ts);

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

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

lastMV = uk;

%基于當前的狀態 xk, 將求解的最優的uk作用在被控對象上（連續時間 ODE）并求解 ODE 得到下一時刻的狀態.

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

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

%將狀態結果保留下來用于下個循環進行測量.

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

end

將規劃的參考軌跡和實際的閉環控制軌跡進行比較。

?

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

第二部分

為被控對象辨識神經網絡狀態

空間模型并用于 MPC

02

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

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

神經網絡 ODE 算法與訓練

Neural ODE 思想介紹

神經網絡狀態空間模型是基于神經網絡 ODE 實現的，這里可以簡單介紹一下實現思路(在前面的文章 數據驅動的動態系統建模-深度學習中簡單介紹了訓練的思想。在視頻[7]

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

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

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

然后計算測量輸出方程

這樣就可以計算損失函數，接下來可以按神經網絡的訓練步驟對參數?[θf，θg]?進行訓練。

Neural ODE 訓練過程實現

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

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

在每次迭代 (each minibatch), 進行預處理.

將輸入數據傳給網絡進行前饋計算?(inference).

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

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

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

%遍歷輪數

for?epoch= 1:numEpochs

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

for i= 1:numIterperEpoch ? ?

iter=iter + 1;

????????%創建每次迭代用于訓練模型參數的 minibatch?

????????數據集，包含控制序列和輸出序列（真值）

[dlU,dlY,timeIds]=

createMiniBatch(neuralOdeTimesteps,?

UtrainCell,YtrainCell);

…

%@modelGradients 函數為了計算損失對模型

參數的梯度，總共完成了三個任務：

· 將輸入數據傳給網絡進行前饋計算(inference).

· 和真值進行比較計算loss

(i.e. compute loss function).

· 計算梯度

[grads,loss]=dlfeval(@modelGradients,

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

parameters,tspan);

%利用計算得到的梯度使用優化求解器 adam?

更新模型權重和偏置

[parameters,averageGrad2,averageSqGrad2]=

adamupdate(parameters,grads,averageGrad,

averageSqGrad,iter,learnRate,gradDecay,

sqGradDecay);

end

end

%modelGradient 函數完成三個任務

function[grad,loss]=modelGradients(dlTarget,

dlX0, dlU, params, tspan)

%neuralStateMode l 函數（接下來列出來的函數）

就是利用初始狀態和控制序列進行網絡推斷計算。

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

%利用模型推斷結果和真值進行損失計算

loss= huber(dlY, dlTarget);

%利用損失進行梯度計算

grad= dlgradient(loss, params);

end

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

%odeModel 函數對應著狀態方程中的F(x,u|θf)，基于初始狀態，和當前神經網絡參數下的狀態方程，利用 dlode45 進行積分求解,得到對應的模型推斷的狀態序列 stateX

%odeModel 對應 F

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

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

stateX0,Params.neuralode);

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

stateControl = cat(1,stateX,controlX);

…(此處省去多層連接)

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

end

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

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

dxdt = cat(1,u,y);

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

… （可以自定義層邏輯）

end

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

工程師更友好的神經網絡狀態空間模型辨識

通過上述腳本實現模型的訓練，可以得到一個訓練好的神經網絡狀態空間模型。對于領域工程師，MATLAB 提供了一種更集成的方式實現上述過程，將上述的六個步驟集成在一個 nlssest 函數中，我們通過這種簡單易用的方式對飛行機器人被控對象模型進行神經網絡狀態空間模型的辨識。?

%創建一個 idNeuralStateSpace 對象，包含6個狀態量和4個控制輸入.

obj=idNeuralStateSpace(6,NumInputs=4);

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

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

LayerSizes=[32 32 32],...

Activations='tanh',...

WeightsInitializer="glorot",...

BiasInitializer="zeros");

% 設置訓練選項

options = nssTrainingOptions('adam');

% 使用nlssest利用實驗數據/仿真數據對構建的進行參數的訓練

% 完成了上述六個步驟

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

訓練完成得到狀態空間模型 nss 后，可以使用 generateMATLABFunction 自動為這個神經網絡狀態空間模型生成狀態方程和雅可比函數解析式。

generateMATLABFunction(nss,'nssStateFcn');

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

% 這是自動生成的神經網絡狀態模型的雅可比函數

? ? %# 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

將訓練的神經網絡狀態空間模型作為MPC中的預測模型

我們把自動生成生成的狀態方程和雅可比函數賦給 MPC 的模型。

nlobj_tracking.Model.StateFcn =

'nssStateFcn';

nlobj_tracking.Jacobian.StateFcn=

'nssStateFcnJacobian';

使用含有神經網絡狀態空間模型的 MPC 再次進行閉環仿真，代碼跟上述使用第一原理的 MPC 一樣，查看追蹤效果。?

將計劃的閉環軌跡與實際的閉環軌跡進行比較。響應與基于第一性原理的預測模型結果相對接近。

總結

本文旨在提供一種將 AI 和模型預測控制結合的方法，可以將訓練的神經網絡狀態空間方程的模型給到模型預測控制器中使用，從而在第一原理模型不存在的情況下也可以實現預測控制。對于模型預測控制本身比較復雜，包括類型選擇和加速等等。

編輯：黃飛

?