# Non-linear differential equation system problem

I have to prepare a system of differential equations to control, and so I need to solve the system for position, velocity, angle and angular velocity of both jaw and pitch angles. The equations are described by (11) (12) and (13) in the link: Two-wheeled balancing LEGO robot

So, my group mates used other softwares to compute the result, and they did get an aswer that is kinda close to what we expected, its just that the calculation fails before the whole pattern could be observed. So, I tried to get the results in Mathematica, but as it seems, the system is complex enough that one simple NDSolve won't help.

I have tried the code bellow:

{g, mb, mw, r, l, w, ip, ij, ixx, iw, jm, ra, kb, kt, bm} = {9.8, 1.5,
0.2, 0.1, 0.2, 1, 10, 10, 5, 5, 10, 10, 10, 10, 0.5};

u = HeavisideTheta[t];
x2[t] := D[x1[t], t];
x4[t] := D[x3[t], t];
x6[t] := D[x5[t], t];

system = {(mb + 2 mw + 2 (iw + jm)/r^2)*D[x2[t], t] +
2/r*(bm + kt*kb/ra)*D[x1[t], t] + (mb*l - 2*jm/r)*D[x6[t], t] -
2/r*(bm + kt*kb/ra)*D[x5[t], t] - mb*l*x5[t]*D[x5[t], t]^2 ==
kt*u/(r*ra), (2 (mw + (iw + jm)/r^2)*w^2 + ixx*x5[t]^2 + ij +
mb*l*x5[t]^2)*D[x4[t], t] +
2*((mb*l^2 + ixx - ij)*x5[t]*D[x5[t], t] +
w^2/r^2*(bm + kt*kb/ra))*D[x3[t], t] == 0,
(mb*l^2 + ip + 2*jm)*D[x6[t], t] + (mb*l - 2*jm/r)*D[x2[t], t] +
2*(bm + kt*kb/ra)*D[x3[t], t] -
2/r*(bm + kt*kb/ra)*D[x1[t], t] - (mb*l^2 + ixx - ij)*x5[t]*
D[x3[t], t]^2 - mb*g*l*x5[t] == -kt*u/ra, x1 == 0,
x2 == 0, x3 == 0, x4 == 0, x5 == 0, x6 == 0};

NDSolve[system, {x1, x2, x3, x4, x5, x6}, t]


Which then returns that this is not an ordinary differential equation. Our teacher has suggested that we use the transformation x'(t)=y(t) so that the equations were all first order differential equations. Any help is appreciated.

• Take a look at ref/FeedbackLinearize#1896001842. (The second example in FeedbackLinearize>Applications>Electromechanical Systems). Also AffineStateSpaceModel will convert the equations to first-order differential equations for you. – Suba Thomas Dec 12 '16 at 20:40
• Thanks, I'll use it from now on; the problem is, our teacher wants that we demonstrate the whole process at this stage. – Lagherta Dec 12 '16 at 21:08

First of all the transformation suggested by your teacher should be part of your system. The derivatives are treated as separate variables. So your system will look like:

{g, mb, mw, r, l, w, ip, ij, ixx, iw, jm, ra, kb, kt, bm} = {9.8, 1.5,
0.2, 0.1, 0.2, 1, 10, 10, 5, 5, 10, 10, 10, 10, 0.5};
u[t_] := HeavisideTheta[t];

system = {(mb + 2 mw + 2 (iw + jm)/r^2)*D[x2[t], t] +
2/r*(bm + kt*kb/ra)*D[x1[t], t] + (mb*l - 2*jm/r)*D[x6[t], t] -
2/r*(bm + kt*kb/ra)*D[x5[t], t] - mb*l*x5[t]*D[x5[t], t]^2 ==
kt*u[t]/(r*ra), (2 (mw + (iw + jm)/r^2)*w^2 + ixx*x5[t]^2 + ij +
mb*l*x5[t]^2)*D[x4[t], t] +
2*((mb*l^2 + ixx - ij)*x5[t]*D[x5[t], t] +
w^2/r^2*(bm + kt*kb/ra))*D[x3[t], t] ==
0, (mb*l^2 + ip + 2*jm)*D[x6[t], t] + (mb*l - 2*jm/r)*
D[x2[t], t] + 2*(bm + kt*kb/ra)*D[x3[t], t] -
2/r*(bm + kt*kb/ra)*D[x1[t], t] - (mb*l^2 + ixx - ij)*x5[t]*
D[x3[t], t]^2 - mb*g*l*x5[t] == -kt*u[t]/ra,
x2[t] == D[x1[t], t], x4[t] == D[x3[t], t], x6[t] == D[x5[t], t] ,
x1 == 0, x2 == 0.0, x3 == 0, x4 == 0.0, x5 == 0,
x6 == 0.0};

NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 10}]


The above gives error because Mathematica tries to differentiate HevisideTheta function at $t=0$

Per @J.M. suggestion I replaced HevisideTheta with u[t_] := UnitStep[t]

And now the NDSolve returns solution, although it complains about singularity. Please check your equations.

• Thanks, the system was so big that I didn't even noticed the problem with the step function. Now, as for the singularity complain, it actually is a stiff system. It is probably the same problem my mates were having when using another softwares. Since it is not a equation problem, is there a way to work around this? – Lagherta Dec 12 '16 at 20:40
• @Lagherta Unfortunately I didn't have time to look deeply enough in the paper you referred to. What do you expect? Some oscillations? Is motor feedback system already included in these equation? – BlacKow Dec 12 '16 at 20:59
• I don't know for sure, but I guess it is as expected, since it is a "two-wheeled balancing lego robot", and although we're planning on applying a PID control later on, for now we're just applying a constant voltage on both motors (and the voltage is the controlled variable). So, i think that when the system fails to calculate any further, it means that the robot has fell over. – Lagherta Dec 12 '16 at 21:18
• Why not replace HeavisideTheta[] with UnitStep[]? The former is really only intended for symbolic calculus. – J. M. is away Dec 13 '16 at 3:00
• @J.M. Good point, I edited the answer. – BlacKow Dec 13 '16 at 14:21