无人机MPC模型预测控制

【本文来自转载】：链接

文章目录

• 前言
• 1. 背景 Background
• • 1.1 滚动时域控制 Receding Horizon Control
  • 1.2 线性模型预测控制 Linear Model Predictive Control
  • 1.3 非线性模型预测控制 Nonlinear Model Predictive Control
  • 1.4 线性鲁棒模型预测控制 Linear Robust Model Predictive Control
  • • 1.4.1 反馈预测 Feedback Predictions
    • 1.4.2 鲁棒不确定不等式 Robust Uncertain Inequalities Satisfaction
• 2. 用于多旋翼系统的基于模型的轨迹跟踪控制器 Model-Based Trajectory Tracking Controller for Multi-rotor System
• • 2.1 多旋翼系统模型 Multirotor System Model
  • 2.2 Linear MPC
  • • 3.2.1 姿态环参数辨识 Attitude Loop Parameters Identification
    • 3.2.2 线性化、解耦和离散化 Linearization, Decoupling and Discretization
    • 2.2.3 ROS集成 ROS Integration
    • 2.2.4 实验结果 Experimental Results
  • 2.3 Nonlinear MPC
  • • 2.3.1 ROS集成 ROS Integration
    • 2.3.2 实验结果 Experimental Results
• 3 结论 Conclusion
• 原文献

前言

本章遵循的一般控制结构是一种级联控制方法，其中一个可靠的适用于系统的低级控制器作为内环，基于模型的轨迹跟踪控制器作为外环运行。

1. 背景 Background

1.1 滚动时域控制 Receding Horizon Control

滚动时域控制(Receding Horizon Control, RHC)对应着控制理论的历史性发展，该理论旨在解决固定时域控制（fixed horizon control）的已知挑战。固定时域优化在时域的目的是计算未来N NN步的控制输入{ u 0 , u 1 , . . . , u N − 1 } {u_0, u_1, …, u_{N-1}}{u0,u1,…,uN−1}。但是他有两个缺点：

• (a) 当控制设计阶段出现意外(未知)扰动，或控制模型与实际系统行为不同时，控制器无法对计算出的控制序列进行解释；
• (b) 在控制序列的最后一步时，由于剩下的时间较少，不能对目标函数进行求解。

为了解决这些限制，RHC提出了计算整个控制序列后，只应用它第一步的结果，然后再次迭代整个过程(滚动时域的方式)求解新的控制序列。RHC策略一般适用于非线性动力学(考虑到状态更新可用):
x ˙ = f ( x , u ) ( 1 ) \dot{x}=f(x,u) \quad \quad \quad (1)x˙=f(x,u)(1)
其中向量f : R n × R m f:\mathbb{R}^n \times \mathbb{R}^mf:Rn×Rm，x ∈ R n × 1 x \in \mathbb{R}^{n\times1}x∈Rn×1表示状态向量，u ∈ R m × 1 u \in \mathbb{R}^{m\times1}u∈Rm×1表示输入向量。通常，RHC的基于状态反馈的优化问题可以定义如下。
min ⁡ z   F ( x t + N ) + ∑ k = 0 N − 1 ∥ x t + k − x t + k r e f ∥ ℓ + ∥ u t + k ∥ ℓ s . t . x t + k + 1 = f ( x t + k , u t + k ) u t + k ∈ U C x t + k ∈ X C x t = x ( t ) ( 2 )

minzF(xt+N)+∑k=0N−1∥∥xt+k−xreft+k∥∥ℓ+∥ut+k∥ℓs.t.min��(��+�)+∑�=0�−1‖��+�−��+����‖ℓ+‖��+�‖ℓ�.�.

\

xt+k+1=f(xt+k,ut+k)ut+k∈UCxt+k∈XCxt=x(t)��+�+1=�(��+�,��+�)��+�∈����+�∈����=�(�)

\quad \quad \quad (2)zminF(xt+N)+k=0∑N−1∥∥∥xt+k−xt+kref∥∥∥ℓ+∥ut+k∥ℓs.t.xt+k+1=f(xt+k,ut+k)ut+k∈UCxt+k∈XCxt=x(t)(2)
其中z = { u t , u t + 1 , . . . . . . u t + N − 1 } z={u_t, u_{t+1}, … u_{t+N-1}}z={ut​,ut+1​,…ut+N−1​}是优化变量，ℓ \ellℓ表示用于每级的权重(惩罚的度量)，F ( x t + N ) F(x_{t+N})F(xt+N​)表示最终状态的权重，x t + k r e f x_{t+k}^{ref}xt+kref​表示参考信号，下标t + k t+kt+k表示一个信号在当前时间t tt之后k kk步处的采样(使用固定的采样时间T s T_sTs​)，而t + k + 1 t+k+1t+k+1则是其下一步的变化，U C \mathfrak{U}CUC​表示输入约束集，X C \mathfrak{X}CXC​表示状态约束集，x ( t ) x(t)x(t)是当前RHC迭代开始时的状态向量的值。这个优化问题的解也是一个最优控制序列{ u t ∗ , u t + 1 ∗ , . . . . . . u t + N − 1 ∗ } {u_t^*, u{t+1}^*, … u{t+N-1}^*}{ut∗​,ut+1∗​,…ut+N−1∗​}，只使用第一步的控制量x t ∗ x_t^*xt∗​作为当前的输入，第二步的输入需要再次迭代后获取（白话：每次求解都可以获得一系列包含未来N NN步的最优控制输入，但是实际使用的时候只取第一个，后面的舍弃）。

在公式1中，F ( x t + N ) F(x_{t+N})F(xt+N)在闭环稳定中起到了至关重要的作用。它迫使系统状态在预测期结束时在一个特定的集合内取值，可以使用李亚普诺夫稳定性来分析、证明每次局部迭代的稳定性。例如对于一个监管问题（x t + k r e f = 0 f o r k = 0 , . . . , N − 1 x_{t+k}^{ref}=0 \quad for \quad k=0, …, N-1xt+kref=0fork=0,…,N−1），和一个“衰减”因子ℓ \ellℓ。上述优化问题的解使得系统在x t = 0 x_t=0xt=0, u t = 0 u_t=0ut=0处稳定，也就意味着系统的终端约束 （terminal constraint）为 x t + k = 0 x_{t+k}=0xt+k=0，如图1所示。然而，全局的稳定性是没有保证的。为此，必须考虑引入terminal cost和terminal constraint。 但是，仅添加terminal constraints可能不可行，而且一般添加约束的优化问题可能会很难解决。 在许多实际情况下，在控制设计过程中不强制执行terminal constraints，而是验证了后验（如果不满足的话，则可以增加预测范围）。

