MPC控制

MPC (basic)

系统模型

$$\begin{array}{rl}{x}_{k+1}=& A{x}_k+B{u}_k+E,\\ y_k=& C{x}_k\end{array}$$

令（预测时域$N_c$，控制时域$N_c$）

$$\begin{array}{r}{X}_k=\left[\begin{array}{c}{x}_{k+1}\\ x_{k+2}\\ ⋮\\ x_{k+N_p}\end{array}\right]\in {\mathbb{R}}^{n{N}_p}，Y_k=\left[\begin{array}{c}{y}_{k+1}\\ y_{k+2}\\ ⋮\\ y_{k+N_p}\end{array}\right],U_k=\left[\begin{array}{c}{u}_k\\ u_{k+1}\\ ⋮\\ u_{k+N_c-1}\end{array}\right]\in {\mathbb{R}}^{m{N}_c}\end{array}$$

则有

$$\begin{array}{rl}{X}_k=& \left[\begin{array}{c}A\\ ⋮\\ A^{N_p}\end{array}\right]x_k+\left[\begin{array}{cccc}B& 0& \cdots & 0\\ AB& B& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ A^{N_c-1}B& A^{N_c-2}B& \cdots & B\\ ⋮& ⋮& \ddots & ⋮\\ A^{N_p-1}B& A^{N_p-2}B& \cdots & A^{N_p-N_c}B\end{array}\right]U_k\\ & +\left[\begin{array}{cccc}I& 0& \cdots & 0\\ A& I& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ A^{N_p-1}& A^{N_p-2}& \cdots & I\end{array}\right](1_{N_p}\otimes E)\\ Y_k=& \left[\begin{array}{c}CA\\ ⋮\\ C{A}^{N_p}\end{array}\right]x_k+\left[\begin{array}{cccc}CB& 0& \cdots & 0\\ CAB& CB& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ C{A}^{N_c-1}B& C{A}^{N_c-2}B& \cdots & CB\\ ⋮& ⋮& \ddots & ⋮\\ C{A}^{N_p-1}B& C{A}^{N_p-2}B& \cdots & C{A}^{N_p-N_c}B\end{array}\right]U_k\\ & +\left[\begin{array}{cccc}C& 0& \cdots & 0\\ CA& C& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ C{A}^{N_p-1}& C{A}^{N_p-2}& \cdots & C\end{array}\right](1_{N_p}\otimes E)\\ =& \Phi x_k+\Theta U_k+\Psi \end{array}$$

上述状态方程和输出方程在$x_{k+N_c}$及以后都是近似的，这是由于控制时域长度小于预测时域导致的。成本函数，如下

$$\begin{array}{rl}cost=& \sum _{i=1}^{N_p}\Vert y_{k+i}-y_{k+i}^{\mathrm{ref}}{\Vert }_Q^2+\sum _{i=0}^{N_c-1}\Vert u_{k+i}{\Vert }_R^2\\ =& {\tilde{Y}}_k^{\top }(Q\otimes I_{N_p}){\tilde{Y}}_k+U_k(R\otimes I_{N_c})U_k\\ =& {\tilde{Y}}_k^{\top }\bar{Q}{\tilde{Y}}_k+U_k^{\top }\bar{R}{U}_k\\ =& (\Phi x_k+\Theta U_k-{\bar{Y}}_k^{\mathrm{ref}}{)}^{\top }\bar{Q}(\Phi x_k+\Theta U_k-{\bar{Y}}_k^{\mathrm{ref}})+U_k^{\top }\bar{R}{U}_k\\ =& (\Phi x_k+\Theta U_k{)}^{\top }\bar{Q}(\Phi x_k+\Theta U_k)+{\bar{Y}}_k^{\mathrm{ref}\top }\bar{Q}{\bar{Y}}_k^{\mathrm{ref}}-2(\Psi x_k+\Theta U_k{)}^{\top }\bar{Q}{\bar{Y}}_k^{\mathrm{ref}}+U_k^{\top }\bar{R}{U}_k\\ =& x_k^{\top }{\Phi }^{\top }\bar{Q}\Phi x_k+{\bar{Y}}_k^{\mathrm{ref}\top }\bar{Q}{\bar{Y}}_k^{\mathrm{ref}}-2{\bar{Y}}_k^{\mathrm{ref}\top }\bar{Q}\Phi x_k-2{\bar{Y}}_k^{\mathrm{ref}\top }\bar{Q}\Theta U_k+2x_k^{\top }{\Phi }^{\top }\bar{Q}\Theta U_k\\ & +U_k^{\top }(\bar{R}+{\Theta }^{\top }\bar{Q}\Theta )U_k\\ =& x_k^{\top }{\Phi }^{\top }\bar{Q}\Phi x_k+{\bar{Y}}_k^{\mathrm{ref}\top }\bar{Q}{\bar{Y}}_k^{\mathrm{ref}}-2{\bar{Y}}_k^{\mathrm{ref}\top }\bar{Q}M{x}_k+2(\Phi x_k-{\bar{Y}}_k^{\mathrm{ref}}{)}^{\top }\bar{Q}\Theta U_k\\ & +U_k^{\top }(\bar{R}+{\Theta }^{\top }\bar{Q}\Theta )U_k\\ =& c+2g{U}_k+U_k^{\top }H{U}_k\end{array}$$

