NDSolve: Using time delays in a WhenEvent

Question

I was wondering whether there is a way to use WhenEvent in a system of delay differential equations, for example:

z[t_] = Piecewise[{{1, t <= 5}, {3/2, 5 < t <= 10}}, 2];
NDSolve[{x'[t] == 1, y'[t] == x'[t] - x'[t - 1], 
 WhenEvent[x[t - 1] == z[t], x[t] -> 0], x[t /; t <= 0] == 0, y[t /; t <= 0] == 0 }, 
 {x[t], y[t]}, {t, 0, 20}];

The issue seems to be that the WhenEvent doesn't recognize x[t-1] as a function and tries to evaluate the condition by throwing in the most recent x and t values. I'm not sure whether this is because WhenEvent can't access the value at x[t - 1] or because I need to input the condition in a way that explicitly tells it to go back and find this value.

---

A more complicated example to illustrate one of the systems I am working with:

n = 3; m = 2;
staff = {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
P = {{{0, 0}, {1, 0}, {0, 1}}, {{0, 0}, {0, 0}, {7/10, 0}}, {{0, 0}, {0, 4/5}, {0, 0}}};
PT = {{75, 360, 720, 0.3}, {75, 420, 930, 0.4}};
PA = {{1, 7/10}, {0, 0}, {0, 3/10}};
v = {{7/2, 5}, {7/2, 9/2}, {15, 20}};

a[t_, i_, k_] := PT[[k, 1]]*PA[[i, k]]*PDF[
    TriangularDistribution[{PT[[k, 2]], PT[[k, 3]]}, ((1 - PT[[k, 4]])*PT[[k, 2]] + PT[[k, 4]]*PT[[k, 3]])]
  , t];

Ns[t_, i_] := Piecewise[Table[{staff[[i, k]], (k - 1)*30 <= t < k*30}, {k, 1, 48}]];

evts = Flatten[Table[With[{i = i}, {
   WhenEvent[Evaluate[Subscript[V, i][t] == 0] && Evaluate[Subscript[S, i][t] <= Ns[t, i]], Subscript[dv, i][t] -> 0]
   , WhenEvent[{Evaluate[Subscript[V, i][t] > 0], Evaluate[Subscript[S, i][t] > Ns[t, i]]}, Subscript[dv, i][t] -> 1]
   , Subscript[dv, i][t /; t <= 0] == 1
  }], {i, 1, n}]];

eqns = Flatten[Table[With[{i = i}, With[{p = p}, {Subscript[U, i, p]'[ t] == (Subscript[dv, i][t])*(Subscript[Q, i, p][t])*Ns[t, i]/(Subscript[V, i][t] + .001) + (1 - Subscript[dv, i][t])*Subscript[L, i, p]'[t - v[[i, p]]], Subscript[U, i, p][t /; t <= 0] == 0}]], {i, 1, n}, {p, 1, m}]];

rls = Flatten[Table[With[{i = i}
   , {Subscript[S, i][t_] -> Sum[With[{p = p}, Subscript[L, i, p]'[t - v[[i, p]]]*v[[i, p]]], {p, 1, m}]
    , Subscript[V, i][t_] -> Sum[With[{p = p}, Subscript[Q, i, p][t - v[[i, p]]]*v[[i, p]]], {p, 1, m}]
    , Table[With[{p = p}
       , {Subscript[Q, i, p][t_] -> Subscript[L, i, p][t] - Subscript[U, i, p][t]
        , Subscript[L, i, p][t_] -> Integrate[a[t, i, p], t] + Sum[With[{jj = jj}, P[[jj, i, p]]*Subscript[U, jj, p][t]], {jj, 1, n}], Subscript[L, i, p]'[t_] -> a[i, t, p] + Sum[With[{jj = jj}, P[[jj, i, p]]*Subscript[U, jj, p]'[t]], {jj, 1, n}]
      }]
      , {p, 1, m}]
     }]
   , {i, 1, n}]];

NDSolve[
  Flatten[{evts, eqns} //. rls]
  , Flatten[Table[ With[{i = i}, {Subscript[dv, i][t], Table[With[{p = p}, {Subscript[U, i, p][t]}], {p, 1, m}]}], {i, 1, n}]]
  , {t, 0, 1000}
  , DiscreteVariables -> Flatten[Table[With[{i = i}, Subscript[dv, i][t]], {i, 1, n}]]];

This example gives fairly decent results without the time delays, but bad things start happening when any element in v gets too large.

Answer 1

Here's something I came up with to get you started. Consider this an extended comment until we can figure out how to make this work..

The idea is to use an auxiliary variable a[t] that is a DiscreteVariable, and use it to find the times at which x crosses z. We set the value of a[t] to be that particular time and then use WhenEvent to reset x one time unit later.

I will use a simplified version of your code that doesn't actually involve the delay piece involving y.

timeLimit = 10;
z[t_] = Piecewise[{{1, t <= 2}, {3/2, 2 < t <= 5}}, 2];
NDSolve[{x'[t] == 1/1.5 t
   , WhenEvent[x[t] == z[t], a[t] -> t]
   , WhenEvent[t == a[t] + 1, {x[t] -> 0, a[t] -> 0}]
   , a[t /; t <= 0] == timeLimit
   , x[t /; t <= 0] == 0}
  , {x[t], a[t]}
  , {t, 0, timeLimit}
  , DiscreteVariables -> {a}];
Show[
 Plot[Evaluate[{x[t], z[t], a[t]} /. %], {t, 0, timeLimit}]
 , Plot[{t, t - 1}, {t, 0, timeLimit}, 
  PlotStyle -> {{Black, Dashed, Opacity[0.5]}}]
 ]

The blue curve is the solution for x[t], the orange curve is z[t], the green curve is a[t], and the dashed lines are the functions t and t - 1 that are acting as aids for the eye to show that a[t] gets set to the current value of t and is then reset to 0 after 1 time unit has elapsed.

There are some important observations.

• We can see that a[t] gets set to t when x[t] crosses z[t]: at these points, x, z, and a all cross.
• We can see that x[t] gets reset to 0 when one time unit has elapsed since a[t] gets set to t.
• Most importantly, we can see that the solution is not perfect: between the second and third resets, x crosses z twice, and a gets set only once. The first crossing is not detected for some reason.

This last point illustrates a second possible weakness of this method: if there is a second crossing of x with z within one unit time step (i.e. before x gets reset), then a might get over-written with the new time, and so x won't reset until one unit time step after the new time.

To illustrate this last problem, change the first two lines in the code above to

timeLimit = 4;
z[t_] = Piecewise[{{1, t <= 2}, {2, 2 < t <= 5}}, 2];

and run the rest of it. This results in:

It's worse than I feared! a[t] gets set by WhenEvent to the current time t when x crosses z at about t == 1.8. a[t] gets set again at t == 2 when z jumps from having a value of 1 to having a value of 2; in other words, this is interpreted as a crossing between x and z even though z exhibits a discontinuity here! Finally, a[t] gets reset one more time when x crosses the new value of z, which occurs before t == 3, and since a[t] was set to 2 at the second "crossing", x won't get reset before t == 3. Finally, x resets at about t == 3.5.

The upshot: this is a good idea, but I don't know exactly how to make it robust.

Answer 2

Almost a solution: I feel that one ought to be able to integrate the delayed values of the variables via equations such as

x1[t] == x[t - 1], y2[t] == y[t - 2]

and use the variables in WhenEvent. It works in the following:

z[t_] = Piecewise[{{1/22, t <= 1}}, 2];
foo = NDSolve[{x'[t] == 2 t, y'[t] == x'[t] - x'[t - 1],
    x1[t] == x[t - 1], y2[t] == y[t - 2],
    WhenEvent[x1[t] + y2[t] > z[t], {x[t] -> 0, y[t] -> 0}], 
    x[t /; t <= 0] == 0, y[t /; t <= 0] == 1/2,
    x1[t /; t <= 0] == 0, y2[t /; t <= 0] == 1/2},
   {x, y, x1, y2}, {t, 0, 3}];

ListLinePlot[{x, y, x1, y2} /. First[foo]]

But that seems to be an accident. Extend the integration to {t, 0, 4}, and you'll get

NDSolve::nderr: Error test failure at t == 3.20710674525188`; unable to continue. >>

Integrate x'[t] == 2 Sqrt[t] instead of x'[t] == 2 t and I get a kernel crash!! Maybe it's a bug??