图1 llustration of the terminal constraint set(Ω)

此外，RHC的最大挑战就是递归可行性。尽管从理论角度绝对推荐，但由于理论或实践的影响，并不能确保可以构建预先保证递归可行性的RHC。一般来说，RHC策略缺乏递归可行性，因此即使有可能找到可行的状态，但最优控制动作将状态向量移动到RHC优化问题不可行的点。虽然一般的可行性分析方法是非常具有挑战性的，但对于特定的情况可能存在强大的工具可以用来求解。例如对于线性系统，可以使用Farkas的引理与双层编程（bilevel programming）结合来搜索缺乏递归可行性的有问题的初始状态。

1.2 线性模型预测控制 Linear Model Predictive Control

在本小节中，我们简要介绍了线性系统MPC背后的理论。 给出了输入变量和状态变量均为线性约束的线性系统的最优控制问题。此外，我们还讨论了线性和二次成本函数情况下的控制输入特性，稳定性和可行性。 为了在考虑干扰的情况下实现无偏移跟踪，我们在系统模型中增加了额外的干扰状态d ( t ) d(t)d(t)来表示模型的误差。 利用观测器（observer）对稳态下的扰动进行估计。本小节将简要讨论观测器的设计和干扰模型。

min ⁡ J 0 ( x 0 , U , X r e f , U r e f ) subject to x k + 1 = A x k + B u k + B d d k ; d k + 1 = d k , k = 0 , … , N − 1 x k ∈ X C , u k ∈ U C x N ∈ X C N x 0 = x ( t 0 ) , d 0 = d ( t 0 ) ( 3 )

minsubject to \quadJ0(x0,U,Xref,Uref)xk+1=Axk+Buk+Bddk;dk+1=dk,k=0,…,N−1xk∈XC,uk∈UCxN∈XCNx0=x(t0),d0=d(t0)min�0(�0,�,����,����)subject to \quad��+1=���+���+����;��+1=��,�=0,…,�−1��∈��,��∈����∈����0=�(�0),�0=�(�0)

\quad \quad \quad (3) minsubject to J0(x0,U,Xref,Uref)xk+1=Axk+Buk+Bddk;dk+1=dk,k=0,…,N−1xk∈XC,uk∈UCxN∈XCNx0=x(t0),d0=d(t0)(3)

在线性状态和输入约束下实现无偏置状态跟踪的最优控制问题如公式3所示，其中J 0 J_0J0是代价函数，X r e f = { x 0 r e f , ⋯   , x N r e f } X^{ref}={x_0{ref},\cdots,x_N^{ref}}Xref={x0ref,⋯,xNref}是参考状态序列，U r e f = { u 0 r e f , ⋯   , u N − 1 } U^{ref}={u_0{ref},\cdots,u_{N-1}}Uref={u0ref,⋯,uN−1}和U r e f = { u 0 r e f , ⋯   , u N − 1 r e f } U^{ref}={u_0{ref},\cdots,u_{N-1}^{ref}}Uref={u0ref,⋯,uN−1ref}分别是控制输入序列和稳态输入序列，B d B_dBd是扰动模型，d k d_kdk是外部干扰，X C 、 U C \mathfrak{X}{C}、\mathfrak{U}{C}XC、UC和X C N \mathfrak{X}_{CN}XCN是正多面体（polyhedra）。扰动模型的选择不是一项简单的任务，它取决于所考虑的系统和预期扰动的类型。优化问题定义为

J 0 ( x 0 , U , X r e f , U r e f ) = ∑ k = 0 N − 1 ( ( x k − x k r e f ) T Q x ( x k − x k r e f ) + ( u k − u k r e f ) T R u ( u k − u k r e f ) + ( u k − u k − 1 ) T R Δ ( u k − u k − 1 ) ) + ( x N − x N r e f ) T P ( x N − x N r e f ) ( 4 )

J0(x0,U,Xref,Uref)=∑k=0N−1⎛⎝⎜⎜⎜(xk−xrefk)TQx(xk−xrefk)+(uk−urefk)TRu(uk−urefk)+(uk−uk−1)TRΔ(uk−uk−1)⎞⎠⎟⎟⎟+(xN−xrefN)TP(xN−xrefN)�0(�0,�,����,����)=∑�=0�−1((��−�����)���(��−�����)+(��−�����)���(��−�����)+(��−��−1)��Δ(��−��−1))+(��−�����)��(��−�����)

\quad \quad \quad (4) J0(x0,U,Xref,Uref)=k=0∑N−1⎝⎜⎜⎜⎛(xk−xkref)TQx(xk−xkref)+(uk−ukref)TRu(uk−ukref)+(uk−uk−1)TRΔ(uk−uk−1)⎠⎟⎟⎟⎞+(xN−xNref)TP(xN−xNref)(4)

其中Q x ⪰ 0 Q_x \succeq 0Qx⪰0是对状态错误的惩罚，R u ⪰ 0 R_u \succeq 0Ru⪰0是对控制输入的惩罚，R Δ ⪰ 0 R_{\Delta} \succeq 0RΔ⪰0是对控制变化率的惩罚，P PP是对终端状态误差的惩罚。

一般而言，除了一些特殊情况而言，其稳定性和可行性难以保证，例如二次调节器的无限时域控制问题(Linear Quadratic Regulator, LQR)。当预测的步长为N NN步时，其稳定性和可行性保证是值得商榷的。理论上，随和预测步长的增加，有利于提高控制器的稳定性和可行性，但计算量会随之增加，而对于空中机器人应用来说，需要在计算能力有限的平台上进行计算，并对机器人进行快速控制。可以通过选择终端代价P PP和终端约束X C N \mathfrak{X}_{CN}XCN来确保闭环的稳定性和可能性。在本章中，我们更多的关注终端权值P的选择，因为它比较容易计算，而终端约束通常比较困难，其在有足够长的预测范围的情况下可以实现实际的稳定性。