其中

$$\begin{array}{rl}& Q\in {\mathbb{S}}_{\ge 0}^n,R\in {\mathbb{S}}_{>0}^m\\ & \bar{Q}=Q\otimes I_{N_p},\bar{R}=R\otimes I_{N_c}\\ & g=(\Phi x_k-{\bar{Y}}_k^{\mathrm{ref}}{)}^{\top }\bar{Q}\Theta ,H=\bar{R}+{\Theta }^{\top }\bar{Q}\Theta ,{\bar{Y}}_k^{\mathrm{ref}}=Y_k^{\mathrm{ref}}-\Psi \end{array}$$

可以对输入$u$或者对其增量施加约束

可以加入松弛因子，增强可解性，如下

$$\begin{array}{rl}cost=& \sum _{i=0}^{N_p}\Vert y_{k+i}-y_{k+i}^{\mathrm{ref}}{\Vert }_Q^2+\sum _{i=0}^{N_c}\Vert u_{k+i}{\Vert }_R^2+w{\lambda }^2\\ =& c+2\bar{g}{\bar{U}}_k+{\bar{U}}_k^{\top }\bar{H}{\bar{U}}_k\end{array}$$

其中

$$\begin{array}{r}\bar{g}=\left[\begin{array}{cc}g& 0\end{array}\right],\bar{H}=\left[\begin{array}{cc}H& 0\\ 0& w\end{array}\right],{\bar{U}}_k=\left[\begin{array}{c}{U}_k\\ \lambda \end{array}\right]\end{array}$$

MPC (incremental)

增量系统模型

$$\begin{array}{rl}{x}_{k+1}=& A{x}_k+B{u}_k+E\\ y_k=& C{x}_k\\ {\bar{x}}_k=& \left[\begin{array}{c}{x}_k\\ u_{k-1}\end{array}\right]\\ {\bar{x}}_{k+1}=& \left[\begin{array}{c}A{x}_k+B{u}_k+E\\ u_k\end{array}\right]=\left[\begin{array}{c}A{x}_k+B{u}_{k-1}+B\mathrm{d}{u}_k+E\\ u_{k-1}+\mathrm{d}{u}_k\end{array}\right]=\left[\begin{array}{cc}A& B\\ O& I\end{array}\right]{\bar{x}}_k+\left[\begin{array}{c}B\\ I\end{array}\right]\mathrm{d}{u}_k+\left[\begin{array}{c}E\\ O\end{array}\right]\\ {\bar{x}}_{k+1}=& \bar{A}{\bar{x}}_k+\bar{B}\mathrm{d}{u}_k+\bar{E}\end{array}$$

进一步，定义增量系统的输出

$$\begin{array}{r}{y}_k=C{x}_k=C\left[\begin{array}{c}I\\ O\end{array}\right]{\bar{x}}_k=\left[\begin{array}{c}C\\ O\end{array}\right]{\bar{x}}_k=\bar{C}{\bar{x}}_k\end{array}$$

运用前面的方法即可

应用于车辆横向控制

