The exponential function is not supported for constraints when using the MTK interface

[rlampel]

Describe the bug 🐞

Trying to use the exponential function exp for constraints leads to an error when using the MTK interface.

Minimal Reproducible Example 👇

using Corleone
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEqTsit5
using ComponentArrays, ForwardDiff
using Optimization
using OptimizationMOI, Ipopt
using LuxCore, Random
using CairoMakie

# We start by defining our system
num_states = 3
num_controls = 1

@variables begin
    r(..) = 1., [tunable = false]
    v(..) = 0., [tunable = false]
    m(..) = 1., [tunable = false]
    u(..) = 0.5, [bounds = (0., 1.0), input = true]
end

@parameters begin
    T = 0.5, [bounds =(1.e-3, Inf), tunable = true]
end

@constants begin
    r_0 = 1., [tunable = false]
    v_0 = 0., [tunable = false]
    m_0 = 1., [tunable = false]
    r_T = 1.01, [tunable = false]
    b = 7., [tunable = false]
    T_max = 3.5, [tunable = false]
    A = 310., [tunable = false]
    k = 500., [tunable = false]
    C = 0.6, [tunable = false]
end

Drag = A * v(t)^2 * exp(-k * (r(t) - r_0))

eqs = [
    D(r(t)) ~ T * v(t)
    D(v(t)) ~ T * (-1 / r(t)^2 + 1 / m(t) * (T_max * u(t) - Drag))
    D(m(t)) ~ T * (-b * u(t))
]

# Define control discretization
tspan = (0.,10.)
dt = 0.1
cgrid = collect(0.0:dt:last(tspan))

# grid_cons_u = [A * v(tᵢ)^2 * exp(-k * (r(tᵢ) - r_0)) ≲ 0.6 for tᵢ in cgrid]   # this leads to a different error
grid_cons_u = [exp(r(tᵢ)) ≲ 0.6 for tᵢ in cgrid]  # this only contains the exponential and leads to the described error

cons = [
    r(last(tspan)) ~ 1.01,
    grid_cons_u...
]

costs = [-m(last(tspan))]

@named oc_problem = System(eqs, t; costs = costs, constraints = cons)

function plot_oc_problem(sol)
    f = Figure()
    ax = CairoMakie.Axis(f[1, 1])
    scatterlines!(ax, sol, vars = collect(1:1:num_states))
    f[1, 2] = Legend(f, ax, "States", framevisible = false)
    ax1 = CairoMakie.Axis(f[2, 1])
    stairs!(ax1, sol, vars = collect(num_states + 1:1:num_states + num_controls))
    f[2, 2] = Legend(f, ax1, "Controls", framevisible = false)
    return f
end

# ## Single Shooting

dynopt = CorleoneDynamicOptProblem(
    oc_problem, [],
    u(t) => cgrid;
    algorithm = Tsit5(),
)

Error & Stacktrace ⚠️