注意，为了保证输入的平滑、避免不必要的震荡，在成本函数(4)中控制输入率Δ u k \Delta u_kΔuk。在成本函数中
(4)，u − 1 u_{−1}u−1​为上一个时间步长应用于系统的实际控制输入。

如前所述，通过对系统模型加入扰动d k d_kdk来捕获建模误差，可以实现无偏置参考跟踪（offset-free reference tracking）。假设我们想跟踪系统输出y k = C x k y_k = Cx_kyk=Cxk，实现稳态偏移无偏置追踪y ∞ = r ∞ y_∞= r_∞y∞=r∞。可以用下面这个简单的观察者来估计这种干扰

[ x ^ k + 1 d ^ k + 1 ] = [ A B d 0 I ] [ x ^ k d ^ k ] + [ B 0 ] u k + [ L x L d ] ( C x ^ k − y m , k ) ( 5 ) \left[

x^k+1dk+1�^�+1��+1

\right]=\left[

A0BdI���0�

\right]\left[

x^kdk�^���

\right]+\left[

B0�0

\right]{{u}_{k}}+\left[

LxLd����

\right](C{{\hat{x}}{k}}-{{y}{m,k}}) \quad \quad \quad (5) [x^k+1dk+1]=[A0BdI][x^kdk]+[B0]uk+LxLd(5)

其中x ^ \hat{x}x^和d ^ \hat{d}d^分别是状态和外部干扰的预估，L x L_xLx和L d L_dLd是观测器增益，y m , k y_{m,k}ym,k是时间k kk时的测量输出。

在稳定观测器的假设下，通过求解以下线性方程组，可以计算出在稳定状态x r e f x^{ref}xref和稳定状态时输入u r e f u^{ref}uref的状态：
[ A − 1 B C 0 ] [ x r e f , k u r e f , k ] = [ − B d d ^ k r k ] ( 6 ) \left[

A−1CB0�−1��0

\right]\left[

xref,kuref,k����,�����,�

\right]=\left[

−Bdd^{krk−���}���

\right] \quad \quad \quad (6) [A−1CB0][xref,kuref,k]=−Bdd^krk

1.3 非线性模型预测控制 Nonlinear Model Predictive Control

飞行器的行为可以用一组非线性微分方程来更好地描述，以便更好的捕捉空气动力和耦合效应。因此，在本小节中，我们将介绍非线性MPC背后的理论，其涉及整个系统动力学，并且在轨迹追踪方面表现出更好的效果。非线性MPC的优化问题如式(7)所示。

min ⁡ ∫ t = 0 T ∥ h ( x ( t ) , u ( t ) ) − y r e f ( t ) ∥ Q 2 d t + ∥ m ( x ( T ) ) − y r e f ( T ) ∥ P 2 s.t. x ˙ = f ( x ( t ) , u ( t ) ) ; u ( t ) ∈ U C x ( t ) ∈ X C x ( 0 ) = x ( t 0 ) . ( 7 )

\begin{aligned} & \min \int_{t=0}^{T}{\left| h\left( x(t),u(t) \right)-{ {y}^{ref}}(t) \right|{Q}^{2}dt+\left| m\left( x(T) \right)-{ {y}^{ref}}(T) \right|}{P}^{2} \ & \begin{aligned} \text{s}\text{.t}\text{.} \quad & \dot{x}=f\left( x(t),u(t) \right); \ {} & u(t)\in {\mathfrak{U}{C}} \ {} & x(t)\in {\mathfrak{X}{C}} \ {} & x(0)=x({ {t}{0}}). \ \end{aligned}\begin{aligned} & \min \int{t=0}^{T}{\left| h\left( x(t),u(t) \right)-{ {y}^{ref}}(t) \right|{Q}^{2}dt+\left| m\left( x(T) \right)-{ {y}^{ref}}(T) \right|}{P}^{2} \ & \begin{aligned} \text{s}\text{.t}\text{.} \quad & \dot{x}=f\left( x(t),u(t) \right); \ {} & u(t)\in {\mathfrak{U}{C}} \ {} & x(t)\in {\mathfrak{X}{C}} \ {} & x(0)=x({ {t}_{0}}). \ \end{aligned}

\ \end{aligned} \quad \quad \quad \quad(7)min∫t=0T∥∥h(x(t),u(t))−yref(t)∥∥Q2dt+∥∥m(x(T))−yref(T)∥∥P2s.t.x˙=f(x(t),u(t));u(t)∈UCx(t)∈XCx(0)=x(t0).(7)

多重映射技术(multiple shooting technique)可以用来解决最优控制问题(Optimal Control Problem, OCP)。在这种方法中，系统动力学通过时间网格t 0 , ⋯   , t N t_0,\cdots,t_Nt0,⋯,tN进行离散化处理。不等式约束和控制动作在同一时间网格上离散化。在每次迭代中求解一个附加的连续性约束的边值问题(Boundary Value Problem, BVP)。由于系统动力学的性质和施加的约束，优化问题成为一个非线性程序(NLP)。该问题的求解采用序列二次规划(Sequential Quadratic Programming, SQP)技术，其中二次规划(Quadratic Programs, QPs)的求解采用主动集方法，使用qpOASES求解器。

1.4 线性鲁棒模型预测控制 Linear Robust Model Predictive Control

尽管MPC公式具有鲁棒性，但是存在存在鲁棒控制变量时需要进一步的鲁棒性保证。线性鲁棒模型预测控制(Robust Model Predictive Control, RMPC)问题可以表述为一个明确求解的极大极小优化问题。作为一个最优度量，我们可以选择最小峰值性能度量(MPPM)的已知鲁棒性。图2概述了相关的构建块。

图2 Overview of the explicit RMPC optimization problem functional components

在这种RMPC方法中，可以考虑以下系统动力学的线性时不变表示:

x k + 1 = A x k + B u k + G w k y k + 1 = C x k ( 8 )

xk+1=Axk+Buk+Gwkyk+1=Cxk��+1=���+���+�����+1=���

