The steady-state branch: reduction, simulation, and optimization

One declared model serves both branches. The dynamic transforms (discretization, the terminal segment, the control profiles) build the horizon problem; the steady branch instead collapses the same declarations to a single point. drto.dynamic_to_steady_state performs the reduction: every derivative reference replaced by zero, time collapsed, the initial condition and terminal pieces removed. drto.steady_state_simulation composes that reduction, fixes the declared controls, and installs the zero objective: the square problem whose solution is the equilibrium under the given inputs.

The model is the Hicks-Ray CSTR from models/hicks.py. Its declared setpoints are a true fixed point of the dynamics, so simulating at the reference inputs must land exactly on the reference steady state: the model-consistency check the reduction makes automatic.

The reduction: the declared model, collapsed

import pyomo.environ as pyo
import pyomo_pounce  # registers pounce with the native factory
from pyomo.contrib.solver.common.factory import SolverFactory
import drto
from models.hicks import hicks

m = hicks(N=5)
ss = pyo.TransformationFactory("drto.dynamic_to_steady_state").create_using(m)
drto.info(ss)

drto registry

states: 2, controls: 2
states	zc (free), zt (free)
dynamics	dzc == (1 - zc)/(u2sf*v2) - k0*zc*exp(- ea/zt)
dynamics	dzt == (ztf - zt)/(u2sf*v2) + k0*zc*exp(- ea/zt) - a0*u1sf*v1*(zt - ztcw)
controls	v1 (free), v2 (free)
tracking stage cost	cost == 10*(zc - zc_ss)**2 + 2*(zt - zt_ss)**2 + (v1 - v1_ss)**2 + 0.5*(v2 - v2_ss)**2
steady-state targets	zc_ss (of zc), zt_ss (of zt)
steady-state control targets	v1_ss (of v1), v2_ss (of v2)

transformations

1. drto.dynamic_to_steady_state: removed=initial condition, tracking terminal cost, collapsed=7 Vars and 3 Constraints to a single point, derivatives=2 fixed at zero

drto.dynamic_to_steady_state also takes the model statement itself and reduces what it returns. hicks is declared hicks(N=5, h=1), so the bare call it makes gives the same grid as above. Reach for this form when the script owns the builder, and for the factory’s create_using when you hold a model you want left alone. A statement that builds its steady form natively, declaring no horizon, comes back untouched, so the dynamic model is never built.

ss2 = drto.dynamic_to_steady_state(hicks)
same = sorted(c.name for c in ss2.component_objects(pyo.Var)) == sorted(
    c.name for c in ss.component_objects(pyo.Var)
)
print(f"same reduction as the factory call above: {same}")
drto.info(ss2)

same reduction as the factory call above: True

drto registry

states: 2, controls: 2
states	zc (free), zt (free)
dynamics	dzc == (1 - zc)/(u2sf*v2) - k0*zc*exp(- ea/zt)
dynamics	dzt == (ztf - zt)/(u2sf*v2) + k0*zc*exp(- ea/zt) - a0*u1sf*v1*(zt - ztcw)
controls	v1 (free), v2 (free)
tracking stage cost	cost == 10*(zc - zc_ss)**2 + 2*(zt - zt_ss)**2 + (v1 - v1_ss)**2 + 0.5*(v2 - v2_ss)**2
steady-state targets	zc_ss (of zc), zt_ss (of zt)
steady-state control targets	v1_ss (of v1), v2_ss (of v2)

transformations

1. drto.dynamic_to_steady_state: removed=initial condition, tracking terminal cost, collapsed=7 Vars and 3 Constraints to a single point, derivatives=2 fixed at zero

Simulation at the reference inputs

sim = pyo.TransformationFactory("drto.steady_state_simulation").create_using(
    m, controls={m.v1: 0.57828, m.v2: 0.49989}
)
SolverFactory("pounce").solve(sim, tee=True)
print(f"equilibrium: zc = {pyo.value(sim.zc):.5f} (target 0.64160), "
      f"zt = {pyo.value(sim.zt):.5f} (target 0.53870)")

