Differential inclusions/Dynamic complementary problems

[SteffenPL]

TLTR
I would like to implement a problem type for differential inclusions (aka differential equations with inequalities).
Is that welcome/fitting to SciML?

Terminology

There are millions of alternative formulations for differential inclusions:
Nonsmooth dynamical systems, Dynamic complementarity problems, Differential equations with inequality constraints, Projective dynamics, set-valued ODEs, 'variational inequalities'.
In the following, I will say differential inclusions, but mostly I mean 'differential equation + inequality constraint', which
is the typical setup for numerical codes.

Current state

At the moment, differential inclusions are only supported through callbacks in SciML. While this allows some examples, such as a bouncing ball, it would be nice to treat differential inclusions similar to differential-algebraic equations, e.g., to provide a framework that allows integration of differential equations with inequality constraints. (I wrote below a bit about the mathematical/numerical motivation).

I'm motivated to implement such an interface (also since I want to do a numerical comparison study with the solver I'm working on), but I wanted to ensure that it fits into SciML and also ask in which form such a project should be developed (external or with the aim of internal?).

Implementation

A precise interface has to be discussed, but the primary idea would be to allow solving for example models like

$$\dot{x} = f(x) + \lambda g(x), \quad g(x) \geq 0, \quad \lambda \geq 0 \quad \text{and} \quad g(x) \lambda = 0.$$

In terms of Julia code, it could look like

# pseudo code:
function rhs!(du, u, p, t)
    x, lambda = u
    du[1] = f(u) + lambda * g'(u)
    du[2] = g(u)
end
 
differential_vars = [true, false]
prob = DIProblem(rhs!, du₀, u₀, tspan; differential_vars = differential_vars)

(There are tons of alternative formulations, going into details is out of scope here.)
In the far distance, it would be amazing to integrate with ModelingToolkit.

Existing solvers

To my knowledge, there are no existing solvers in Julia. However, there are some projects which implemented
solvers for specific differential inclusions:

• https://en.wikipedia.org/wiki/Siconos ( https://github.com/siconos/siconos )

Focusing on contact problems, there are many solvers who at their core, also solve differential inclusions:

• https://projectchrono.org/
• http://mimetics-engineering.fr/index.php/en/lmgc90-2/

The list is not very extensive, I have to do my homework better. But the field is also very cluttered...

My question:

• Can I work on defining such a DIProblem type with the aim of being compatible with SciML?

If you are not planning to add differential inclusions, I just make a small independent Julia package and safe some time with both bothering to interface ;)

---

Bonus: Application and math background

A good reference on differential inclusions is Numerical Methods for Nonsmooth Dynamical Systems.

To cite some typical applications taken from the book:

• Electrical Circuits with Diodes,
• Mechanical Contact Systems (with Coulomb Friction or Impact Laws).

(My own research is in mathematical biology, where one can simplify cell migration models to get faster numerics.
Generally, I think many agent-based models would benefit from fast solvers for differential inclusions.)

Why differential inclusions/non-smooth dynamical systems?

A good motivation for differential inclusions is to look at the numerical challenge of the alternatives.

For example, let's say we want to solve $\dot x = f(x)$ but such that $x \in S \subset \mathbb{R}^n$.
A typical approach would be to add penalty terms, e.g.

$$\dot x = f(x) - \varepsilon^{-1} \nabla_x \mathrm{dist}(S, x).$$

This works, but, as $\varepsilon \to 0$, the equations become very stiff. Which typically leads to the necessity to
find the ideal trade-off for the choice of this hyperparameter $\varepsilon$.

With some convex analysis (subdifferentials), one can compute the equations which occurs in the limit $\varepsilon \to 0$:

$$\dot x \in f(x) - N(S,x)$$

where $N(S,x)$ is the normal cone of the set $S$. (That is actually nothing else as the subdifferential $\nabla_x \lim_{\varepsilon \to 0} \frac{\mathrm{dist}(S,x)}{\varepsilon}$).

If one can compute $P_S(x)$ well enough, one can obtain unconditionally stable numerical methods. Moreover, gets ride of the annoying hyperparameter. (Numerical benefit + modelling benefit.)

Math outlook:

A lot of ideas from differential-algebraic equations generalize to differential inclusions. For example, for numerical methods which use the Lagrangian multiplier formulation as starting point, index-reduction plays a big role. Here, similar structural simplifications as done with ModelingToolkit could be very interesting, also to explore new numerics.

[SteffenPL]

Some updates. (Please feel free to close if irrelevant here!)

I'm starting now to implement solvers via DiscreteCallbacks. That seems more economical for me, as it doesn't need to be part of SciML then.

A very simple example is given here: https://gist.github.com/SteffenPL/84cc7a07a1b12a3a6156cf04c7957961
which allows to do:

cons = (
    u -> u[1],            # u[1] ≥ 0
    u -> 2.0 - u[2],   # u[2] ≤ 2
    u -> abs(u[1] - u[2]) - 1.0   # |u[1] - u[2]| ≥ 1
)

rhs!(du,u,p,t) = du .= -p[1] * u
u0 = [0.0, 2.0]
p = [10.0] 

ode = ODEProblem( rhs!, u0, (0.0,1.0), p)
prob = DIProblem(ode, cons)

alg = (ode = Euler(), lincompl = OneStepPJ())     # combine the Euler method with one step of the projected Jacobi method
sol = solve(prob, alg, dt = 0.001)

In contrast to DAEs, most solvers for differential inclusions require the solution of linear complementarity problems,
e.g. one tries to find $\lambda \in \mathbb{R}^m$ such that

$$A \lambda + q \geq 0, \quad \lambda \geq 0, \quad \lambda^T (A \lambda + q) \geq 0$$

where all equations are component-wise and with given vectors $q \in \mathbb{R}^m$ and $A \in \mathbb{A}^{m \times m}$

I am working on a small package that has a few more tools than the above, most imporantly

• some types to exploit the sparsity of the constraints (e.g. do not compute full gradients but only non-zero parts, that is quite relevant in applications) [that is essentially done, but not very general yet]
• collecting more typical solvers for linear complementarity problems

[ChrisRackauckas]

This looks great. I've love to see this all in the SciML interface. @avik-pal and crew have started working on LCPs so you may want to be in touch with them and collaborate.

I think what we want to do is make problem types for complementarity problems first and build out those solvers (potentially into nonlinear solve, depending on how the infra is done), then do the DIProblem stuff.

[SteffenPL]

Yes, I think LCP and then DIProblem makes perfect sense! I would be happy to collaborate or contribute to @avik-pal project!

I will probably still continue my small repository for experimenting with possible interfaces.

[avik-pal]

@SteffenPL, this looks quite interesting. We have been developing tools for sensitivity analysis of complementarity problems, so we only have basic solvers for LCPs and MCPs (like smooth nonlinear reformulation). These are already compatible with the SciML interface so integrating them should be trivial.

I can add you to that project, and we can brainstorm potential next steps for DIs there if you are interested.

Also tagging @yonatanwesen, who is working on complementarity problems with me.

[SteffenPL]

Cool, I would be happy to join the project!