Cold start on the IDAES CSTR

The saponification CSTR carries a full algebraic layer: the outlet stream, the reaction-rate cascade, and the indexed holdup’s water member. drto.cold_start_dynamic seeds the states on straight lines to their declared targets and solves that algebra pointwise, so the guess is consistent everywhere except the declared dynamics.

The plots show the model before the initialization and after it; the report shows the per-point solves working through the algebra.

%matplotlib inline
import pyomo.environ as pyo
import matplotlib.pyplot as plt
import drto
from drto import plot_controls, plot_stage_cost, plot_states
from idaes.core.util.model_statistics import degrees_of_freedom
from models.idaes_cstr import DC_START, F_IN, VOLUME, build

declared = build()
drto.info(declared)

drto registry

states: 5, controls: 12
horizon	fs._time (ContinuousSet, 11 points)
states	fs.cstr.control_volume.material_holdup_Liq_NaOH (free, mol), fs.cstr.control_volume.material_holdup_Liq_EthylAcetate (free, mol), fs.cstr.control_volume.material_holdup_Liq_SodiumAcetate (free, mol), fs.cstr.control_volume.material_holdup_Liq_Ethanol (free, mol), fs.cstr.control_volume.energy_holdup (free, J)
dynamics	fs.cstr.control_volume.material_accumulation[t,p,j] == fs.cstr.control_volume.properties_in[t].flow_vol*fs.cstr.control_volume.properties_in[t].conc_mol_comp[j] - fs.cstr.control_volume.properties_out[t].flow_vol*fs.cstr.control_volume.properties_out[t].conc_mol_comp[j] + fs.cstr.control_volume.rate_reaction_generation[t,p,j] for t in fs.cstr.control_volume.SetProduct_OrderedSet (mol/s)
dynamics	fs.cstr.control_volume.energy_accumulation[t,Liq] == fs.props.dens_mol*fs.props.cp_mol*fs.cstr.control_volume.properties_in[t].flow_vol*(fs.cstr.control_volume.properties_in[t].temperature - fs.props.temperature_ref) - fs.props.dens_mol*fs.props.cp_mol*fs.cstr.control_volume.properties_out[t].flow_vol*(fs.cstr.control_volume.properties_out[t].temperature - fs.props.temperature_ref) + fs.cstr.control_volume.heat[t] + (- fs.rxn.dh_rxn[R1]*fs.cstr.control_volume.rate_reaction_extent[t,R1]) for t in fs._time (W)
controls	fs.cstr.control_volume.heat (piecewise_constant, free, W), fs.cstr._flow_vol_inlet_ref (piecewise_constant, free, m**3/s)
tracking stage cost	cost[t] == ((fs.cstr.control_volume.energy_holdup[t,'Liq'] - eng_ss[Liq])/scale_E)**2 + ((fs.cstr.control_volume.heat[t] - duty_ss)/scale_Q)**2 + ((fs.cstr._flow_vol_inlet_ref[t] - flow_ss)/scale_F)**2 + ((fs.cstr.control_volume.material_holdup[t,'Liq','NaOH'] - ss_naoh)/scale_M)**2 + ((fs.cstr.control_volume.material_holdup[t,'Liq','EthylAcetate'] - ss_ea)/scale_M)**2 + ((fs.cstr.control_volume.material_holdup[t,'Liq','SodiumAcetate'] - ss_sa)/scale_M)**2 + ((fs.cstr.control_volume.material_holdup[t,'Liq','Ethanol'] - ss_etoh)/scale_M)**2 for t in sorted(fs._time)[:-1]
terminal cost	term == ((fs.cstr.control_volume.material_holdup[10.0,Liq,NaOH] - ss_naoh)/scale_M)**2 + ((fs.cstr.control_volume.material_holdup[10.0,Liq,EthylAcetate] - ss_ea)/scale_M)**2 + ((fs.cstr.control_volume.material_holdup[10.0,Liq,SodiumAcetate] - ss_sa)/scale_M)**2 + ((fs.cstr.control_volume.material_holdup[10.0,Liq,Ethanol] - ss_etoh)/scale_M)**2 + ((fs.cstr.control_volume.energy_holdup[10.0,Liq] - eng_ss[Liq])/scale_E)**2
initial conditions	mat0[j] == fs.cstr.control_volume.material_holdup[0.0,'Liq',j] for j in OrderedScalarSet (mol)
initial conditions	eng0[p] == fs.cstr.control_volume.energy_holdup[0.0,p] for p in fs.props.phase_list (J)
steady-state targets	ss_naoh (of fs.cstr.control_volume.material_holdup_Liq_NaOH, mol), ss_ea (of fs.cstr.control_volume.material_holdup_Liq_EthylAcetate, mol), ss_sa (of fs.cstr.control_volume.material_holdup_Liq_SodiumAcetate, mol), ss_etoh (of fs.cstr.control_volume.material_holdup_Liq_Ethanol, mol), eng_ss (of fs.cstr.control_volume.energy_holdup, J)
steady-state control targets	duty_ss (of fs.cstr.control_volume.heat, W), flow_ss (of fs.cstr._flow_vol_inlet_ref, m**3/s)

