为了淡化Oleg Stepanovich Kozlov的演讲“技术系统中的控制”的数学运算，我们在此发布一个示例，将这些演讲的知识应用于实践。

在本文中，Oleg Stepanovich的学生Alexander Shchekaturov描述了一种八旋翼飞机模型的创建，同时揭示了在物理系统的结构建模环境中工作的秘密和示范技术。

本文提供了逻辑上完成的描述，但仅是建模的全部工作的一部分，但这部分足以介绍无人驾驶型无人机的动态建模主题。

, , : , . ( ), , , — , , , / , , .. .., — .

1

: 10…25 , 4, 6 8- (), , , 8, 12 16 T-motors Antigravity 6007 KV320 ( 2 50% ). .

: , , .

2

, - . , — ( , , ). [1] [2].

, ( ), , :

1. ( ) , , – () . , - . ( , 0), :

   $$

   $$F_M(t)=C_T \cdot \omega^2(t), \\ M_M(t)=C_Q \cdot \omega^2(t),$$
   
   
   , ( ):
   
   $$$C_T=2.02268 \cdot 10^{-4}H \cdot c^2$,
   
   $$$M_M ≈0 \ \ \cdot \cdot ^2$ — ,
   
   $$$ω(t)$ – , /.
2. () , , , . (4 , 4 ) . , - – , . , $$$F_M (t)$ [], , , . , , ( ) .
3. : ) , ) , . I B ( inertial – body – ). : $$$x_I$-, $$$y_I$- , $$$z_I$-, $$$x_B$ — , $$$y_B$- , $$$z_B$- (. ):
   
   
   
   , B :
   
   

   $$

   $$\vec {r_{M1}} = (l_1,0,0)^T, \ \ \ \ \vec {r_{M2}}=\frac{1}{\sqrt 2} (l_2,l_2,0)^T \\ \vec {r_{M3}} = (0,l_1,0)^T, \ \ \ \ \vec {r_{M4}}=\frac{1}{\sqrt 2} (-l_2,l_2,0)^T \\ \vec {r_{M5}} = (l_1,0,0)^T, \ \ \ \ \vec {r_{M6}}=\frac{1}{\sqrt 2} (-l_2,-l_2,0)^T \\ \vec {r_{M7}} = (0,-l_1,0)^T, \ \ \ \ \vec {r_{M8}}=\frac{1}{\sqrt 2} (l_2,-l_2,0)^T$$
   
   
   $$$l_1$ – 1, 3, 4 5- , $$$l_2$ – 2, 4, 6 8- .
   
   

   ( ) , , ( 1…5 , – ). $$$\gamma$, :
   
   

   $$

   $$\vec{e_{M1}}=(0,-sin(\gamma),-cos(\gamma))^T, \ \ \ \vec{e_{M2}}=(-\frac{sin(\gamma)}{\sqrt2},\frac{sin(\gamma)}{\sqrt2},-cos(\gamma))^T,\\ \vec{e_{M3}}=(sin(\gamma),0,-cos(\gamma))^T, \ \ \ \vec{e_{M4}}=(-\frac{sin(\gamma)}{\sqrt2},-\frac{sin(\gamma)}{\sqrt2},-cos(\gamma))^T,\\ \vec{e_{M5}}=(0,sin(\gamma),-cos(\gamma))^T, \ \ \ \vec{e_{M6}}=(\frac{sin(\gamma)}{\sqrt2},-\frac{sin(\gamma)}{\sqrt2},-cos(\gamma))^T,\\ \vec{e_{M7}}=(-sin(\gamma),0,-cos(\gamma))^T, \ \ \ \vec{e_{M8}}=(\frac{sin(\gamma)}{\sqrt2},\frac{sin(\gamma)}{\sqrt2},-cos(\gamma))^T.$$
   
   
   (, B, ) , . , , .
   
   [1].
4. , . :
   
   
   
   4.1) . $$$z_I$. – , . ( , , ).
   
   
   
   4.2) – 8, , .
   
   
   
   4.3) (/ ) – . , ( ). – , « » ( ).
   
   
   
   4.4) – . .
   
   
5. , :
   
   
   
   5.1) . , - – . .
   
   
   
   5.2) – , .
   
   
   
   5.3) – , , .
   
   
   
   5.4) .
   
   
   
   5.5) – .
   
   
   
   5.6）VMG推力产生的力矩。由于螺旋桨不在直升机的重心上，因此每个螺旋桨都会产生自己的转弯力矩。这也许是决定飞机在太空中方位的主要因素。

