The equations and code are below.
The short version is:
I know the error is because NDSolve
is doing something in the background to handle the discontinuous ODE. Adding the following to NDSolve
produces the correct result but it's incredibly ineffecient.
Method -> "DiscontinuityProcessing" -> False, AccuracyGoal -> 6,
PrecisionGoal -> 6, MaxSteps -> 10^7
Solving it this was is incredibly slow for a second order system. What's the right way to handle discontinuous ODEs?
The ODE: $$ \theta 1''(t) = \frac{-\theta 1'(t)^2 \sin(\theta 1(t)) \cos(\theta 1(t))+2 \sin(\theta 1(t))-u(t) \cos(\theta 1(t))}{2-\cos^2(\theta 1(t))} $$ where $$ u(t) = 2 \tan(\theta 1(t))-\theta 1'(t)^2 \sin(\theta 1(t)) + u_1(\theta1(t),\theta1'(t))(2 - \cos^2(\theta1(t))) $$ and $$ u_1(z_1,z_2) = \begin{cases}\begin{cases}1 & z_2 \leq -\sqrt{2z_1} \newline-1 & z_2 >-\sqrt{2z_1} \end{cases} & z_1 > 0 \newline\begin{cases}1 & z_2 < \sqrt{-2z_1}\newline-1 & z_2 \geq \sqrt{-2z_1} \end{cases} & z_1 \leq 0 \newline\end{cases} $$ For reference, the first two terms of the control eliminate the right hand side of $\theta1''$, and the last term is a time optimal control that drives the state to $(\theta1 = 0, \theta1'=0)$, which comes from the Pontryagin Maximum Principle.
I tried to use the following code to solve it
sol = NDSolve[{eqs, θ1[0] == -Pi/6, θ1'[0] == -1},θ1[t],{t, 0, 10}];
The reason the control has such an ugly piecewise definition is to check that we're never taking the square root of a negative number. So, it's not clear to me what is happening that's skipping all these cases and evaluating the square root anyway.
I see in the reference page for NDSolve
there is a WhenEvent
option that can be implemented, but it's not clear how to specify that the action should only be taken once (when the state reaches the switching curve, and not continue to cycle through the discrete variables for each step that it remains on the switching curve).
ClearAll[x1, x, y, s, u, eqs, u1, sol, uop];
l1 = 1;
m1 = 1;
mCart = 1;
g = 1;
x = s[t] + l1 Sin[θ1[t]];
y = l1 Cos[θ1[t]];
T = (1/2) m1 (D[x, t]^2 + D[y, t]^2) + (1/2) mCart (s'[t])^2;
V = m1 g ( y - 1);
L = T - V;
ELE = Flatten[
Solve[
FullSimplify[D[D[L, {{θ1'[t], s'[t]}}], t]] -
FullSimplify[D[L, {{θ1[t], s[t]}}]] == {0, u[t]},
{θ1''[t], s''[t]}
] /. {Rule -> Equal}
]
u1[t] = (u[t] /. Solve[ELE[[1]] /. {θ1''[t] -> 0}, u[t]] //
FullSimplify)[[1]]
optimalcontrol[z1_, z2_] =
Piecewise[{
{Piecewise[{
{1, z2 <= -Sqrt[2 z1]},
{-1, z2 > -Sqrt[2 z1]}}], z1 > 0},
{Piecewise[{
{1, z2 < Sqrt[-2 z1]},
{-1, z2 >= Sqrt[-2 z1]},
{0, z1 == 0 && z2 == 0}}], z1 <= 0}
}]
eqs = ELE[[1]] /. {u[t] ->
u1[t] + optimalcontrol[θ1[t], θ1'[
t]] (-2 + Cos[θ1[t]]^2)/Cos[θ1[t]]} //
FullSimplify
sol = NDSolve[
{eqs, θ1[0] == -Pi/6, θ1'[0] == -1},
θ1[t],
{t, 0, 10}
];
eqs
written with Mathematica code rather than $\LaTeX$. If you're having difficulty in making it look nice, press Ctrl+Shift+I to convert it toInputForm
before copying. For more info, read this: mathematica.meta.stackexchange.com/a/1585/1871 $\endgroup$