以后轴中心为参考点
$$\begin{array}{r}\left[\begin{array}{c}\dot{x}\\ \dot{y}\\ \dot{\varphi }\end{array}\right]=\left[\begin{array}{c}{v}_r\cos \varphi \\ v_r\sin \varphi \\ \frac{v_r}{L}\tan {\delta }_f\end{array}\right],\left[\begin{array}{c}y\\ \varphi \end{array}\right]=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ \varphi \end{array}\right]\\ \bar{x}={\left[\begin{array}{ccc}x& y& \varphi \end{array}\right]}^{\top },\bar{y}={\left[\begin{array}{cc}y& \varphi \end{array}\right]}^{\top },u={\delta }_f,{\bar{y}}_k=C{\bar{x}}_k\end{array}$$

一阶泰勒展开，有

$$\begin{array}{rl}\dot{\bar{x}}=& {\dot{\bar{x}}}_k+\left[\begin{array}{ccc}0& 0& -v_{rk}\sin {\varphi }_k\\ 0& 0& v_{rk}\cos {\varphi }_k\\ 0& 0& 0\end{array}\right](\bar{x}-{\bar{x}}_k)+\left[\begin{array}{c}0\\ 0\\ \frac{v_{rk}}{L{\cos }^2{\delta }_{fk}}\end{array}\right](u-({\varphi }_k+{\delta }_{fk})u_k)\\ =& \left[\begin{array}{ccc}0& 0& -v_{rk}\sin {\varphi }_k\\ 0& 0& v_{rk}\cos {\varphi }_k\\ 0& 0& 0\end{array}\right]\bar{x}+\left[\begin{array}{c}0\\ 0\\ \frac{v_{rk}}{L{\cos }^2{\delta }_{fk}}\end{array}\right]u+\left[\begin{array}{c}{v}_{rk}\cos {\varphi }_k+{\varphi }_k{v}_{rk}\sin {\varphi }_k\\ v_{rk}\sin {\varphi }_k-{\varphi }_k{v}_{rk}\cos {\varphi }_k\\ \frac{v_{rk}}{L}\tan {\delta }_{fk}-\frac{{\delta }_{fk}{v}_{rk}}{L{\cos }^2{\delta }_{fk}}\end{array}\right]\\ =& A\bar{x}+Bu+E\end{array}$$

以前轴中心为参考点
$$\begin{gathered} \left[\begin{array}{c}{\dot{x}}_f\\ {\dot{y}}_f\\ \dot{\varphi }\end{array}\right]=\left[\begin{array}{c}{v}_f\cos (\varphi +{\delta }_f)\\ v_f\sin (\varphi +{\delta }_f)\\ \frac{v_f\sin {\delta }_f}{L}\end{array}\right],\left[\begin{array}{c}y\\ \varphi \end{array}\right]=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ \varphi \end{array}\right] \\ \bar{x}={\left[\begin{array}{ccc}x& y& \varphi \end{array}\right]}^{\top },\bar{y}={\left[\begin{array}{cc}y& \varphi \end{array}\right]}^{\top },u=\delta ,{\bar{y}}_k=C{\bar{x}}_k \end{gathered}$$
一阶泰勒展开
$$\begin{array}{rl}\dot{\bar{x}}=& {\dot{\bar{x}}}_k+\left[\begin{array}{ccc}0& 0& -v_{fk}\sin ({\varphi }_k+{\delta }_{fk})\\ 0& 0& v_{fk}\cos ({\varphi }_k+{\delta }_{fk})\\ 0& 0& 0\end{array}\right](\bar{x}-{\bar{x}}_k)+\left[\begin{array}{c}-v_{fk}\sin ({\varphi }_k+{\delta }_{fk})\\ v_{fk}\cos ({\varphi }_k+{\delta }_{fk})\\ \frac{v_{fk}\cos {\delta }_{fk}}{L}\end{array}\right](u-u_k)\\ =& A\bar{x}+Bu+\left[\begin{array}{c}{v}_{fk}\cos ({\varphi }_k+{\delta }_{fk})+({\varphi }_k+{\delta }_{fk})v_{fk}\sin ({\varphi }_k+{\delta }_{fk})\\ v_{fk}\sin ({\varphi }_k+{\delta }_{fk})-({\varphi }_k+{\delta }_{fk})v_{fk}\cos ({\varphi }_k+{\delta }_{fk})\\ \frac{v_k\sin {\delta }_{fk}-{\delta }_{fk}{v}_{fk}\cos {\delta }_k}{L}\end{array}\right]\\ =& A\bar{x}+Bu+E\end{array}$$

