Different initial conditions / trajectories when updating parameters with setp+remake vs re-solving steady state (ModelingToolkit ODEProblem)

Specific Domains Modelling & Simulations
1

Hi all,

I’m confused about how remake interacts with initialization / consistent initial conditions in a ModelingToolkit workflow where I:

  1. Solve a “static” system using a SteadyStateProblem to get an equilibrium.
  2. Use that equilibrium to initialize a “dynamic” (pulsatile) system.
  3. Then I change parameters and re-solve the dynamic system.

Depending on whether I update parameters via SymbolicIndexingInterface (SII) setp + remake(prob, p=...) or via re-solving the steady state and remaking with new u0, I get different results, mainly in the initial value of a state variable and the early transient.

I’m attaching two plots showing the discrepancy.

SII_approach
SII_approach599×399 22.8 KB

Remake_u0_approach
Remake_u0_approach600×398 24.5 KB

Minimal reproducible example

using ModelingToolkit
using DifferentialEquations
using SciCompDSL
using ModelingToolkit: t_nounits as t, D_nounits as D
using SymbolicIndexingInterface
using SciMLStructures: Tunable, canonicalize, replace
using PreallocationTools
using Plots

@connector Pin begin
  p(t)
  q(t), [connect = Flow]
end

@mtkmodel OnePort begin
  @components begin
      in = Pin()
      out = Pin()
  end
  @variables begin
      Δp(t)
      q(t)
  end
  @equations begin
      Δp ~ in.p - out.p
      0 ~ in.q + out.q
      q ~ in.q
  end
end

@mtkmodel ConstantPressure begin
  @components begin
      node = Pin()
  end
  @parameters begin
      P, [description = "Pressure"]
  end
  @equations begin
      node.p ~ P
  end
end

@mtkmodel ConstantFlow begin
  @components begin
      node = Pin()
  end
  @parameters begin
      Q, [description = "Flow"]
  end
  @equations begin
      node.q ~ -Q
  end
end

@component function DrivenFlow(; name, Factor=1.0, fun)
    @named node = Pin()
    sts = []
    ps = @parameters Factor = Factor
    eqs = [ node.q ~ -Factor * fun(t) ]
    compose(ODESystem(eqs, t, sts, ps; name=name), node)
end

@mtkmodel Resistor begin
  @extend Δp, q = oneport = OnePort()
  @parameters begin
      R, [description = "Resistance"]
  end
  @equations begin
      Δp ~ q * R
  end
end

@mtkmodel Compliance begin
  @components begin
      in = Pin()
      out = Pin()
  end
  @parameters begin
      k, [description = "Elastance Coefficient"]
      p0, [description = "Offset pressure"]
  end
  @variables begin
      V(t)
      C(t)
  end
  @equations begin
      out.q ~ -in.q
      D(V) ~ in.q
      in.q * (in.p - out.p - p0) * k ~ (D(in.p)-D(out.p))
      C ~ 1/(k*(in.p - out.p - p0))
  end
end

# -------------------------
# Numerical values
vals = Dict{Symbol, Float64}()
vals[:Pbaseline_v] = 10*1330
vals[:Pv_v] = 5*1330
vals[:Routflow_v] = 8*60*1330
vals[:Qinjection_v] = 0.01
vals[:E_v] = 0.6
vals[:p0_v] = 0
vals[:V_v] = 100
vals[:Qprod_v] = (vals[:Pbaseline_v] - vals[:Pv_v]) / vals[:Routflow_v]
vals[:Q_mean_v] = 0.01
vals[:hr_v] = 72.0
vals[:amp_v] = 2
vals[:f_v] = vals[:hr_v] / 60

# -------------------------
# Static system
@named Q_Inj = ConstantFlow(Q = 0.0)
@named Q_Prod = ConstantFlow(Q = vals[:Qprod_v])
@named Ro = Resistor(R = vals[:Routflow_v])
@named Cc = Compliance(k = vals[:E_v], p0 = vals[:p0_v])
@named Psinus = ConstantPressure(P = vals[:Pv_v])

elements = [Q_Inj, Q_Prod, Ro, Cc, Psinus]
assembly = [
    connect(Q_Inj.node, Q_Prod.node, Ro.in, Cc.in),
    connect(Ro.out, Cc.out, Psinus.node)
]
@mtkcompile sys_static = System(assembly, t; systems = elements)

pmap = Dict(
    Q_Inj.Q => 0.0,
    Q_Prod.Q => vals[:Qprod_v],
    Ro.R => vals[:Routflow_v],
    Cc.k => vals[:E_v],
    Cc.p0 => vals[:p0_v],
    Psinus.P => vals[:Pv_v],
)

u0 = [
    Cc.V => vals[:V_v],
    Cc.in.p => vals[:Pbaseline_v]
]

prob_steady = SteadyStateProblem(sys_static, merge(Dict(u0), pmap))
sol_steady = solve(prob_steady, SSRootfind())

