FLEXPEN.mws

Dynamic Equations for the
Flexible Pendulum (FLEXPEN) Experiment

©  2012 Quanser Consulting Inc.

URL: http://www.quanser.com

NOTE: This worksheet presents the general dynamic equations modelling a Flexible Pendulum mounted on a linear cart.

Description

The Quanser_Tools  Package

Worksheet Initialization

>    restart: interface( imaginaryunit = j ):

>    with( LinearAlgebra ):

>    libname := "C:/Program Files/Quanser/Maple Repository", ".", libname:

>    with( Quanser_Tools );

HTM deriveA deriveB deriveF kinetic_energy lagrange_equations moment_of_inertia n_norm potential_energy write_ABCD_to_Mfile

environment variable representing the order of series calculations

>    Order := 2:

Notations

Generalized Coordinates: q[i] 's

The generalized coordinates are also called Lagrangian coordinates.

x[c](t)  =  driving cart (i.e. IP02 ) linear position

α t  = Position angle of the solid pendulum about the vertical: the zero angle [mod 2*pi] is defined when the IP is in the perfect upright position.

  γ t  = Deflection angle of the flexible pendulum.

>    q := [ x[c](t), alpha(t), gamma(t) ];

q := [x[c](t), alpha(t), gamma(t)]

Nq = number of Lagrangian coordinates

Nq is also the number of position states.

>    Nq := nops( q ):

Np = number of pendulums

>    Np := Nq - 1:

qd = first-order time derivative of the generalized coordinates, e.g. position and angular velocities

>    qd := map( diff, q, t );

qd := [diff(x[c](t),t), diff(alpha(t),t), diff(gamma(t),t)]

Cartesian Coordinates of the Flexible Pendulum's Centre Of Gravity

conventions:

1) α = 0  corresponds to the stiff pendulum being perfectly vertical, pointing upwards.

2) γ = 0  corresponds to the flexible pendulum being perfectly stiff. Having both γ = 0  and α = 0  means the flexible pendulum is vertical, point upwards, and stiff.

3) Positive rotation is CCW when facing the cart.

x[sp]  = x-coordinate of the solid pendulum (i.e. subscript "sp") COG absolute Cartesian position

y[sp]  = y-coordinate of the solid pendulum (i.e. subscript "sp") COG absolute Cartesian position

>    x[sp] := q[1] - l[sp] * sin( q[2] );
y[sp] := l[sp] * cos( q[2] );

x[sp] := x[c](t)-l[sp]*sin(alpha(t))

y[sp] := l[sp]*cos(alpha(t))

x[fp]  = x-coordinate of the flexible pendulum (i.e. subscript "fp") COG absolute Cartesian position

y[fp]  = y-coordinate of the flexible pendulum (i.e. subscript "fp") COG absolute Cartesian position

>    x[fp] :=  q[1] - l[sp] * sin( q[2] ) - l[fp] * sin(q[3]);
y[fp] := l[sp] * cos( q[2] ) + l[fp] * cos(q[3]);

x[fp] := x[c](t)-l[sp]*sin(alpha(t))-l[fp]*sin(gamma(t))

y[fp] := l[sp]*cos(alpha(t))+l[fp]*cos(gamma(t))

xd[sp]  = x-coordinate of the stiff IP's COG absolute Cartesian velocity

yd[sp]  = y-coordinate of the stiff IP's COG absolute Cartesian velocity

>    xd[sp] := diff( x[sp], t );
yd[sp] := diff( y[sp], t );

xd[sp] := diff(x[c](t),t)-l[sp]*cos(alpha(t))*diff(alpha(t),t)

yd[sp] := -l[sp]*sin(alpha(t))*diff(alpha(t),t)

xd[fp]  = x-coordinate of the stiff IP's COG absolute Cartesian velocity

yd[fp]  = y-coordinate of the stiff IP's COG absolute Cartesian velocity

>    xd[fp] := diff( x[fp], t );
yd[fp] := diff( y[fp], t );

xd[fp] := diff(x[c](t),t)-l[sp]*cos(alpha(t))*diff(alpha(t),t)-l[fp]*cos(gamma(t))*diff(gamma(t),t)

yd[fp] := -l[sp]*sin(alpha(t))*diff(alpha(t),t)-l[fp]*sin(gamma(t))*diff(gamma(t),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 ) };

subs_Xq := {alpha(t) = X[2], gamma(t) = X[3], x[c](t) = X[1]}

subs_Xqd := {diff(x[c](t),t) = X[4], diff(alpha(t),t) = X[5], diff(gamma(t),t) = X[6]}

Substitution set for the input(s).

set the input to be the motor voltage:

>    #subs_U := { V[m] = U[1] }:

set the input to be the cart's driving force: F[c]  (if not expressed as a function of the motor voltage):

>    subs_U := { F[c] = U[1] }:

Nu = number of inputs; U = input (row) vector (e.g. U = [ V[m] ] )

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