以质心为参考点
$$\begin{gathered} \left[\begin{array}{c}{\dot{x}}_c\\ {\dot{y}}_c\\ \dot{\varphi }\end{array}\right]=\left[\begin{array}{c}{v}_c\cos (\varphi +\beta )\\ v_c\sin (\varphi +\beta )\\ \frac{v_c\cos \beta (\tan {\delta }_f-\tan {\delta }_r)}{L}\end{array}\right],\beta =\arctan \frac{l_f\tan {\delta }_r+l_r\tan {\delta }_f}{L} \\ \bar{x}={\left[\begin{array}{ccc}x& y& \varphi \end{array}\right]}^{\top },\bar{y}={\left[\begin{array}{cc}y& \varphi \end{array}\right]}^{\top },u={\left[\begin{array}{cc}{\delta }_f& {\delta }_r\end{array}\right]}^{\top },{\bar{y}}_k=C{\bar{x}}_k \end{gathered}$$
一阶泰勒展开
$$\begin{array}{rl}\dot{\bar{x}}=& {\dot{\bar{x}}}_k+\left[\begin{array}{ccc}0& 0& -v_{ck}\sin ({\varphi }_k+{\beta }_k)\\ 0& 0& v_{ck}\cos ({\varphi }_k+{\beta }_k)\\ 0& 0& 0\end{array}\right](\bar{x}-{\bar{x}}_k)\\ & +\left[\begin{array}{ccc}-v_{ck}\sin ({\varphi }_k+{\beta }_k){\beta }_k^{\prime f}& -v_{ck}\sin ({\varphi }_k+{\beta }_k){\beta }_k^{\prime r}& \\ v_{ck}\cos ({\varphi }_k+{\beta }_k){\beta }_k^{\prime f}& v_{ck}\cos ({\varphi }_k+{\beta }_k){\beta }_k^{\prime r}& \\ -\frac{v_{ck}({\beta }_k^{\prime f}\sin {\beta }_k(\tan {\delta }_{kf}-\tan {\delta }_{kr})+\cos {\beta }_k{\sec }^2{\delta }_{kf})}{L}& -\frac{v_{ck}({\beta }_k^{\prime r}\sin {\beta }_k(\tan {\delta }_{kf}-\tan {\delta }_{kr})-\cos {\beta }_k{\sec }^2{\delta }_{kr})}{L}& \end{array}\right](u-u_k)\\ =& A\bar{x}+Bu+\left[\begin{array}{c}{v}_{ck}\cos ({\varphi }_k+{\beta }_k)+({\varphi }_k+{\beta }_k^{\prime f}{\delta }_{fk}+{\beta }_k^{\prime r}{\delta }_{rk})v_{ck}\sin ({\varphi }_k+{\beta }_k)\\ v_{ck}\sin ({\varphi }_k+{\beta }_k)-({\varphi }_k+{\beta }_k^{\prime f}{\delta }_{fk}+{\beta }_k^{\prime r}{\delta }_{rk})v_{ck}\cos ({\varphi }_k+{\beta }_{fk})\\ v_{ck}\frac{(\cos {\beta }_k+({\beta }_k^{\prime f}{\delta }_{fk}+{\beta }_k^{\prime r}{\delta }_{rk})\sin {\beta }_k)(\tan {\delta }_{fk}-\tan {\delta }_{rk})+\cos {\beta }_k({\sec }^2{\delta }_{kf}{\delta }_{rk}-{\sec }^2{\delta }_{rk}{\delta }_{rk})}{L}\end{array}\right]\\ =& A\bar{x}+Bu+E\\ {\beta }_k^{\prime f}=& \frac{l_r{\sec }^2{\delta }_f}{L(1+{\beta }_0^2)}\\ {\beta }_k^{\prime r}=& \frac{l_f{\sec }^2{\delta }_r}{L(1/+{\beta }_0^2)}\\ {\beta }_0=& \frac{l_f\tan {\delta }_r+l_r\tan {\delta }_f}{L}\end{array}$$

注释：实际应用中可用${\delta }_{k-1}\approx {\delta }_k$，但由于转角一般很小，一般直接将其忽略。同时对于车体作弊坐标系有${\bar{x}}_k=0$。故模型可以被简化为：
以后轴中心为参考点