八旋翼动力学的三个非线性方程

考虑到公认的假设以及将直升飞机建模为单个刚体这一事实，直升飞机动力学模型的基础非常简单-这是牛顿第二定律，在矢量形式下看起来像这样，只有两个简单方程式：

$$

$$\frac{d \vec{p}(t)}{dt}=\vec{F}(t), \\ \frac{d \vec{L}(t)}{dt}=\vec{M}(t),$$

哪里 $$$\vec {p}(t))=m \cdot \vec{v}(t)$是直升机的冲动，并且 $$$\vec{L} (t)=I(t)\cdot \vec{\omega}(t)$ -直升机的角动量， $$$m$ -它的质量， $$$I(t)$ -惯性张量（惯性矩的线性算子）。 $$$\vec{F}(t)$ 和 $$$\vec{M}(t)$ 直升机上所有力量和所有瞬间之和的本质。

, , , ( ). - , - , , , , . – , , B, , ( ) . [1], .

将直升机的动量和角动量的值代入牛顿第二定律的给定方程式中，对于惯性坐标系，我们得到：

$$

$$\frac{d \vec{v_I}(t)}{dt}=\frac{1}{m}\vec{F_I}(t),\\ \frac{d \vec{\omega_I}(t)}{dt}= I^{-1}(\vec{M_I}(t)-\frac{dI_I}{dt}\cdot \vec{\omega_I}).$$

在这些等式中，数学上，右手侧非常复杂，并且由于与直升机相关联的坐标系旋转，因此无法轻松获得相关联的坐标系B（其中右手侧更简单）中的动力学方程。

– . , – . . : φ (roll), θ (pitch) / ψ (yaw), I B $$$R_{IB}$ $$$R_BI=R_{IB}^T$, I, B, , , : $$$\vec{F_I }(t)=R_{BI} \vec{F_B}(t)$, $$$\vec{F_B }(t)=R_{IB} \vec{F_I}(t)$.

, : I , , B. , , $$$\vec F(t)$ ( ), . – . – , , . , , . .

- , , , B. , $$$\vec{ω_B}(t)$ $$$\vec{v_B}(t)$, ( – ) I , , .

I B, cos() = c() sin() = s():

$$

$$R_{IB}=\\ =\left( \begin{array}{ccc} c(\theta)\cdot s(\psi) & c(\theta)\cdot s(\psi) & -s(\theta)\\ s(\varphi) \cdot s(\theta) \cdot c(\psi)- c(\varphi) \cdot s(\psi)&s(\varphi) \cdot s(\theta) \cdot s(\psi) + c(\varphi) \cdot C(\psi) & s(\varphi) \cdot c(\theta)\\ c(\varphi) \cdot s(\theta) \cdot c(\psi) + s(\varphi) \cdot s(\psi)&s(\varphi) \cdot s(\theta) \cdot s(\psi) - s(\varphi) \cdot c(\psi) & c(\varphi) \cdot c(\theta) \end{array} \right)$$

再次注意，该矩阵在每个时刻显然都不同，因为 直升机的方位角发生变化。但是这些只是代数计算，没有任何微分方程。矩阵本身不是常数，它具有时间导数。

通过结构建模方法实现的SimInTech环境中的矩阵如图1所示。

图1.从系统I到坐标系B的旋转矩阵

然后，对于线速度的矢量，您可以编写：

$$

$$\vec {v_B} (t)=R_{IB} \cdot \vec {v_I}(t),$$

对于角动量

$$

$$\vec {L_B} (t)=R_{IB} \cdot \vec{L_I }(t).$$

省略中间计算（包括旋转矩阵的导数，该导数最终等于对象的角速度与旋转矩阵本身的矢量乘积，用相反的符号表示），对于所考虑的第一个动力学方程式，我们在与直升机相关的坐标系B中获得以下表达式：

$$

$$\frac{d\vec{v_B}}{dt} =- \vec{\omega_B}(t) \times \vec{v_B}(t)+ \frac{1}{m}\vec{F_B}(t),$$

对于第二个方程，考虑到 $$$\vec{L_B } (t)=I_B\cdot \vec{\omega_B} (t)$，

并且由于在旋转的边界坐标系中，惯性张量是一个常数，并且其时间导数等于零，并且$$$\frac{d\vec{L_B } (t)}{dt}=I_B\cdot \frac{d\vec{\omega_B} (t)}{dt}$ 我们得到：