subs_Xqdd := {diff(x[c](t),`$`(t,2)) = Xd[4], diff(alpha(t),`$`(t,2)) = Xd[5], diff(gamma(t),`$`(t,2)) = Xd[6]}

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 ) ];

Xqdd := [Xd[4], Xd[5], Xd[6]]

substitution set to linearize the state-space matrices (i.e. A and B)

about the quiescent null state vector (small-displacement theory)

>    subs_XUzero := { 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 Lagragian of the system.

Total Potential Energy: V[T]  

The total potential energy can be expressed in terms of the generalized coordinates alone.

V[e] = Total Elastic Potential Energy of the system

>    V[e] := potential_energy( 'spring', K[s], q[3] );

V[e] := 1/2*K[s]*gamma(t)^2

V[g] = Total Gravitational Potential Energy of the system

initialization:

>    V[g] := potential_energy( 'gravity', M[sp], g, y[sp] ) + potential_energy( 'gravity', M[fp], g, y[fp] );

V[g] := M[sp]*g*l[sp]*cos(alpha(t))+M[fp]*g*(l[sp]*cos(alpha(t))+l[fp]*cos(gamma(t)))

V[T] = Total Potential Energy of the system

>    V[T] := simplify( V[g] + V[e] );

V[T] := M[sp]*g*l[sp]*cos(alpha(t))+M[fp]*g*l[sp]*cos(alpha(t))+M[fp]*g*l[fp]*cos(gamma(t))+1/2*K[s]*gamma(t)^2

Total Kinetic Energy: T[T]

The total kinetic energy can be expressed in terms of the generalized coordinates and their first-time derivatives.

Tt[c]  = translational kinetic energy of the motorized cart (e.g. IP02)

M[c]  = cart total mass

>    Tt[c] := kinetic_energy( 'translation', M[c], qd[1] );

Tt[c] := 1/2*M[c]*diff(x[c](t),t)^2

The cart's directions of translation and rotation are orthogonal.

Tr[c]  = kinetic energy due to rotation tied to the motorized cart

J[m]  = motor armature inertia; K[g]  = gear ratio; r[mp]  = motor pinion radius

>    Tr[c] := kinetic_energy( 'rotation', J[m], omega[m](t) ):
'Tr[c]' = Tr[c];

Tr[c] = 1/2*J[m]*omega[m](t)^2

>    omega[m](t) := K[g] * qd[1] / r[mp];

omega[m](t) := K[g]*diff(x[c](t),t)/r[mp]

Since J[m]*K[g]^2/(r[mp]^2)  < 0.14 kg < M[c] , the rotational kinetic energy on the IP02 can be neglected:

>    Tr[c] := 0:

T[T_c]  = Total Kinetic Energy of the cart system

initialization:

>    T[T_c] := Tt[c] + Tr[c];

T[T_c] := 1/2*M[c]*diff(x[c](t),t)^2

The rigid pendulum, or STIFFPEN, motion is composed of one rotation and one translation, both directions are orthogonal.

T[r,sp] = stiff pendulum's rotational kinetic kinergy;

  T[t,sp]  = stiff pendulum translational kinetic kinergy

Rotational energy of STIFFPEN

>    Tr[sp] := kinetic_energy( 'rotation', J[sp], qd[2] );

Tr[sp] := 1/2*J[sp]*diff(alpha(t),t)^2

Translational energy of STIFFPEN

>    v[sp] := n_norm( [ xd[sp], yd[sp] ], 2 ):
Tt[sp] := kinetic_energy( 'translation', M[sp], v[sp] );

Tt[sp] := 1/2*M[sp]*((diff(x[c](t),t)-l[sp]*cos(alpha(t))*diff(alpha(t),t))^2+l[sp]^2*sin(alpha(t))^2*diff(alpha(t),t)^2)

Total kinetic energy of STIFFPEN

>    T[T_sp] := Tr[sp] + Tt[sp];

T[T_sp] := 1/2*J[sp]*diff(alpha(t),t)^2+1/2*M[sp]*((diff(x[c](t),t)-l[sp]*cos(alpha(t))*diff(alpha(t),t))^2+l[sp]^2*sin(alpha(t))^2*diff(alpha(t),t)^2)

The flexible pendulum, or FLEXPEN, motion is composed of one rotation and one translation, both directions are orthogonal.

T[r,fp] = FLEXPEN rotational kinetic kinergy;

  T[t,fp]  = FLEXPEN translational kinetic kinergy

Rotational energy of flexpen.

>    Tr[fp] := kinetic_energy( 'rotation', J[fp], qd[3] );

Tr[fp] := 1/2*J[fp]*diff(gamma(t),t)^2

Translational energy of flexpen

>    v[fp] := n_norm( [ xd[fp], yd[fp] ], 2 ):
Tt[fp] := kinetic_energy( 'translation', M[fp], v[fp] );

Tt[fp] := 1/2*M[fp]*((diff(x[c](t),t)-l[sp]*cos(alpha(t))*diff(alpha(t),t)-l[fp]*cos(gamma(t))*diff(gamma(t),t))^2+(-l[sp]*sin(alpha(t))*diff(alpha(t),t)-l[fp]*sin(gamma(t))*diff(gamma(t),t))^2)
Tt[fp] := 1/2*M[fp]*((diff(x[c](t),t)-l[sp]*cos(alpha(t))*diff(alpha(t),t)-l[fp]*cos(gamma(t))*diff(gamma(t),t))^2+(-l[sp]*sin(alpha(t))*diff(alpha(t),t)-l[fp]*sin(gamma(t))*diff(gamma(t),t))^2)

Total kinetic energy of FLEXPEN.

>    T[T_fp] := Tr[fp] + Tt[fp];

T[T_fp] := 1/2*J[fp]*diff(gamma(t),t)^2+1/2*M[fp]*((diff(x[c](t),t)-l[sp]*cos(alpha(t))*diff(alpha(t),t)-l[fp]*cos(gamma(t))*diff(gamma(t),t))^2+(-l[sp]*sin(alpha(t))*diff(alpha(t),t)-l[fp]*sin(gamma(t...
T[T_fp] := 1/2*J[fp]*diff(gamma(t),t)^2+1/2*M[fp]*((diff(x[c](t),t)-l[sp]*cos(alpha(t))*diff(alpha(t),t)-l[fp]*cos(gamma(t))*diff(gamma(t),t))^2+(-l[sp]*sin(alpha(t))*diff(alpha(t),t)-l[fp]*sin(gamma(t...

Total kinetic energy of system

>    T_temp := T[T_c] + T[T_sp] + T[T_fp]:
T[T] := collect( T_temp, diff );

T[T] := (1/2*M[c]+1/2*M[sp]+1/2*M[fp])*diff(x[c](t),t)^2+((-M[sp]*l[sp]*cos(alpha(t))-M[fp]*l[sp]*cos(alpha(t)))*diff(alpha(t),t)-M[fp]*l[fp]*cos(gamma(t))*diff(gamma(t),t))*diff(x[c](t),t)+(1/2*M[sp]*...
T[T] := (1/2*M[c]+1/2*M[sp]+1/2*M[fp])*diff(x[c](t),t)^2+((-M[sp]*l[sp]*cos(alpha(t))-M[fp]*l[sp]*cos(alpha(t)))*diff(alpha(t),t)-M[fp]*l[fp]*cos(gamma(t))*diff(gamma(t),t))*diff(x[c](t),t)+(1/2*M[sp]*...
T[T] := (1/2*M[c]+1/2*M[sp]+1/2*M[fp])*diff(x[c](t),t)^2+((-M[sp]*l[sp]*cos(alpha(t))-M[fp]*l[sp]*cos(alpha(t)))*diff(alpha(t),t)-M[fp]*l[fp]*cos(gamma(t))*diff(gamma(t),t))*diff(x[c](t),t)+(1/2*M[sp]*...
T[T] := (1/2*M[c]+1/2*M[sp]+1/2*M[fp])*diff(x[c](t),t)^2+((-M[sp]*l[sp]*cos(alpha(t))-M[fp]*l[sp]*cos(alpha(t)))*diff(alpha(t),t)-M[fp]*l[fp]*cos(gamma(t))*diff(gamma(t),t))*diff(x[c](t),t)+(1/2*M[sp]*...
T[T] := (1/2*M[c]+1/2*M[sp]+1/2*M[fp])*diff(x[c](t),t)^2+((-M[sp]*l[sp]*cos(alpha(t))-M[fp]*l[sp]*cos(alpha(t)))*diff(alpha(t),t)-M[fp]*l[fp]*cos(gamma(t))*diff(gamma(t),t))*diff(x[c](t),t)+(1/2*M[sp]*...

Generalized Forces: Q[i] 's

The non-conservative forces corresponding to the generalized coordinates are: F[c]  and the viscous damping forces, where

F[c]  = linear force generated by the motorized cart (e.g.  IP02)

B[sp] = STIFFPEN viscous friction torque coefficient (a.k.a. viscous damping)

B[fp] = FLEXPEN viscous friction torque coefficient (a.k.a. viscous damping)

B[eq] = cart viscous damping force coefficient

B[eq]  = 0 and B[sp]  = 0 and B[fp]  = 0 if all viscous dampings are neglected

>    B[eq] := 0: B[sp] := 0: B[fp] := 0:

Q[i] = generalized force applied on generalized coordinate q[i]
the Coulomb friction is neglected, as well as the force due to the pendulum and acting on the linear cart

>    Q[1] := F[c] - B[eq] * qd[1];

Q[1] := F[c]

>    Q[2] := - B[sp] * qd[2];

Q[2] := 0

>    Q[3] := - B[fp] * qd[3];

Q[3] := 0

F[c]  = Linear force produced by the motor at the motor pinion: cart driving force

rm: comment the following line out if u = F[c] , uncomment it if u = V[m]

>    #F[c] := eta[g] * K[g] * eta[m] * K[t] * ( V[m] * r[mp] - K[g] * K[m] * qd[1] ) / R[m] / r[mp]^2;

>    Q := [ seq( Q[i], i=1..Nq ) ];

Q := [F[c], 0, 0]

Euler-Lagrange's Equations

For a N -DOF system, the Lagrange's equations can be written:

Diff(Diff(L,qdot[i]),t)-Diff(L,q[i]) = Q[i]   for i = 1  through N

where:

  Q[i] 's are special combinations of external forces and called the generalized forces,

  q[1] , ..., q[N] , are N independent coordinates chosen to describe the system and called the generalized coordinates ,

 and L  is the Lagrangian  of the system.

L  is defined by:

L = T - U

where T  is the total kinetic energy of the system and U  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 ) } );

EOM_orig := [(M[sp]*l[sp]*sin(alpha(t))+M[fp]*l[sp]*sin(alpha(t)))*diff(alpha(t),t)^2+M[fp]*l[fp]*sin(gamma(t))*diff(gamma(t),t)^2+(M[fp]+M[c]+M[sp])*diff(x[c](t),`$`(t,2))+(-M[sp]*l[sp]*cos(alpha(t))-...
EOM_orig := [(M[sp]*l[sp]*sin(alpha(t))+M[fp]*l[sp]*sin(alpha(t)))*diff(alpha(t),t)^2+M[fp]*l[fp]*sin(gamma(t))*diff(gamma(t),t)^2+(M[fp]+M[c]+M[sp])*diff(x[c](t),`$`(t,2))+(-M[sp]*l[sp]*cos(alpha(t))-...
EOM_orig := [(M[sp]*l[sp]*sin(alpha(t))+M[fp]*l[sp]*sin(alpha(t)))*diff(alpha(t),t)^2+M[fp]*l[fp]*sin(gamma(t))*diff(gamma(t),t)^2+(M[fp]+M[c]+M[sp])*diff(x[c](t),`$`(t,2))+(-M[sp]*l[sp]*cos(alpha(t))-...
EOM_orig := [(M[sp]*l[sp]*sin(alpha(t))+M[fp]*l[sp]*sin(alpha(t)))*diff(alpha(t),t)^2+M[fp]*l[fp]*sin(gamma(t))*diff(gamma(t),t)^2+(M[fp]+M[c]+M[sp])*diff(x[c](t),`$`(t,2))+(-M[sp]*l[sp]*cos(alpha(t))-...
EOM_orig := [(M[sp]*l[sp]*sin(alpha(t))+M[fp]*l[sp]*sin(alpha(t)))*diff(alpha(t),t)^2+M[fp]*l[fp]*sin(gamma(t))*diff(gamma(t),t)^2+(M[fp]+M[c]+M[sp])*diff(x[c](t),`$`(t,2))+(-M[sp]*l[sp]*cos(alpha(t))-...
EOM_orig := [(M[sp]*l[sp]*sin(alpha(t))+M[fp]*l[sp]*sin(alpha(t)))*diff(alpha(t),t)^2+M[fp]*l[fp]*sin(gamma(t))*diff(gamma(t),t)^2+(M[fp]+M[c]+M[sp])*diff(x[c](t),`$`(t,2))+(-M[sp]*l[sp]*cos(alpha(t))-...
EOM_orig := [(M[sp]*l[sp]*sin(alpha(t))+M[fp]*l[sp]*sin(alpha(t)))*diff(alpha(t),t)^2+M[fp]*l[fp]*sin(gamma(t))*diff(gamma(t),t)^2+(M[fp]+M[c]+M[sp])*diff(x[c](t),`$`(t,2))+(-M[sp]*l[sp]*cos(alpha(t))-...
EOM_orig := [(M[sp]*l[sp]*sin(alpha(t))+M[fp]*l[sp]*sin(alpha(t)))*diff(alpha(t),t)^2+M[fp]*l[fp]*sin(gamma(t))*diff(gamma(t),t)^2+(M[fp]+M[c]+M[sp])*diff(x[c](t),`$`(t,2))+(-M[sp]*l[sp]*cos(alpha(t))-...

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 );

EOM_states := [(M[sp]*l[sp]*sin(X[2])+M[fp]*l[sp]*sin(X[2]))*X[5]^2+M[fp]*l[fp]*sin(X[3])*X[6]^2+(M[fp]+M[c]+M[sp])*Xd[4]+(-M[sp]*l[sp]*cos(X[2])-M[fp]*l[sp]*cos(X[2]))*Xd[5]-M[fp]*l[fp]*cos(X[3])*Xd[6...
EOM_states := [(M[sp]*l[sp]*sin(X[2])+M[fp]*l[sp]*sin(X[2]))*X[5]^2+M[fp]*l[fp]*sin(X[3])*X[6]^2+(M[fp]+M[c]+M[sp])*Xd[4]+(-M[sp]*l[sp]*cos(X[2])-M[fp]*l[sp]*cos(X[2]))*Xd[5]-M[fp]*l[fp]*cos(X[3])*Xd[6...
EOM_states := [(M[sp]*l[sp]*sin(X[2])+M[fp]*l[sp]*sin(X[2]))*X[5]^2+M[fp]*l[fp]*sin(X[3])*X[6]^2+(M[fp]+M[c]+M[sp])*Xd[4]+(-M[sp]*l[sp]*cos(X[2])-M[fp]*l[sp]*cos(X[2]))*Xd[5]-M[fp]*l[fp]*cos(X[3])*Xd[6...
EOM_states := [(M[sp]*l[sp]*sin(X[2])+M[fp]*l[sp]*sin(X[2]))*X[5]^2+M[fp]*l[fp]*sin(X[3])*X[6]^2+(M[fp]+M[c]+M[sp])*Xd[4]+(-M[sp]*l[sp]*cos(X[2])-M[fp]*l[sp]*cos(X[2]))*Xd[5]-M[fp]*l[fp]*cos(X[3])*Xd[6...
EOM_states := [(M[sp]*l[sp]*sin(X[2])+M[fp]*l[sp]*sin(X[2]))*X[5]^2+M[fp]*l[fp]*sin(X[3])*X[6]^2+(M[fp]+M[c]+M[sp])*Xd[4]+(-M[sp]*l[sp]*cos(X[2])-M[fp]*l[sp]*cos(X[2]))*Xd[5]-M[fp]*l[fp]*cos(X[3])*Xd[6...
EOM_states := [(M[sp]*l[sp]*sin(X[2])+M[fp]*l[sp]*sin(X[2]))*X[5]^2+M[fp]*l[fp]*sin(X[3])*X[6]^2+(M[fp]+M[c]+M[sp])*Xd[4]+(-M[sp]*l[sp]*cos(X[2])-M[fp]*l[sp]*cos(X[2]))*Xd[5]-M[fp]*l[fp]*cos(X[3])*Xd[6...

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
  for k from 1 to Np do
    EOM_ser[i] := subsop( 1 = convert( series( op( 1, EOM_ser[i] ), X[ k+1 ] ), polynom ), EOM_ser[i] );
  end do:
  EOM_ser[i] := simplify( EOM_ser[i] );
end do:

>    EOM_ser;

[-Xd[5]*M[sp]*l[sp]-Xd[5]*M[fp]*l[sp]+X[5]^2*X[2]*M[sp]*l[sp]+X[5]^2*X[2]*M[fp]*l[sp]+Xd[4]*M[fp]+Xd[4]*M[c]+Xd[4]*M[sp]-M[fp]*l[fp]*Xd[6]+M[fp]*l[fp]*X[6]^2*X[3] = U[1], -Xd[4]*M[sp]*l[sp]-Xd[4]*M[fp]...
[-Xd[5]*M[sp]*l[sp]-Xd[5]*M[fp]*l[sp]+X[5]^2*X[2]*M[sp]*l[sp]+X[5]^2*X[2]*M[fp]*l[sp]+Xd[4]*M[fp]+Xd[4]*M[c]+Xd[4]*M[sp]-M[fp]*l[fp]*Xd[6]+M[fp]*l[fp]*X[6]^2*X[3] = U[1], -Xd[4]*M[sp]*l[sp]-Xd[4]*M[fp]...
[-Xd[5]*M[sp]*l[sp]-Xd[5]*M[fp]*l[sp]+X[5]^2*X[2]*M[sp]*l[sp]+X[5]^2*X[2]*M[fp]*l[sp]+Xd[4]*M[fp]+Xd[4]*M[c]+Xd[4]*M[sp]-M[fp]*l[fp]*Xd[6]+M[fp]*l[fp]*X[6]^2*X[3] = U[1], -Xd[4]*M[sp]*l[sp]-Xd[4]*M[fp]...
[-Xd[5]*M[sp]*l[sp]-Xd[5]*M[fp]*l[sp]+X[5]^2*X[2]*M[sp]*l[sp]+X[5]^2*X[2]*M[fp]*l[sp]+Xd[4]*M[fp]+Xd[4]*M[c]+Xd[4]*M[sp]-M[fp]*l[fp]*Xd[6]+M[fp]*l[fp]*X[6]^2*X[3] = U[1], -Xd[4]*M[sp]*l[sp]-Xd[4]*M[fp]...
[-Xd[5]*M[sp]*l[sp]-Xd[5]*M[fp]*l[sp]+X[5]^2*X[2]*M[sp]*l[sp]+X[5]^2*X[2]*M[fp]*l[sp]+Xd[4]*M[fp]+Xd[4]*M[c]+Xd[4]*M[sp]-M[fp]*l[fp]*Xd[6]+M[fp]*l[fp]*X[6]^2*X[3] = U[1], -Xd[4]*M[sp]*l[sp]-Xd[4]*M[fp]...

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;

F[i] = Matrix(%id = 51054820)

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_XUzero, Fi_lin[ i, k ] );
  end do;
end do:

>    'F_lin' = Fi_lin;

F_lin = Matrix(%id = 51056196)

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[m] }:
subs_U_rev := { U[1] = F[c] }:

>    eom_collect_list := { seq( diff( q[i], t ), i=1..Nq ), seq( q[i], i=1..Nq ) };

eom_collect_list := {x[c](t), alpha(t), gamma(t), diff(x[c](t),t), diff(alpha(t),t), diff(gamma(t),t)}

>    if not assigned( J[sp] ) then
  eom_collect_list := eom_collect_list union { J[sp] };
end if;

eom_collect_list := {x[c](t), alpha(t), gamma(t), diff(x[c](t),t), diff(alpha(t),t), diff(gamma(t),t), J[sp]}

>    if not assigned( J[fp] ) then
  eom_collect_list := eom_collect_list union { J[fp] };
end if;

eom_collect_list := {J[fp], x[c](t), alpha(t), gamma(t), diff(x[c](t),t), diff(alpha(t),t), diff(gamma(t),t), J[sp]}

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 );

>    #for i from 1 to Nq do
#  Xd_nl[i+Nq] := simplify( Xd[i+Nq] ):
#  unassign( 'Xd[i+Nq]' ):
#end do:

pretty display w.r.t. the named system states

>    #for i from 1 to Nq do
#  Xd_nl[i+Nq] := simplify( subs( subs_U_rev, subs_Xq_rev, subs_Xqd_rev, #Xd_nl[i+Nq] ) ):
#end do:

>    #for i from 1 to Nq do
#  diff( qd[i], t ) = collect( Xd_nl[i+Nq], eom_collect_list ):
#end do:

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_avsq_list := { X[Nq+2]^2 = 0, X[Nq+3]^2 = 0 }:

>    for i from 1 to Nq do
  Xd[i+Nq] := subs( subs_avsq_list, Xd[i+Nq] );
end:

pretty display w.r.t. the named system states

>    for i from 1 to Nq do
  diff( qd[i], t ) = collect( subs( subs_U_rev, subs_Xq_rev, subs_Xqd_rev, Xd[i+Nq] ), eom_collect_list );
end do;

diff(x[c](t),`$`(t,2)) = -(((2*M[sp]*l[sp]^2*M[fp]*g+M[fp]^2*l[sp]^2*g+M[sp]^2*l[sp]^2*g)*alpha(t)+M[sp]*l[sp]^2*F[c]+M[fp]*l[sp]^2*F[c]+J[sp]*F[c])*J[fp]+(-M[fp]^2*l[fp]^2*gamma(t)^2*l[sp]^2*F[c]+(-l[...
diff(x[c](t),`$`(t,2)) = -(((2*M[sp]*l[sp]^2*M[fp]*g+M[fp]^2*l[sp]^2*g+M[sp]^2*l[sp]^2*g)*alpha(t)+M[sp]*l[sp]^2*F[c]+M[fp]*l[sp]^2*F[c]+J[sp]*F[c])*J[fp]+(-M[fp]^2*l[fp]^2*gamma(t)^2*l[sp]^2*F[c]+(-l[...
diff(x[c](t),`$`(t,2)) = -(((2*M[sp]*l[sp]^2*M[fp]*g+M[fp]^2*l[sp]^2*g+M[sp]^2*l[sp]^2*g)*alpha(t)+M[sp]*l[sp]^2*F[c]+M[fp]*l[sp]^2*F[c]+J[sp]*F[c])*J[fp]+(-M[fp]^2*l[fp]^2*gamma(t)^2*l[sp]^2*F[c]+(-l[...
diff(x[c](t),`$`(t,2)) = -(((2*M[sp]*l[sp]^2*M[fp]*g+M[fp]^2*l[sp]^2*g+M[sp]^2*l[sp]^2*g)*alpha(t)+M[sp]*l[sp]^2*F[c]+M[fp]*l[sp]^2*F[c]+J[sp]*F[c])*J[fp]+(-M[fp]^2*l[fp]^2*gamma(t)^2*l[sp]^2*F[c]+(-l[...
diff(x[c](t),`$`(t,2)) = -(((2*M[sp]*l[sp]^2*M[fp]*g+M[fp]^2*l[sp]^2*g+M[sp]^2*l[sp]^2*g)*alpha(t)+M[sp]*l[sp]^2*F[c]+M[fp]*l[sp]^2*F[c]+J[sp]*F[c])*J[fp]+(-M[fp]^2*l[fp]^2*gamma(t)^2*l[sp]^2*F[c]+(-l[...
diff(x[c](t),`$`(t,2)) = -(((2*M[sp]*l[sp]^2*M[fp]*g+M[fp]^2*l[sp]^2*g+M[sp]^2*l[sp]^2*g)*alpha(t)+M[sp]*l[sp]^2*F[c]+M[fp]*l[sp]^2*F[c]+J[sp]*F[c])*J[fp]+(-M[fp]^2*l[fp]^2*gamma(t)^2*l[sp]^2*F[c]+(-l[...
diff(x[c](t),`$`(t,2)) = -(((2*M[sp]*l[sp]^2*M[fp]*g+M[fp]^2*l[sp]^2*g+M[sp]^2*l[sp]^2*g)*alpha(t)+M[sp]*l[sp]^2*F[c]+M[fp]*l[sp]^2*F[c]+J[sp]*F[c])*J[fp]+(-M[fp]^2*l[fp]^2*gamma(t)^2*l[sp]^2*F[c]+(-l[...

diff(alpha(t),`$`(t,2)) = -(((M[fp]^2*g+M[sp]*g*M[c]+M[sp]^2*g+2*M[fp]*g*M[sp]+M[fp]*g*M[c])*alpha(t)+F[c]*M[sp]+F[c]*M[fp])*J[fp]+((-M[fp]^3*l[fp]^2*g+M[fp]^2*l[fp]*K[s]-M[c]*M[fp]^2*l[fp]^2*g+M[sp]*M...
diff(alpha(t),`$`(t,2)) = -(((M[fp]^2*g+M[sp]*g*M[c]+M[sp]^2*g+2*M[fp]*g*M[sp]+M[fp]*g*M[c])*alpha(t)+F[c]*M[sp]+F[c]*M[fp])*J[fp]+((-M[fp]^3*l[fp]^2*g+M[fp]^2*l[fp]*K[s]-M[c]*M[fp]^2*l[fp]^2*g+M[sp]*M...
diff(alpha(t),`$`(t,2)) = -(((M[fp]^2*g+M[sp]*g*M[c]+M[sp]^2*g+2*M[fp]*g*M[sp]+M[fp]*g*M[c])*alpha(t)+F[c]*M[sp]+F[c]*M[fp])*J[fp]+((-M[fp]^3*l[fp]^2*g+M[fp]^2*l[fp]*K[s]-M[c]*M[fp]^2*l[fp]^2*g+M[sp]*M...
diff(alpha(t),`$`(t,2)) = -(((M[fp]^2*g+M[sp]*g*M[c]+M[sp]^2*g+2*M[fp]*g*M[sp]+M[fp]*g*M[c])*alpha(t)+F[c]*M[sp]+F[c]*M[fp])*J[fp]+((-M[fp]^3*l[fp]^2*g+M[fp]^2*l[fp]*K[s]-M[c]*M[fp]^2*l[fp]^2*g+M[sp]*M...
diff(alpha(t),`$`(t,2)) = -(((M[fp]^2*g+M[sp]*g*M[c]+M[sp]^2*g+2*M[fp]*g*M[sp]+M[fp]*g*M[c])*alpha(t)+F[c]*M[sp]+F[c]*M[fp])*J[fp]+((-M[fp]^3*l[fp]^2*g+M[fp]^2*l[fp]*K[s]-M[c]*M[fp]^2*l[fp]^2*g+M[sp]*M...
diff(alpha(t),`$`(t,2)) = -(((M[fp]^2*g+M[sp]*g*M[c]+M[sp]^2*g+2*M[fp]*g*M[sp]+M[fp]*g*M[c])*alpha(t)+F[c]*M[sp]+F[c]*M[fp])*J[fp]+((-M[fp]^3*l[fp]^2*g+M[fp]^2*l[fp]*K[s]-M[c]*M[fp]^2*l[fp]^2*g+M[sp]*M...

diff(gamma(t),`$`(t,2)) = ((l[sp]^2*M[fp]^2*g*l[fp]*M[c]+M[fp]*l[fp]*M[sp]^2*l[sp]^2*g+M[c]*M[fp]*l[fp]*M[sp]*l[sp]^2*g+M[fp]^3*l[fp]*l[sp]^2*g+2*M[fp]^2*l[fp]*M[sp]*l[sp]^2*g)*gamma(t)*alpha(t)^2+((M[...
diff(gamma(t),`$`(t,2)) = ((l[sp]^2*M[fp]^2*g*l[fp]*M[c]+M[fp]*l[fp]*M[sp]^2*l[sp]^2*g+M[c]*M[fp]*l[fp]*M[sp]*l[sp]^2*g+M[fp]^3*l[fp]*l[sp]^2*g+2*M[fp]^2*l[fp]*M[sp]*l[sp]^2*g)*gamma(t)*alpha(t)^2+((M[...
diff(gamma(t),`$`(t,2)) = ((l[sp]^2*M[fp]^2*g*l[fp]*M[c]+M[fp]*l[fp]*M[sp]^2*l[sp]^2*g+M[c]*M[fp]*l[fp]*M[sp]*l[sp]^2*g+M[fp]^3*l[fp]*l[sp]^2*g+2*M[fp]^2*l[fp]*M[sp]*l[sp]^2*g)*gamma(t)*alpha(t)^2+((M[...
diff(gamma(t),`$`(t,2)) = ((l[sp]^2*M[fp]^2*g*l[fp]*M[c]+M[fp]*l[fp]*M[sp]^2*l[sp]^2*g+M[c]*M[fp]*l[fp]*M[sp]*l[sp]^2*g+M[fp]^3*l[fp]*l[sp]^2*g+2*M[fp]^2*l[fp]*M[sp]*l[sp]^2*g)*gamma(t)*alpha(t)^2+((M[...
diff(gamma(t),`$`(t,2)) = ((l[sp]^2*M[fp]^2*g*l[fp]*M[c]+M[fp]*l[fp]*M[sp]^2*l[sp]^2*g+M[c]*M[fp]*l[fp]*M[sp]*l[sp]^2*g+M[fp]^3*l[fp]*l[sp]^2*g+2*M[fp]^2*l[fp]*M[sp]*l[sp]^2*g)*gamma(t)*alpha(t)^2+((M[...
diff(gamma(t),`$`(t,2)) = ((l[sp]^2*M[fp]^2*g*l[fp]*M[c]+M[fp]*l[fp]*M[sp]^2*l[sp]^2*g+M[c]*M[fp]*l[fp]*M[sp]*l[sp]^2*g+M[fp]^3*l[fp]*l[sp]^2*g+2*M[fp]^2*l[fp]*M[sp]*l[sp]^2*g)*gamma(t)*alpha(t)^2+((M[...
diff(gamma(t),`$`(t,2)) = ((l[sp]^2*M[fp]^2*g*l[fp]*M[c]+M[fp]*l[fp]*M[sp]^2*l[sp]^2*g+M[c]*M[fp]*l[fp]*M[sp]*l[sp]^2*g+M[fp]^3*l[fp]*l[sp]^2*g+2*M[fp]^2*l[fp]*M[sp]*l[sp]^2*g)*gamma(t)*alpha(t)^2+((M[...

Determine the System State-Space Matrices: A, B, C, and D

>    A_ss := Matrix( Nx, Nx ):

>    A_ss := deriveA( Xqdd, A_ss, Nq, subs_XUzero ): 'A' = A_ss;

A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)
A = Matrix(%id = 51056612)

>    B_ss := Matrix( Nx, Nu ):

>    B_ss := deriveB( Xqdd, B_ss, Nq, subs_XUzero ): 'B' = B_ss;

B = Matrix(%id = 135183400)
B = Matrix(%id = 135183400)
B = Matrix(%id = 135183400)
B = Matrix(%id = 135183400)
B = Matrix(%id = 135183400)
B = Matrix(%id = 135183400)
B = Matrix(%id = 135183400)
B = Matrix(%id = 135183400)
B = Matrix(%id = 135183400)
B = Matrix(%id = 135183400)

>    C_ss := IdentityMatrix( Ny, Nx );

C_ss := Matrix(%id = 135183832)

>    D_ss := Matrix( Ny, Nu, 0 );

D_ss := Matrix(%id = 135184248)

>   

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 := "FLEXPEN_ABCD_eqns.m":

unassign variables, when necessary, for those present in the "Matlab_Notations" substitution set

>    if assigned( B[eq] ) then
  unassign( 'B[eq]' );
end if;

>    if assigned( B[sp] ) then
  unassign( 'B[sp]' );
end if;

>    if assigned( J[sp] ) then
  unassign( 'J[sp]' );
end if;

>    if assigned( B[fp] ) then
  unassign( 'B[fp]' );
end if;

>    if assigned( J[fp] ) then
  unassign( 'J[fp]' );
end if;

substitution set containing a notation consistent with that used in the MATLAB design script(s)

>    Matlab_Notations := { M[c] = Mc, M[sp] = Msp, l[sp] = lsp, J[sp] = Jsp, B[sp] = Bsp,  M[fp] = Mfp, l[fp] = lfp, J[fp] = Jfp, B[fp] = Bfp, B[eq] = Beq, K[t] = Kt, K[m]=Km, K[g]=Kg, r[mp] = r_mp, eta[m] = Eff_m, eta[g] = Eff_g, R[m]=Rm, K[s] = Ks }:

>    Experiment_Name := "IP02+FLEXPEN System":

>    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 ):

>   

>   

>   

>   

>   

>   

>   

>   

>   

>   

Click here to go back to top: Description Section.