ModelingToolKit: Cannot `convert` an object of type SymbolicUtils.BasicSymbolic{Float64}

Specific Domains Modelling & Simulations
1 
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations: solve

f(t) = t
@register_symbolic f(t)

@mtkmodel FOL begin
    @variables begin
        y(t)
        z(t)
    end
    @equations begin
        y ~ f(t)
        D(z) ~ D(y)
    end
end

@mtkbuild fol = FOL()
prob = ODEProblem(fol, [fol.z => 1], (0.0, 18.0), []);
sol = solve(prob);

yields

MethodError: Cannot `convert` an object of type SymbolicUtils.BasicSymbolic{Float64} to an object of type Float64
Stacktrace:
  [1] setindex!(A::Vector{Float64}, x::SymbolicUtils.BasicSymbolic{Float64}, i1::Int64)
    @ Base .\array.jl:1021
  [2] macro expansion
    @ C:\Users\a1058035\.julia\packages\SymbolicUtils\EGhOJ\src\code.jl:418 [inlined]
  [3] macro expansion
    @ C:\Users\a1058035\.julia\packages\Symbolics\Nk48t\src\build_function.jl:546 [inlined]
  [4] macro expansion
    @ C:\Users\a1058035\.julia\packages\SymbolicUtils\EGhOJ\src\code.jl:375 [inlined]
  [5] macro expansion
    @ C:\Users\a1058035\.julia\packages\RuntimeGeneratedFunctions\M9ZX8\src\RuntimeGeneratedFunctions.jl:163 [inlined]
  [6] macro expansion
    @ .\none:0 [inlined]
  [7] generated_callfunc
    @ .\none:0 [inlined]
  [8] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…}, ::Vector{…}, ::Float64)
    @ RuntimeGeneratedFunctions C:\Users\a1058035\.julia\packages\RuntimeGeneratedFunctions\M9ZX8\src\RuntimeGeneratedFunctions.jl:150
  [9] (::ModelingToolkit.var"#f#717"{…})(du::Vector{…}, u::Vector{…}, p::ModelingToolkit.MTKParameters{…}, t::Float64)
    @ ModelingToolkit C:\Users\a1058035\.julia\packages\ModelingToolkit\U2JaS\src\systems\diffeqs\abstractodesystem.jl:342
 [10] (::SciMLBase.Void{ModelingToolkit.var"#f#717"{…}})(::Vector{Float64}, ::Vararg{Any})
    @ SciMLBase C:\Users\a1058035\.julia\packages\SciMLBase\rR75x\src\utils.jl:482
 [11] (::FunctionWrappers.CallWrapper{…})(f::SciMLBase.Void{…}, arg1::Vector{…}, arg2::Vector{…}, arg3::ModelingToolkit.MTKParameters{…}, arg4::Float64)
    @ FunctionWrappers C:\Users\a1058035\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:65
 [12] macro expansion
    @ C:\Users\a1058035\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:137 [inlined]
 [13] do_ccall
    @ C:\Users\a1058035\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:125 [inlined]
 [14] FunctionWrapper
    @ C:\Users\a1058035\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:144 [inlined]
 [15] _call
    @ C:\Users\a1058035\.julia\packages\FunctionWrappersWrappers\9XR0m\src\FunctionWrappersWrappers.jl:12 [inlined]
 [16] FunctionWrappersWrapper
    @ C:\Users\a1058035\.julia\packages\FunctionWrappersWrappers\9XR0m\src\FunctionWrappersWrappers.jl:10 [inlined]
 [17] Void
    @ C:\Users\a1058035\.julia\packages\SciMLBase\rR75x\src\utils.jl:482 [inlined]
 [18] (::FunctionWrappers.CallWrapper{…})(f::SciMLBase.Void{…}, arg1::Vector{…}, arg2::Vector{…}, arg3::ModelingToolkit.MTKParameters{…}, arg4::Float64)
    @ FunctionWrappers C:\Users\a1058035\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:65
 [19] macro expansion
    @ C:\Users\a1058035\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:137 [inlined]
 [20] do_ccall
    @ C:\Users\a1058035\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:125 [inlined]
 [21] FunctionWrapper
    @ C:\Users\a1058035\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:144 [inlined]
 [22] _call
    @ C:\Users\a1058035\.julia\packages\FunctionWrappersWrappers\9XR0m\src\FunctionWrappersWrappers.jl:12 [inlined]
 [23] FunctionWrappersWrapper
    @ C:\Users\a1058035\.julia\packages\FunctionWrappersWrappers\9XR0m\src\FunctionWrappersWrappers.jl:10 [inlined]
 [24] ODEFunction
    @ C:\Users\a1058035\.julia\packages\SciMLBase\rR75x\src\scimlfunctions.jl:2297 [inlined]
 [25] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.Tsit5Cache{…})
    @ OrdinaryDiffEq C:\Users\a1058035\.julia\packages\OrdinaryDiffEq\kMujN\src\perform_step\low_order_rk_perform_step.jl:799
 [26] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.DefaultCache{…})
    @ OrdinaryDiffEq C:\Users\a1058035\.julia\packages\OrdinaryDiffEq\kMujN\src\perform_step\composite_perform_step.jl:34
 [27] __init(prob::ODEProblem{…}, alg::OrdinaryDiffEq.CompositeAlgorithm{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::@Kwargs{…})
    @ OrdinaryDiffEq C:\Users\a1058035\.julia\packages\OrdinaryDiffEq\kMujN\src\solve.jl:524
 [28] __init (repeats 4 times)
    @ C:\Users\a1058035\.julia\packages\OrdinaryDiffEq\kMujN\src\solve.jl:11 [inlined]
 [29] #__solve#500
    @ C:\Users\a1058035\.julia\packages\OrdinaryDiffEq\kMujN\src\solve.jl:6 [inlined]
 [30] __solve
    @ C:\Users\a1058035\.julia\packages\OrdinaryDiffEq\kMujN\src\solve.jl:1 [inlined]
 [31] solve_call(_prob::ODEProblem{…}, args::OrdinaryDiffEq.CompositeAlgorithm{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\a1058035\.julia\packages\DiffEqBase\c8MAQ\src\solve.jl:612
 [32] solve_call
    @ C:\Users\a1058035\.julia\packages\DiffEqBase\c8MAQ\src\solve.jl:569 [inlined]
 [33] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::ModelingToolkit.MTKParameters{…}, args::OrdinaryDiffEq.CompositeAlgorithm{…}; kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\a1058035\.julia\packages\DiffEqBase\c8MAQ\src\solve.jl:1080
 [34] solve_up
    @ C:\Users\a1058035\.julia\packages\DiffEqBase\c8MAQ\src\solve.jl:1066 [inlined]
 [35] #solve#51
    @ C:\Users\a1058035\.julia\packages\DiffEqBase\c8MAQ\src\solve.jl:1003 [inlined]
 [36] solve
    @ C:\Users\a1058035\.julia\packages\DiffEqBase\c8MAQ\src\solve.jl:993 [inlined]

Did I do anything wrong? If I change f(t) to sin(t), it works.

2

I’m not an expert here, but do you need to register f(t) as symbolic? I think it is necessary to use @register_symbolic f(t) if you define your own function, typically one that Julia does not know how to differentiate [e.g., an interpolation function of experimental data, a step input with if statements, etc.]. But not if you use built-in functions such as sin(t), t, etc.

Again, I’m not an expert.

I also notice that you differentiate your input function f(t) in your model. It may be that by registering the function via @register_symbolic that you tell MTK that it should not attempt to differentiate the function, and then you explicitly differentiate it in the equations part.

Again, I’m just guessing here.

3

That seems to be true… But in this system taking the derivative isn’t necessary. The solution is just z = 1 + f(t). I guess if I write z = initial_z + y, and make initial_z a parameter, I can achieve the same affect, but it’s not elegant.

To be more concrete, the engineers I work with have an equation of the form
(z2 - z1) = 54.213(y2 - y1) (1 and 2 are different time points). I figured that the D(z) ~ 54.213D(y) was the natural way to write that in differential equation form. Since y contains input data, I followed Getting Started with ModelingToolkit.jl · ModelingToolkit.jl (sciml.ai) and hit the problem in this post. Am I missing something?

4

Hm. So how do the engineers that you work with use this equation:
z_2 - z_1 = 54.213(y_2 - y_1)
Specifically,

  • are y_1 and y_2 known?
  • is z_1 known?
  • is the constant always 54.213?

A similar problem could be the relationship between position s and velocity v:

v = \frac{\mathrm{d}s}{\mathrm{d}t}

If velocity v is constant, the solution is (of course):

s_2 - s_1 = v\cdot (t_2 - t_1)

Here, the solution is simple, so one would not pose this as a differential equation.
I probably miss something in your problem.

5

As written the model starts as index 2 IIRC, so you do need to differentiate it in order to get an index 1 in order to run the standard processes. That’s a limitation of Pantelides.

But the error (which we should improve BTW, I wouldn’t say it’s a good error message) is then saying the derivative is not defined on the registered function.