********************************************************************************

                    ####    ###   /   # /#   #/  ####  #####
                    #   #  #   # /#   #/ ##  /  #      #
                    ####   #   #/ #   /  # #/#  #      ####
                    #      #   /  #  /#  # /##  #      #
                    #       ##/    #/#   #/  #   ####  #####

********************************************************************************
This program contains POUNCE, a pure-Rust interior-point optimization solver
for nonlinear, conic, and global problems (its NLP core is ported from Ipopt).
Released under the Eclipse Public License (EPL) — drop-in compatible with Ipopt.
         For more information visit https://github.com/jkitchin/pounce
********************************************************************************

This is POUNCE version 0.10.0, running with linear solver FERAL.

Reading C:\Users\Devin\AppData\Local\Temp\tmpd4ukcycm\unknown007ejgvo.nl...
pounce: wrote C:\Users\Devin\AppData\Local\Temp\tmpd4ukcycm\unknown007ejgvo.sol
Parsed 2 vars, 2 cons, jac_nnz=4, h_nnz=2 in 0.01s
Problem class: NLP. Selected solver: NLP filter line-search interior-point (pounce-nlp) [solver_selection=auto].

Number of nonzeros in equality constraint Jacobian...:        4
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        2

Total number of variables............................:        2
                     variables with only lower bounds:        1
                variables with lower and upper bounds:        1
                     variables with only upper bounds:        0
Total number of equality constraints.................:        2
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter      objective   inf_pr   inf_du lg(mu)    ||d|| lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 1.12e-07 1.00e+00   -1.0 0.00e+00      - 0.00e+00 0.00e+00   0
   1  0.0000000e+00 9.30e-14 4.83e-07   -1.7 1.99e-06      - 1.00e+00 1.00e+00h  1
   2  0.0000000e+00 6.94e-18 1.51e-09   -3.8 1.75e-12      - 1.00e+00 1.00e+00h  1
   3  0.0000000e+00 2.08e-17 1.84e-11   -5.7 2.11e-16      - 1.00e+00 1.00e+00t  0
   4  0.0000000e+00 4.86e-17 2.51e-14   -8.6 7.97e-17      - 1.00e+00 1.00e+00t  0


Number of Iterations....: 4

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   2.5058834496277814e-14    2.5058834496277814e-14
Constraint violation....:   4.8572257327350599e-17    4.8572257327350599e-17
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5059035596778805e-09    2.5059035596778805e-09
Overall NLP error.......:   2.5059035596778805e-09    2.5059035596778805e-09


Number of objective function evaluations             = 5
Number of objective gradient evaluations             = 5
Number of equality constraint evaluations            = 5
Number of inequality constraint evaluations          = 5
Number of equality constraint Jacobian evaluations   = 5
Number of inequality constraint Jacobian evaluations = 5
Number of Lagrangian Hessian evaluations             = 4
Total seconds in POUNCE                              = 0.002

EXIT: Optimal Solution Found.

POUNCE 0.10.0: Optimal Solution Found.
Status: Solve_Succeeded
equilibrium: zc = 0.64160 (target 0.64160), zt = 0.53870 (target 0.53870)

drto.info(sim)

drto registry

states: 2, controls: 2
states	zc (free), zt (free)
dynamics	dzc == (1 - zc)/(u2sf*v2) - k0*zc*exp(- ea/zt)
dynamics	dzt == (ztf - zt)/(u2sf*v2) + k0*zc*exp(- ea/zt) - a0*u1sf*v1*(zt - ztcw)
controls	v1 (fixed), v2 (fixed)
steady-state targets	zc_ss (of zc), zt_ss (of zt)
steady-state control targets	v1_ss (of v1), v2_ss (of v2)

transformations

1. drto.dynamic_to_steady_state: removed=initial condition, tracking terminal cost, collapsed=7 Vars and 3 Constraints to a single point, derivatives=2 fixed at zero
2. drto.build_objective: objective=zero
3. drto.steady_state_simulation: controls=v1=0.57828, v2=0.49989, dropped=tracking stage cost

