MTK and JuMP

Specific Domains Modelling & Simulations
, ,
1

I would like to use MTK to build a physical model and than use the equations and convert them into a JuMP model. Any ideas how to do that?

My thoughs (based on Automatic Conversion of Julia Code to C Functions · Symbolics.jl):

  • This could be a method that takes an ODESystem and a JuMP.Model and attaches the needed variables and constraints to it.
  • Apply build_function and afterwards feed JuMP.VariableRefs through it and obtain expressions
  • Use this RHS expressions and do a simple forward Euler discretization.
using ModelingToolkit

@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

eqs = [
    D(x) ~ σ*(y-x),
    D(y) ~ x*(ρ-z)-y,
    D(z) ~ x*y - β*z,
]

sys = ODESystem(eqs, t, [x,y,z], [σ, ρ, β])

fnc = eval(build_function(
    equations(sys), 
    states(sys), 
    parameters(sys), 
    t, 
    target=Symbolics.JuliaTarget(), expression=Val{false},
)[1])

However, I still get Equations back.

julia> fnc(zeros(3), [10, 28, 8/3], 0.0)
3-element Array{Equation,1}:
 Differential(t)(x(t)) ~ σ*(y(t) - x(t))
 Differential(t)(y(t)) ~ x(t)*(ρ - z(t)) - y(t)
 Differential(t)(z(t)) ~ x(t)*y(t) - (β*z(t))

import JuMP
m = JuMP.Model()
JuMP.@variable(m, x[1:3])

julia> fnc(x, [10, 28, 8/3], 0.0)
3-element Array{Equation,1}:
 Differential(t)(x(t)) ~ σ*(y(t) - x(t))
 Differential(t)(y(t)) ~ x(t)*(ρ - z(t)) - y(t)
 Differential(t)(z(t)) ~ x(t)*y(t) - (β*z(t))
2

There’s a way to do this with an arbitrary differential equation stepper:

https://mtk.sciml.ai/dev/tutorials/nonlinear_optimal_control/

It’s on the ControlSystem but easy to map over the ODESystem.

3

Thx!
Do you mean the runge_kutta_discretize(sys::ControlSystem?

4

Yes. Needs a bit more to do everything but it’s close

5

@ChrisRackauckas @ohmsweetohm1 Dear all, not sure this is the best place to ask but: what is the current state of MTK in terms of re-usability / translation towards other interfaces (such as JuMP) for optimisation?

For instance: how can I use MTK and Ipopt?

6

MTK hooks into Optimization.jl.

http://optimization.sciml.ai/stable/tutorials/symbolic/

Optimization.jl has an MOI interface to hit the JuMP solvers directly.

http://optimization.sciml.ai/stable/optimization_packages/mathoptinterface/

7

Thanks for the prompt answer: having a look at it.