Um die wütende Mathematik von Oleg Stepanovich Kozlovs Vorlesungen " Kontrolle in technischen Systemen " zu verwässern , veröffentlichen wir hier ein Beispiel für die praktische Anwendung des Wissens aus diesen Vorlesungen.

In diesem Artikel beschreibt Alexander Shchekaturov, ein Schüler von Oleg Stepanovich, die Erstellung eines Modells eines Oktokopters, der gleichzeitig die Geheimnisse enthüllt und die Techniken des Arbeitens in der Umgebung der strukturellen Modellierung physikalischer Systeme demonstriert.

Dieser Artikel enthält eine Beschreibung eines logisch abgeschlossenen, aber nur einen Teil der gesamten Modellierungsarbeit. Dieser Teil reicht jedoch für eine Einführung in das Thema der dynamischen Modellierung eines UAV vom Typ Drohne aus.

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

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) Momente aus den Schubkräften der VMG. Da sich die Propeller nicht im Schwerpunkt des Hubschraubers befinden, erzeugt jeder Propeller sein eigenes Drehmoment. Dies ist vielleicht der Hauptfaktor, der die Ausrichtung der Drohne im Weltraum bestimmt.

3 Nichtlineare Gleichungen der Oktokopterdynamik

Unter Berücksichtigung der getroffenen Annahmen und der Tatsache, dass der Copter als ein einziger starrer Körper modelliert wird, ist die Grundlage des Copter-Dynamikmodells sehr einfach - dies ist Newtons zweites Gesetz, das in Vektorform so aussieht, dass es nur zwei einfache Gleichungen gibt:

$$

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

Wo $$$\vec {p}(t))=m \cdot \vec{v}(t)$Ist der Impuls des Hubschraubers, und $$$\vec{L} (t)=I(t)\cdot \vec{\omega}(t)$ - der Drehimpuls des Hubschraubers, $$$m$ - seine Masse, $$$I(t)$ - Trägheitstensor (linearer Operator des Trägheitsmoments). $$$\vec{F}(t)$ und $$$\vec{M}(t)$ die Essenz der Summe aller Kräfte und aller Momente, die auf den Hubschrauber wirken.

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

Wenn wir die Werte für den Impuls und den Drehimpuls des Copters in die gegebenen Gleichungen des zweiten Newtonschen Gesetzes einsetzen, erhalten wir für das Trägheitskoordinatensystem:

$$

$$\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}).$$

In diesen Gleichungen sind die rechten Seiten sehr mathematisch komplex, und die Dynamikgleichungen im zugehörigen Koordinatensystem B (wo die rechten Seiten einfacher sind) können nicht so einfach erhalten werden, da sich das dem Copter zugeordnete Koordinatensystem dreht.

– . , – . . : φ (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)$$

Beachten Sie erneut, dass diese Matrix zu jedem Zeitpunkt offensichtlich anders sein wird Die Orientierungswinkel des Copters ändern sich. Dies sind jedoch nur algebraische Berechnungen ohne Differentialgleichungen. Die Matrix selbst ist keine Konstante und hat eine Zeitableitung.

Durch die Methode der Strukturmodellierung implementiert, sieht die Matrix in der SimInTech-Umgebung wie in Abbildung 1 dargestellt aus.

Abbildung 1. Rotationsmatrix von System I zu Koordinatensystem B.

Dann können Sie für den Vektor der linearen Geschwindigkeit schreiben:

$$

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

und für den Drehimpuls

$$

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

Ohne Zwischenberechnungen (einschließlich der Ableitung der Rotationsmatrix, die letztendlich gleich dem Vektorprodukt der Winkelgeschwindigkeit des Objekts durch die Rotationsmatrix selbst ist, genommen mit dem entgegengesetzten Vorzeichen) erhalten wir für die erste der betrachteten Dynamikgleichungen den folgenden Ausdruck im kopterbezogenen Koordinatensystem B :

$$

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

und für die zweite Gleichung unter Berücksichtigung dessen $$$\vec{L_B } (t)=I_B\cdot \vec{\omega_B} (t)$,

Und da in einem rotierenden Koordinatensystem , das gebundene Trägheitstensor eine Konstante ist und dessen zeitliche Ableitung gleich Null ist , und$$$\frac{d\vec{L_B } (t)}{dt}=I_B\cdot \frac{d\vec{\omega_B} (t)}{dt}$ wir bekommen:

$$

$$\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))).$$

Insgesamt wurde die Aufzeichnung II des Newtonschen Gesetzes für den rotierenden Rahmen B im Vergleich zu den ursprünglichen Gleichungen durch zwei Vektorprodukte ergänzt.

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

( B), I, , , , .

Das einzige, was wir mit den Gleichungen noch nicht gemacht haben, ist, dass wir keine Ausdrücke für die Kräfte und Momente geschrieben haben, die auf den Hubschrauber wirken. Lassen Sie es uns unten im B- Koordinatensystem tun . Gemäß den Annahmen berücksichtigen und bezeichnen wir durch Indizes: M - die Arbeit der Motoren, nur in Bezug auf die erzeugte Schubkraft und die daraus resultierenden Momente, D - die Kraft des Luftwiderstands (zusammen mit dem Wind), O - äußere Störung, im Allgemeinen Null, und vom Benutzer willkürlich die Schwerkraft eingestellt - Es entsteht kein Drehmoment der Kräfte:

