Dynamic Equations for the  

Quarter Car Model 

© 2003 Quanser Consulting Inc. 

URL: http://www.quanser.com 

NOTE: This worksheet presents the general dynamic equations modelling a quarter car suspension. 

 

Description 

This worksheet presents the general dynamic equations modelling the  quarter car suspension. 

Specifically, this worksheet is used to derive the state-space matrices for the Quanser's Quarter Car Suspension experiment. 

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

[HTM, deriveA, deriveB, deriveF, kinetic_energy, lagrange_equations, moment_of_inertia, n_norm, potential_energy, write_ABCD_to_Mfile]
[HTM, deriveA, deriveB, deriveF, kinetic_energy, lagrange_equations, moment_of_inertia, n_norm, potential_energy, write_ABCD_to_Mfile]
(2.1)
 

environment variable representing the order of series calculations 

> Order := 2:
 

Notations 

Generalized Coordinates: 's 

Figure 1: System Schematic 

Image 

The generalized coordinates are also called Lagrangian coordinates. 

=  the position of mass 1, the ground plane 

=  the position of mass 1 presenting a tire center of mass measured from is resting position 

position of mass 2, representing the car body measured from its resting position 

 

> q1 := [ x[1](t), x[2](t) ,z[r](t) , x[1](t)-z[r](t), x[2](t)-x[1](t)];
 

[x[1](t), x[2](t), z[r](t), `+`(x[1](t), `-`(z[r](t))), `+`(x[2](t), `-`(x[1](t)))] (3.1.1)
 

> q := [ x[1](t), x[2](t)  ];
 

[x[1](t), x[2](t)] (3.1.2)
 

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, e.g. position and angular velocities 

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

[diff(x[1](t), t), diff(x[2](t), t)] (3.1.3)
 

> Nq := nops( q ):
 

> u := z[r](t);
 

z[r](t) (3.1.4)
 

> ud:= diff(u,t);
 

diff(z[r](t), t) (3.1.5)
 

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

 

{x[1](t) = X[1], x[2](t) = X[2]}
{diff(x[1](t), t) = X[3], diff(x[2](t), t) = X[4]} (3.2.1)
 

Substitution set for the input(s). 

set the input to be the current cmd to the two motors and gravity 

> subs_U := { ud = U[1],F[c]=U[2],g=U[3]};
 

{g = U[3], F[c] = U[2], diff(z[r](t), t) = U[1]} (3.2.2)
 

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

> Nu := nops( subs_U );
 

3 (3.2.3)
 

substitution set for the position states' second time derivatives 

> subs_Xqdd := { seq( diff( q[i], t$2 ) =  Xd[i+Nq], i=1..Nq ) };
 

{diff(diff(x[1](t), t), t) = Xd[3], diff(diff(x[2](t), t), t) = Xd[4]} (3.2.4)
 

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

[Xd[3], Xd[4]] (3.2.5)
 

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

{U[1] = 0, U[2] = 0, U[3] = 0, X[1] = 0, X[2] = 0, X[3] = 0, X[4] = 0} (3.2.6)
 

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. 

Vg[T] = Total Gravitational Potential Energy of the system 

> Vg[1] := potential_energy( 'gravity', M[us], g, q[1]);
 

`*`(M[us], `*`(g, `*`(x[1](t)))) (4.1.1)
 

> Vg[2] := potential_energy( 'gravity', M[s], g, q[2]);
 

`*`(M[s], `*`(g, `*`(x[2](t)))) (4.1.2)
 

> Vg[T] := Vg[1] + Vg[2] ;
 

`+`(`*`(M[us], `*`(g, `*`(x[1](t)))), `*`(M[s], `*`(g, `*`(x[2](t))))) (4.1.3)
 

K[s] = spring torsional stiffness 

Ve[1] = first spring potetnial energy 

> Ve[1] := potential_energy( 'spring', K[us], ( q[1] - u ) );
 

`+`(`*`(`/`(1, 2), `*`(K[us], `*`(`^`(`+`(x[1](t), `-`(z[r](t))), 2))))) (4.1.4)
 

Ve[2] = second spring potetnial energy 

> Ve[2] := potential_energy( 'spring', K[s], ( q[2] - q[1] ) );
 

`+`(`*`(`/`(1, 2), `*`(K[s], `*`(`^`(`+`(x[2](t), `-`(x[1](t))), 2))))) (4.1.5)
 

Ve[T] = Total Elastic Potential Energy of the system 

> Ve[T] := Ve[1] + Ve[2];
 

`+`(`*`(`/`(1, 2), `*`(K[us], `*`(`^`(`+`(x[1](t), `-`(z[r](t))), 2)))), `*`(`/`(1, 2), `*`(K[s], `*`(`^`(`+`(x[2](t), `-`(x[1](t))), 2))))) (4.1.6)
 