$$\begin{array}{rl}\dot{\bar{x}}=& \left[\begin{array}{ccc}0& 0& 0\\ 0& 0& v_{rk}\\ 0& 0& 0\end{array}\right]\bar{x}+\left[\begin{array}{c}0\\ 0\\ \frac{v_{rk}}{L{\cos }^2{\delta }_{fk}}\end{array}\right]u+\left[\begin{array}{c}{v}_{rk}\\ 0\\ \frac{v_{rk}}{L}\tan {\delta }_{fk}-\frac{{\delta }_{fk}{v}_{rk}}{L{\cos }^2{\delta }_{fk}}\end{array}\right]\\ =& A\bar{x}+Bu+E\end{array}$$

以前轴中心为参考点

$$\begin{array}{rl}\dot{\bar{x}}=& \left[\begin{array}{ccc}0& 0& -v_{fk}\sin {\delta }_{fk}\\ 0& 0& v_{fk}\cos {\delta }_{fk}\\ 0& 0& 0\end{array}\right]\bar{x}+\left[\begin{array}{c}-v_{fk}\sin {\delta }_{fk}\\ v_{fk}\cos {\delta }_{fk}\\ \frac{v_{fk}\cos {\delta }_{fk}}{L}\end{array}\right]u+\left[\begin{array}{c}{v}_{fk}\cos {\delta }_{fk}+{\delta }_{fk}{v}_{fk}\sin {\delta }_{fk}\\ v_{fk}\sin {\delta }_{fk}-{\delta }_{fk}{v}_{fk}\cos {\delta }_{fk}\\ \frac{v_k\sin {\delta }_{fk}-{\delta }_{fk}{v}_{fk}\cos {\delta }_k}{L}\end{array}\right]\\ =& A\bar{x}+Bu+E\end{array}$$

以质心为参考点

$$\begin{array}{rl}\dot{\bar{x}}=& \left[\begin{array}{ccc}0& 0& -v_{ck}\sin {\beta }_k\\ 0& 0& v_{ck}\cos {\beta }_k\\ 0& 0& 0\end{array}\right]\bar{x}\\ & +\left[\begin{array}{ccc}-v_{ck}{\beta }_k^{\prime f}\sin {\beta }_k& -v_{ck}{\beta }_k^{\prime r}\sin {\beta }_k& \\ v_{ck}{\beta }_k^{\prime f}\cos {\beta }_k& v_{ck}{\beta }_k^{\prime r}\cos {\beta }_k& \\ -\frac{v_{ck}({\beta }_k^{\prime f}\sin {\beta }_k(\tan {\delta }_{kf}-\tan {\delta }_{kr})+\cos {\beta }_k{\sec }^2{\delta }_{kf})}{L}& -\frac{v_{ck}({\beta }_k^{\prime r}\sin {\beta }_k(\tan {\delta }_{kf}-\tan {\delta }_{kr})-\cos {\beta }_k{\sec }^2{\delta }_{kr})}{L}& \end{array}\right]u\\ & +\left[\begin{array}{c}{v}_{ck}\cos {\beta }_k+({\beta }_k^{\prime f}{\delta }_{fk}+{\beta }_k^{\prime r}{\delta }_{rk})v_{ck}\sin {\beta }_k\\ v_{ck}\sin {\beta }_k-({\beta }_k^{\prime f}{\delta }_{fk}+{\beta }_k^{\prime r}{\delta }_{rk})v_{ck}\cos {\beta }_{fk}\\ v_{ck}\frac{(\cos {\beta }_k+({\beta }_k^{\prime f}{\delta }_{fk}+{\beta }_k^{\prime r}{\delta }_{rk})\sin {\beta }_k)(\tan {\delta }_{fk}-\tan {\delta }_{rk})+\cos {\beta }_k({\sec }^2{\delta }_{kf}{\delta }_{rk}-{\sec }^2{\delta }_{rk}{\delta }_{rk})}{L}\end{array}\right]\\ =& A\bar{x}+Bu+E\end{array}$$

注释：笔者认为在高速场景下，前后轮转角对横向控制的影响较大，因此这里简化中没有直接将前后轮转角近似为0，但到底影响有多大需要测试
低速时，再次近似
以后轴中心为参考点