transformations (none)

The setpoint and the start

The steady solve fills the targets, and the initial condition is the setpoint composition knocked back along the reaction, hot.

def plot_rates(m):
    """The undeclared algebraic pair: the rate constant and the rate."""
    rxn = m.fs.cstr.control_volume.reactions
    grid = sorted(m.fs.time)
    fig, axes = plt.subplots(1, 2, figsize=(9, 3), constrained_layout=True)
    for ax, name, vals in (
        (axes[0], "k_rxn [m3/mol/s]",
         [pyo.value(rxn[t].k_rxn) for t in grid]),
        (axes[1], "reaction_rate[R1] [mol/m3/s]",
         [pyo.value(rxn[t].reaction_rate["R1"]) for t in grid]),
    ):
        ax.plot(grid, vals, ".-")
        ax.set_title(name)
        ax.set_xlabel("time [s]")

ss = pyo.TransformationFactory("drto.steady_state_simulation").create_using(
    declared, controls={declared.fs.cstr.control_volume.heat.name: 0.0,
                        declared.fs.cstr.inlet.flow_vol.name: F_IN})
assert degrees_of_freedom(ss) == 0
print(drto.initialize_steady_state(ss))
res = drto.scaled_solve(ss, tee=True)

cvs = ss.fs.cstr.control_volume
print(f"setpoint  T = {pyo.value(cvs.properties_out[0.0].temperature):7.3f} K"
      f"   NaOH = {pyo.value(cvs.properties_out[0.0].conc_mol_comp['NaOH']):6.3f} mol/m3")

# the controller models below are their own builds, so the numbers are
# collected here and written onto each of them
target = {name: pyo.value(cvs.material_holdup["Liq", j])
          for j, name in (("NaOH", "ss_naoh"), ("EthylAcetate", "ss_ea"),
                          ("SodiumAcetate", "ss_sa"), ("Ethanol", "ss_etoh"))}
# the start: the setpoint composition with conversion knocked back by
# DC_START along the reaction stoichiometry, hot by T_START - T_ss
start = {j: pyo.value(cvs.material_holdup["Liq", j]) + sgn * DC_START * VOLUME
         for j, sgn in (("NaOH", 1), ("EthylAcetate", 1),
                        ("SodiumAcetate", -1), ("Ethanol", -1))}
energy = {k: pyo.value(cvs.energy_holdup[k]) for k in declared.eng_ss}


def setpoint(model):
    """Write the steady solve's numbers onto a freshly built controller."""
    for name, value in target.items():
        model.component(name).set_value(value)
    for j, value in start.items():
        model.mat0[j] = value
    for k, value in energy.items():
        model.eng_ss[k] = value


drto.info(ss)

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

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

********************************************************************************
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.11.0, running with linear solver FERAL.

Reading C:\Users\Devin\AppData\Local\Temp\tmpazq_8uei\unknownanzxbsyi.nl...
Parsed 8 vars, 8 cons, jac_nnz=22, h_nnz=9 in 0.00s
Problem class: NLP. Selected solver: NLP filter line-search interior-point (pounce-nlp) [solver_selection=auto].

Number of nonzeros in equality constraint Jacobian...:       22
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        9

Total number of variables............................:        8
                     variables with only lower bounds:        5
                variables with lower and upper bounds:        1
                     variables with only upper bounds:        0
