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

**URL:** <https://discourse.julialang.org/t/different-initial-conditions-trajectories-when-updating-parameters-with-setp-remake-vs-re-solving-steady-state-modelingtoolkit-odeproblem/135676>\
**Category:** Modelling & Simulations\
**Tags:** sciml\
**Created:** [February 16, 2026, 10:47am UTC](https://discourse.julialang.org/t/different-initial-conditions-trajectories-when-updating-parameters-with-setp-remake-vs-re-solving-steady-state-modelingtoolkit-odeproblem/135676 "2026-02-16T10:47:44Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![celia.bat](https://avatars.discourse-cdn.com/v4/letter/c/ee59a6/32.png) [@celia.bat](https://discourse.julialang.org/u/celia.bat)\
**Post date:** [February 16, 2026, 10:47am UTC](https://discourse.julialang.org/t/different-initial-conditions-trajectories-when-updating-parameters-with-setp-remake-vs-re-solving-steady-state-modelingtoolkit-odeproblem/135676/1 "2026-02-16T10:47:44Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/5/a/5a8b196570b4e90c035d57ed7fbe88991d570a4c.png)  
 ![Remake_u0_approach](https://global.discourse-cdn.com/julialang/original/3X/b/a/bad4e3bb5107729a6eb33ad23e6b33ee9624598e.png)

### Minimal reproducible example

```julia-auto
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 !

---

<div class="post-metadata">

**Author:** ![Vitor\_Patricio\_Canta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vitor_patricio_canta/32/50861_2.png) [@Vitor\_Patricio\_Canta](https://discourse.julialang.org/u/Vitor_Patricio_Canta)\
**Post date:** [February 16, 2026, 12:22pm UTC](https://discourse.julialang.org/t/different-initial-conditions-trajectories-when-updating-parameters-with-setp-remake-vs-re-solving-steady-state-modelingtoolkit-odeproblem/135676/2 "2026-02-16T12:22:03Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![celia.bat](https://avatars.discourse-cdn.com/v4/letter/c/ee59a6/32.png) [@celia.bat](https://discourse.julialang.org/u/celia.bat)\
**Post date:** [February 16, 2026, 4:20pm UTC](https://discourse.julialang.org/t/different-initial-conditions-trajectories-when-updating-parameters-with-setp-remake-vs-re-solving-steady-state-modelingtoolkit-odeproblem/135676/3 "2026-02-16T16:20:50Z")

</div>

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.