Initializing from the steady state

drto.initialize_steady_state seeds a model from its own equilibrium: one function, both declared shapes. A dynamic, discretized model is initialized from a throwaway reduced clone, the equilibrium solved by pyomo-pounce’s fill, project, block-solve pipeline with the declared controls as the decisions, then broadcast flat across the horizon. A steady-state model (authored steady, or a drto.dynamic_to_steady_state reduction) is initialized in place by the same pipeline. Both return a report worth printing.

The model is the Hicks-Ray CSTR from models/hicks.py. Requires the pounce extra (pip install drto[pounce]).

The dynamic path: a flat starting trajectory

import pyomo.environ as pyo
import drto
from models.hicks import hicks

m = hicks(N=5)
pyo.TransformationFactory("dae.collocation").apply_to(m, wrt=m.t, nfe=5, ncp=3, scheme="LAGRANGE-RADAU")
report = drto.initialize_steady_state(m, controls={m.v1: 0.57828, m.v2: 0.49989})
print(report)
print()
print(f"zc across the grid: {min(pyo.value(m.zc[t]) for t in m.t):.5f} "
      f"to {max(pyo.value(m.zc[t]) for t in m.t):.5f} (flat at the equilibrium)")

drto initialize_steady_state (reduce a clone -> solve -> broadcast)
  broadcast: 5 variables across 16 grid points, 32 derivatives zeroed
  pyomo-pounce initialize (fill -> repair -> block-solve)
    decisions held: 2
    values filled : 2
    projection    : optimal
    pyomo-pounce block_initialize
      decisions fixed   : 0
      system square     : True
      blocks solved     : 2 (1 by Newton 1x1, 1 subsystem solves)
      vars initialized  : 3
      left untouched    : 0 underdetermined, 0 overdetermined

zc across the grid: 0.64160 to 0.64160 (flat at the equilibrium)

The steady path: a reduced model, initialized in place

The same call on a model with no horizon runs the pipeline directly and the solved equilibrium lands in the variable values.

ss = pyo.TransformationFactory("drto.dynamic_to_steady_state").create_using(hicks(N=5))
report = drto.initialize_steady_state(ss)
print(report)
print()
print(f"equilibrium: zc = {pyo.value(ss.zc):.5f}, zt = {pyo.value(ss.zt):.5f}")

pyomo-pounce initialize (fill -> repair -> block-solve)
  decisions held: 2
  values filled : 2
  projection    : optimal
  pyomo-pounce block_initialize
    decisions fixed   : 0
    system square     : True
    blocks solved     : 2 (1 by Newton 1x1, 1 subsystem solves)
    vars initialized  : 3
    left untouched    : 0 underdetermined, 0 overdetermined

equilibrium: zc = 0.64160, zt = 0.53870