Total number of equality constraints.................:        8
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 2.60e-16 1.00e+00   -1.0 0.00e+00      - 0.00e+00 0.00e+00   0
   1  0.0000000e+00 2.08e-17 0.00e+00   -1.0 3.64e-14      - 9.91e-01 1.00e+00h  1


Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   2.0816681711721685e-17    2.0816681711721685e-17
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.0816681711721685e-17    2.0816681711721685e-17
  ...above the per-row floating-point noise floor:   0.0000000000000000e+00
  (the strict convergence test judges this value; the residual below it is finer than the row's own arithmetic can resolve. Set primal_noise_floor_kappa = 0 to disable.)


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

EXIT: Optimal Solution Found.

POUNCE 0.11.0: Optimal Solution Found.
Status: Solve_Succeeded
pounce: wrote C:\Users\Devin\AppData\Local\Temp\tmpazq_8uei\unknownanzxbsyi.sol
setpoint  T = 304.039 K   NaOH = 24.286 mol/m3

drto registry

states: 5, controls: 2
states	fs.cstr.control_volume.material_holdup_Liq_NaOH (free), fs.cstr.control_volume.material_holdup_Liq_EthylAcetate (free), fs.cstr.control_volume.material_holdup_Liq_SodiumAcetate (free), fs.cstr.control_volume.material_holdup_Liq_Ethanol (free), fs.cstr.control_volume.energy_holdup (free)
dynamics	fs.cstr.control_volume.material_accumulation[o,_r] == fs.cstr.control_volume.properties_in[0.0].flow_vol*fs.cstr.control_volume.properties_in[0.0].conc_mol_comp[_r] - fs.cstr.control_volume.properties_out[0.0].flow_vol*fs.cstr.control_volume.properties_out[0.0].conc_mol_comp[_r] + fs.cstr.control_volume.rate_reaction_generation[o,_r] for o in fs.props._phase_component_set (inc)
dynamics	fs.cstr.control_volume.energy_accumulation[Liq] == fs.props.dens_mol*fs.props.cp_mol*fs.cstr.control_volume.properties_in[0.0].flow_vol*(fs.cstr.control_volume.properties_in[0.0].temperature - fs.props.temperature_ref) - fs.props.dens_mol*fs.props.cp_mol*fs.cstr.control_volume.properties_out[0.0].flow_vol*(fs.cstr.control_volume.properties_out[0.0].temperature - fs.props.temperature_ref) + fs.cstr.control_volume.heat - fs.rxn.dh_rxn[R1]*fs.cstr.control_volume.rate_reaction_extent[R1] (inc)
controls	fs.cstr.control_volume.heat (fixed), fs.cstr._flow_vol_inlet_ref (fixed, m**3/s)
steady-state targets	ss_naoh (of fs.cstr.control_volume.material_holdup_Liq_NaOH, mol), ss_ea (of fs.cstr.control_volume.material_holdup_Liq_EthylAcetate, mol), ss_sa (of fs.cstr.control_volume.material_holdup_Liq_SodiumAcetate, mol), ss_etoh (of fs.cstr.control_volume.material_holdup_Liq_Ethanol, mol), eng_ss (of fs.cstr.control_volume.energy_holdup, J)
steady-state control targets	duty_ss (of fs.cstr.control_volume.heat, W), flow_ss (of fs.cstr._flow_vol_inlet_ref, m**3/s)

transformations

1. drto.dynamic_to_steady_state: removed=initial condition, tracking terminal cost, collapsed=8 Vars and 8 Constraints to a single point, derivatives=2 fixed at zero, blocks=3 time-indexed Block(s) collapsed to the steady member, 30 member(s) removed
2. drto.build_objective: objective=zero
3. drto.steady_state_simulation: controls=fs.cstr.control_volume.heat=0.0, fs.cstr._flow_vol_inlet_ref=1.0, dropped=tracking stage cost

The controller model

drto.dynamic_optimization takes the builder rather than a model and reaches the controller in one call: it discretizes the horizon at one finite element per declared sample, appends the terminal segment, and assembles the objective. The builder fixes the feed at the declared sample points, which is all an undiscretized model has, and the call fixes the members collocation adds to it. The cold start below runs on the assembled model and covers the tail with it. The steady reduction above branches from the plain declarations, which is why it runs on its own build.

m = drto.dynamic_optimization(build, N=10, infinite_horizon=True)
setpoint(m)
drto.info(m)

drto registry

states: 5, controls: 2
horizon	fs._time (ContinuousSet, 31 points)
states	fs.cstr.control_volume.material_holdup_Liq_NaOH (free, mol), fs.cstr.control_volume.material_holdup_Liq_EthylAcetate (free, mol), fs.cstr.control_volume.material_holdup_Liq_SodiumAcetate (free, mol), fs.cstr.control_volume.material_holdup_Liq_Ethanol (free, mol), fs.cstr.control_volume.energy_holdup (free, J)
dynamics	fs.cstr.control_volume.material_accumulation[t,p,j] == fs.cstr._flow_vol_inlet_ref[0.0]*fs.cstr.control_volume.properties_in[t].conc_mol_comp[j] - fs.cstr.control_volume.properties_out[t].flow_vol*fs.cstr.control_volume.properties_out[t].conc_mol_comp[j] + fs.cstr.control_volume.rate_reaction_generation[t,p,j] for t in fs.cstr.control_volume.SetProduct_OrderedSet (mol/s)
dynamics	fs.cstr.control_volume.energy_accumulation[t,Liq] == fs.props.dens_mol*fs.props.cp_mol*fs.cstr._flow_vol_inlet_ref[0.0]*(fs.cstr.control_volume.properties_in[t].temperature - fs.props.temperature_ref) - fs.props.dens_mol*fs.props.cp_mol*fs.cstr.control_volume.properties_out[t].flow_vol*(fs.cstr.control_volume.properties_out[t].temperature - fs.props.temperature_ref) + fs.cstr.control_volume.heat[0.0] - fs.rxn.dh_rxn[R1]*fs.cstr.control_volume.rate_reaction_extent[t,R1] for t in fs._time (W)
controls	fs.cstr.control_volume.heat (piecewise_constant, free, W), fs.cstr._flow_vol_inlet_ref (piecewise_constant, free, m**3/s)
tracking stage cost	cost[t] == ((fs.cstr.control_volume.energy_holdup[t,'Liq'] - eng_ss[Liq])/scale_E)**2 + ((fs.cstr.control_volume.heat[t] - duty_ss)/scale_Q)**2 + ((_flow_vol_inlet_ref[t] - flow_ss)/scale_F)**2 + ((fs.cstr.control_volume.material_holdup[t,'Liq','NaOH'] - ss_naoh)/scale_M)**2 + ((fs.cstr.control_volume.material_holdup[t,'Liq','EthylAcetate'] - ss_ea)/scale_M)**2 + ((fs.cstr.control_volume.material_holdup[t,'Liq','SodiumAcetate'] - ss_sa)/scale_M)**2 + ((fs.cstr.control_volume.material_holdup[t,'Liq','Ethanol'] - ss_etoh)/scale_M)**2 for t in sorted(fs._time)[:-1]
terminal cost	term == ((fs.cstr.control_volume.material_holdup[10.0,Liq,NaOH] - ss_naoh)/scale_M)**2 + ((fs.cstr.control_volume.material_holdup[10.0,Liq,EthylAcetate] - ss_ea)/scale_M)**2 + ((fs.cstr.control_volume.material_holdup[10.0,Liq,SodiumAcetate] - ss_sa)/scale_M)**2 + ((fs.cstr.control_volume.material_holdup[10.0,Liq,Ethanol] - ss_etoh)/scale_M)**2 + ((fs.cstr.control_volume.energy_holdup[10.0,Liq] - eng_ss[Liq])/scale_E)**2
initial conditions	mat0[j] == fs.cstr.control_volume.material_holdup[0.0,'Liq',j] for j in OrderedScalarSet (mol)
initial conditions	eng0[p] == fs.cstr.control_volume.energy_holdup[0.0,p] for p in fs.props.phase_list (J)
steady-state targets	ss_naoh (of fs.cstr.control_volume.material_holdup_Liq_NaOH, mol), ss_ea (of fs.cstr.control_volume.material_holdup_Liq_EthylAcetate, mol), ss_sa (of fs.cstr.control_volume.material_holdup_Liq_SodiumAcetate, mol), ss_etoh (of fs.cstr.control_volume.material_holdup_Liq_Ethanol, mol), eng_ss (of fs.cstr.control_volume.energy_holdup, J)
steady-state control targets	duty_ss (of fs.cstr.control_volume.heat, W), flow_ss (of fs.cstr._flow_vol_inlet_ref, m**3/s)

transformations

1. drto.infinite_horizon: segment=3 elements x 5 Legendre points, beta=1.2, gamma=0.01563827, profile=collocation, horizon=kept, infinite tail appended, algebraic=3 components replicated, blocks=3 time-indexed Block(s): 9 components, 3 equation families replicated, partial=algebraic entries of 2 indexed Var(s) copied per member, terminal_cost=terminal deactivated (the tail owns it), terminal=soft pin z(tau=1)=z_s on 5 states, mu=1000.0
2. drto.parameterize: profiles=fs.cstr.control_volume.heat (piecewise_constant), fs.cstr._flow_vol_inlet_ref (piecewise_constant)
3. drto.build_objective: objective=sum of 45 weighted cost terms
4. drto.dynamic_optimization: horizon=kept, tracking_weight=(one stage cost declared), sensitivity=5 initial-condition Params declared

Before

The builder’s values: flat, and unaware of the initial condition. The rate constant and the reaction rate sit at the property package’s defaults.

plot_states(m);
plt.show()

plot_rates(m)
plt.show()

plot_controls(m)
plt.show()

The cold start

States on their lines, controls at their targets, the tail at rest on the targets with the pin slacks at zero, and one block solve working the outlet stream, the reaction cascade, and the water member consistent at every grid point, the segment’s included. The solves run in the model’s own units: each is a square block solved in calculation order, so no scaling is involved.

print(drto.cold_start_dynamic(m, scale={"J": 1e7, "W": 1e6}))

drto cold_start_dynamic (interpolate to the targets)
  states        : 5 on a line across 31 grid points
  derivatives   : 155 at the profile's slope
  controls      : 2 at their targets
  segment       : 12 copies at the targets, slacks at zero
  point solves  : run (pyomo-pounce block solve)
  pyomo-pounce initialize (fill -> repair -> block-solve)
    decisions held: 300
    values filled : 0
    projection    : skipped
    pyomo-pounce block_initialize
      decisions fixed   : 0
      system square     : True
      blocks solved     : 1032 (1032 by Newton 1x1, 0 subsystem solves)
      vars initialized  : 1032
      left untouched    : 0 underdetermined, 0 overdetermined

After

States on their lines, controls at their targets, and the algebra, including the rate constant and the reaction rate, consistent with them at every grid point.

plot_states(m);
plt.show()

plot_rates(m)
plt.show()

plot_controls(m);
plt.show()

The first solve

The controller from the cold-started guess. The cold start above ran with the flowsheet’s operating magnitudes stated as its scale source, energy in joules near 1e7 and duties in watts near 1e6, since the factors have to exist before any values do. With every block solved, the initialized point is worth measuring: drto.scaled_solve reads the values the cold start left and hands the factors to pounce with the model. The log shows the iteration count from the linear start.

res = drto.scaled_solve(m, tee=True)
print(res.termination_condition.name)

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

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

********************************************************************************
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.11.0, running with linear solver FERAL.

Number of nonzeros in equality constraint Jacobian...:     4461
Number of nonzeros in inequality constraint Jacobian.:       20
Number of nonzeros in Lagrangian Hessian.............:      645

Total number of variables............................:     1408
                     variables with only lower bounds:      342
                variables with lower and upper bounds:       92
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1338
Total number of inequality constraints...............:        4
        inequality constraints with only lower bounds:        4
   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  1.3163721e+04 3.01e+07 1.00e+04   -1.0 0.00e+00      - 0.00e+00 0.00e+00   0
   1 -6.9232067e+03 5.60e+07 9.96e+03   -1.0 1.98e+10      - 3.00e-03 4.29e-03f  1
   2 -5.1980760e+04 1.03e+08 1.18e+04   -1.0 6.97e+09      - 9.04e-05 1.66e-02f  1
   3 -6.0362959e+04 1.02e+08 1.18e+04   -1.0 9.39e+09      - 1.46e-03 2.77e-03f  1
   4 -7.3847158e+04 1.02e+08 1.13e+04   -1.0 1.76e+09      - 1.63e-02 1.66e-02f  1
   5 -7.5323425e+04 1.00e+08 1.10e+04   -1.0 2.53e+08      - 8.51e-02 2.47e-02f  1
   6 -4.0688532e+04 5.77e+07 7.92e+03   -1.0 1.10e+08      - 1.27e-01 4.43e-01H  1
   7 -4.8483811e+04 4.16e+07 8.88e+03   -1.0 9.03e+07      - 4.36e-01 2.75e-01f  1
   8 -1.4877995e+04 1.29e+07 1.31e+04   -1.0 2.44e+07      - 6.05e-01 7.35e-01h  1
   9  2.1761144e+03 2.30e+07 2.86e+03   -1.0 4.96e+06      - 7.12e-01 1.00e+00h  1
iter      objective   inf_pr   inf_du lg(mu)    ||d|| lg(rg) alpha_du alpha_pr  ls
  10  3.2743356e+03 5.26e+06 1.37e+03   -1.0 3.85e+06      - 1.00e+00 1.00e+00h  1
  11  3.4997005e+03 2.24e+05 3.38e+02   -1.0 3.32e+06      - 9.67e-01 1.00e+00h  1
  12  3.5381984e+03 6.23e+03 2.20e+00   -1.0 2.90e+05      - 1.00e+00 1.00e+00h  1
  13  3.5386484e+03 2.07e+00 5.00e-04   -1.0 1.13e+03      - 1.00e+00 1.00e+00h  1
  14  3.5376444e+03 2.87e+01 8.41e-02   -2.5 9.32e+03      - 1.00e+00 1.00e+00f  1
  15  3.5376596e+03 3.53e-04 1.38e-06   -2.5 4.25e+01      - 1.00e+00 1.00e+00h  1
  16  3.5376328e+03 2.20e-02 6.50e-05   -3.8 2.59e+02      - 1.00e+00 1.00e+00f  1
  17  3.5376313e+03 6.80e-05 2.01e-07   -5.7 1.44e+01      - 1.00e+00 1.00e+00f  1
  18  3.5376313e+03 3.87e-07 3.14e-11   -8.6 1.78e-01      - 1.00e+00 1.00e+00h  1


Number of Iterations....: 18

                                   (scaled)                 (unscaled)
Objective...............:   3.5376313199950287e+03    3.5376313199950287e+03
Dual infeasibility......:   3.1411961817469705e-11    3.1411961817469705e-11
Constraint violation....:   8.0597204257792932e-13    3.8743019104003906e-07
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5131323994155104e-09    2.5131323994155104e-09
Overall NLP error.......:   2.5131323994155104e-09    3.8743019104003906e-07


Number of objective function evaluations             = 20
Number of objective gradient evaluations             = 19
Number of equality constraint evaluations            = 20
Number of inequality constraint evaluations          = 20
Number of equality constraint Jacobian evaluations   = 19
Number of inequality constraint Jacobian evaluations = 19
Number of Lagrangian Hessian evaluations             = 19
Number of linear solver quality escalations          = 1
Total seconds in POUNCE                              = 0.129

EXIT: Optimal Solution Found.

POUNCE 0.11.0: Optimal Solution Found.
convergenceCriteriaSatisfied

plot_states(m);
plt.show()

plot_controls(m);
plt.show()

From scratch, with an exponential profile

The same model rebuilt and refilled, cold-started on the normalized exponential decay instead of the line (time constant 3 s), and the same solve again: compare the iteration count with the linear start above.

m2 = drto.dynamic_optimization(build, N=10, infinite_horizon=True)
setpoint(m2)

print(drto.cold_start_dynamic(m2, profile="exponential", time_constant=.3,
                              scale={"J": 1e7, "W": 1e6}))

drto cold_start_dynamic (interpolate to the targets)
  states        : 5 on an exponential decay (time constant 0.3) across 31 grid points
  derivatives   : 155 at the profile's slope
  controls      : 2 at their targets
  segment       : 12 copies at the targets, slacks at zero
  point solves  : run (pyomo-pounce block solve)
  pyomo-pounce initialize (fill -> repair -> block-solve)
    decisions held: 300
    values filled : 0
    projection    : skipped
    pyomo-pounce block_initialize
      decisions fixed   : 0
      system square     : True
      blocks solved     : 1032 (1032 by Newton 1x1, 0 subsystem solves)
      vars initialized  : 1032
      left untouched    : 0 underdetermined, 0 overdetermined

plot_states(m2)
plt.show()

res = drto.scaled_solve(m2, tee=True)
print(res.termination_condition.name)

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

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

********************************************************************************
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.11.0, running with linear solver FERAL.

Number of nonzeros in equality constraint Jacobian...:     4461
Number of nonzeros in inequality constraint Jacobian.:       20
Number of nonzeros in Lagrangian Hessian.............:      645

Total number of variables............................:     1408
                     variables with only lower bounds:      342
                variables with lower and upper bounds:       92
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1338
Total number of inequality constraints...............:        4
        inequality constraints with only lower bounds:        4
   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  3.4982381e+03 1.58e+08 1.00e+04   -1.0 0.00e+00      - 0.00e+00 0.00e+00   0
   1  3.3263197e+03 1.56e+08 9.87e+03   -1.0 3.58e+09      - 3.62e-02 1.27e-02f  1
   2  1.3595565e+02 1.55e+08 9.77e+03   -1.0 2.18e+09      - 1.69e-04 1.10e-02f  1
   3 -6.3067936e+03 1.28e+08 8.07e+03   -1.0 7.78e+08      - 4.41e-02 1.73e-01f  1
   4 -1.5827093e+04 1.21e+08 7.62e+03   -1.0 6.17e+08      - 1.05e-01 5.62e-02f  1
   5 -3.8786105e+04 6.66e+07 4.30e+03   -1.0 3.57e+08      - 9.98e-02 4.47e-01f  1
   6 -2.8984780e+04 4.22e+07 2.99e+03   -1.0 1.29e+08      - 2.15e-01 3.67e-01h  1
   7  8.3912648e+02 4.77e+07 6.91e+03   -1.0 4.30e+07      - 3.49e-01 1.00e+00h  1
   8  3.2175734e+03 3.75e+06 8.63e+02   -1.0 1.31e+07      - 5.66e-01 9.87e-01h  1
   9  3.5306518e+03 4.73e+05 1.10e+02   -1.0 5.68e+05      - 9.29e-01 1.00e+00h  1
iter      objective   inf_pr   inf_du lg(mu)    ||d|| lg(rg) alpha_du alpha_pr  ls
  10  3.5386059e+03 2.68e+02 6.88e-01   -1.0 1.33e+05      - 1.00e+00 1.00e+00h  1
  11  3.5378217e+03 1.92e+01 9.69e-03   -1.7 7.64e+03      - 1.00e+00 1.00e+00f  1
  12  3.5376322e+03 1.23e+00 6.24e-04   -3.8 1.95e+03      - 1.00e+00 1.00e+00f  1
  13  3.5376313e+03 7.81e-05 4.18e-08   -5.7 1.62e+01      - 1.00e+00 1.00e+00h  1
  14  3.5376313e+03 4.02e-07 5.48e-12   -8.6 1.79e-01      - 1.00e+00 1.00e+00h  1


Number of Iterations....: 14

                                   (scaled)                 (unscaled)
Objective...............:   3.5376313199950309e+03    3.5376313199950309e+03
Dual infeasibility......:   5.4805221775741653e-12    5.4805221775741653e-12
Constraint violation....:   1.0301564289105780e-12    4.0206442675216891e-07
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5076045313673449e-09    2.5076045313673449e-09
Overall NLP error.......:   2.5076045313673449e-09    4.0206442675216891e-07


Number of objective function evaluations             = 15
Number of objective gradient evaluations             = 15
Number of equality constraint evaluations            = 15
Number of inequality constraint evaluations          = 15
Number of equality constraint Jacobian evaluations   = 15
Number of inequality constraint Jacobian evaluations = 15
Number of Lagrangian Hessian evaluations             = 15
Number of linear solver quality escalations          = 1
Total seconds in POUNCE                              = 0.161

EXIT: Optimal Solution Found.

POUNCE 0.11.0: Optimal Solution Found.
convergenceCriteriaSatisfied

plot_states(m2);
plt.show()