# \[pre-ANN\] DifferentialInclusions.jl

**URL:** <https://discourse.julialang.org/t/pre-ann-differentialinclusions-jl/99604>\
**Category:** Package Announcements\
**Tags:** ordinarydiffeq\
**Created:** [May 30, 2023, 12:51pm UTC](https://discourse.julialang.org/t/pre-ann-differentialinclusions-jl/99604 "2023-05-30T12:51:52Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [May 30, 2023, 12:51pm UTC](https://discourse.julialang.org/t/pre-ann-differentialinclusions-jl/99604/1 "2023-05-30T12:51:52Z")

</div>

This is a very very early announcement 🙂_(please let me know if I’m doing bullshit 😃 )_

I’m working on a package for non-smooth dynamical systems, _aka differential inclusions, aka dynamical complementarity problems._

## Current state (almost trivial):

> **[GitHub - SteffenPL/DifferentialInclusions.jl](https://github.com/SteffenPL/DifferentialInclusions.jl)**
>
> Contribute to SteffenPL/DifferentialInclusions.jl development by creating an account on GitHub.

Right now, the package is just a few lines of code and has only one feature:

- It can solve a first-order dynamic complementarity problem of this form

\dot x = f(x) + \sum\_{j} \lambda\_j \nabla g\_j(x) \\ g\_j(x) \geq 0, \quad \lambda\_j \geq 0, \quad g\_j(x) \lambda\_j = 0.

with a numerical method that _approximates_ Moreau’s catch-up scheme

x\_{k+1} = P\_S( x\_k + f(x\_k) ) \\ \text{where } S = \{ x \in \mathbb{R}^d \mid g\_j(x) \geq 0 \quad \text{for all } 1 \leq j \leq m \}.

With a suitable approximation of the projections P\_S. That is related to solving **linear complementarity problems**.

## Next steps:

I’m planning on extending the package to provide:

- first: more solvers (first: classical linear complementarity solvers, then, interfaces to other Julia packages)
- sparse interfaces for constraints (e.g., for non-overlap between objects in contact mechanics)
- later: better, more general interface

## Small demo:

```julia
using DifferentialInclusions, OrdinaryDiffEq

cons = (
    u -> u[1], 
    u -> 2.0 - u[2], 
    u -> abs(u[1] - u[2]) - 1.0
)

ode = ODEProblem( (u, p, t) -> -p[1] * u, [0.0, 2.0], (0.0,1.0), [10.0])
prob = DIProblem(ode, cons)

# OSPJ: one-step projected Jacobi method
alg = ProjectiveMethod(OSPJ(),Euler())  
sol = solve(prob, alg, dt = 0.001)    

```

## Questions:

> If there are related Julia packages, or, if you are interested in such problems, let me know 🙂

Related github issue: [Differential inclusions/Dynamic complementary problems · Issue #958 · SciML/DifferentialEquations.jl · GitHub](https://github.com/SciML/DifferentialEquations.jl/issues/958)

* * *

Note: I’m preparing the package mainly for a publication where I want to compare a few numerical methods.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 31, 2023, 1:51pm UTC](https://discourse.julialang.org/t/pre-ann-differentialinclusions-jl/99604/2 "2023-05-31T13:51:05Z")

</div>

I commented on the issue. But there is a big push going right now on SciML for complementary problems with @avikpal and a few others in the lab, and it would be good to collaborate on this and get this all fully into the interface.

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [May 31, 2023, 1:53pm UTC](https://discourse.julialang.org/t/pre-ann-differentialinclusions-jl/99604/3 "2023-05-31T13:53:14Z")

</div>

Amazing! I will try to get in touch with @avikpal, and I can discuss more details with him. Of course, I would be very happy to collaborate.

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [June 23, 2023, 3:46pm UTC](https://discourse.julialang.org/t/pre-ann-differentialinclusions-jl/99604/4 "2023-06-23T15:46:44Z")

</div>

## [Update & Questions]

The package seems to be usable for first-order DI.

- I added the classical projected Gauss-Seidel methods and `PSOR`.
- I implemented **sparse constraints**.
- Convergence tests a la `DiffEqDevTools.jl` is possible (see images for the test case of non-overlapping spheres)

_There are a few questions I had during coding (I can also create separate posts if that’s better…)_:

- Is there a data type for `OnceDifferentiable` functions which keeps a `DiffResult` and is type stable? Somehow [NLSolversBase.OnceDifferentiable](https://github.com/JuliaNLSolvers/NLSolversBase.jl/blob/master/src/objective_types/oncedifferentiable.jl) is not type stable with respect to the function `f`.

- Going further, I’m kind of searching for a way to lazily define `~ N^2` many constraints but with \mathcal{O}(1) memory, e.g., without creating an `OnceDifferentiable` for each constraint individually. I’m not sure how such an interface could look like… Below is my attempt:

The code below shows my current attempt to implement the constraints

\Vert X\_i - X\_j \Vert \geq 2R \quad \text{for all } 1 \leq i \< j \leq N.

in a way that allows the solver to extract `u[:,i]` and `u[:,j]` as static vectors/matrices, such that AD works nice and fast. However, the interface feels a bit complex.

### Example

```julia
N = 40
u0 = 2 * ( rand(2, N) .- 0.5 )

cs = let
    R = 0.2
    f = (u) -> sum(x -> x^2, u[:,1] - u[:,2]) - R^2
    
    con = TSOnceDifferentiable(f, MMatrix{2,2}(zeros(2,2)))  
    l_inds = LinearIndices(u0)

    pairs = ((i,j) for i in 1:N for j in 1:i-1)
    pairs_indices = ( SMatrix{2,2,Int64,4}[[l_inds[:,i] l_inds[:,j]] for (i,j) in pairs ] )

    # one constraint function which will be evaluated for each of the indices
    SparseConstraints(pairs_indices, con)
end 

ode = ODEProblem((du,u,p,t) -> (@. du = -p.gamma * u), u0, (0.0,1.0), (gamma = 1.0,))
prob = DIProblem(ode, cs)
alg = ProjectiveMethod(PBD(), Euler())

sol_ = solve(ode, Euler(), dt = 1e-4); # 3.369 ms (40044 allocations: 22.54 MiB)
sol = solve(prob, alg, dt = 1e-4); # 85.476 ms (514091 allocations: 30.74 MiB)

```

Note that solving `prob::DIProblem` is of course slower than the `ode::ODEProblem` since one needs to compute `~ N^2` constraints.

### Plots

 ![spheres_terminal](https://global.discourse-cdn.com/julialang/original/3X/2/f/2f6306a4367f5a4c964a947bae19cc99c2b361e4.png)  
 ![spheres_error_work](https://global.discourse-cdn.com/julialang/original/3X/4/4/44310ddc7c4a384b475b5bb522e6646df6851529.png)