4
$\begingroup$

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

which returns the warning enter image description here

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}
   ];
$\endgroup$
5
  • $\begingroup$ There are examples in the docs. If those approaches do not work for you, then the details of your problem are needed to figure out a better approach, if in fact it is possible. $\endgroup$
    – Michael E2
    Dec 4, 2022 at 0:46
  • $\begingroup$ Sorry, SE was struggling with the latex in the post and only the first paragraph posted. The full post is now up. $\endgroup$
    – Nolan King
    Dec 4, 2022 at 0:52
  • $\begingroup$ @NolanKing Please, pay attention, that $\theta 1=\pm \pi/2, \theta 1'=0$ are also equilibrium states. $\endgroup$ Dec 4, 2022 at 6:01
  • 1
    $\begingroup$ Please show us definition of 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 to InputForm before copying. For more info, read this: mathematica.meta.stackexchange.com/a/1585/1871 $\endgroup$
    – xzczd
    Dec 4, 2022 at 8:21
  • $\begingroup$ The code is attached now. $\endgroup$
    – Nolan King
    Dec 4, 2022 at 13:52

1 Answer 1

4
$\begingroup$

Solution 1

Simplifying definition of optimalcontrol to

optimalcontrol[z1_, z2_] = 
 With[{pm1 = Piecewise[{{1, z1 > 0}}, -1]}, 
  Piecewise[{{-pm1, pm1 z2 + Sqrt[2 z1 // Abs] > 0}}, pm1]]

and make use of Simplify`PWToUnitStep and introduce a smooth approximate UnitStep:

appro = With[{k = 10^5}, ArcTan[k #]/Pi + 1/2 &];

sol = 
   NDSolveValue[
    Simplify`PWToUnitStep@
      Simplify`PWToUnitStep@{eqs, θ1[0] == -π/6, 
        θ1'[0] == -1} /. UnitStep -> appro, θ1, {t, 0, 10}]; // AbsoluteTiming
(* {0.0235535, Null} *)

Plot[sol[t], {t, 0, 10}]

enter image description here

Solution 2

It turns out that, after turning off DiscontinuityProcessing, discard the step control of NDSolve by setting FixedStep as the Method of NDSolve works wonders:

sol = 
   NDSolveValue[{eqs, θ1[0] == -Pi/6, θ1'[0] == -1}, θ1, {t, 0, 10}, 
     Method -> {FixedStep, Method -> Automatic, DiscontinuityProcessing -> False}, 
     StartingStepSize -> 0.001]; // AbsoluteTiming
(* {0.10727, Null} *)

ListPlot[sol, Frame -> True]

enter image description here

I've used the hidden syntax of ListPlot, see this post for more info.

This is surprising, I should say. It's rare to see example in which adaptive ODE solver is inferior to fixed-step ODE solver.

Solution 3

Though the default adaptive ODE solver (it's LSODA, I believe) doesn't work well on the problem, StiffnessSwitching does if we adjust MaxStepSize a bit:

sol = 
   NDSolveValue[{eqs, θ1[0] == -Pi/6, θ1'[0] == -1}, θ1, {t, 0, 10}, 
     Method -> {StiffnessSwitching, DiscontinuityProcessing -> False}, 
     MaxStepSize -> 0.1]; // AbsoluteTiming
(* {0.0140824, Null} *)

ListPlot[sol, Frame -> True, PlotRange -> All]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.