NDSolve for optimal control theory

Question

For an optimal control problem, I need to solve a differential equation in which (a) the system is formulated according to a function and (b) initial/terminal conditions are discontinuous.

Here is an example of the code I would like to implement

fo[{u_, v_}] := If[u > v, {3 u - v^2, u^2}, {3 u - v^2, v^2}]

system = {{u'[t], v'[t]} == fo[{u[t], v[t]}], u[0] == 3, v[10] == 4}

NDSolve[system, {u[t],v[t]}, {t, 0, 0.3}]

Of course, there are more dimensions and function is more complicated. I have no idea if this is possible, and if so on how to proceed. Any suggestion is welcome.

Answer 1

You can get NDSolve to try to answer by setting up the equations individually. In the OP's case, the message suggests that this type of problem, in which the constraints on u and v are at different times and there is discrete switch (by the If statement), is not solvable.

fo[{u_, v_}] := {3 u - v^2, If[u > v, u^2, v^2]};

system = {{u'[t], v'[t]} == fo[{u[t], v[t]}], u[0] == 3, v[10] == 4};

NDSolve[system, {u[t], v[t]}, {t, 0, 0.3}]

NDSolve::bvdisc: NDSolve is not currently able to solve boundary value problems with discrete variables. >>

Update 2: Considerations explained in the first update below lead to this alternative:

NDSolve[system, {u[t], v[t]}, {t, 0, 0.3},
 Method -> {Automatic, "DiscontinuityProcessing" -> False}]

The shooting method NDSolve uses fails with the boundary condition v[10] == 4 due to a singularity that occurs with the automatically chosen initial conditions, but it succeeds with the updated BC v[0.3] == 4. (A singularity will occur before t = t0 + 1/v[t0], if v[t0] > 0, since v'[t] >= v[t]^2.) NDSolve complained about an unknown method if I replaced Automatic with "Shooting", with or without the "StartingInitialConditions" option.

---

Update: Extending the idea above. one could give separate NumericQ-protected components. The use of NumericQ protects the function(s) from being analyzed symbolically. This prevents symbolic analysis of potential discontinuities, which are treated as discrete events and the source of the error above. This generally can lead to a more difficult time-integration, but not always.

Clear[fov];
fov[{u_?NumericQ, v_?NumericQ}] :=
  If[u > v, u^2, v^2];

system = {{u'[t], v'[t]} == {3 u[t] - v[t]^2, fov[{u[t], v[t]}]}, 
   u[0] == 3, v[0.3] == 4};

{sol} = NDSolve[system, {u, v}, {t, 0, 0.3}]

Alternatively, one can try to implement the shooting method, with NDSolve is using here, oneself. There are several answer on the site which do this. The following ones seem related to the type of problem in this question:

How to use NDSolve with discontinuities at internal boundaries?

Can NDSolve handle discontinuous data?

Implementing one's own shooting method is sometimes useful. For instance, when I tweak the parameters in this problem, NDSolve sometimes runs into singularities with some of its chosen starting initial conditions for the shooting method. In that case, the shooting method seems to always fail. One can set manual "StartingInitialConditions", but sometimes they're hard to find. One approach in the following answer implements a shooting method that searches for starting initial conditions that extend the range of integration to push back a singularity as far as possible.

Nonlinear differential equation: numerical solution