parameter as a function of other parameters

[tpdsantos]

Question❓

I have a model of a continuously stirred tank (CST):

@mtkmodel CST begin
    @structural_parameters begin
        NC = 1
    end
    @parameters begin
        Q = 1.0
        V
    end
    @variables begin
        cin(t)[1:NC] = nothing
        cout(t)[1:NC] = 0.0
    end
    @equations begin
        D(cout) ~ Q/V * (cin - cout)
    end
end

I also have a model of a Plugflow Device (PF) but it's a bit large and is not that relevant for the question. The important part is that also has a parameter V
I wanted to simulate the hydrodynamic behaviour of a fluid in the tubes of an equipment by assuming that behaviour is a mix of CST and PF, in the sense that a part of the equipment volume behaves as CST and the rest as PF, so I did the following:

@mtkmodel CST_ED begin
    @structural_parameters begin
        NC = 1 # number of components
    end
    @parameters begin
        α  = 0.25
        Ve = 2.77
        Q  = 1.0
        cin[1:NC] = [0.5i for i in 1:NC]
    end
    @components begin
        T = CST(NC = NC)
        C = PF(NC = NC)
    end
    @equations begin
        T.V ~ α * Ve
        C.V ~ (1-α) * Ve
        T.cin ~ cin
        C.cin ~ T.cout
        C.Q ~ Q
    end
end

This approach does not work because the equations involving the parameters do not go to the observed section, which leaves me with $N$ variables and $N+2$ equations. How would I solve this without changing C.V and T.V to variables?

[baggepinnen]

do you need to do T = CST(NC = NC, cin = cin), and not add the T.cin ~ cin equation?

[tpdsantos]

To add, the main purpose is to perform parameter estimation to determine $\alpha$. I can also solve this problem by having $\alpha$ outside the model and then use setp to change both C.V and T.V according to $\alpha$, but that's not as elegant as changing only $\alpha$ as being an internal parameter of the @mtkmodel.

[tpdsantos]

do you need to do T = CST(NC = NC, cin = cin), and not add the T.cin ~ cin equation?

In this case, T.cin is a variable and I wanted to keep it that way, since I use it for other models and I want to leave it as general as possible. When I use it as components, then I decide if T.cin is fixed by some parameter or is the outlet of another equipment.

[baggepinnen]

T.V is a variable, so T.V ~ α * Ve should not be included and instead isntantiate CST like this. CST(V = α * Ve).

The same for PF.V with the equation C.V ~ (1-α) * Ve.

[SebastianM-C]

To add, the main purpose is to perform parameter estimation to determine α . I can also solve this problem by having α outside the model and then use setp to change both C.V and T.V according to α , but that's not as elegant as changing only α as being an internal parameter of the @mtkmodel.

I think that parameter dependencies or equations involving parameters are what would you need to express that. See https://docs.sciml.ai/ModelingToolkit/stable/tutorials/initialization/#Initialization-of-parameters for more about this.

The solution proposed above should be equivalent to parameter initialization if I understand correctly.

[tpdsantos]

T.V is a variable, so T.V ~ α * Ve should not be included and instead isntantiate CST like this. CST(V = α * Ve).

The same for PF.V with the equation C.V ~ (1-α) * Ve.

this does work but only in initialization, then T.V and C.V will be branded as an actual value instead of an expression involving $\alpha$ and Ve.
Consider the example:

@mtkmodel CST_2 begin
    @structural_parameters begin
        NC = 1
    end
    @parameters begin
        α = 0.5
        Ve = 1.0
        cin[1:NC] = [1.0 for i in 1:NC]
    end
    @components begin
        T1 = CST(NC = NC, V = α*Ve)
        T2 = CST(NC = NC, V = (1-α)*Ve)
    end
    @equations begin
        T1.cin ~ cin
        T2.cin ~ cin
    end
end

@mtkbuild m = CST_2()

u0 = unknowns .=> 0.0
prob = ODEProblem(m, u0, (0.0, 5.0))

getter = getp(m, [m.T.V, m.C.V])
setter! = setp(m, [m.α, m.Ve])

solve(prob, Tsit5())
getter(prob) # [0.5, 0.5]

setter!(prob, [0.25, 1.0])
solve(prob, Tsit5())
getter!(prob) # [0.5, 0.5]

This is because the expression provided is only valid in the initialization step, the parameter itself does not remain an expression dependent of other parameters, that's why it's not updated after several solves. The question is: can we create a parameter expression that is retained throughout the entire solution of the problem?

[ChrisRackauckas]

No, that is a bug because the initialization is not ran in this case. Parameter initialization is exactly what fixes that. #3347 should fix this case.

[AayushSabharwal]

This works:

julia> @mtkmodel CST begin
           @structural_parameters begin
               NC = 1
           end
           @parameters begin
               Q = 1.0
               V, [guess = 1.0] # this guess makes V solvable during initialization
           end
           @variables begin
               cin(t)[1:NC] = nothing
               cout(t)[1:NC] = 0.0
           end
           @equations begin
               D(cout) ~ Q/V * (cin - cout)
           end
       end
julia> @mtkmodel CST_2 begin
           @structural_parameters begin
               NC = 1
           end
           @parameters begin
               α = 0.5
               Ve = 1.0
               cin[1:NC] = [1.0 for i in 1:NC]
           end
           @components begin
               T1 = CST(NC = NC, V = α*Ve)
               T2 = CST(NC = NC, V = (1-α)*Ve)
           end
           @equations begin
               T1.cin ~ cin
               T2.cin ~ cin
           end
       end
julia> @mtkbuild m = CST_2()
julia> u0 = unknowns(m) .=> 0.0
julia> prob = ODEProblem(m, u0, (0.0, 5.0))
julia> getter = getp(m, [m.T1.V, m.T2.V])
julia> setter! = setp(m, [m.α, m.Ve])
julia> getter(solve(prob))
julia> getter(solve(prob))
2-element Vector{Float64}:
 0.5
 0.5
julia> setter!(prob, [0.25, 1.0])
julia> getter(solve(prob))
2-element Vector{Float64}:
 0.25
 0.75

The problem won't be updated with new values of m.T1.V, m.T2.V because initialization isn't solved until solve runs.

[AayushSabharwal]

This doesn't require #3347, it works right now

[tpdsantos]

so the secret is to add the guess to the parameter?

[AayushSabharwal]

It's not much of a secret :D https://docs.sciml.ai/ModelingToolkit/stable/tutorials/initialization/#Initialization-of-parameters documents how parameter initialization works.

[tpdsantos]

Always skimmed through this topic and never actually read this... Thanks, we can close this now