# -------------------------
# Dynamic system (pulsatile injection)
@named Q_Inj_dynamic = DrivenFlow(fun = t -> vals[:Q_mean_v] + vals[:amp_v]*cos(2π*vals[:f_v]*t))
elements_dynamic = [Q_Inj_dynamic, Q_Prod, Ro, Cc, Psinus]
assembly_dynamic = [
    connect(Q_Inj_dynamic.node, Q_Prod.node, Ro.in, Cc.in),
    connect(Ro.out, Cc.out, Psinus.node)
]
@mtkcompile sys_dynamic = System(assembly_dynamic, t; systems = elements_dynamic)

pmap_dynamic = Dict(
    Q_Inj_dynamic.Factor => 1.0,
    Q_Prod.Q => vals[:Qprod_v],
    Ro.R => vals[:Routflow_v],
    Cc.k => vals[:E_v],
    Cc.p0 => vals[:p0_v],
    Psinus.P => vals[:Pv_v],
)

vars_to_keep = first.(u0)
u0_ss = [var => sol_steady[var] for var in vars_to_keep]

prob_dynamic = ODEProblem(sys_dynamic, merge(Dict(u0_ss), pmap_dynamic), (0, 10*60), saveat=1e-2)
sol_dynamic = solve(prob_dynamic, Rodas5P(), maxiters=1e9)

# -------------------------
# Parameter change
p_change = Dict(
    Psinus.P => 8645.0,
    Ro.R => 638400.0,
    Cc.k => 0.6
)

# ==============
# Approach A: SII setp + remake(prob, p=ps_new)
params_to_update = collect(keys(pmap_dynamic))
setter_dynamic = setp(prob_dynamic, params_to_update)

ps = parameter_values(prob_dynamic)
param_buffer = copy(canonicalize(Tunable(), ps)[1])
diffcache = DiffCache(param_buffer)

function update_parameters_and_solve_SII(p_change)
    ps = parameter_values(prob_dynamic)

    buffer = get_tmp(diffcache, Float64[])
    copyto!(buffer, canonicalize(Tunable(), ps)[1])
    ps_new = replace(Tunable(), ps, buffer)

    new_param_values = [get(p_change, param, pmap_dynamic[param]) for param in params_to_update]
    setter_dynamic(ps_new, new_param_values)

    prob_new = remake(prob_dynamic, p=ps_new)
    return solve(prob_new, Rodas5P(), maxiters=1e9)
end

sol_dynamic_new_SII = update_parameters_and_solve_SII(p_change)

# ==============
# Approach B: re-solve steady state with new params -> remake dynamic with new u0
function update_parameters_and_solve_redoSteady(p_change)
    pmap_steady_new = merge(pmap, p_change)
    prob_steady_new = remake(prob_steady, p=pmap_steady_new)
    sol_steady_new = solve(prob_steady_new, SSRootfind())

    vars_to_keep = first.(u0)
    u0_ss_new = [var => sol_steady_new[var] for var in vars_to_keep]

    pmap_dynamic_new = merge(pmap_dynamic, p_change)
    prob_dynamic_new = remake(prob_dynamic, p=pmap_dynamic_new, u0=u0_ss_new)
    return solve(prob_dynamic_new, Rodas5P(), maxiters=1e9)
end

sol_dynamic_new_redoSteady = update_parameters_and_solve_redoSteady(p_change)

plot(sol_dynamic, idxs=[Cc.in.p/1330], label="Dynamic")
plot!(sol_dynamic_new_SII, idxs=[Cc.in.p/1330], label="Dynamic new", title="Replacing parameters with the SII approach")

What I observe

  • sol_dynamic_new_SII and sol_dynamic_new_redoSteady do not match.
  • The main mismatch appears right at t = 0, as if the initial condition is not being treated consistently (or re-initialization behaves differently).

Questions

  1. In a workflow with a static system and a dynamic system where the dynamic system changes an element (here ConstantFlow → DrivenFlow), is it considered good practice to:
  • solve a steady state on the static system
  • then use its equilibrium as u0 for the dynamic system?
  1. In Approach A (SII setp + remake), am I missing an initialization / consistent IC step?
  • Does remake(prob, p=...) ever recompute consistent initial conditions for MTK-generated problems?
  • If not, what is the recommended way to ensure u0 is consistent after parameter changes?
  1. (If relevant) Should I be using ODEProblem(sys_dynamic, ...) with an explicit initialization system / initializealg, or something like remake(prob; u0=...) every time parameters change?

Thanks a lot for any guidance !

2

Sorry I am not in your field. but my take is:

  1. In my area (groundwater modelling) it is considered good practice for us to have a initial s.s. model as initial value (u0).
  2. Your second approach seems more correct because the steady-state condition may change, given the new parameters. Without studying your system in detail, I would reccomend that your initialization also solves the steady-state problem to update u0 as well.
  3. I think you can use remake(...). But make sure your initialization is consistent and solves the s.s. problem as well.
3

Thanks, that’s very helpful !

I think you’re right: if the parameters change, the steady-state may change as well, so reusing the old u0 can be inconsistent.