MathematicalSystems

[blegat]

There is a common system representation that is shared by a few control projects such as JuliaReach and Dionysos.jl in https://github.com/JuliaReach/MathematicalSystems.jl/. As discussed after your talk in Paris, it would be useful if OptimalControl.jl was supporting the system to be representing in this form.

cc @adrienbanse @schillic @mforets

[ocots]

Thanks for the comment, we will take a look at it. This seems effectively a nice idea.

[jbcaillau]

@blegat Hi Benoît, many thanks for the feedback. Sure, it definitely makes sense to be able to exchange between various high level representations / modelling formats. We were aware of MTK.jl, not of MathematicalSystems.jl which looks very nice and general. Two points here:

(i) The ConstrainedBlackBoxControlContinuousSystem type seems tailored to describe general (deterministic and explicit) systems

$$ \dot{x}(t) = f(x(t), u(t)) $$

plus state constraints, $x(t) \in X$, and control constraints, $u(t) \in U$. Is it also possible to add

• finite dimensional parameters to be optimised (such as the final time, $t_f$, whenever it's free)
• boundary constraints, $(x(0), x(t_f)) \in B \subset \mathbf{R^{2n}}$ (where $n$ is the state dim)
• mixed constraints, $(x(t), u(t)) \in Z \subset \mathbf{R}^n \times \mathbf{R}^m$ (where $m$ is the control dim)?

For instance, how would you describe this example (I hope the description is transparent enough 🤞🏾) with MathematicalSystems.jl?

@def ocp begin
    tf ∈ R, variable
    t ∈ [ 0, tf ], time
    x ∈ R², state
    u ∈ R, control
    tf ≥ 0
    -1 ≤ u(t) ≤ 1
    q = x₁
    v = x₂
    q(0) == 1
    v(0) == 2
    q(tf) == 0
    v(tf) == 0
    0 ≤ q(t) ≤ 5
    -2 ≤ v(t) ≤ 3
    (q^2 + v^2 + u^2)(t) ≤ 100
    ẋ(t) == [ v(t), u(t) ]
    tf → min
end

(ii) A strong motivation for us to map OptimalControl.jl models towards other descriptions is the availability of associated optimal control solvers; or, more specifically, of discretisation tools from an abstract optimal control description towards a mathematical program (that can eventually feed an optimisation solver). Are such tools available in MathematicalSystems.jl ecosystem?

[mforets]

The aim in MathematicalSystems.jl has been to make it as thin or lightweight as possible, yet with expressive power to specialize on system types if the developer needs so (hence the ugly struct names, though there is a macro meant to simplify the process of defining the system).

Addressing some of your questions:

• The package has some discretize features meant to go from *Continuous to *Discrete systems.

• Additional functionalities are a bit scattered around such as in ReachabilityAnalysis.jl which has methods to apply automatic differentiation leveraging MTK, or even to solve associated ODEs. Again, this is a bit scattered so if someone is using this in their package it would probably be nicer to just outsource some functionalities to a separate Julia pkg (or integrate them into MathematicalSystems as optional deps).

[jbcaillau]

Hi @mforets ; thanks for the input. I see that there are indeed some discretisation features for affine control systems. Going to general nonlinear systems is key to us (most of our applications are indeed from nonlinear optimal control), plus having a good blend of differential geometric tools to address both direct and indirect methods. There has been very nice progress on top of JuMP with InfiniteOpt.jl, though the description remains a bit far from what a mathematician would write (*).

That the DSL is smoothly operating with efficient solvers (optimisation, in particular) is also very important for us. At the moment, this is not completely true for MTK, e.g. Convergence will occur 🤞🏾

(*) Compare, e.g., this example with this:

using OptimalControl

ρ = 0.025              # discount rate
k = 100.0              # utility bliss point
T = 10.0               # life horizon
r = 0.05               # interest rate
B₀ = 100.0             # endowment
u(c; k=k) = -(c - k)^2 # utility function

@def ocp begin
    t ∈ [ 0, T ], time
    B ∈ R, state
    c ∈ R, control 
    B(0) == B₀
    B(T) == 0
    Ḃ(t) == r * B(t) - c(t)
    ∫( exp(-ρ * t) * u(c(t)) ) → max 
end

sol = solve(ocp; grid_size=1_000)

plot(sol)