How to Create Global Parameters in ModelingToolkit.jl

Specific Domains Modelling & Simulations
1

Is it possible to create global parameters? For example, when modeling a gas or fluid, it would be convenient to define them globally, similar to how is done in SimScape with a domain definition.

I attempted to create an ODESystem with my fluid parameters and then pass this named system to other systems that rely on the fluid properties. This works, but my “fluid” ODESystem then becomes a subsystem of each component ODESystem, and now if I want to change density (i.e. fluid.rho_0), I need to set it in many places.

For example, if I look at the parameters of my system sys, then fluid is defined multiple times…

julia> parameters(sys)
12-element Vector{Any}:
 i4sim₊f₊vSA₊fluid₊rho_0
 i4sim₊f₊vSA₊fluid₊bulk
 i4sim₊f₊vBR₊fluid₊rho_0
 i4sim₊f₊vBR₊fluid₊bulk
 i4sim₊piston₊v1₊fluid₊rho_0
 i4sim₊piston₊v1₊fluid₊bulk
 i4sim₊piston₊v2₊fluid₊rho_0
 i4sim₊piston₊v2₊fluid₊bulk
 i4sim₊piston₊m₊g
 i4sim₊piston₊m₊m
 i4sim₊p_s₊p_source
 i4sim₊p_r₊p_source

Is there a way to define a global parameter that doesn’t need to be added as a sub-system or parameter to each ODESystem component?

2

Hello.

Have you tried to define rho as a variable in your subsytems and then define a connection between them to propagate the parameters ?

This would allow you to set rho_0 only in one source subsystem.
For example :

function SourceFluid(;name, rho_0 = 1.0)
    val = rho_0
    @named p = Pin()
    @parameters rho_0
    @variables rho 
    eqs = [
           0 ~ rho - rho_0
           0 ~ p.rho - rho
          ]
    ODESystem(eqs, t, [rho], [rho_0], systems=[p], defaults=Dict(rho_0 => val, rho => rho_0), name=name)
end

function vSA(;name)
    @named p = Pin()
    @named n = Pin()
    @variables rho bulk
    eqs = [
           0 ~ p.rho - rho
           0 ~ n.rho - rho
           0 ~ bulk - ...
          ]
    ODESystem(eqs, t, [rho, bulk], [], systems=[p, n], defaults=Dict(rho => rho), name=name)
end

It is an adaptation of the example in the doc

3

There’s also the GlobalScope feature that was recently added. No docs, yet, but see the discussion here.

4

Thanks for this tip. I gave GlobalScope() a try, and technically it works well, but I can see a duplication of the global parameter in parent components. For example, the following code produces comp with 2 parameters named rho

using ModelingToolkit

@variables t
@variables p(t)
@parameters rho
D = Differential(t)

rho = GlobalScope(rho)

fluid(;name) = ODESystem(Equation[], t, [p], [rho], defaults=[rho=>1000, p=>0], name=name)

@named A = fluid()
@named B = fluid()

eqs = [
    D(B.p) ~ A.p
    0 ~ (A.p) - (A.rho)
]

@named comp = ODESystem(eqs, t, [], [], systems=[A, B])

Here is the result for comp

julia> comp
Model comp with 2 equations
States (2):
  A₊p(t) [defaults to 0]
  B₊p(t) [defaults to 0]
Parameters (2):
  rho [defaults to 1000]
  rho [defaults to 1000]

The solution works well, if I run the model, it does seem to set rho globally

ps = [
    rho=>900
]
prob = ODEProblem(comp, [], (0,0.1), ps)
solve(prob, ImplicitEuler())

With result…

retcode: Success
Interpolation: 3rd order Hermite
t: 9-element Vector{Float64}:
 0.0
 1.0e-6
 1.2e-6
 3.2e-6
 2.3199999999999998e-5
 0.00022319999999999998
 0.0022231999999999994
 0.02222319999999999
 0.1