$$

$$\frac{d\vec{\omega_B}}{dt} =I_B^{-1}( \vec{M_B}(t) - \vec{\omega_B}(t) \times (I_B \cdot \vec{\omega_B}(t))).$$

总之，与原始方程相比，旋转框架B的牛顿定律的记录II补充了两个向量积。

( 6 ) – . 6 – . 6 , SimInTech Simulink Scilab .

( B), I, , , , .

我们尚未根据方程式完成的唯一事情是，我们还没有写出作用在直升机上的力和力矩的表达式。下面在B坐标系中进行操作。根据这些假设，我们考虑并用指数表示：M-发动机的工作，仅就产生的推力和力矩而言，D-空气阻力（与风共同作用），O-外部干扰，通常为零，并且用户可以任意设置重力-它不会产生力量的转折时刻：

$$

$$\vec {F_B}(t) = \vec {F_M}(t)+\vec {F_D}(t)+\vec {F_O}(t)+m\cdot g \cdot R_{IB}\cdot\vec{e_{IZ}}\\ \vec {M_B}(t) = \vec {M_M}(t)+\vec {M_D}(t)+\vec {M_O}(t)$$

链接坐标系中的重力将根据飞机的方向“旋转”。

让我们更详细地写出这些术语等于什么：

$$

$$\vec {F_M}(t) = C_T\cdot\omega_{M1}^2(t)\cdot \vec {e_{M1}}+ C_T\cdot\omega_{M2}^2(t)\cdot \vec {e_{M2}}+ C_T\cdot\omega_{M3}^2(t)\cdot \vec {e_{M3}}+\\ +C_T\cdot\omega_{M4}^2(t)\cdot \vec {e_{M4}}+ C_T\cdot\omega_{M5}^2(t)\cdot \vec {e_{M5}}+ C_T\cdot\omega_{M6}^2(t)\cdot \vec {e_{M6}}+\\ +C_T\cdot\omega_{M7}^2(t)\cdot \vec {e_{M7}}+ C_T\cdot\omega_{M8}^2(t)\cdot \vec {e_{M8}}$$

这里已经提到了VMG而不是直升飞机的旋转角速度。回想一下，牵引力的单位矢量是为坐标系B编写的，并且它们是常数-很明显，它们可以计算1次，总系数在这里被排除在外，可以大大简化计算。如果我们要写下系统I中的动力学方程，则该表达式将通过每个ort的旋转矩阵增加8个以上的乘法，并且在计算的每个步骤中都必须考虑该表达式。

空气阻力（无风时）：

$$

$$\vec{F_{D}} =-0.5 \rho \cdot C_D \cdot \left( \begin{array}{c} A_{yz}\cdot v_x \cdot |v_x|\\ A_{xz}\cdot v_y \cdot |v_y| \\ A_{xy}\cdot v_z \cdot |v_z| \end{array} \right)$$

外部干扰-零，应用户要求，他可以稍后在计算之前或在模拟过程中设置该值或该值。

我们将发动机推力的力矩写为：

$$

$$\vec {M_M}(t) =\vec{r_{M1}} \times \vec {e_{M1}} \cdot C_T\cdot\omega_{M1}^2(t) + \vec{r_{M2}} \times \vec {e_{M2}} \cdot C_T\cdot\omega_{M2}^2(t)+\\ +\vec{r_{M3}} \times \vec {e_{M3}} \cdot C_T\cdot\omega_{M3}^2(t)+ \vec{r_{M4}} \times \vec {e_{M4}} \cdot C_T\cdot\omega_{M4}^2(t)+\\ +\vec{r_{M5}} \times \vec {e_{M5}} \cdot C_T\cdot\omega_{M5}^2(t)+ \vec{r_{M6}} \times \vec {e_{M6}} \cdot C_T\cdot\omega_{M6}^2(t)+\\ \vec{r_{M7}} \times \vec {e_{M7}} \cdot C_T\cdot\omega_{M7}^2(t)+ \vec{r_{M8}} \times \vec {e_{M8}} \cdot C_T\cdot\omega_{M8}^2(t)$$

空气阻力力矩（有关更多详细信息，请参见[1]）：

$$

$$\vec{M_{D}} =-0.5 \rho \cdot C_D \cdot \left( \begin{array}{c} A_{yz}\cdot v_x \cdot |v_x| \cdot l_x\\ A_{xz}\cdot v_y \cdot |v_y| \cdot l_y\\ 8A_{xy}\cdot v_z \cdot |v_z| \cdot l_z \end{array} \right)$$

