Mupad to setup state-space equations of vibrations using Lagrange’s equation

If you want to solve a dynamic system’s response (like a vibrating system) to arbitrary inputs using matlab’s control system toolbox it helps to have a symbolic math application setup the equations for you. There are many ways to skin a cat but here is the overall approach followed by a technique to use mupad / matlab’s symbolic math engine to setup the state space system.

Overall approach:

1.       Reduce the real system to a schematic of lumped parameters (springs, masses, dampers).

a.       This is the conceptual reduction of a continuous system to a finite dimensional system and you would probably like to do this manually.

2.       Set up the Lagrange equations. Many vibrations textbooks will show you how, eg. Rao, S.S., Mechanical Vibrations.

a.       This is easy to setup (the kinetic energy, potential energy and Rayleigh dissipation equations).

3.       Convert to equations of motion (mx’’ + cx’ + kx = F)

a.       This is tedious for complicated systems and you could spend a lot of time debugging your algebra if you do this manually. – this is the step this post attempts to help you automate.

4.       Convert to state space

a.       Somewhat tedious and this sample script below does a little to address it.

Mupad code: this code will setup equations of motion in state space form y’ = Ay+Bu for a m1 attached to ground with spring/mass and m2 and m3 are parallelly or independently attached to m1 with spring/damper. The system is not all that important. The way to use mupad / matlab to get the SS equations is the point here.

T:= 1/2 * m1 * x1md^2 + 1/2*m2*x2md^2 + 1/2*m3*x3md^2

V:= 1/2*k1*x1^2+1/2*k2*(x2-x1)^2+1/2*k3*(x3-x1)^2

R := 1/2*c1*x1d^2+1/2*c2*(x2d-x1d)^2+1/2*c3*(x3d-x1d)^2

dof1m:=diff(T, x1md);
dof1m:=subs(dof1m, x1md=x1'(t));
dof1m:=rewrite(diff(dof1m, t), D);
dof1:=collect(dof1m - diff(T, x1) + diff(R, x1d) + diff(V, x1), [x1, x2, x3, x1d, x2d, x3d]);
dof1:=collect(subs(dof1, x1=y1, x1d=y2, x2 = y3, x2d = y4, x3=y5, x3d=y6), [y1, y2, y3, y4, y5, y6]);
dof1:=subs(dof1, y1''(t) = y2d);
dof1:=(dof1=F1);
solve(dof1, y2d, IgnoreAnalyticConstraints)

dof2m:=diff(T, x2md);
dof2m:=subs(dof2m, x2md=x2'(t));
dof2m:=rewrite(diff(dof2m, t), D);
dof2:=collect(dof2m - diff(T, x2) + diff(R, x2d) + diff(V, x2), [x1, x2, x3, x1d, x2d, x3d]);
dof2:=collect(subs(dof2, x1=y1, x1d=y2, x2 = y3, x2d = y4, x3=y5, x3d=y6), [y1, y2, y3, y4, y5, y6]);
dof2:=subs(dof2, y3''(t) = y4d);
dof2:=(dof2=F2);
solve(dof2, y4d, IgnoreAnalyticConstraints)

dof3m:=diff(T, x3md);
dof3m:=subs(dof3m, x3md=x3'(t));
dof3m:=rewrite(diff(dof3m, t), D);
dof3:=collect(dof3m - diff(T, x3) + diff(R, x3d) + diff(V, x3), [x1, x2, x3, x1d, x2d, x3d]);
dof3:=collect(subs(dof3, x1=y1, x1d=y2, x2 = y3, x2d = y4, x3=y5, x3d=y6), [y1, y2, y3, y4, y5, y6]);
dof3:=subs(dof3, y5''(t) = y6d);
dof3:=(dof3=F3);
solve(dof3, y6d, IgnoreAnalyticConstraints)

Why would I choose such a solution approach?

1.       Why not a finite element model: Lumped parameter models are much faster and are helpful to have to try a large number of scenarios out quickly.

2.       Matlab is way cool.

3.       Why not differential equations solution: sure, why not? But I like thinking in laplace transforms / transfer functions / state space.

4.       Control system toolbox with functions like lsim, bode etc is very neat.

5.       You could potentially integrate your system into a feedback loop and design controllers very conveniently once you have your dynamics setup.

6.       “Ok, but why the Lagrange approach? I can write the equations of motion by inspection.” More power to you! I get confused with signs and subscripts with large number of DOF.

7.       “So why bother to post this. What have you really added?” I couldn’t find a good answer on the web on how to use mupad to setup a lagrangian approach to equations of motion. I spent a good number of hours on the wknd and this am figuring this out and want to share my hard-won wisdom. And I’ve been telling myself for years to become a more active blogger.

8.       “Is this the best way to do this in mupad?” Not likely. This is the attempt of a newbie with symbolic math (and I don’t plan to get a whole lot better in the symbolic world). But its more substantial that anything I found on the internets.