7

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.

m_goldberg
  • 107,779
  • 16
  • 103
  • 257
Christian Mullon
  • 71
  • 2
  • 2
    You could try to model this with a WhenEvent, perhaps. – user21 Feb 12 '16 at 00:52
  • Thanks, I have discovered the WhenEvent function. I have tried this:

    T = 0.1; fo[{u_, v_}] := If[u > v, {3 u - v^2, u^2}, {3 u - v^2, v^2}] { Derivative[1][X][t] == fo[X[t]],
    WhenEvent[Evaluate[t] == 0, X[t] -> ReplacePart[X[t], 1 -> 3]], WhenEvent[Evaluate[t] == T/2, X[t] -> ReplacePart[X[t], 2 -> 10]] } sol = NDSolve[system, X[t], {t, 0, T}] Plot[X[t] /. sol, {t, 0, T}]

    But it does not give the right solution

     – Christian Mullon Feb 12 '16 at 08:26
  • what is the control and state variables ? $u$ or $v$ ? in this case, the initial condition for state is given but how can you give a value for optimal initial condition for control variable ? – optimal control Feb 12 '16 at 13:33
  • Why do you not edit you question to contain the WhenEvent version. Also explains what system is and you might mention any error messages NDSolve gives you. – user21 Feb 12 '16 at 17:16
  • 'system' is a huge simplification of my optimal control problem, with two features (a) a vector expression, (b) initial condition on a variable, terminal condition on the other. – Christian Mullon Feb 13 '16 at 09:51
  • When I run the above program, It works, I get a plot but function v does not satisfy the terminal condition – Christian Mullon Feb 13 '16 at 09:53

1 Answers1

8

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

Community
  • 1
Michael E2
  • 235,386
  • 17
  • 334
  • 747