\quad \quad \quad \quad(8)xk+1=Axk+Buk+Gwkyk+1=Cxk(8)
其中x k ∈ X x_k \in Xxk​∈X、u k ∈ U u_k \in Uuk​∈U干扰w k ∈ W w_k \in Wwk​∈W是已知且有界的。本文中考虑有界的干扰(W ∞ = { w : ∣ ∣ x ∣ ∣ ∞ ≤ 1 } W_∞={w:||x||∞\leq1 }W∞​={w:∣∣x∣∣∞​≤1})。因此，针对上述的系统和外加干扰，可以定义成RMPC问题。下面表示所预测的输出、状态、输入和干扰的串联版本，其中[ k + i ∣ k ] [k + i |k][k+i∣k]表示从时间k kk开始的时间k + i k + ik+i的值。
Y = ( y k ∣ k T y k + 1 ∣ k T ⋯ y k + N − 1 ∣ k T ) ( 9 ) \mathcal{Y}=(y{k|k}^T \quad y_{k+1|k}^T\quad \cdots \quad y_{k+N-1|k}^T) \quad \quad \quad \quad(9)Y=(yk∣kT​yk+1∣kT​⋯yk+N−1∣kT​)(9)
X = ( x k ∣ k T x k + 1 ∣ k T ⋯ x k + N − 1 ∣ k T ) ( 10 ) \mathcal{X}=(x_{k|k}^T \quad x_{k+1|k}^T \quad \cdots \quad x_{k+N-1|k}^T) \quad \quad \quad \quad(10)X=(xk∣kT​xk+1∣kT​⋯xk+N−1∣kT​)(10)
U = ( u k ∣ k T u k + 1 ∣ k T ⋯ u k + N − 1 ∣ k T ) ( 11 ) \mathcal{U}=(u_{k|k}^T \quad u_{k+1|k}^T\quad \cdots \quad u_{k+N-1|k}^T) \quad \quad \quad \quad(11)U=(uk∣kT​uk+1∣kT​⋯uk+N−1∣kT​)(11)
W = ( w k ∣ k T w k + 1 ∣ k T ⋯ w k + N − 1 ∣ k T ) ( 12 ) \mathcal{W}=(w_{k|k}^T \quad w_{k+1|k}^T\quad \cdots \quad w_{k+N-1|k}^T) \quad \quad \quad \quad(12)W=(wk∣kT​wk+1∣kT​⋯wk+N−1∣kT​)(12)
其中X ∈ X N = X × X ⋯ × X \mathcal{X} \in \mathbb{X}^N=X \times X \cdots \times XX∈XN=X×X⋯×X，U ∈ U N = U × U ⋯ × U \mathcal{U} \in \mathbb{U}^N=U \times U \cdots \times UU∈UN=U×U⋯×U，W ∈ W N = W × W ⋯ × W \mathcal{W} \in \mathbb{W}^N=W \times W \cdots \times WW∈WN=W×W⋯×W。预测的状态和输出与当前状态、未来的控制输入和干扰呈线性相关，因此如下所示:

X = A x k ∣ k + B U + G W Y = C X ( 13 )

X=Axk|k+BU+GWY=CX�=���|�+��+���=��

\quad \quad \quad \quad(13)

X=Axk∣k+BU+GWY=CX(13)

其中A \mathcal{A}A、B \mathcal{B}B、C \mathcal{C}C和G \mathcal{G}G是堆叠的状态向量矩阵。因此，基于MPPM的RMPC问题(MPPM - RMPC)可表述为:

min ⁡ u   max ⁡ w   ∥ Y ∥ ∞ , ∥ Y ∥ ∞ = max ⁡ ∥ y k + j ∣ k ∥ ∞ subject to u k + j ∣ k ∈ U , ∀ w ∈ W x k + j ∣ k ∈ X , ∀ w ∈ W w k + j ∣ k ∈ W ( 14 )

minumaxwsubject to \quad∥Y∥∞,∥Y∥∞=max∥∥yk+j|k∥∥∞uk+j|k∈U,∀w∈Wxk+j|k∈X,∀w∈Wwk+j|k∈Wmin�max�‖�‖∞,‖�‖∞=max‖��+�|�‖∞subject to \quad��+�|�∈�,∀�∈���+�|�∈�,∀�∈���+�|�∈�

\quad \quad \quad \quad(14)

uminwmaxsubject to ∥Y∥∞,∥Y∥∞=max∥∥yk+j∣k∥∥∞uk+j∣k∈U,∀w∈Wxk+j∣k∈X,∀w∈Wwk+j∣k∈W(14)

1.4.1 反馈预测 Feedback Predictions

按照上述公式，优化问题将趋于保守，因为优化本质上计算的是开环控制序列。反馈预测是一种对遵循后退时域方法的知识进行编码的方法。为了将它们合并，必须假定一种反馈控制结构。在已知导致凸问题空间的反馈参数化中，选择如下:

1.4.2 鲁棒不确定不等式 Robust Uncertain Inequalities Satisfaction

2. 用于多旋翼系统的基于模型的轨迹跟踪控制器 Model-Based Trajectory Tracking Controller for Multi-rotor System

在本节中，我们提出一个多转子系统的简化模型，可用于基于模型的控制来实现轨迹跟踪，并提出一个线性和非线性的模型预测控制器来实现轨迹跟踪。

2.1 多旋翼系统模型 Multirotor System Model

多旋翼系统的6DoF位姿可以通过固定惯性坐标系W WW（世界坐标系）和机身坐标系B BB来定义，如图3所示。我们用p pp表示坐标系W WW中坐标系B BB的原点位置，用R RR表示坐标系W WW中坐标系B BB的旋转矩阵。此外，我们用φ φφ、θ θθ和ψ \psiψ来表示飞行器的横摇、俯仰和偏航角。在该模型中，我们假设一个低级的姿态控制器能够跟踪期望的横摇和俯仰φ d φ_dφd、θ d θ_dθd角，并具有一阶行为。一阶内环近似为MPC提供了充分的信息来考虑低电平控制器的行为。通过经典的系统辨识技术，可以辨识出回路的一阶参数。用于多转子系统轨迹跟踪的非线性模型如式(37)所示。

图3 Illustration of the Firefly hexacopter from Ascending Technologies with attached body fixed frame B and inertial frame W

