[pre-ANN] DifferentialInclusions.jl

Package Announcements
1

This is a very very early announcement :slight_smile: (please let me know if I’m doing bullshit :smiley: )

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

Current state (almost trivial):

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:

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 :slight_smile:

Related github issue: Differential inclusions/Dynamic complementary problems · Issue #958 · SciML/DifferentialEquations.jl · GitHub

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

10 Likes
2

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.

1 Like
3

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.

4

[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 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

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
spheres_error_work

4 Likes