Dynamic Equations for the 3-DOF Helicopter
© 2004 Quanser Consulting Inc.
Description
The Lagrange's method is used to obtain the dynamic model of the system.
The Quanser_Tools Package
The Quanser_Tools module defines generic procedures and data in relation to determining the state-space representation of all the Quanser experiments. Specifically, this means deriving and solving the Lagrange's equations of the Quanser systems.
The quanser repository containing the Quanser_Tools package is implemented in the 2 following files: quanser.ind and quanser.lib. If these two files are not readily available, they can be generated by executing the Maple worksheet titled: quanser_tools.mws.
To install the Quanser_Tools package, copy the two files quanser.ind and quanser.lib into a directory of your choice, like for example: "C:\Program Files\Quanser\Maple Repository".
To use the Quanser_Tools package in a Maple worksheet, add the path to its disk location to the Maple global variable libname. For example, this can be achieved by the following Maple command:
libname := "C:/Program Files/Quanser/Maple Repository", libname:
Worksheet Initialization
> | restart: interface( imaginaryunit = j ): |
> | with( LinearAlgebra ): |
> | libname := "C:/Program Files/Quanser/Maple Repository", ".", libname: |
> | with( Quanser_Tools ); |
environment variable representing the order of series calculations
> | Order := 2: |
Notations
Generalized Coordinates:
's
The generalized coordinates are also called Lagrangian coordinates.
> | q := [ epsilon(t), p(t), lambda(t) ]; |
Nq = number of Lagrangian coordinates
Nq is also the number of position states.
> | Nq := nops( q ): |
qd = first-order time derivative of the generalized coordinates
> | qd := map( diff, q, t ); |
Cartesian Coordinates of the Moving Bodies
> | HTM_BASE_TO_TRAVEL := HTM( 'rot', 'Z', -lambda(t) ); |
> | HTM_TRAVEL_TO_CW := Multiply( HTM( 'rot', 'X', epsilon(t) ), HTM( 'trans', 0, -L[w], 0 ) ); |
> | HTM_TRAVEL_TO_HB := Multiply( HTM( 'rot', 'X', epsilon(t) ), HTM( 'trans', 0, L[a], 0 ) ); |
> | HTM_HB_TO_FM := Multiply( HTM( 'rot', 'Y', -p(t) ), HTM( 'trans', L[h], 0, 0 ) ); |
> | HTM_HB_TO_BM := Multiply( HTM( 'rot', 'Y', -p(t) ), HTM( 'trans', -L[h], 0, 0 ) ); |
> | HTM_BASE_TO_CW := Multiply( HTM_BASE_TO_TRAVEL, HTM_TRAVEL_TO_CW ); |
> | HTM_BASE_TO_HB := Multiply( HTM_BASE_TO_TRAVEL, HTM_TRAVEL_TO_HB ); |
> | HTM_BASE_TO_FM := Multiply( HTM_BASE_TO_HB, HTM_HB_TO_FM ); |
> | HTM_BASE_TO_BM := Multiply( HTM_BASE_TO_HB, HTM_HB_TO_BM ); |
> | x[cw] := HTM_BASE_TO_CW[ 1, 4 ];
y[cw] := HTM_BASE_TO_CW[ 2, 4 ]; z[cw] := HTM_BASE_TO_CW[ 3, 4 ]; |
> | x[fm] := HTM_BASE_TO_FM[ 1, 4 ];
y[fm] := HTM_BASE_TO_FM[ 2, 4 ]; z[fm] := HTM_BASE_TO_FM[ 3, 4 ]; |
> | x[bm] := HTM_BASE_TO_BM[ 1, 4 ];
y[bm] := HTM_BASE_TO_BM[ 2, 4 ]; z[bm] := HTM_BASE_TO_BM[ 3, 4 ]; |
> | xd[cw] := diff( x[cw], t );
yd[cw] := diff( y[cw], t ); zd[cw] := diff( z[cw], t ); |
> | xd[fm] := diff( x[fm], t );
yd[fm] := diff( y[fm], t ); zd[fm] := diff( z[fm], t ); |
> | xd[bm] := diff( x[bm], t );
yd[bm] := diff( y[bm], t ); zd[bm] := diff( z[bm], t ); |
State-Space Variables
The chosen states should at least include the generalized coordinates and their first-time derivatives.
X is the state vector.
In the state vector X: Lagrangian coordinates are first, followed by their first-time derivatives, and finally any other states, as required.
Substitution sets for the states (to obtain time-independent state equations).
> | subs_Xq := { seq( q[i] = X[i], i=1..Nq ) };
subs_Xqd := { seq( qd[i] = X[i+Nq], i=1..Nq ) }; |
Substitution set for the input(s).
> | subs_U := { V[f] = U[1], V[b] = U[2] }; |
Nu = number of inputs; U = input (row) vector
> | Nu := nops( subs_U ): |
substitution set for the position states' second time derivatives
> | subs_Xqdd := { seq( diff( q[i], t$2 ) = Xd[i+Nq], i=1..Nq ) }; |
second time derivatives of the position states (written as time-independent variables).
The set of unknowns is obtained from this list to solve the Lagrange's equations of motion.
> | Xqdd := [ seq( Xd[i+Nq], i=1..Nq ) ]; |
substitution set to linearize the state-space matrices (i.e. A and B)
about the quiescent null state vector (small-displacement theory)
> | subs_XU_op := { seq( X[i] = 0, i=1..2*Nq ), seq( U[i] = 0, i=1..Nu ) }: |
Nx = dim( X ) = total number of states (should be greater than or equal to: 2 * Nq)
Ny = chosen number of outputs
> | Nx := 2 * Nq + 0:
Ny := Nq: |
Total Potential and Kinetic Energies of the System
The total potential and kinetic energies are needed to calculate the Lagrangian of the system.
Total Potential Energy:
The total potential energy can be expressed in terms of the generalized coordinates alone.
Ve[T] = Total Elastic Potential Energy of the system
> | Ve[T] := 0: |
Vg[T] = Total Gravitational Potential Energy of the system
initialization:
> | Vg[cw] := potential_energy( 'gravity', M[c], g, z[cw] ); |
> | Vg[fm] := potential_energy( 'gravity', M[f], g, z[fm] ); |
> | Vg[bm] := potential_energy( 'gravity', M[b], g, z[bm] ); |
> | Vg[T] := Vg[cw] + Vg[fm] + Vg[bm]: |
V[T] = Total Potential Energy of the system
> | V[T] := simplify( Ve[T] + Vg[T] ); |
Total Kinetic Energy:
The total kinetic energy can be expressed in terms of the generalized coordinates and their first-time derivatives.
> | v[cw] := n_norm( [ xd[cw], yd[cw], zd[cw] ], 2 ): |
> | Tt[cw] := kinetic_energy( 'translation', M[c], v[cw] ); |
> | v[fm] := n_norm( [ xd[fm], yd[fm], zd[fm] ], 2 ): |
> | Tt[fm] := kinetic_energy( 'translation', M[f], v[fm] ); |
> | v[bm] := n_norm( [ xd[bm], yd[bm], zd[bm] ], 2 ): |
> | Tt[bm] := kinetic_energy( 'translation', M[b], v[bm] ); |
= Total Kinetic Energy of the system
> | T[T] := Tt[cw] + Tt[fm] + Tt[bm]:
T[T] := collect( T[T], diff ); |
Generalized Forces:
's
Qe: on epsilon: elevation
> | Q[1] := L[a] * K[m] * ( V[f] + V[b] ); |
Qp: on pitch, p
> | Q[2] := K[m] * ( V[f] - V[b] ) * L[h]; |
Qt: on travel, lambda
> | Q[3] := L[a] * K[m] * ( V[f] + V[b] ) * sin( p(t) ); |
> | Q := [ seq( Q[i], i=1..Nq ) ]; |
quiescent voltage: operating point
> | EOM_OP_PT := L[w] * M[c] * g = L[a] * g * ( M[f] + M[b] ) + L[a] * ( V[op] + V[op] ) * K[m]; |
> | V[op] := solve( EOM_OP_PT, V[op] ); |
> | subs_U_op := { V[f] = V[f] - V[op], V[b] = V[b] - V[op] }: |
> | Q := subs( subs_U_op, Q ); |
Euler-Lagrange's Equations
For a N-DOF system, the Lagrange's equations can be written:
for
through
where:
's are special combinations of external forces and called the generalized forces,
, ...,
, are N independent coordinates chosen to describe the system and called the generalized coordinates,
and
is the Lagrangian of the system.
is defined by:
where
is the total kinetic energy of the system and
the total potential energy of the system.
> | EOM_orig := lagrange_equations( T[T], V[T], q, Q ): |
this is to display the EOM's
> | EOM_orig := collect( EOM_orig, { seq( diff( q[i], t$2 ), i=1..Nq ), seq( diff( q[i], t ), i=1..Nq ), seq( q[i], i=1..Nq ) } ); |
Express the Euler-Lagrange equations of motion as functions of the states:
1) substitute (i.e. name) the acceleration states first!
2) then substitute the velocity states!
3) and only after, the position states, and the inputs!
> | EOM_states := subs( subs_Xqd, subs( subs_Xqdd, EOM_orig ) ): |
> | EOM_states := subs( subs_Xq, subs_U, EOM_states ); |
Linearization in the EOM's of the Trigonometric Functions
Linearization of the equations of motion around the quiescent point of operation (in order to solve them).
Here, linearization around the zero angles, i.e. for small-amplitude oscillations.
Linearization around: alpha0 = 0, and alpha_dot0 = 0
> | EOM_ser := EOM_states: |
Generalized series expansions of the trigonometric functions is used (for small angles).
> | for i from 1 to Nq do
EOM_ser[i] := subsop( 1 = convert( series( op( 1, EOM_ser[i] ), X[1] ), polynom ), EOM_ser[i] ); EOM_ser[i] := subsop( 1 = convert( series( op( 1, EOM_ser[i] ), X[2] ), polynom ), EOM_ser[i] ); EOM_ser[i] := subsop( 1 = convert( series( op( 1, EOM_ser[i] ), X[3] ), polynom ), EOM_ser[i] ); EOM_ser[i] := simplify( EOM_ser[i] ); end do: |
> | EOM_ser; |
Additional Insight: Inertia (or mass) Matrix: Fi
The nonlinear system of equations resulting from the Lagrangian mechanics can be written in the following matrix form:
F( q ) . qdd + G( q, qd ) . qd + H( q ) . q = L( q, qd, u )
F, G, and H are called, respectively, the mass, damping, and stiffness matrices.
They are symmetric in form.
The inertia (a.k.a. mass) matrix, F, gives indications regarding the coupling existing in the system.
> | Fi := Matrix( Nq, Nq ): |
> | for i from 1 to Nq do
for k from 1 to Nq do Fi[ i, k ] := simplify( diff( op( 1, EOM_states[i] ), Xd[k+Nq] ) ); Fi[ i, k ] := collect( combine( Fi[ i, k ], trig ), cos ); end do; end do: |
> | 'F[i]' = Fi; |
Linearization of the inertia matrix for small-displacements
> | Fi_lin := Matrix( Nq, Nq ): |
> | for i from 1 to Nq do
for k from 1 to Nq do Fi_lin[ i, k ] := Fi[ i, k ]; Fi_lin[ i, k ] := convert( series( Fi_lin[ i, k ], X[ 2 ] ), polynom ); Fi_lin[ i, k ] := subs( subs_XU_op, Fi_lin[ i, k ] ); end do; end do: |
> | 'F_lin' = Fi_lin; |
Solving the Euler-Lagrange's Equations
To solve the Euler-Lagrange's equations, they need to be linear.
Reverse State Substitution for Pretty Display of the Solved EOM's
only for pretty print
> | subs_Xq_rev := { seq( X[i] = q[i], i=1..Nq ) }:
subs_Xqd_rev := { seq( X[i+Nq] = qd[i], i=1..Nq ) }: |
> | subs_U_rev := { U[1] = V[f], U[2] = V[b] }; |
> | eom_collect_list := { seq( diff( q[i], t ), i=1..Nq ), seq( q[i], i=1..Nq ) }; |
Solution to the Non-Linear Equations of Motion
Solve the non-linear form of the equations of motion for the states' second time derivatives
> | Xqdd_solset_nl := solve( convert( EOM_states, set ), convert( Xqdd, set ) ): |
> | assign( Xqdd_solset_nl ); |
> | Xd_nl[Nq+1] := simplify( Xd[Nq+1] ):
Xd_nl[Nq+2] := simplify( Xd[Nq+2] ): Xd_nl[Nq+3] := simplify( Xd[Nq+3] ): unassign( 'Xd[Nq+1]', 'Xd[Nq+2]', 'Xd[Nq+3]' ): |
pretty display w.r.t. the named system states
> | Xd_nl[Nq+1] := simplify( subs( subs_U_rev, subs_Xq_rev, subs_Xqd_rev, Xd_nl[Nq+1] ) ):
diff( qd[1], t ) = collect( Xd_nl[Nq+1], eom_collect_list ); |
> | Xd_nl[Nq+2] := simplify( subs( subs_U_rev, subs_Xq_rev, subs_Xqd_rev, Xd_nl[Nq+2] ) ):
diff( qd[2], t ) = collect( Xd_nl[Nq+2], eom_collect_list ); |
> | Xd_nl[Nq+3] := simplify( subs( subs_U_rev, subs_Xq_rev, subs_Xqd_rev, Xd_nl[Nq+3] ) ):
diff( qd[3], t ) = collect( Xd_nl[Nq+3], eom_collect_list ); |
Solution to the Linearized EOM's
Solve the linear form of the equations of motion for the states' second time derivatives
> | Xqdd_solset_ser := solve( convert( EOM_ser, set ), convert( Xqdd, set ) ); |
> | assign( Xqdd_solset_ser ); |
Moreover, for small angles
> | subs_small_angles_list := { X[1]^2 = 0, X[2]^2 = 0, X[3]^2 = 0, X[Nq+1]^2 = 0, X[Nq+2]^2 = 0, X[Nq+3]^2 = 0, M[b] = M[f] }: |
> | Xd[Nq+1] := algsubs( sin(X[2]) = X[2], Xd[Nq+1] ): |
> | Xd[Nq+1] := subs( subs_small_angles_list, Xd[Nq+1] ): |
> | Xd[Nq+1] := algsubs( X[Nq+1] * X[Nq+2] = 0, Xd[Nq+1] ): |
> | Xd[Nq+1] := algsubs( X[Nq+1] * X[Nq+3] = 0, Xd[Nq+1] ): |
> | Xd[Nq+1] := algsubs( X[Nq+2] * X[Nq+3] = 0, Xd[Nq+1] ): |
> | Xd[Nq+1] := algsubs( X[1] * X[2] = 0, Xd[Nq+1] ): |
> | Xd[Nq+1] := algsubs( X[1] * X[3] = 0, Xd[Nq+1] ): |
> | Xd[Nq+1] := algsubs( X[2] * X[3] = 0, Xd[Nq+1] ); |
> | Xd[Nq+2] := algsubs( sin(X[2]) = X[2], Xd[Nq+2] ): |
> | Xd[Nq+2] := subs( subs_small_angles_list, Xd[Nq+2] ): |
> | Xd[Nq+2] := algsubs( X[Nq+1] * X[Nq+2] = 0, Xd[Nq+2] ): |
> | Xd[Nq+2] := algsubs( X[Nq+1] * X[Nq+3] = 0, Xd[Nq+2] ): |
> | Xd[Nq+2] := algsubs( X[Nq+2] * X[Nq+3] = 0, Xd[Nq+2] ): |
> | Xd[Nq+2] := algsubs( X[1] * X[2] = 0, Xd[Nq+2] ): |
> | Xd[Nq+2] := algsubs( X[1] * X[3] = 0, Xd[Nq+2] ): |
> | Xd[Nq+2] := algsubs( X[2] * X[3] = 0, Xd[Nq+2] ); |
> | Xd[Nq+3] := algsubs( sin(X[2]) = X[2], Xd[Nq+3] ): |
> | Xd[Nq+3] := subs( subs_small_angles_list, Xd[Nq+3] ): |
> | Xd[Nq+3] := algsubs( X[Nq+1] * X[Nq+2] = 0, Xd[Nq+3] ): |
> | Xd[Nq+3] := algsubs( X[Nq+1] * X[Nq+3] = 0, Xd[Nq+3] ): |
> | Xd[Nq+3] := algsubs( X[Nq+2] * X[Nq+3] = 0, Xd[Nq+3] ): |
> | Xd[Nq+3] := algsubs( X[1] * X[2] = 0, Xd[Nq+3] ): |
> | Xd[Nq+3] := algsubs( X[1] * X[3] = 0, Xd[Nq+3] ): |
> | Xd[Nq+3] := algsubs( X[2] * X[3] = 0, Xd[Nq+3] ); |
pretty display w.r.t. the named system states
1st Linearized Equation Of Motion in the time domain
> | L_EOM_1_DT := diff( qd[1], t ) = collect( subs( subs_U_rev, subs_Xq_rev, subs_Xqd_rev, Xd[Nq+1] ), eom_collect_list ):
L_EOM_1_DT; |
2nd Linearized Equation Of Motion in the time domain
> | L_EOM_2_DT := diff( qd[2], t ) = collect( subs( subs_U_rev, subs_Xq_rev, subs_Xqd_rev, Xd[Nq+2] ), eom_collect_list ):
L_EOM_2_DT; |
2nd Linearized Equation Of Motion (i.e., re. alpha) in the time domain
> | L_EOM_3_DT := diff( qd[3], t ) = collect( subs( subs_U_rev, subs_Xq_rev, subs_Xqd_rev, Xd[Nq+3] ), eom_collect_list ):
L_EOM_3_DT; |
Determine the System State-Space Matrices: A, B, C, and D
> | A_ss := Matrix( Nx, Nx ): |
> | A_ss := deriveA( Xqdd, A_ss, Nq, subs_XU_op ): 'A' = A_ss; |
> | #lprint( A_ss ); |
> | B_ss := Matrix( Nx, Nu ): |
> | B_ss := deriveB( Xqdd, B_ss, Nq, subs_XU_op ): 'B' = B_ss; |
> | #lprint( B_ss ); |
> | C_ss := IdentityMatrix( Nx, Nx ); |
> | #lprint( C_ss ); |
> | D_ss := Matrix( Nx, Nu, 0 ); |
> | #lprint( D_ss ); |
Write A, B, C, and D to a Matlab file
Save the state-space matrices A, B, C and D to a MATLAB file.
> | Matlab_File_Name := "HELI3D": |
> | Matlab_File_Name := cat( Matlab_File_Name, "_ABCD_eqns_VfVb.m" ); |
substitution set containing a notation consistent with that used in the MATLAB design script(s)
> | Matlab_Notations := { L[w] = Lw, L[h] = Lh, M[f] = mf, L[a] = La, M[b] = mb, M[c] = mc, K[m] = Km, g = g }: |
> | Experiment_Name := "3-DOF Helicopter": |
> | write_ABCD_to_Mfile( Matlab_File_Name, Experiment_Name, Matlab_Notations, A_ss, B_ss, C_ss, D_ss ); |
Procedure Printing
default:
> | #interface( verboseproc = 1 ); |
> | #eval( lagrange_equations ); |
> | #eval( deriveB ); |
> |
> |
> |