$$

$$\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)$$

Die Schwerkraft im verknüpften Koordinatensystem "dreht" sich basierend auf der Ausrichtung des Flugzeugs.

Lassen Sie uns genauer schreiben, was die Begriffe bedeuten:

$$

$$\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}}$$

Die Winkeldrehgeschwindigkeiten des VMG, nicht des Copters, sind hier bereits erwähnt. Denken Sie daran, dass die Einheitsvektoren der Zugkräfte für das Koordinatensystem B geschrieben sind und dort Konstanten sind - es ist offensichtlich, dass sie einmal berechnet werden können, der Gesamtkoeffizient hier aus der Klammer genommen wird und die Berechnungen stark vereinfacht werden können. Wenn wir die Gleichungen der Dynamik in System I aufschreiben würden , würde dieser Ausdruck über 8 weitere Multiplikationen mit der Rotationsmatrix jedes Ortes wachsen und dies müsste bei jedem Schritt der Berechnung berücksichtigt werden.

Luftwiderstandskraft (bei Windstille):

$$

$$\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)$$

Externe Störung - Null, auf Wunsch des Benutzers kann er diesen oder jenen Wert später, vor der Berechnung oder während der Simulation einstellen.

Wir schreiben den Moment der Schubkräfte der Motoren wie folgt:

$$

$$\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)$$

Luftwiderstandsmoment (weitere Einzelheiten siehe [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)$$

Beachten Sie noch einmal, dass die Berechnung der Präzession und der reaktiven Momente der HMG in diesem Artikel der Kürze halber weggelassen wird.

Um Fehler beim Übergang von Vektor- zu Skalargleichungen zu vermeiden, die auf die Achsen geschrieben sind, ist es einfacher, das Paket MathCAD oder Maple zu nutzen, bei dem die meisten Konvertierungen in symbolischer Form automatisiert durchgeführt werden können, und die erforderlichen 6 Gleichungen der Dynamik zu erhalten, die von den Achsen des sich bewegenden Koordinatensystems Bed und aufgezeichnet werden .

In der kompaktesten Form sehen die erhaltenen und gelösten Gleichungen der Dynamik folgendermaßen aus:

$$

$$\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).$$

Indem man sie integriert und die Werte der Geschwindigkeiten im B-System erhält, kann man die Geschwindigkeits- und Orientierungswinkel im Trägheitskoordinatensystem I berechnen:

$$

$$\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).$$

Über die Transformationsmatrix von der Winkelgeschwindigkeit des Copters zur Drehgeschwindigkeit um die Eulerwinkel $$$W_{BI}$ siehe Details in [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$und zeichnen Sie das Ergebnis in die SimInTech-Schaltung. Ohne die Berechnungen (der Leser kann sie unabhängig auf einem Blatt Papier oder mit Hilfe einer symbolischen mathematischen Software ausführen) präsentieren wir die endgültige Form der Gleichungen. Aufgrund der Umständlichkeit geben wir auf der rechten Seite Begriffe an.

Zugkraft der Propellergruppen $$$\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),$$

Wo $$$F_{Mi}(t)$ Ist die Zugkraft des i-ten VMG zum aktuellen Zeitpunkt.

Luftwiderstandskraft $$$\vec{F_D}(t)$ Bei Windstille - die oben angegebenen Formeln sind die Projektionen gleich:

$$

$$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|.\\$$

Die Schwerkraft, in der Projektion auf die Achse des sich bewegenden Koordinatensystems, der Begriff $$$g\cdot R_{IB}\cdot \vec {e_{Iz}}$ wird offensichtlich gleich sein:

$$

$$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})$$

Die Summe der Momente der Zugkräfte aller VMGs:

$$

$$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).$$

Wo $$$F_{Mi} (t)$ - Zugkraft des i-ten VMG zum aktuellen Zeitpunkt, $$$l_1,l_2$ - Schulterkräfte (Länge der Balken des Copter-Rahmens für ungerade und gerade VMGs).

Luftwiderstandsmoment (bei Windstille):

$$

$$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.$$

Das Vektorprodukt der Winkelgeschwindigkeit und das Produkt des Trägheitstensors und der Winkelgeschwindigkeit:

$$

$$-\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).$$

Nachdem wir diese Terme in die Differentialgleichung für die Winkelgeschwindigkeit eingesetzt und die in der SimInTech-Umgebung erhaltene implementiert haben, erhalten wir die folgenden Strukturdiagramme:

Abbildung 8. Winkelbeschleunigung des Copters um die Achse $$$x_{B}$...

Abbildung 9. Moment der Schubkräfte von VMG um die Achse $$$x_{B}$...

Abbildung 10. Widerstandskraft der Luft beim Drehen $$$x_{B}$...

3.3 Gleichungssystem für die Dynamik von Oktokoptern in struktureller Form

, , ( 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.