V[T] = Total Potential Energy of the system 

> V[T] := simplify( Ve[T] + Vg[T] );
 

`+`(`*`(`/`(1, 2), `*`(K[us], `*`(`^`(x[1](t), 2)))), `-`(`*`(K[us], `*`(x[1](t), `*`(z[r](t))))), `*`(`/`(1, 2), `*`(K[us], `*`(`^`(z[r](t), 2)))), `*`(`/`(1, 2), `*`(K[s], `*`(`^`(x[2](t), 2)))), `-...
`+`(`*`(`/`(1, 2), `*`(K[us], `*`(`^`(x[1](t), 2)))), `-`(`*`(K[us], `*`(x[1](t), `*`(z[r](t))))), `*`(`/`(1, 2), `*`(K[us], `*`(`^`(z[r](t), 2)))), `*`(`/`(1, 2), `*`(K[s], `*`(`^`(x[2](t), 2)))), `-...
(4.1.7)
 

 

Total Kinetic Energy:  

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

= kinetic energy of the the Mass 0, the lower road plane 

qd[1] = road plane linear velocity 

> Tr[1] := kinetic_energy( 'translation', M[us], qd[1] );
 

`+`(`*`(`/`(1, 2), `*`(M[us], `*`(`^`(diff(x[1](t), t), 2))))) (4.2.1)
 

= kinetic energy of the the Mass 1, the middle tire platform 

qd[2] = Tire plane linear velocity 

> Tr[2] := kinetic_energy( 'translation', M[s], qd[2] );
 

`+`(`*`(`/`(1, 2), `*`(M[s], `*`(`^`(diff(x[2](t), t), 2))))) (4.2.2)
 

=  kinetic energy of the the Mass 2, the top car platform 

qd[3] = Tire plane linear velocity 

= Total Kinetic Energy of the system 

> T[T] := Tr[1] + Tr[2] ;
 

`+`(`*`(`/`(1, 2), `*`(M[us], `*`(`^`(diff(x[1](t), t), 2)))), `*`(`/`(1, 2), `*`(M[s], `*`(`^`(diff(x[2](t), t), 2))))) (4.2.3)
 

> T[T] := collect( T[T], diff );
 

`+`(`*`(`/`(1, 2), `*`(M[us], `*`(`^`(diff(x[1](t), t), 2)))), `*`(`/`(1, 2), `*`(M[s], `*`(`^`(diff(x[2](t), t), 2))))) (4.2.4)
 

Generalized Forces: 's 

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

= Current to motor 1, which generates a force on the bottom platform,  

= Current to motor 2, which generates a force between the top and middle platform.,  

= tire dashpot viscous damping coefficient 

= suspension strut dashpot viscous damping coefficient 

Q[i] = generalized force applied on generalized coordinate q[i] 

> Q[1] := -F[c] + B[s]* (qd[2]-qd[1]) - B[us]*(qd[1] - ud);
 

>
 

`+`(`-`(F[c]), `*`(B[s], `*`(`+`(diff(x[2](t), t), `-`(diff(x[1](t), t))))), `-`(`*`(B[us], `*`(`+`(diff(x[1](t), t), `-`(diff(z[r](t), t))))))) (5.1)
 

> Q[2] := F[c] - B[s] * (qd[2]-qd[1]);
 

`+`(F[c], `-`(`*`(B[s], `*`(`+`(diff(x[2](t), t), `-`(diff(x[1](t), t))))))) (5.2)
 

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

[`+`(`-`(F[c]), `*`(B[s], `*`(`+`(diff(x[2](t), t), `-`(diff(x[1](t), t))))), `-`(`*`(B[us], `*`(`+`(diff(x[1](t), t), `-`(diff(z[r](t), t))))))), `+`(F[c], `-`(`*`(B[s], `*`(`+`(diff(x[2](t), t), `-`...
[`+`(`-`(F[c]), `*`(B[s], `*`(`+`(diff(x[2](t), t), `-`(diff(x[1](t), t))))), `-`(`*`(B[us], `*`(`+`(diff(x[1](t), t), `-`(diff(z[r](t), t))))))), `+`(F[c], `-`(`*`(B[s], `*`(`+`(diff(x[2](t), t), `-`...
(5.3)
 

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

[`+`(`*`(M[us], `*`(diff(diff(x[1](t), t), t))), `*`(K[us], `*`(x[1](t))), `-`(`*`(K[us], `*`(z[r](t)))), `-`(`*`(K[s], `*`(x[2](t)))), `*`(K[s], `*`(x[1](t))), `*`(M[us], `*`(g))) = `+`(`-`(F[c]), `*...
[`+`(`*`(M[us], `*`(diff(diff(x[1](t), t), t))), `*`(K[us], `*`(x[1](t))), `-`(`*`(K[us], `*`(z[r](t)))), `-`(`*`(K[s], `*`(x[2](t)))), `*`(K[s], `*`(x[1](t))), `*`(M[us], `*`(g))) = `+`(`-`(F[c]), `*...
[`+`(`*`(M[us], `*`(diff(diff(x[1](t), t), t))), `*`(K[us], `*`(x[1](t))), `-`(`*`(K[us], `*`(z[r](t)))), `-`(`*`(K[s], `*`(x[2](t)))), `*`(K[s], `*`(x[1](t))), `*`(M[us], `*`(g))) = `+`(`-`(F[c]), `*...
(6.1)
 

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

[`+`(`*`(M[us], `*`(diff(diff(x[1](t), t), t))), `*`(`+`(K[us], K[s]), `*`(x[1](t))), `*`(M[us], `*`(g)), `-`(`*`(K[us], `*`(z[r](t)))), `-`(`*`(K[s], `*`(x[2](t))))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[us...
[`+`(`*`(M[us], `*`(diff(diff(x[1](t), t), t))), `*`(`+`(K[us], K[s]), `*`(x[1](t))), `*`(M[us], `*`(g)), `-`(`*`(K[us], `*`(z[r](t)))), `-`(`*`(K[s], `*`(x[2](t))))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[us...
[`+`(`*`(M[us], `*`(diff(diff(x[1](t), t), t))), `*`(`+`(K[us], K[s]), `*`(x[1](t))), `*`(M[us], `*`(g)), `-`(`*`(K[us], `*`(z[r](t)))), `-`(`*`(K[s], `*`(x[2](t))))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[us...
(6.2)
 

> EOM_New := collect( EOM_orig,{seq( diff( q[i], t$2 ), i=1..Nq ), seq( diff( q[i], t ), i=1..Nq ), K[s],K[us]},distributed,factor);
 

[`+`(`*`(`+`(x[1](t), `-`(x[2](t))), `*`(K[s])), `*`(M[us], `*`(g)), `*`(`+`(x[1](t), `-`(z[r](t))), `*`(K[us])), `*`(M[us], `*`(diff(diff(x[1](t), t), t)))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[us])), `*`(...
[`+`(`*`(`+`(x[1](t), `-`(x[2](t))), `*`(K[s])), `*`(M[us], `*`(g)), `*`(`+`(x[1](t), `-`(z[r](t))), `*`(K[us])), `*`(M[us], `*`(diff(diff(x[1](t), t), t)))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[us])), `*`(...
[`+`(`*`(`+`(x[1](t), `-`(x[2](t))), `*`(K[s])), `*`(M[us], `*`(g)), `*`(`+`(x[1](t), `-`(z[r](t))), `*`(K[us])), `*`(M[us], `*`(diff(diff(x[1](t), t), t)))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[us])), `*`(...
(6.3)
 

> subs_Xq := { seq( q1[i] = X[i-3], i=4..5 ), x[1](t)-x[2](t)=-X[2]   };
 

{`+`(x[1](t), `-`(x[2](t))) = `+`(`-`(X[2])), `+`(x[1](t), `-`(z[r](t))) = X[1], `+`(x[2](t), `-`(x[1](t))) = X[2]} (6.4)
 

> EOM_states := subs( subs_Xq, subs( subs_Xqdd, EOM_New ) );
 

[`+`(`-`(`*`(X[2], `*`(K[s]))), `*`(M[us], `*`(g)), `*`(X[1], `*`(K[us])), `*`(M[us], `*`(Xd[3]))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[us])), `*`(diff(x[1](t), t))), `-`(F[c]), `*`(B[s], `*`(diff(x[2](t), ...
[`+`(`-`(`*`(X[2], `*`(K[s]))), `*`(M[us], `*`(g)), `*`(X[1], `*`(K[us])), `*`(M[us], `*`(Xd[3]))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[us])), `*`(diff(x[1](t), t))), `-`(F[c]), `*`(B[s], `*`(diff(x[2](t), ...
(6.5)
 

> EOM_states := subs( subs_Xqd, subs_U, EOM_states );
 

[`+`(`-`(`*`(X[2], `*`(K[s]))), `*`(M[us], `*`(U[3])), `*`(X[1], `*`(K[us])), `*`(M[us], `*`(Xd[3]))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[us])), `*`(X[3])), `-`(U[2]), `*`(B[s], `*`(X[4])), `*`(B[us], `*`(...
[`+`(`-`(`*`(X[2], `*`(K[s]))), `*`(M[us], `*`(U[3])), `*`(X[1], `*`(K[us])), `*`(M[us], `*`(Xd[3]))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[us])), `*`(X[3])), `-`(U[2]), `*`(B[s], `*`(X[4])), `*`(B[us], `*`(...
(6.6)
 

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

Additional Insight: Eliminating Gravity Terms 

> xeq := [x[eq1](t), x[eq2](t)];
 

[x[eq1](t), x[eq2](t)] (6.1.1)
 

> seq(  q[i] = xeq[i] , i=1..Nq );
 

x[1](t) = x[eq1](t), x[2](t) = x[eq2](t) (6.1.2)
 

> EOM_NoGrav := simplify( subs( [seq(  q[i] = xeq[i] , i=1..Nq), seq( diff( q[i], t) = 0 , i=1..Nq),F[c]=0,z[r](t)=0] ,  EOM_New  ) );
 

[`+`(`*`(K[s], `*`(x[eq1](t))), `-`(`*`(K[s], `*`(x[eq2](t)))), `*`(M[us], `*`(g)), `*`(x[eq1](t), `*`(K[us])), `*`(M[us], `*`(diff(0, t)))) = `*`(B[us], `*`(diff(0, t))), `+`(`*`(M[s], `*`(g)), `*`(K...
[`+`(`*`(K[s], `*`(x[eq1](t))), `-`(`*`(K[s], `*`(x[eq2](t)))), `*`(M[us], `*`(g)), `*`(x[eq1](t), `*`(K[us])), `*`(M[us], `*`(diff(0, t)))) = `*`(B[us], `*`(diff(0, t))), `+`(`*`(M[s], `*`(g)), `*`(K...
(6.1.3)
 

> xeq := solve( EOM_NoGrav, xeq );
 

[[x[eq1](t) = `+`(`-`(`/`(`*`(g, `*`(`+`(M[us], M[s]))), `*`(K[us])))), x[eq2](t) = `+`(`-`(`/`(`*`(g, `*`(`+`(`*`(M[s], `*`(K[us])), `*`(K[s], `*`(M[us])), `*`(K[s], `*`(M[s]))))), `*`(K[us], `*`(K[s... (6.1.4)
 

> x[eq1] := -(M[s]+M[us])*g/K[us] ;                                                                                                     x[eq2] := -M[s]*g/K[s]-(M[s]+M[us])*g/K[us] ;
 

 

`+`(`-`(`/`(`*`(g, `*`(`+`(M[us], M[s]))), `*`(K[us]))))
`+`(`-`(`/`(`*`(M[s], `*`(g)), `*`(K[s]))), `-`(`/`(`*`(g, `*`(`+`(M[us], M[s]))), `*`(K[us])))) (6.1.5)
 

> z_var := [ z[us](t), z[s](t) ];                                                                                              subs_zx := { x[1](t) = z_var[1]+x[eq1] ,     x[2](t) = z_var[2]+x[eq2] };
subs_zxd := { seq( diff( q[i], t) = diff( z_var[i], t) , i=1..Nq ) };                                                       subs_zxdd := { seq( diff( q[i], t$2 ) =  diff( z_var[i], t$2), i=1..Nq ) };
 

 

 

 

[z[us](t), z[s](t)]
{x[1](t) = `+`(z[us](t), `-`(`/`(`*`(g, `*`(`+`(M[us], M[s]))), `*`(K[us])))), x[2](t) = `+`(z[s](t), `-`(`/`(`*`(M[s], `*`(g)), `*`(K[s]))), `-`(`/`(`*`(g, `*`(`+`(M[us], M[s]))), `*`(K[us]))))}
{diff(x[1](t), t) = diff(z[us](t), t), diff(x[2](t), t) = diff(z[s](t), t)}
{diff(diff(x[1](t), t), t) = diff(diff(z[us](t), t), t), diff(diff(x[2](t), t), t) = diff(diff(z[s](t), t), t)} (6.1.6)
 

> EOM_noG := subs( subs_zxd, subs( subs_zxdd, EOM_orig ) );
 

[`+`(`*`(M[us], `*`(diff(diff(z[us](t), t), t))), `*`(`+`(K[us], K[s]), `*`(x[1](t))), `*`(M[us], `*`(g)), `-`(`*`(K[us], `*`(z[r](t)))), `-`(`*`(K[s], `*`(x[2](t))))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[u...
[`+`(`*`(M[us], `*`(diff(diff(z[us](t), t), t))), `*`(`+`(K[us], K[s]), `*`(x[1](t))), `*`(M[us], `*`(g)), `-`(`*`(K[us], `*`(z[r](t)))), `-`(`*`(K[s], `*`(x[2](t))))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[u...
[`+`(`*`(M[us], `*`(diff(diff(z[us](t), t), t))), `*`(`+`(K[us], K[s]), `*`(x[1](t))), `*`(M[us], `*`(g)), `-`(`*`(K[us], `*`(z[r](t)))), `-`(`*`(K[s], `*`(x[2](t))))) = `+`(`*`(`+`(`-`(B[s]), `-`(B[u...
(6.1.7)
 

> EOM_noG := subs( subs_zx, EOM_noG ) ;
 

[`+`(`*`(M[us], `*`(diff(diff(z[us](t), t), t))), `*`(`+`(K[us], K[s]), `*`(`+`(z[us](t), `-`(`/`(`*`(g, `*`(`+`(M[us], M[s]))), `*`(K[us])))))), `*`(M[us], `*`(g)), `-`(`*`(K[us], `*`(z[r](t)))), `-`...
[`+`(`*`(M[us], `*`(diff(diff(z[us](t), t), t))), `*`(`+`(K[us], K[s]), `*`(`+`(z[us](t), `-`(`/`(`*`(g, `*`(`+`(M[us], M[s]))), `*`(K[us])))))), `*`(M[us], `*`(g)), `-`(`*`(K[us], `*`(z[r](t)))), `-`...
[`+`(`*`(M[us], `*`(diff(diff(z[us](t), t), t))), `*`(`+`(K[us], K[s]), `*`(`+`(z[us](t), `-`(`/`(`*`(g, `*`(`+`(M[us], M[s]))), `*`(K[us])))))), `*`(M[us], `*`(g)), `-`(`*`(K[us], `*`(z[r](t)))), `-`...
[`+`(`*`(M[us], `*`(diff(diff(z[us](t), t), t))), `*`(`+`(K[us], K[s]), `*`(`+`(z[us](t), `-`(`/`(`*`(g, `*`(`+`(M[us], M[s]))), `*`(K[us])))))), `*`(M[us], `*`(g)), `-`(`*`(K[us], `*`(z[r](t)))), `-`...
(6.1.8)
 

> EOM_noG_S := collect( EOM_noG,{  K[s], K[us]},distributed,factor);
 

[`+`(`*`(`+`(z[us](t), `-`(z[s](t))), `*`(K[s])), `*`(`+`(z[us](t), `-`(z[r](t))), `*`(K[us])), `*`(M[us], `*`(diff(diff(z[us](t), t), t)))) = `+`(`-`(`*`(B[s], `*`(diff(z[us](t), t)))), `-`(`*`(diff(...
[`+`(`*`(`+`(z[us](t), `-`(z[s](t))), `*`(K[s])), `*`(`+`(z[us](t), `-`(z[r](t))), `*`(K[us])), `*`(M[us], `*`(diff(diff(z[us](t), t), t)))) = `+`(`-`(`*`(B[s], `*`(diff(z[us](t), t)))), `-`(`*`(diff(...
[`+`(`*`(`+`(z[us](t), `-`(z[s](t))), `*`(K[s])), `*`(`+`(z[us](t), `-`(z[r](t))), `*`(K[us])), `*`(M[us], `*`(diff(diff(z[us](t), t), t)))) = `+`(`-`(`*`(B[s], `*`(diff(z[us](t), t)))), `-`(`*`(diff(...
(6.1.9)
 

> subs_Zq := [ z_var[1]-z[r](t)=Z[1], z_var[2]-z_var[1]=Z[2], z_var[1]-z_var[2]=-Z[2]   ];                                          subs_Zqd := { seq( diff ( z_var, t )[i] = Z[i+Nq], i=1..Nq ) };                                                          subs_Zqdd := { seq( diff( z_var[i], t$2 ) =  Zd[i+Nq], i=1..Nq ) };
 

 

 

[`+`(z[us](t), `-`(z[r](t))) = Z[1], `+`(z[s](t), `-`(z[us](t))) = Z[2], `+`(z[us](t), `-`(z[s](t))) = `+`(`-`(Z[2]))]
{diff(z[s](t), t) = Z[4], diff(z[us](t), t) = Z[3]}
{diff(diff(z[s](t), t), t) = Zd[4], diff(diff(z[us](t), t), t) = Zd[3]} (6.1.10)
 

> EOM_states_noG := subs( subs_Zq, subs( subs_Zqdd, EOM_noG_S ) );
 

[`+`(`-`(`*`(Z[2], `*`(K[s]))), `*`(Z[1], `*`(K[us])), `*`(M[us], `*`(Zd[3]))) = `+`(`-`(`*`(B[s], `*`(diff(z[us](t), t)))), `-`(`*`(diff(z[us](t), t), `*`(B[us]))), `-`(F[c]), `*`(B[s], `*`(diff(z[s]...
[`+`(`-`(`*`(Z[2], `*`(K[s]))), `*`(Z[1], `*`(K[us])), `*`(M[us], `*`(Zd[3]))) = `+`(`-`(`*`(B[s], `*`(diff(z[us](t), t)))), `-`(`*`(diff(z[us](t), t), `*`(B[us]))), `-`(F[c]), `*`(B[s], `*`(diff(z[s]...
(6.1.11)
 

> EOM_states_noG := subs( subs_Zqd, subs_U, EOM_states_noG );
 

[`+`(`-`(`*`(Z[2], `*`(K[s]))), `*`(Z[1], `*`(K[us])), `*`(M[us], `*`(Zd[3]))) = `+`(`-`(`*`(B[s], `*`(Z[3]))), `-`(`*`(Z[3], `*`(B[us]))), `-`(U[2]), `*`(B[s], `*`(Z[4])), `*`(B[us], `*`(U[1]))), `+`...
[`+`(`-`(`*`(Z[2], `*`(K[s]))), `*`(Z[1], `*`(K[us])), `*`(M[us], `*`(Zd[3]))) = `+`(`-`(`*`(B[s], `*`(Z[3]))), `-`(`*`(Z[3], `*`(B[us]))), `-`(U[2]), `*`(B[s], `*`(Z[4])), `*`(B[us], `*`(U[1]))), `+`...
(6.1.12)
 

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:
 

Here, Fi shows a linear system 

> 'F[i]' = Fi;
 

F[i] = Matrix(%id = 38443588) (6.2.1)
 

Solving the Euler-Lagrange's Equations 

Reverse State Substitution for Pretty Display of the Solved EOM's 

only for pretty print 

> subs_Xq_rev := { seq( X[i] = q1[i+3], i=1..Nq ) };
subs_Xqd_rev := { seq( X[i+Nq] = qd[i], i=1..Nq ) };
 

 

{X[1] = `+`(x[1](t), `-`(z[r](t))), X[2] = `+`(x[2](t), `-`(x[1](t)))}
{X[3] = diff(x[1](t), t), X[4] = diff(x[2](t), t)} (7.1.1)
 

> #subs_U_rev := { U[1] = V[m] }:
subs_U_rev := { U[1] = ud,U[2] = F[c],U[3]=g};
 

{U[1] = diff(z[r](t), t), U[2] = F[c], U[3] = g} (7.1.2)
 

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

{g, F[c], diff(x[1](t), t), diff(x[2](t), t), diff(z[r](t), t), x[1](t), x[2](t)} (7.1.3)
 

Solution to the (Linear) EOM's 

Solve the (linear) equations of motion for the states' second time derivatives 

> Xqdd_solset_ser := solve( convert( EOM_states, set ), convert( Xqdd, set ) );
 

{Xd[3] = `+`(`-`(`/`(`*`(`+`(`-`(`*`(X[2], `*`(K[s]))), `*`(M[us], `*`(U[3])), `*`(X[1], `*`(K[us])), `-`(`*`(B[s], `*`(X[4]))), `*`(B[s], `*`(X[3])), `*`(X[3], `*`(B[us])), U[2], `-`(`*`(B[us], `*`(U...
{Xd[3] = `+`(`-`(`/`(`*`(`+`(`-`(`*`(X[2], `*`(K[s]))), `*`(M[us], `*`(U[3])), `*`(X[1], `*`(K[us])), `-`(`*`(B[s], `*`(X[4]))), `*`(B[s], `*`(X[3])), `*`(X[3], `*`(B[us])), U[2], `-`(`*`(B[us], `*`(U...
(7.2.1)
 

> assign( Xqdd_solset_ser );
 

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

> diff( qd[1], t ) = collect( subs( subs_U_rev, subs_Xq_rev, subs_Xqd_rev, Xd[Nq+1] ), eom_collect_list );
 

diff(diff(x[1](t), t), t) = `+`(`-`(g), `-`(`/`(`*`(F[c]), `*`(M[us]))), `-`(`/`(`*`(`+`(B[s], B[us]), `*`(diff(x[1](t), t))), `*`(M[us]))), `/`(`*`(B[s], `*`(diff(x[2](t), t))), `*`(M[us])), `/`(`*`(...
diff(diff(x[1](t), t), t) = `+`(`-`(g), `-`(`/`(`*`(F[c]), `*`(M[us]))), `-`(`/`(`*`(`+`(B[s], B[us]), `*`(diff(x[1](t), t))), `*`(M[us]))), `/`(`*`(B[s], `*`(diff(x[2](t), t))), `*`(M[us])), `/`(`*`(...
(7.2.2)
 

> diff( qd[2], t ) = collect( subs( subs_U_rev, subs_Xq_rev, subs_Xqd_rev, Xd[Nq+2] ), eom_collect_list);
 

diff(diff(x[2](t), t), t) = `+`(`-`(g), `/`(`*`(F[c]), `*`(M[s])), `/`(`*`(B[s], `*`(diff(x[1](t), t))), `*`(M[s])), `-`(`/`(`*`(B[s], `*`(diff(x[2](t), t))), `*`(M[s]))), `/`(`*`(K[s], `*`(x[1](t))),... (7.2.3)
 

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_ss[2,3]:=-1: 'A' = A_ss;
 

A = Matrix(%id = 150073664) (8.1)
 

> #lprint( A_ss );
 

The gravity has been eliminated from equations of motion. and Nu should be reduced by one. 

> Nu := Nu-1:
 

> B_ss := Matrix( Nx, Nu ):
 

> B_ss := deriveB( Xqdd, B_ss, Nq, subs_XUzero ): B_ss[1,1]:=-1: 'B' = B_ss;
 

B = Matrix(%id = 150073920) (8.2)
 

> #lprint( B_ss );
 

> #C_ss := IdentityMatrix( Ny, Nx ):
C_ss := IdentityMatrix( 4, 4 ):
'C' = C_ss;
 

C = Matrix(%id = 150090344) (8.3)
 

> #lprint( C_ss );
 

 

> #D_ss := Matrix( Ny, Nu, 0 ):
D_ss := Matrix( Nx, Nu, 0 ):
'D' = D_ss;
 

D = Matrix(%id = 150090496) (8.4)
 

> #lprint( D_ss );
 

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

> X_bar := Vector( [ z[s](t)-z[us](t),diff(z_var[2],t), z[us](t)-z[r](t), diff( z_var[1], t ) ] ):                                   U_bar := Vector( [ diff( z[r](t),t ) , F[c]  ] ):                                                                             Y_bar := Vector( [ z[s](t)-z[us](t) , diff( z[s](t),t$2) ] ):  'x'=X_bar, 'u'=U_bar, 'y'=Y_bar;
 

x = Vector[column](%id = 150107012), u = Vector[column](%id = 150162452), y = Vector[column](%id = 150178864) (9.1)
 

> T_ss := Matrix( Nx, Nx, [ [0, 1, 0, 0], [0, 0, 0, 1], [1, 0, 0, 0], [0, 0, 1, 0] ] );
 

`assign`(T_ss, Matrix(%id = 150211652)) (9.2)
 

> A_bar := Multiply( Multiply( T_ss, A_ss ), MatrixInverse( T_ss ) ): 'A'=A_bar;
 

A = Matrix(%id = 150212920) (9.3)
 

> B_bar := Multiply( T_ss, B_ss ): 'B'=B_bar;
 

B = Matrix(%id = 150213196) (9.4)
 

> C_bar := Matrix( 2, Nx , [ [ 1 , 0 , 0 , 0 ] ,  [ A_bar[2,1] , A_bar[2,2] , A_bar[2,3] , A_bar[2,4] ] ] ):
 

> D_bar := Matrix( 2, Nu, [ [0 ,0 ] , [ 0 , B_bar[2,2] ] ] ): 'C'=C_bar,'D'=D_bar;
 

C = Matrix(%id = 150213276), D = Matrix(%id = 150213376) (9.5)
 

> 'diff(x(t),t)=A*x(t)+B*u'; 'y=C*x(t)+Du';
 

> Characteristic := simplify(Determinant(S*IdentityMatrix( Nx, Nx )-A_bar)): 'det(s*I-A)'=Characteristic;
 

 

 

diff(x(t), t) = `+`(`*`(A, `*`(x(t))), `*`(B, `*`(u)))
y = `+`(`*`(C, `*`(x(t))), Du)
det(`+`(`*`(s, `*`(I)), `-`(A))) = `/`(`*`(`+`(`*`(`^`(S, 4), `*`(M[s], `*`(M[us]))), `*`(`^`(S, 3), `*`(M[s], `*`(B[s]))), `*`(`^`(S, 3), `*`(M[s], `*`(B[us]))), `*`(`^`(S, 2), `*`(M[s], `*`(K[us])))...
det(`+`(`*`(s, `*`(I)), `-`(A))) = `/`(`*`(`+`(`*`(`^`(S, 4), `*`(M[s], `*`(M[us]))), `*`(`^`(S, 3), `*`(M[s], `*`(B[s]))), `*`(`^`(S, 3), `*`(M[s], `*`(B[us]))), `*`(`^`(S, 2), `*`(M[s], `*`(K[us])))...
(9.6)
 

Determine the Transtere Functions 

> Tf := simplify(Multiply(Multiply(C_bar,MatrixInverse(S*IdentityMatrix( Nx, Nx )-A_bar)),B_bar)+D_bar):                                    Tf := collect(Tf,S,factor): '(z[s](t)-z[us](t))/diff(z[r](t),t), (z[s](t)-z[us](t))/F[c] , diff(z[s](t),t$2)/diff(z[r](t),t),diff(z[s](t),t$2)/F[c] ';  'C*(s*I-A_bar)^(-1),B+D';
 

 

`/`(`*`(`+`(z[s](t), `-`(z[us](t)))), `*`(diff(z[r](t), t))), `/`(`*`(`+`(z[s](t), `-`(z[us](t)))), `*`(F[c])), `/`(`*`(diff(z[s](t), `$`(t, 2))), `*`(diff(z[r](t), t))), `/`(`*`(diff(z[s](t), `$`(t, ...
`/`(`*`(C), `*`(`+`(`*`(s, `*`(I)), `-`(A_bar)))), `+`(B, D) (10.1)
 

> '(z[s](t)-z[us](t))/diff(z[r](t),t)'=Tf[1,1];                                                                                              '(z[s](t)-z[us](t))/F[c]'=Tf[1,2];                                                                                                      'diff(z[s](t),t$2)/diff(z[r](t),t)'=Tf[2,1];                                                                                            'diff(z[s](t),t$2)/F[c]'=Tf[2,2];
 

>
 

 

 

 

`/`(`*`(`+`(z[s](t), `-`(z[us](t)))), `*`(diff(z[r](t), t))) = `+`(`-`(`/`(`*`(S, `*`(M[s], `*`(`+`(K[us], `*`(S, `*`(B[us])))))), `*`(`+`(`*`(`^`(S, 4), `*`(M[s], `*`(M[us]))), `*`(`+`(`*`(M[s], `*`(...
`/`(`*`(`+`(z[s](t), `-`(z[us](t)))), `*`(diff(z[r](t), t))) = `+`(`-`(`/`(`*`(S, `*`(M[s], `*`(`+`(K[us], `*`(S, `*`(B[us])))))), `*`(`+`(`*`(`^`(S, 4), `*`(M[s], `*`(M[us]))), `*`(`+`(`*`(M[s], `*`(...
`/`(`*`(`+`(z[s](t), `-`(z[us](t)))), `*`(F[c])) = `/`(`*`(`+`(`*`(`+`(M[us], M[s]), `*`(`^`(S, 2))), `*`(S, `*`(B[us])), K[us])), `*`(`+`(`*`(`^`(S, 4), `*`(M[s], `*`(M[us]))), `*`(`+`(`*`(M[s], `*`(...
`/`(`*`(`+`(z[s](t), `-`(z[us](t)))), `*`(F[c])) = `/`(`*`(`+`(`*`(`+`(M[us], M[s]), `*`(`^`(S, 2))), `*`(S, `*`(B[us])), K[us])), `*`(`+`(`*`(`^`(S, 4), `*`(M[s], `*`(M[us]))), `*`(`+`(`*`(M[s], `*`(...
`/`(`*`(diff(z[s](t), `$`(t, 2))), `*`(diff(z[r](t), t))) = `/`(`*`(S, `*`(`+`(`*`(B[s], `*`(`^`(S, 2), `*`(B[us]))), `*`(`+`(`*`(K[s], `*`(B[us])), `*`(B[s], `*`(K[us]))), `*`(S)), `*`(K[s], `*`(K[us...
`/`(`*`(diff(z[s](t), `$`(t, 2))), `*`(diff(z[r](t), t))) = `/`(`*`(S, `*`(`+`(`*`(B[s], `*`(`^`(S, 2), `*`(B[us]))), `*`(`+`(`*`(K[s], `*`(B[us])), `*`(B[s], `*`(K[us]))), `*`(S)), `*`(K[s], `*`(K[us...
`/`(`*`(diff(z[s](t), `$`(t, 2))), `*`(F[c])) = `/`(`*`(`^`(S, 2), `*`(`+`(`*`(`^`(S, 2), `*`(M[us])), `*`(S, `*`(B[us])), K[us]))), `*`(`+`(`*`(`^`(S, 4), `*`(M[s], `*`(M[us]))), `*`(`+`(`*`(M[s], `*...
`/`(`*`(diff(z[s](t), `$`(t, 2))), `*`(F[c])) = `/`(`*`(`^`(S, 2), `*`(`+`(`*`(`^`(S, 2), `*`(M[us])), `*`(S, `*`(B[us])), K[us]))), `*`(`+`(`*`(`^`(S, 4), `*`(M[s], `*`(M[us]))), `*`(`+`(`*`(M[s], `*...
(10.2)
 

> 'R[0]=[ B  , A*B , A^2*B , A^3*B ]';
 

R[0] = [B, `*`(A, `*`(B)), `*`(`^`(A, 2), `*`(B)), `*`(`^`(A, 3), `*`(B))] (10.3)
 

>
 

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

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

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

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

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

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

> Matlab_Notations := { M[0]=M0, M[1] = M1, M[2] = M2, M[3] = M3, K[1] = K1, K[2] = K2, K[s] = Ks, B[1] = B1, B[2] = B2, B[3] = B3, K[1] = K1, K[2] = K2, Km[1] = Km1, Km[2]=Km2, N[1]=N1,N[2]=N2}:
 

> Experiment_Name := "Active_suspension_Experiment";
 

Active_suspension_Experiment (11.1)
 

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