# 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] == 0,
x2[0] == 0, x3[0] == 0, x4[0] == 0, x5[0] == 0, x6[0] == 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. Dec 12, 2016 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. Dec 12, 2016 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] == 0, x2[0] == 0.0, x3[0] == 0, x4[0] == 0.0, x5[0] == 0,
x6[0] == 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? Dec 12, 2016 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? Dec 12, 2016 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. Dec 12, 2016 at 21:18
• Why not replace HeavisideTheta[] with UnitStep[]? The former is really only intended for symbolic calculus. Dec 13, 2016 at 3:00
• @J.M. Good point, I edited the answer. Dec 13, 2016 at 14:21