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))
```