p ˙ ( t ) = v ( t ) v ˙ ( t ) = R ( ψ , θ , ϕ ) ( 0 0 T ) + ( 0 0 − g ) − ( A x 0 0 0 A y 0 0 0 A z ) v ( t ) + d ( t ) ϕ ˙ ( t ) = 1 τ ϕ ( K ϕ ϕ d ( t ) − ϕ ( t ) ) θ ˙ ( t ) = 1 τ θ ( K θ θ d ( t ) − θ ( t ) ) ( 37 )

p˙(t)=v(t)v˙(t)=R(ψ,θ,ϕ)⎛⎝⎜00T⎞⎠⎟+⎛⎝⎜00−g⎞⎠⎟−⎛⎝⎜Ax000Ay000Az⎞⎠⎟v(t)+d(t)ϕ˙(t)=1τϕ(Kϕϕd(t)−ϕ(t))θ˙(t)=1τθ(Kθθd(t)−θ(t))�˙(�)=�(�)�˙(�)=�(�,�,�)(00�)+(00−�)−(��000��000��)�(�)+�(�)�˙(�)=1��(����(�)−�(�))�˙(�)=1��(����(�)−�(�))

\quad \quad \quad \quad(37)p˙(t)=v(t)v˙(t)=R(ψ,θ,ϕ)⎝⎛00T⎠⎞+⎝⎛00−g⎠⎞−⎝⎛Ax000Ay000Az⎠⎞v(t)+d(t)ϕ˙(t)=τϕ1(Kϕϕd(t)−ϕ(t))θ˙(t)=τθ1(Kθθd(t)−θ(t))(37)
其中v vv代表无人机的速度，g gg是重力加速度，T TT是质量标准化推力，A x A_xAx​、A y A_yAy​和A z A_zAz​表示质量归一化阻力系数，d dd是外部干扰，τ φ τ_φτφ​、K φ K_φKφ​和τ θ τ_θτθ​、K θ K_θKθ​分别为滚转角和俯仰角内环行为的时间常数和增益。

本章假设的级联控制器结构如图4所示。

图4 Controller scheme for multi-rotor system. n1 . . . nm are the i − th rotor speed and y is the measured vehicle state

2.2 Linear MPC

在本小节中，我们将展示如何制定一个线性MPC来实现多旋翼系统的轨迹跟踪，并将其集成到ROS中。利用CVXGEN框架生成c代码求解器，求解式(3)中的优化问题。CVXGEN利用问题结构生成了一个求解凸优化问题的高速求解器。为了清晰起见，我们在这里重写了优化问题，并展示了如何使用CVXGEN生成自定义求解器。优化问题定义如下：

min ⁡ u   ∑ k = 0 N − 1 ( ( x k − x k r e f ) T Q x ( x k − x k r e f ) + ( u k − u k r e f ) T R u ( u k − u k r e f ) + ( u k − u k − 1 ) T R Δ ( u k − u k − 1 ) ) + ( x N − x N r e f ) T P ( x N − x N r e f ) subject to x k + 1 = A x k + B u k + B d d k ； d k + 1 = d k , k = 0 , … , N − 1 u k ∈ U C x 0 = x ( t 0 ) , d 0 = d ( t 0 ) ( 38 )

minusubject to\quad∑k=0N−1⎛⎝⎜⎜⎜(xk−xrefk)TQx(xk−xrefk)+(uk−urefk)TRu(uk−urefk)+(uk−uk−1)TRΔ(uk−uk−1)⎞⎠⎟⎟⎟+(xN−xrefN)TP(xN−xrefN)xk+1=Axk+Buk+Bddk；dk+1=dk,k=0,…,N−1uk∈UCx0=x(t0),d0=d(t0)min�∑�=0�−1((��−�����)���(��−�����)+(��−�����)���(��−�����)+(��−��−1)��Δ(��−��−1))+(��−�����)��(��−�����)subject to\quad��+1=���+���+����；��+1=��,�=0,…,�−1��∈���0=�(�0),�0=�(�0)

\quad \quad \quad \quad(38) uminsubject tok=0∑N−1⎝⎜⎜⎜⎛(xk−xkref)TQx(xk−xkref)+(uk−ukref)TRu(uk−ukref)+(uk−uk−1)TRΔ(uk−uk−1)⎠⎟⎟⎟⎞+(xN−xNref)TP(xN−xNref)xk+1=Axk+Buk+Bddk；dk+1=dk,k=0,…,N−1uk∈UCx0=x(t0),d0=d(t0)(38)

为生成上述优化问题的求解器，CVXGEN中使用了以下问题描述。

dimensions
  m = 3 # dimension of i n p u t s .
  nd = 3 # dimension of d i s t u r b a n c e s .
  nx = 8 # dimension of s t a t e ve c t o r .
  T = 18 # horizon − 1 .
  end
 
  parameters
  A (nx,nx) # dynamics mat r ix .
  B (nx,m) # t r a n s f e r matrix .
  Bd ( nx , nd ) # di s t u r b a n c e t r a n s f e r mat r ix
  Q_x ( nx , nx ) psd # s t a t e cost , p o s i t i v e semidi f ined .
  P (nx,nx) psd # f i n a l s t a t e penal ty , p o s i t i v e semidi f ined .
  R_u (m,m) psd # input penal ty , p o s i t i v e semidi f ined .
  R_delta (m,m) psd # de l t a input penal ty , p o s i t i v e semidi f ined .
  x [0] ( nx ) # i n i t i a l s t a t e .
  d (nd) # di s t u r b a n c e s .
  u_prev (m) # previous input applied to the system .
  u_max (m) # input ampl i tude l i m i t .
  u_min (m) # input ampl i tude l i m i t .
  x_ss [ t ] ( nx ) , t =0. .T+1 # r e f e r e n c e s t a t e .
  u_ss [ t ] (m) , t =0. .T # r e f e r enc e input .
end

variables
  x [ t ] ( nx ) , t =1. .T+1 # s t a t e .
  u [ t ] (m) , t =0. .T # input .
  end
 
  minimize
  quad ( x[0]−x_ss [0] , Q_x) + quad ( u[0]−u_ss [0] , R_u) + quad ( u [0] −
u_prev , R_del ta ) + sum[ t =1 . .T] ( quad ( x [ t ]−x_ss [ t ] , Q_x) + quad (
u [ t ]−u_ss [ t ] , R_u) + quad ( u [ t ] − u [ t −1] , R_delta ) )+quad ( x [T+1]−
x_ss [T+1] , P)
  s u b j e c t to
  x [ t +1] == A∗x [ t ] + B∗u [ t ] + Bd∗d , t=0..T # dynamics
  u_min <= u [ t ] <= u_max , t = 0 . .T # input c o n s t r a i n t .