再次注意，为简洁起见，本文中省略了HMG的进动力矩和反应力矩的计算。

为避免在轴上写出从向量到标量方程的过渡时出错，可以更轻松地利用MathCAD或Maple类型的软件包，在其中可以以符号形式自动执行大多数转换，并获得由Bed和的移动系统的轴记录的所需的6个动力学方程。

以最紧凑的形式，获得和求解的动力学方程式如下：

$$

$$\frac{d\vec{v_B}}{dt}=\frac{1}{m}\left(\vec{F_M}(t) +\vec{F_D}(t) + \vec{F_O}(t)\right)+ gR_{IB}\vec{e_Iz} -\vec{\omega_B}(t)\times\vec{v_B}(t),\\ \frac{d\vec{\omega_B}}{dt} = I_B^{-1}\left (\vec{M_M}(t) +\vec{M_D}(t)- \vec{\omega_B}(t)\times(I_B\cdot\vec{\omega_B}(t) \right).$$

通过对它们进行积分并获得B系统中的速度值，可以计算出惯性坐标系I中的速度和方向角：

$$

$$\vec{v_I}(t) = R_{BI}\cdot \vec{v_B}(t),\\ \left( \begin{array}{c} \dot{\varphi} \\ \dot{\theta} \\ \dot{\psi} \end{array} \right) = W_{BI} \cdot\vec{\omega_B}(t).$$

关于从直升机的角速度到欧拉角的转弯速率的转换矩阵 $$$W_{BI}$ 请参阅[1]中的详细信息。

3

– SimInTech. , , – . , , – . , , .

, – B 2 , ( ) 6 , 6 : . 6DOF , .. . , 6 , 6 ( ). . , ( ) – . , 12 . , , 18 . , 12 6 ( ) .

3.1

– , «», , . – 6- , . 2.

2 «» , 6+3+3 «». , ( ) – $$$a_{Bx},a_{By},a_{Bz},ω_{Bx},ω_{By},ω_{Bz}$, – (, ).

I ( B ) , ,

, «» , ( ) I, B $$$W_{BI}$. .

2.

« » « » – , 1 .

– ( a = F/m), x, 3:

3.

, .. - , – . , SimInTech ( , ).

W(B->I) B , . 4 [1] .

4. $$$W_{BI}$

3.2

, 2 $$$x_B,y_B,z_B$并在SimInTech电路中绘制结果。省略计算（读者可以在一张纸上或借助某些符号数学软件独立地执行计算），我们给出了方程的最终形式。由于麻烦，我们将在右侧用术语引用。

螺旋桨组的牵引力 $$$\vec{F_M }(t):$

$$

$$F_{}(t)=\frac{1}{\sqrt{2}} \left(-F_{2}(t)-F_{M4}(t)+F_{M6}(t)+F_{M8}(t)\right) \cdot sin(\gamma)+\\ +\left(F_{M3}(t)-F_{M7}(t) \right)\cdot sin(\gamma),\\ F_{}(t)=\frac{1}{\sqrt{2}} \left(F_{2}(t)-F_{M4}(t)-F_{M6}(t)+F_{M8}(t)\right) \cdot sin(\gamma)+\\ + \left( F_{M5}(t)-F_{M1}(t) \right ) \cdot sin(\gamma),\\ F_{z}(t)= \left(F_{M1}(t)+F_{2}(t)+F_{M3}(t)+F_{M4}(t)+\\ +F_{M5}(t)+F_{M6}(t)+F_{M8}(t)+F_{M8}(t)\right) \cdot cos(\gamma),$$

哪里 $$$F_{Mi}(t)$ 是当前第i个VMG的牵引力。

空气阻力 $$$\vec{F_D}(t)$ 在没有风的情况下-上面给出了公式，预测将相等：

$$

$$F_{Dx}(t) = 0.5 \cdot \rho \cdot C_D\cdot A_{yz}\cdot v_x\cdot |v_x|,\\ F_{D}(t) = 0.5 \cdot \rho \cdot C_D\cdot A_{xz}\cdot v_\cdot |v_|,\\ F_{Dz}(t) = 0.5 \cdot \rho \cdot C_D\cdot A_{xy}\cdot v_z\cdot |v_z|.\\$$

重力，在移动坐标系的轴上的投影中，术语 $$$g\cdot R_{IB}\cdot \vec {e_{Iz}}$ 显然等于：

$$

$$g\cdot R_{IB}\cdot \vec {e_{Iz}} =\left( \begin{array}{c} -sin(\theta)\\ sin(\phi)\cdot cos(\theta) \\ cos(\phi)\cdot cos(\theta) \end{array} \right)$$

– :

$$

$$-\vec{\omega_B}(t)\times \vec{v_B}(t) =\left( \begin{array}{c} v_{Bz}\cdot \omega_{By} - v_{By} \cdot \omega_{Bz} \\ v_{Bx}\cdot \omega_{Bz} - v_{Bz} \cdot \omega_{Bx} \\ v_{By}\cdot \omega_{Bx} - v_{Bx} \cdot \omega_{By} \end{array} \right).$$

, SimInTech, , , 5, 6 7 ( $$$x_B$, ):

5. $$$x_{B}$.

6. $$$x_{B}$.

7. $$$x_{B}$.

, ( 2 – ). . , – :

$$

$$\vec{M_m}(t)+\vec{M_D}(t)-\vec{\omega_B}(t) \times(I_B \cdot \vec{\omega_B})$$

所有VMG的牵引力矩之和：

$$

$$M_{Mx}(t)=\left[ \frac{l_2}{\sqrt2} \left( -F_{M2}(t)-F_{M4}(t)+F_{M6}(t)+F_{M8}(t) \right) +\\ +l_1(F_{M7}(t)-F_{M3}(t))\right] \cdot cos(\gamma), \\ M_{M}(t)=\left[ \frac{l_2}{\sqrt2} \left( F_{M2}(t)-F_{M4}(t)-F_{M6}(t)+F_{M8}(t) \right) +\\ +l_1(F_{M1}(t)-F_{M5}(t))\right] \cdot cos(\gamma),\\ M_{Mz}(t)= \left[ l_2 ( F_{M2}(t)+F_{M4}(t)+F_{M6}(t)+F_{M8}(t)) +\\ -l_1(F_{M1}(t)+F_{M3}(t)+F_{M5}(t)+F_{M7}(t)) \right] \cdot sin(\gamma).$$

哪里 $$$F_{Mi} (t)$ -当前第i个VMG的牵引力， $$$l_1,l_2$ -肩力（奇数和偶数VMG的直升飞机镜架的光束长度）。

空气阻力力矩（无风时）：

$$

$$M_{Dx}(t) = - 0.5 \cdot \rho \cdot C_D\cdot A_{yz}\cdot \omega_{Bx}\cdot |\omega_{Bx}|\cdot l_x,\\ M_{D}(t) = - 0.5 \cdot \rho \cdot C_D\cdot A_{xz}\cdot \omega_{By}\cdot |\omega_{By}|\cdot l_y,\\ M_{Dz}(t) = - 0.5 \cdot \rho \cdot C_D\cdot 8\cdot A_{xy}\cdot \cdot \omega_{Bz}\cdot |\omega_{Bz}|\cdot l_z.$$

角速度的矢量乘积与惯性张量和角速度的乘积：

$$

$$-\vec{\omega_B}(t)\times (I_B \cdot\vec{\omega_B}(t)) =\left( \begin{array}{c} (I_{yy}-I_{zz})\cdot \omega_{By}\cdot \omega_{Bz} \\ (I_{zz}-I_{xx})\cdot \omega_{Bx}\cdot \omega_{Bz} \\ (I_{xx}-I_{yy})\cdot \omega_{Bx}\cdot \omega_{By} \end{array} \right).$$

因此，将这些项代入角速度的微分方程，并实现在SimInTech环境中获得的项后，我们得到以下结构图：

图8.直升机绕轴的角加速度 $$$x_{B}$...

图9. VMG绕轴的推力力矩 $$$x_{B}$...

图10.旋转时的空气阻力 $$$x_{B}$...

3.3结构形式的八旋翼动力学方程组

, , ( 11). , , . (, , ), 6+ , , , . 12.

11. .

12. , .

, SimInTech, .

( , ), ( ) . , // - — , .

, () – B I, – I ( B .. – ).

, .

, , .

1. Design, Modeling and Control of an Octocopter, Oscar Oscarson, Royal Institute of Technology 2015
2. Development, Modelling and Control of a Multirotor Vehicle, Markus Mikkelsen, Umeå University 2015.