drto.steady_state_simulation also takes the model statement itself and prepares what it returns. Nothing is discretized on this path: the transformation composes the reduction, and a statement that constructs its steady form natively takes that reduction’s skip. Reach for this form when the script owns the builder, and for the factory’s create_using when you hold a model you want left alone.

sim2 = drto.steady_state_simulation(
    hicks, controls={"v1": 0.57828, "v2": 0.49989}
)
SolverFactory("pounce").solve(sim2)
print(f"equilibrium: zc = {pyo.value(sim2.zc):.5f} (target 0.64160), "
      f"zt = {pyo.value(sim2.zt):.5f} (target 0.53870)")

equilibrium: zc = 0.64160 (target 0.64160), zt = 0.53870 (target 0.53870)

Optimization with the controls free

drto.steady_state_optimization is the third mode on the steady branch: the same reduction, the controls left free, and the declared cost kept as the objective. It also takes the model statement itself. No example builder declares an economic stage cost, so on hicks this is the tracking-only path feature 009 allows, the steady point nearest the declared targets. Those targets are an equilibrium, so the optimum is the point the simulation above found, reached here with nothing held.

rto = drto.steady_state_optimization(hicks)
SolverFactory("pounce").solve(rto)
print(f"optimum: zc = {pyo.value(rto.zc):.5f}, zt = {pyo.value(rto.zt):.5f}, "
      f"v1 = {pyo.value(rto.v1):.5f}, v2 = {pyo.value(rto.v2):.5f}")
drto.info(rto)

optimum: zc = 0.64160, zt = 0.53870, v1 = 0.57828, v2 = 0.49989

drto registry

states: 2, controls: 2
states	zc (free), zt (free)
dynamics	dzc == (1 - zc)/(u2sf*v2) - k0*zc*exp(- ea/zt)
dynamics	dzt == (ztf - zt)/(u2sf*v2) + k0*zc*exp(- ea/zt) - a0*u1sf*v1*(zt - ztcw)
controls	v1 (free), v2 (free)
tracking stage cost	cost == 10*(zc - zc_ss)**2 + 2*(zt - zt_ss)**2 + (v1 - v1_ss)**2 + 0.5*(v2 - v2_ss)**2
steady-state targets	zc_ss (of zc), zt_ss (of zt)
steady-state control targets	v1_ss (of v1), v2_ss (of v2)

transformations

1. drto.dynamic_to_steady_state: removed=initial condition, tracking terminal cost, collapsed=7 Vars and 3 Constraints to a single point, derivatives=2 fixed at zero
2. drto.build_objective: objective=sum of 1 weighted cost terms
3. drto.steady_state_optimization: controls=free, tracking_weight=(one stage cost declared)

The equilibrium locus

The same one-line simulation, swept over the coolant flow at the reference residence time: the map of resting points the controller moves between.

import matplotlib.pyplot as plt
import numpy as np

# below v1 of about 0.34 no equilibrium exists within the variable bounds
# (too little cooling), so only converged solves make the locus
v1_eq, zc_eq, zt_eq = [], [], []
solver = SolverFactory("pounce")
for v1 in np.linspace(0.25, 1.0, 31):
    sim = pyo.TransformationFactory("drto.steady_state_simulation").create_using(
        m, controls={m.v1: float(v1), m.v2: 0.49989}
    )
    res = solver.solve(sim, load_solutions=False)
    if drto.scaling.solved_to_optimality(res):
        res.solution_loader.load_vars()
        v1_eq.append(v1)
        zc_eq.append(pyo.value(sim.zc))
        zt_eq.append(pyo.value(sim.zt))

fig, axes = plt.subplots(1, 2, figsize=(10.0, 3.2), sharex=True)
for ax, vals, name, ref in zip(axes, (zc_eq, zt_eq), ("zc", "zt"), (0.6416, 0.5387)):
    ax.plot(v1_eq, vals, "o-", color="C0", markersize=3)
    ax.axhline(ref, color="C0", linestyle=":")
    ax.axvline(0.57828, color="grey", linewidth=0.8)
    ax.set_title(name)
    ax.set_xlabel("v1")
fig.tight_layout()