u: 9-element Vector{Vector{Float64}}:
 [900.0, 0.0]
 [900.0, 0.0009]
 [900.0, 0.00108]
 [900.0, 0.0028799999999999997]
 [900.0, 0.02088]
 [900.0, 0.20088]
 [900.0, 2.0008799999999995]
 [900.0, 20.00087999999999]
 [900.0, 90.0]
5

Open an issue.

6

Follow-up from April 2021 discussion…

I am trying to move two quantities out of my component constructors: (i) a constant (gravity; should never be changed), and (ii) a shared liquid density (parameter, in principle, I would like to be able to change this).

Here is what I try to do:

# 1. Define parameter and constant:
@parameters rho = 1 [unit=u"kg/L", description="water density"]
@constants g = 98.1 [unit=u"dm/s^2", description="gravity"]

# 2. Give parameter and constant global scope
rho = GlobalScope(rho)
g = GlobalScope(g)


# Component constructor for open tank water level model
#
@component function Tank(; name, A=5, md_c=25, p_s=1e4, eps_v=1e-9)
    #
    # Symbolic model parameters
    params = @parameters begin 
        A=A,            [unit=u"m^2", description = "Cross sectional tank area"]
        md_c=md_c,      [unit=u"kg/s", description = "Effluent valve capacity"]
        p_s=p_s,        [unit=u"dPa", description = "Scaling pressure in valve model"]
        eps_v=eps_v,    [unit=u"1", description = "Softening parameter in valve model"]
    end
    # Symbolic model variables, with initial values needed
    vars = @variables begin
        m(t)=1.5*rho*A,    [unit=u"kg", description = "Liquid mass"]
        md_i(t),    [unit=u"kg/s", description = "Influent mass flow rate"]
        md_e(t),    [unit=u"kg/s", description = "Effluent mass flow rate"]
        V(t),       [unit=u"L", description = "Liquid volume"]
        h(t),       [unit=u"dm", description = "level"]
        Dp(t),      [unit=u"dPa", description = "Pressure drop across valve"]
    end
    # Model equations
    eqs = [
        D(m) ~ md_i-md_e
        m ~ rho*V
        V ~ A*h
        Dp ~ rho*g*h
        md_e ~ md_c*s_ssqrt(Dp/p_s; eps=eps_v)
    ]
    #
    return ODESystem(eqs, t, vars, params; name = name)
end

I then build a system of two tanks…

@component function Sys2Tank(; name)
    # Components used
    subsys = @named begin
        tank_1 = Tank()
        tank_2 = Tank()
    end
    # Equations for connecting components
    eqs = [
        tank_1.md_i ~ md(t)
        tank_2.md_i ~ tank_1.md_e
    ]
    sys = ODESystem(eqs, t; name = name)
    return compose(sys, subsys...)
end

I then simplify/compile the symbolic model:

@mtkcompile csys2tank = Sys2Tank()

and everything works so far.

But when I try to create a numeric problem, I get an error message that g and rho does not have a numeric value:

prob_sys = ODEProblem(csys2tank, [], tspan)

[Observe: I have set tspan = (0,10). Also, I need to specify function md(t) some place, and register it symbolically. Function s_ssqrt() is just a smoothed signed square root function.]

Some Questions:

  1. Does the above make sense?
  2. If I do parameters(csys2tank), both g and rho are listed among the parameters. I am a little surprised that constant g is listed as a parameter?
  3. Clearly, setting values to g and rho has not worked?
  4. How do I specify the values of g and rho in the ODEProblem constructor?
7 
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D

@parameters g rho
g   = GlobalScope(g)
rho = GlobalScope(rho)

function Tank(; name, A = 1.0)
    @parameters A = A
    @variables h(t) = 1.0 p(t)
    System([p ~ rho*g*h, D(h) ~ -A*p/(rho*g)], t; name)
end

@named tank1 = Tank(A = 1.0)
@named tank2 = Tank(A = 2.0)

@mtkcompile sys = System(Equation[], t; systems = [tank1, tank2],
                         initial_conditions = Dict(g => 9.81, rho => 1000.0))

prob = ODEProblem(sys, [], (0.0, 1.0))
solve(prob, Tsit5())   # Success