end
1234567891011121314151617181920212223242526272829303132333435363738

在将生成的求解器集成到ROS中之前，我们讨论了如何估计姿态环参数τ φ τ_φτφ，K φ K_φKφ，τ θ τ_θτθ，K θ K_θKθ，以及离散系统模型A,B,Bd的推导。

3.2.1 姿态环参数辨识 Attitude Loop Parameters Identification

如果不了解无人机上上使用的姿态控制器，则建议进行姿态环识别。 许多市售平台就是这种情况。 姿态命令和实际无人机姿态（使用运动捕捉系统（如果可用）或惯性测量单元（IMU）估算）均带有准确的时间戳记。 通常收集两个数据集，一个用于参数估计，另一个用于验证目的。

要使用可用脚本进行参数识别，请遵循以下步骤:

1. 准备系统，并确保您能够记录带时间戳的姿态命令和实际无人机姿态。
2. 执行飞行，将命令和无人机姿态记录在包文件中。 在飞行过程中，尽可能的激活无人机的每个轴。
3. 重复飞行测试以收集验证数据集。
4. 在提供的脚本中设置正确的包文件名和主题名称，然后运行它。
5. 控制器参数将以确认百分比显示在屏幕上，以确认标识的有效性。

3.2.2 线性化、解耦和离散化 Linearization, Decoupling and Discretization

在设计控制器时，可以假定无人机以小姿态角、飞行器航向与惯性坐标系x轴对齐(即ψ = 0 ψ= 0ψ=0)的情况下，这时可以近似计算悬停状态下的飞行器动力学。围绕悬停状态的线性化方程组可以写成

x ˙ ( t ) = ( 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 − A x 0 0 g 0 0 0 0 0 − A y 0 0 − g 0 0 0 0 0 − A z 0 0 0 0 0 0 0 0 − 1 τ ϕ 0 0 0 0 0 0 0 0 − 1 τ θ ) x ( t ) + ( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 K ϕ τ ϕ 0 0 0 K θ τ θ 0 ) u ( t ) + ( 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 ) d ( t ) ( 39 ) \dot{x}(t)=\left(

000000000000000000000000100−Ax00000100−Ay00000100−Az00000g00−1τϕ00000−g00−1τθ000100000000100000000100000−��00�00000−��00−�00000−��00000000−1��00000000−1��

\right)x(t)+\left(

000000Kϕτϕ00000000Kθτθ00000000000000000000000000����000����0

\right)u(t)+\left(

000100000000100000000100000000000100010001000000

\right)d(t) \quad \quad \quad \quad(39) x˙(t)=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛000000000000000000000000100−Ax00000100−Ay00000100−Az00000g00−τϕ100000−g00−τθ1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞x(t)+⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛000000τϕKϕ00000000τθKθ00000000⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞u(t)+⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛000100000000100000000100⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞d(t)(39)
其中状态向量x = { p T , v T , W ϕ , W θ } T x={p^T,vT,{}^{W}\phi, {}^{W}\theta}Tx={pT,vT,Wϕ,Wθ}T，输入向量u = { W ϕ d , W θ d , T } T u={{}^{W}\phi_d, {}^{{W}\theta_d,T}}Tu={Wϕd​,Wθd​,T}T，干扰向量d = ( d x , d y , d z ) T d=(d_x,d_y,d_z)^Td=(dx​,dy​,dz​)T。公式（38）中三个矩阵分别对应连续时间模型中的A c A_cAc​、B c B_cBc​和B d , c B_{d,c}Bd,c​。注意，在这个线性化中，我们将姿态标记为在惯性坐标系W WW，以消除模型中的偏航角ψ ψψ。在惯性坐标系中计算姿态控制动作W φ d W_{φ_d}Wφd​​、W θ d W_{θ_d}Wθd​​，然后绕z zz轴旋转将姿态控制动作转化为机体框架。机体框架B BB中的控制作用为

( ϕ d θ d ) = ( cos ⁡ ψ sin ⁡ ψ − sin ⁡ ψ cos ⁡ ψ ) ( W ϕ d W θ d ) ( 40 ) \left(

ϕdθd����

\right)\text{=}\left(

cosψ−sinψsinψcosψcos⁡�sin⁡�−sin⁡�cos⁡�

\right)\left(

WϕdWθd������

\right) \quad \quad \quad \quad(40) (ϕdθd)=(cosψ−sinψsinψcosψ)(WϕdWθd)(40)

在计算控制输入后，建议加入前馈项来补偿耦合效应，以获得更好的跟踪性能。为了补偿非零姿态对推力的影响，实现更好的动态轨迹跟踪，采用了以下补偿方案。

T ~ = T + g cos ⁡ ϕ cos ⁡ θ + B z ¨ d ϕ ~ d = g ϕ d − B y ¨ d T ~ θ ~ d = g θ d + B x ¨ d T ~ ( 41 )

T_{=T+gcosϕcosθ+Bz¨dϕ}d=gϕd−By¨dT_θd=gθd+Bx¨dT_�=�+�cos⁡�cos⁡�+��¨��_{�=���−��¨��}�_{�=���+��¨��}

\quad \quad \quad \quad(41)T_{=cosϕcosθT+g+Bz¨dϕ}d=T_{gϕd−By¨dθ}d=T~gθd+Bx¨d(41)
其中B x ¨ d {}^{B}\ddot{x}{d}Bx¨d​、B y ¨ d {}^{B}\ddot{y}{d}By¨​d​和B z ¨ d {}^{B}\ddot{z}_{d}Bz¨d​是机体坐标系下期望轨迹加速度，带波浪符号的量是实际应用的控制输入。

由于控制器是在离散时间内实现的，因此有必要对系统动力学进行离散化(39)。这可以如下所示
A = e A c T s B = ∫ 0 T s e A c τ d τ B c B d = ∫ 0 T s e A c τ d τ B d , c ( 42 )