ERROR: exp is not defined for Expr.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] error_f_symbolic(f::Function, T::Type)
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/eVL77/src/methods.jl:363
  [3] promote_symtype(::typeof(exp), T::DataType)
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/eVL77/src/methods.jl:335
  [4] _promote_symtype(f::Function, args::SymbolicUtils.SmallVec{…})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/eVL77/src/types.jl:1075
  [5] maketerm(::Type{…}, f::Any, args::SymbolicUtils.SmallVec{…}, metadata::Nothing)
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/eVL77/src/terminterface.jl:259
  [6] (::SymbolicUtils.DefaultSubstituter{…})(ex::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/eVL77/src/substitute.jl:198
  [7] (::SymbolicUtils.DefaultSubstituter{…})(ex::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/eVL77/src/substitute.jl:180
  [8] #substitute#378
    @ ~/.julia/packages/SymbolicUtils/eVL77/src/substitute.jl:263 [inlined]
  [9] substitute
    @ ~/.julia/packages/SymbolicUtils/eVL77/src/substitute.jl:261 [inlined]
 [10] (::CorleoneModelingToolkitExtension.var"#25#42"{Symbol, Dict{Any, Any}})(::Tuple{Int64, Inequality})
    @ CorleoneModelingToolkitExtension ~/.../NLP_Solvers/Corleone.jl/ext/MTKExtension/optimal_control.jl:131
 [11] iterate
    @ ./generator.jl:47 [inlined]
 [12] collect_to!(dest::Vector{…}, itr::Base.Generator{…}, offs::Int64, st::Tuple{…})
    @ Base ./array.jl:892
 [13] collect_to_with_first!
    @ ./array.jl:870 [inlined]
 [14] collect(itr::Base.Generator{Base.Iterators.Enumerate{…}, CorleoneModelingToolkitExtension.var"#25#42"{…}})
    @ Base ./array.jl:844
 [15] map(f::Function, A::Base.Iterators.Enumerate{Vector{Any}})
    @ Base ./abstractarray.jl:3314
 [16] CorleoneDynamicOptProblem(sys::System, inits::Vector{…}, controls::Pair{…}; algorithm::Tsit5{…}, shooting::Nothing, tspan::Nothing, sensealg::SciMLBase.NoAD, kwargs::@Kwargs{})
    @ CorleoneModelingToolkitExtension ~/.../NLP_Solvers/Corleone.jl/ext/MTKExtension/optimal_control.jl:130
 [17] top-level scope
    @ ~/.../NLP_Solvers/Corleone.jl/examples/julia-mtk-benchmarks/oc_goddarts_rocket.jl:80

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

⌃ [336ed68f] CSV v0.10.15
  [13f3f980] CairoMakie v0.15.8
  [b0b7db55] ComponentArrays v0.15.32
  [2751e0db] Corleone v0.0.3 `~/Dokumente/Promotion/NLP_Solvers/Corleone.jl`
  [6590a4b1] CorleoneOED v0.0.2 `~/Dokumente/Promotion/NLP_Solvers/Corleone.jl/lib/CorleoneOED`
  [a93c6f00] DataFrames v1.8.1
  [f6369f11] ForwardDiff v1.3.2
  [b6b21f68] Ipopt v1.14.0
  [bb33d45b] LuxCore v1.5.3
  [961ee093] ModelingToolkit v11.10.0
  [7f7a1694] Optimization v5.5.0
  [22f7324a] OptimizationLBFGSB v1.4.0
  [fd9f6733] OptimizationMOI v1.2.0
  [1dea7af3] OrdinaryDiffEq v6.108.0
  [b1df2697] OrdinaryDiffEqTsit5 v1.9.0
⌃ [1ed8b502] SciMLSensitivity v7.95.0
  [2efcf032] SymbolicIndexingInterface v0.3.46
  [3a884ed6] UnPack v1.0.2

• Output of versioninfo()

Julia Version 1.10.10
Commit 95f30e51f41 (2025-06-27 09:51 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × AMD Ryzen 7 4800H with Radeon Graphics
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver2)
Threads: 1 default, 0 interactive, 1 GC (on 16 virtual cores)
Environment:
  LD_LIBRARY_PATH = /opt/gurobi1103/linux64/lib
  JULIA_EDITOR = code
  JULIA_VSCODE_REPL = 1

Additional context

Add any other context about the problem here.

[AlCap23]

This is due to k and not exp :D

Currently we did not support parameters or constants in the constraints. I am fixing this right now 😄

[rlampel]

Thank you for the response! Sorry for the confusion, my main problem is that the line
grid_cons_u = [exp(r(tᵢ)) ≲ 0.6 for tᵢ in cgrid]
does not work, but
grid_cons_u = [2.7182818^(r(tᵢ)) ≲ 0.6 for tᵢ in cgrid]
does. (So it's not a huge problem after all 😅.)