$$\begin{array}{rl}\dot{\bar{x}}=& \left[\begin{array}{ccc}0& 0& 0\\ 0& 0& v_{rk}\\ 0& 0& 0\end{array}\right]\bar{x}+\left[\begin{array}{c}0\\ 0\\ \frac{v_{rk}}{L}\end{array}\right]u+\left[\begin{array}{c}{v}_{rk}\\ 0\\ 0\end{array}\right]\\ =& A\bar{x}+Bu+E\end{array}$$

以前轴中心为参考点

$$\begin{array}{rl}\dot{\bar{x}}=& \left[\begin{array}{ccc}0& 0& 0\\ 0& 0& v_{fk}\\ 0& 0& 0\end{array}\right]\bar{x}+\left[\begin{array}{c}0\\ v_{fk}\\ \frac{v_{fk}}{L}\end{array}\right]u+\left[\begin{array}{c}{v}_{fk}\\ 0\\ 0\end{array}\right]\\ =& A\bar{x}+Bu+E\end{array}$$

以质心为参考点

$$\begin{array}{rl}\dot{\bar{x}}=& \left[\begin{array}{ccc}0& 0& 0\\ 0& 0& v_{ck}\\ 0& 0& 0\end{array}\right]\bar{x}+\left[\begin{array}{ccc}0& 0& \\ v_{ck}\frac{l_r}{L}& v_{ck}\frac{l_f}{L}& \\ -\frac{v_{ck}}{L}& \frac{v_{ck}}{L}& \end{array}\right]u+\left[\begin{array}{c}{v}_{ck}\\ 0\\ 0\end{array}\right]\\ =& A\bar{x}+Bu+E\\ {\beta }_k^{\prime f}=& \frac{l_r}{L},{\beta }_k^{\prime r}=\frac{l_f}{L},{\beta }_0=1\end{array}$$

差分离散

$$\begin{array}{r}{\bar{x}}_{k+1}=(I+\mathrm{d}tA){\bar{x}}_k+\mathrm{d}tB{u}_k+\mathrm{d}tE=\bar{A}{\bar{x}}_k+\bar{B}{u}_k+\bar{E},\end{array}$$

中点欧拉离散

$$\begin{array}{r}{y}_{n+1}=y(t+\mathrm{d}t)=y_k+\mathrm{d}tf(t+\frac{1}{2}\mathrm{d}t,{\bar{y}}_k+\frac{\mathrm{d}t}{2}f(t,{\bar{y}}_k))=y_k+\mathrm{d}tf(t+\frac{1}{2}\mathrm{d}t,\frac{1}{2}({\bar{y}}_k+{\bar{y}}_{k+1}))\end{array}$$

用$t+0.5\mathrm{d}t$的导数值来代替整个周期的导数（平均值）

$$\begin{array}{rl}{\bar{x}}_{k+1}=& {\bar{x}}_k+\mathrm{d}t(\frac{1}{2}A({\bar{x}}_{k+1}+{\bar{x}}_k)+B{u}_k+E)\\ =& (I-\frac{\mathrm{d}t}{2}A{)}^{-1}\left((I+\frac{\mathrm{d}t}{2}A){\bar{x}}_k+\mathrm{d}tB{u}_k+\mathrm{d}tE\right)\end{array}$$

4阶龙格-库塔

$$\begin{array}{r}\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{\mathrm{d}t}{6}(k_1+2k_2+2k_3+k_4)\\ k_1& =f(t_n,y_n)\\ k_2& =f(t_n+\frac{\mathrm{d}t}{2},\frac{\mathrm{d}t}{2}{k}_1)\\ k_3& =f(t_n+\frac{\mathrm{d}t}{2},\frac{\mathrm{d}t}{2}{k}_2)\\ k_4& =f(t_n+\mathrm{d}t,\mathrm{d}t{k}_3)\end{array}\end{array}$$

4阶龙格-库塔法即近似值与真值之间的误差是四阶导数的高阶无穷小（欧拉中值法为二阶）。计算复杂，略。

【注释】无论中点欧拉还是龙格-库塔，都需要带入半采样周期，并依赖其预估采样周期内某点的导数。但是对于扰动或者控制输入，半周期处的值是不可知的（$u_{k+\frac{1}{2}}$），因此估计出来的值不一定是准确的，因此，非状态自治系统采用这两种方式得到离散模型的精度是否更高？高多少？都是未可知的。

应用上面的额MPC方法即可。