A=eAcTsB=∫Ts0eAcτdτBcBd=∫Ts0eAcτdτBd,c�=������=∫0������������=∫0����������,�

\quad \quad \quad \quad(42)A=eAcTsB=∫0TseAcτdτBcBd=∫0TseAcτdτBd,c(42)
其中T s T_sTs​为采样时间。终端代价P PP矩阵的计算是通过迭代求解代数Riccati方程来完成的。

2.2.3 ROS集成 ROS Integration

在ROS中集成求解器时，通常需要创建一个ROS节点来将控制器连接到ROS环境，而控制器、估计器和其他组件实现为c++共享库，在编译时链接到节点。控制器节点将无人机状态作为nav里程计消息，并封装成自定义消息类型RollPitchYawRateThrust。规划出的轨迹可以发送到ROS主题MultiDOFJointTrajectory上，也可以作为单个PoseStamped类型的轨迹点。将整个预期轨迹传递给单个点的优势在于，MPC可以考虑未来的预期轨迹并做出相应的反应。图5通过ros图概述了节点和话题之间的通信关系。

图5 ROS graph diagram showing various nodes and topics to control multi-rotor system

控制器每次接收到新的里程计消息时，就计算并发布一个控制动作。因此，这就需要状态估计器以期望的控制速率发布里程测量信息。我们推荐控制器使用多传感器融合的模块化框架(Modular Framework for MultiSensor Fusion)。在实现这种控制器时，需要尽可能地减少循环中的延迟。例如可以通过ROS:: transporthents()和ros::TransportHints()来发送数据给里程计话题的订阅者。

控制器节点发送以下消息：

1. Current desired reference as tf.

图6 Through dynamic reconfiguration it is possible to change the controller parameters online for tuning purposes

1. Desired trajectory as rviz marker.
2. Predicted vehicle state as rviz marker.

可以通过修改图6中的参数来动态修改控制器的参数。控制器的开源代码连接为https://github.com/ethz-asl/mav_control_rw.

图7 Aggressive trajectory tracking performance using linear MPC controller running onboard of Firefly hexacopter from Ascending Technologies

2.2.4 实验结果 Experimental Results

为了验证控制器的性能，我们追踪一段轨迹。控制器运行在上Ascending Technologies (AscTec)的一架Firefly六轴无人机上，机载一台NUC i7计算机。采用外部运动捕捉系统Vicon作为姿态传感器，利用多传感器融合框架(MSF)将其与车载仿真系统融合。控制器运行在100Hz，预测的时长为20步。图7展示了控制器的追踪性能。

2.3 Nonlinear MPC

在本小节中，我们利用无人机的非线性模型设计连续时间非线性模型预测控制器。用于生成非线性求解器的工具包是ACADO，它能够为一般的最优控制问题OCP（optimal control problems）快速生成C语言求解器。优化问题可以写成

min ⁡ U   ∫ t = 0 T ( x ( t ) − x r e f ( t ) ) T Q x ( x ( t ) − x r e f ( t ) ) + ( u ( t ) − u r e f ( t ) ) R u ( u ( t ) − u r e f ( t ) ) d t + ( x ( T ) − x r e f ( T ) ) P ( x ( T ) − x r e f ( T ) ) s.t. x ˙ = f ( x , u ) ; u ( t ) ∈ U C x ( 0 ) = x ( t 0 ) . ( 43 )

\begin{aligned} \underset{U}{\mathop{\min }}, & \int_{t=0}^{T}{ { {\left( x(t)-{ {x}^{ref}}(t) \right)}^{T}}}{ {Q}{x}}\left( x(t)-{ {x}^{ref}}(t) \right)+\left( u(t)-{ {u}^{ref}}(t) \right){ {R}{u}}\left( u(t)-{ {u}^{ref}}(t) \right)dt \ & +\left( x(T)-{ {x}^{ref}}(T) \right)P\left( x(T)-{ {x}^{ref}}(T) \right) \ & \begin{aligned} \text{s}\text{.t}\text{.}\quad & \dot{x}=f\left( x,u \right); \ {} & u(t)\in {\mathfrak{U}{C}} \ {} & x(0)=x({ {t}{0}}). \ \end{aligned}\begin{aligned} \underset{U}{\mathop{\min }}, & \int_{t=0}^{T}{ { {\left( x(t)-{ {x}^{ref}}(t) \right)}^{T}}}{ {Q}{x}}\left( x(t)-{ {x}^{ref}}(t) \right)+\left( u(t)-{ {u}^{ref}}(t) \right){ {R}{u}}\left( u(t)-{ {u}^{ref}}(t) \right)dt \ & +\left( x(T)-{ {x}^{ref}}(T) \right)P\left( x(T)-{ {x}^{ref}}(T) \right) \ & \begin{aligned} \text{s}\text{.t}\text{.}\quad & \dot{x}=f\left( x,u \right); \ {} & u(t)\in {\mathfrak{U}{C}} \ {} & x(0)=x({ {t}{0}}). \ \end{aligned}

\ \end{aligned} \quad \quad \quad \quad(43) Umin∫t=0T(x(t)−xref(t))TQx(x(t)−xref(t))+(u(t)−uref(t))Ru(u(t)−uref(t))dt+(x(T)−xref(T))P(x(T)−xref(T))s.t.x˙=f(x,u);u(t)∈UCx(0)=x(t0).(43)

为了生成上述优化问题的求解器，在ACADO的Matlab界面中使用下面的问题描述。

1 c l c ;
2 c l e a r a l l ;
3 close a l l ;
4
5 Ts = 0 . 1 ; %sampling time
6 EXPORT = 1 ;
7
8 D i f f e r e n t i a l S t a t e p o s i t i o n (3) ve l o c i t y (3) r o l l p i t c h yaw;
9 Control r o l l _ r e f p i t c h _ r e f t h r u s t ;
10
11 %online data r epr e s ent data t h a t can be pa s s ed to the s o l v e r online .
12 OnlineData r o l l _ t a u ;
13 OnlineData r o l l _ g a i n ;
14 OnlineData p i t c h _ t a u ;
15 OnlineData pi t c h _ g a i n ;
16 Onl ineData yaw_rate_command ;
17 OnlineData dr a g _ c o e f f i c i e n t (3) ;
18 OnlineData e x t e r n a l _ d i s t u r b a n c e s (3) ;
19
20 n_XD = length ( d i f f S t a t e s ) ;
21 n_U = length ( c o n t r o l s ) ;
22
23 g = [ 0 ; 0 ; 9 . 8 0 6 6 ] ;
24
25 %% Di f f e r e n t i a l Equation
26
27 %de f ine v e h i c l e body z−axi s expres sed in i n e r t i a l frame .
28 z_axis = [ ( cos (yaw) ∗ s i n ( p i t c h ) ∗cos ( r o l l )+s i n (yaw) ∗ s i n ( rol l ) ) ; . . .
29 ( s i n (yaw) ∗ s i n ( p i t c h ) ∗cos ( r o l l )−cos (yaw) ∗ s i n ( rol l ) ) ; . . .
30 cos ( p i t c h ) ∗cos ( r o l l ) ] ;
31
32 d r o l l = (1/ r o l l _ t a u ) ∗( r o l l _ g a i n ∗ r o l l _ r e f − r o l l ) ;
33 dpitch = (1/ p i t c h _ t a u ) ∗( pi t c h _ g a i n ∗ p i t c h _ r e f − p i t c h ) ;
34
35 f = dot ( [ p o s i t i o n , ve loc i ty ; r o l l ; p i t c h ; yaw] ) == . . .
36 [ ve l o c i t y . . . ;
37 z_axis ∗ t h r u s t − g − diag ( dr a g _ c o e f f i c i e n t ) ∗ v e l o c i t y +
e x t e r n a l _ d i s t u r b a n c e s ; . . .
38 d r o l l ; . . .
39 dpitch ; . . .
40 yaw_rate_command ] ;
41
42 h = [ p o s i t i o n ; ve l o c i t y ; r o l l ; p i t c h ; r o l l _ r e f ; p i t c h _ r e f ; z_axis (3) ∗
t h r u s t −g (3) ] ;
43
44 hN = [ p o s i t i o n ; ve l o c i t y ] ;
45
46
47 %% NMPCexport
48 acadoSet ( ’ problemname ’ , ’ nmpc_trajectory_tracking ’ ) ;
49
50 N = 20;
51 ocp = acado .OCP( 0.0 , N∗Ts , N ) ;
52
53 W_mat = eye ( length ( h ) ) ;
54 WN_mat = eye ( length (hN) ) ;
55 W = acado . BMatrix (W_mat ) ;
56 WN = acado . BMatrix (WN_mat) ;
57
58 ocp . minimizeLSQ ( W, h ) ;
59 ocp . minimizeLSQEndTerm ( WN, hN ) ;
60 ocp . subjectTo(−deg2rad (45) <= [ r o l l _ r e f ; p i t c h _ r e f ] <= deg2rad (45) ) ;
61 ocp . subjectTo ( g (3) /2.0 <= t h r u s t <= g (3) ∗1.5) ;
62 ocp . setModel ( f ) ;
63
64 mpc = acado . OCPexport ( ocp ) ;
65 mpc . s e t ( ’HESSIAN_APPROXIMATION’ , ’GAUSS_NEWTON’ ) ;
66 mpc . s e t ( ’DISCRETIZATION_TYPE ’ , ’MULTIPLE_SHOOTING’ ) ;
67 mpc . s e t ( ’SPARSE_QP_SOLUTION’ , ’FULL_CONDENSING_N2’ ) ;
68 mpc . s e t ( ’INTEGRATOR_TYPE’ , ’INT_IRK_GL4 ’ ) ;
69 mpc . s e t ( ’NUM_INTEGRATOR_STEPS’, N ) ;
70 mpc . s e t ( ’QP_SOLVER’ , ’QP_QPOASES’ ) ;
71 mpc . s e t ( ’HOTSTART_QP’ , ’NO’ ) ;
72 mpc . s e t ( ’LEVENBERG_MARQUARDT’, 1e−10 ) ;
73 mpc . s e t ( ’LINEAR_ALGEBRA_SOLVER’ , ’GAUSS_LU’ ) ;
74 mpc . s e t ( ’IMPLICIT_INTEGRATOR_NUM_ITS’, 4 );
75 mpc . s e t ( ’CG_USE_OPENMP’ , ’YES’ ) ;
76 mpc . s e t ( ’CG_HARDCODE_CONSTRAINT_VALUES’ , ’NO’ ) ;
77 mpc . s e t ( ’CG_USE_VARIABLE_WEIGHTING_MATRIX’ , ’NO’ ) ;
78
79
80
81 i f EXPORT
82 mpc . exportCode ( ’mav_NMPC_trajectory_tracking ’ ) ;
83
84 cd mav_NMPC_trajectory_tracking
85 make_acado_solver ( ’ . . / mav_NMPC_trajectory_tracking ’ )
86 cd . .
87 end
1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889

2.3.1 ROS集成 ROS Integration

非线性控制器的ROS集成遵循先前介绍的线性MPC的相同准则。 可以在此处找到先前介绍的控制器的开源实现 https://github.com/ethz-asl/mav_control_rw。

2.3.2 实验结果 Experimental Results

非线性MPC的设置方式与线性MPC相同， 唯一的不同是，我们目前正在以更具波动性的（aggressive）轨迹评估控制器，如图8所示。

图8 Aggressive trajectory tracking performance using nonlinear MPC controller running onboard of Firefly hexacopter from Ascending Technologies

3 结论 Conclusion

在这一章中，我们概述了多种类型的无人机基于模型的控制策略。提出了多旋翼系统的鲁棒MPC、多旋翼系统的线性MPC和多旋翼系统及固定翼无人机的非线性MPC策略。这些控制策略在实际实验中进行了评价，本章给出了性能评价。所提出的控制器可作为开放源码ROS包在https://github.com/ethz-asl/mav_control_rw和https:// github.com/ethz-asl/mav_control_fw，前者是多旋翼，后者是固定翼。

原文献

• Kamel, M. , Stastny, T. , Alexis, K. , & Siegwart, R. . (2017). Model Predictive Control for Trajectory Tracking of Unmanned Aerial Vehicles Using Robot Operating System. Robot Operating System (ROS) The Complete Reference, Volume 2. Springer International Publishing.

【注】文章来自转载，如有侵权，请联系作者。