Hicks-Ray CSTR: forward integration#
drto.dynamic_simulation frees nothing. It applies the declared
control profiles, fixes the controls, and installs the constant-zero objective,
so the solver is handed the square forward integration of the model exactly as
declared. The horizon is kept. This page calls it as a function taking the
model statement, which builds the model, discretizes it at one finite element
per declared sample, and applies the transformation, so one call reaches a
system ready to solve.
The declared model is the Hicks-Ray CSTR from models/hicks.py.
Holding its controls at their declared steady-state values and integrating from
different starting states shows why this reactor is the standard NMPC example:
the declared setpoint is a real equilibrium, but an unstable one, sitting
between two stable branches. Case 1 starts on it and stays. Case 2 starts just
below it and falls away. The closing cell perturbs both directions.
The mode installs no cost, so the stage and terminal cost equations leave the model and there is no stage-cost panel to plot.
Case 1: started on the declared steady state#
With the controls held at v1_ss and v2_ss and the state started at the
declared target, nothing should move. This is the check that the declared
equilibrium really is one.
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
from drto import plot_states, plot_controls
# the declared control targets, read once off a one-interval build. The one
# call below returns the model already prepared, so the values it fixes the
# controls at have to be in hand before it runs
V1_SS = pyo.value(hicks(N=1).v1_ss)
V2_SS = pyo.value(hicks(N=1).v2_ss)
def simulate(zc0, zt0, v1=None, N=50):
"""Integrate the CSTR forward from (zc0, zt0) with the controls held."""
m = drto.dynamic_simulation(
hicks, N=N,
controls={"v1": V1_SS if v1 is None else v1, "v2": V2_SS})
m.zc_hat.set_value(zc0)
m.zt_hat.set_value(zt0)
return m
m1 = simulate(0.6416, 0.5387) # exactly the declared steady state
SolverFactory("pounce").solve(m1, tee=True)
********************************************************************************
#### ### / # /# #/ #### #####
# # # # /# #/ ## / # #
#### # #/ # / # #/# # ####
# # / # /# # /## # #
# ##/ #/# #/ # #### #####
********************************************************************************
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\tmpwuuqyuw7\unknownr12nwloo.nl...
Parsed 602 vars, 602 cons, jac_nnz=2396, h_nnz=300 in 0.02s
Problem class: NLP. Selected solver: NLP filter line-search interior-point (pounce-nlp) [solver_selection=auto].
Number of nonzeros in equality constraint Jacobian...: 2396
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 300
Total number of variables............................: 602
variables with only lower bounds: 150
variables with lower and upper bounds: 150
variables with only upper bounds: 0
Total number of equality constraints.................: 602
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
pounce: wrote C:\Users\Devin\AppData\Local\Temp\tmpwuuqyuw7\unknownr12nwloo.sol
1 0.0000000e+00 9.89e-08 6.43e-04 -1.7 3.50e-04 - 1.00e+00 1.00e+00h 1
2 0.0000000e+00 6.54e-13 3.86e-07 -3.8 7.80e-07 - 1.00e+00 1.00e+00h 1
3 0.0000000e+00 1.79e-15 1.91e-11 -5.7 2.05e-12 - 1.00e+00 1.00e+00h 1
4 0.0000000e+00 1.46e-15 3.35e-14 -8.6 8.09e-14 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 4
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 3.3470486581962944e-14 3.3470486581962944e-14
Constraint violation....: 1.4628199706775459e-15 1.4628199706775459e-15
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5059035599189102e-09 2.5059035599189102e-09
Overall NLP error.......: 2.5059035599189102e-09 2.5059035599189102e-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.012
EXIT: Optimal Solution Found.
POUNCE 0.10.0: Optimal Solution Found.
Status: Solve_Succeeded
<pyomo.contrib.solver.common.results.Results at 0x25b5a4f2030>
drto.info(m1)
drto registry
transformations
| states: 2, controls: 2 | |
| horizon | t (ContinuousSet, 151 points) |
| states | zc (free), zt (free) |
| dynamics | dzc[t] == (1 - zc[t])/(u2sf*v2[t]) - k0*zc[t]*exp(- ea/zt[t]) for t in t |
| dynamics | dzt[t] == (ztf - zt[t])/(u2sf*v2[t]) + k0*zc[t]*exp(- ea/zt[t]) - a0*u1sf*v1[t]*(zt[t] - ztcw) for t in t |
| controls | v1 (piecewise_constant, fixed), v2 (piecewise_constant, fixed) |
| initial conditions | zc[0] == zc_hat |
| initial conditions | zt[0] == zt_hat |
| steady-state targets | zc_ss (of zc), zt_ss (of zt) |
| steady-state control targets | v1_ss (of v1), v2_ss (of v2) |
drto.parameterize: profiles=v1 (piecewise_constant), v2 (piecewise_constant)drto.build_objective: objective=zerodrto.dynamic_simulation: horizon=kept, controls=v1=0.57828, v2=0.49989, dropped=tracking stage cost, tracking terminal cost
tN = m1.t.last()
print(f"start zc {pyo.value(m1.zc[0]):.5f} zt {pyo.value(m1.zt[0]):.5f}")
print(f"end zc {pyo.value(m1.zc[tN]):.5f} zt {pyo.value(m1.zt[tN]):.5f}")
print(f"declared zc {pyo.value(m1.zc_ss):.5f} zt {pyo.value(m1.zt_ss):.5f}")
start zc 0.64160 zt 0.53870
end zc 0.64125 zt 0.53896
declared zc 0.64160 zt 0.53870
plot_states(m1);
Case 2: started just below it#
The model’s own initial state (zc_hat 0.625, zt_hat 0.525) sits a little
below the target. With the same held controls the reactor does not return to
the setpoint, it cools and runs out to the low-conversion branch. The setpoint
is unstable in open loop, which is exactly what drto.dynamic_optimization is for.
m2 = simulate(0.625, 0.525) # the model's declared initial state
SolverFactory("pounce").solve(m2, tee=True)
********************************************************************************
#### ### / # /# #/ #### #####
# # # # /# #/ ## / # #
#### # #/ # / # #/# # ####
# # / # /# # /## # #
# ##/ #/# #/ # #### #####
********************************************************************************
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\tmpwiz3h3u9\unknown1v0fxnhr.nl...
Parsed 602 vars, 602 cons, jac_nnz=2396, h_nnz=300 in 0.02s
Problem class: NLP. Selected solver: NLP filter line-search interior-point (pounce-nlp) [solver_selection=auto].
Number of nonzeros in equality constraint Jacobian...: 2396
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 300
Total number of variables............................: 602
variables with only lower bounds: 150
variables with lower and upper bounds: 150
variables with only upper bounds: 0
Total number of equality constraints.................: 602
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 6.87e-02 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1r 0.0000000e+00 6.87e-02 9.99e+02 -1.2 0.00e+00 - 0.00e+00 7.90e-07R 14
2r 0.0000000e+00 3.76e-02 9.87e+02 -1.2 3.66e+01 - 3.88e-02 1.44e-03f 1
3r 0.0000000e+00 4.12e-02 9.73e+02 -1.2 6.31e-01 - 2.08e-01 2.49e-02f 1
4r 0.0000000e+00 4.21e-02 7.39e+02 -1.2 5.78e-02 - 7.13e-01 2.40e-01f 1
5r 0.0000000e+00 2.02e-02 4.83e+02 -1.2 1.32e-01 - 2.47e-02 1.00e+00f 1
6r 0.0000000e+00 1.81e-02 1.99e+02 -1.2 2.27e-02 4.0 5.94e-01 5.94e-01h 1
7r 0.0000000e+00 2.48e-02 1.55e+02 -1.2 3.47e-02 3.5 1.38e-01 4.04e-01f 1
8r 0.0000000e+00 1.59e-02 7.02e+02 -1.2 8.44e-02 - 1.00e+00 1.06e-01f 1
9r 0.0000000e+00 1.61e-02 5.21e+00 -1.2 8.51e-03 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10r 0.0000000e+00 1.53e-02 3.06e-01 -1.2 1.71e-02 - 1.00e+00 1.00e+00h 1
11r 0.0000000e+00 1.67e-02 3.30e-01 -2.8 2.55e-02 - 9.38e-01 1.00e+00f 1
12r 0.0000000e+00 1.67e-02 1.18e+00 -2.8 1.06e-03 3.0 1.00e+00 1.00e+00h 1
13r 0.0000000e+00 1.67e-02 8.81e-02 -2.8 7.63e-03 - 1.00e+00 1.00e+00h 1
14r 0.0000000e+00 1.68e-02 2.52e-02 -4.2 6.79e-05 2.6 1.00e+00 1.00e+00f 1
15r 0.0000000e+00 1.68e-02 1.95e-02 -4.2 1.58e-04 2.1 1.00e+00 1.00e+00h 1
16r 0.0000000e+00 1.68e-02 1.96e-02 -4.2 4.75e-04 1.6 1.00e+00 1.00e+00h 1
17r 0.0000000e+00 1.68e-02 1.97e-02 -4.2 1.44e-03 1.1 1.00e+00 1.00e+00h 1
18r 0.0000000e+00 1.68e-02 2.02e-02 -4.2 4.42e-03 0.7 1.00e+00 1.00e+00h 1
19r 0.0000000e+00 1.68e-02 2.19e-02 -4.2 1.44e-02 0.2 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20r 0.0000000e+00 1.68e-02 3.00e-02 -4.2 5.90e-02 -0.3 1.00e+00 1.00e+00h 1
21r 0.0000000e+00 1.68e-02 2.93e-02 -4.2 2.14e-02 0.1 1.00e+00 1.00e+00h 1
22r 0.0000000e+00 2.06e-02 7.71e-02 -4.2 1.59e-01 -0.3 1.00e+00 1.00e+00h 1
23r 0.0000000e+00 1.66e-02 2.16e-01 -4.2 1.29e-01 0.1 1.00e+00 1.00e+00h 1
24r 0.0000000e+00 1.56e-02 5.67e+00 -4.2 3.22e-01 0.5 1.00e+00 6.50e-01f 1
25r 0.0000000e+00 1.32e-02 7.11e+01 -4.2 8.67e-02 1.8 2.48e-03 1.00e+00f 1
26r 0.0000000e+00 8.53e-03 8.10e+00 -4.2 4.30e-02 2.3 1.00e+00 9.23e-01f 1
27r 0.0000000e+00 1.71e-03 4.44e+01 -4.2 2.25e-02 2.7 1.00e+00 3.91e-01f 1
pounce: wrote C:\Users\Devin\AppData\Local\Temp\tmpwiz3h3u9\unknown1v0fxnhr.sol
28r 0.0000000e+00 1.71e-03 6.85e+02 -4.2 8.06e-02 2.2 1.00e+00 4.23e-05f 1
29r 0.0000000e+00 1.87e-04 3.20e+02 -4.2 6.61e-03 - 1.00e+00 6.80e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r 0.0000000e+00 3.05e-06 2.79e-01 -4.2 3.30e-03 - 1.00e+00 1.00e+00h 1
31r 0.0000000e+00 3.87e-10 2.23e-05 -4.2 4.12e-05 - 1.00e+00 1.00e+00h 1
32r 0.0000000e+00 3.49e-11 3.63e-10 -6.3 6.42e-08 - 1.00e+00 1.00e+00f 1
33r 0.0000000e+00 2.16e-14 1.33e-11 -9.0 8.40e-10 - 1.00e+00 1.00e+00f 1
33 0.0000000e+00 2.16e-14 8.94e-01 -1.7 5.49e+01 - 0.00e+00 7.90e-07R 14
34 0.0000000e+00 1.71e-15 2.00e-07 -1.7 5.21e-13 - 1.00e+00 1.00e+00h 1
35 0.0000000e+00 2.12e-15 1.50e-09 -3.8 6.75e-16 - 1.00e+00 1.00e+00t 0
36 0.0000000e+00 1.60e-15 1.85e-11 -5.7 1.37e-15 - 1.00e+00 1.00e+00t 0
37 0.0000000e+00 1.81e-15 2.51e-14 -8.6 1.21e-15 - 1.00e+00 1.00e+00t 0
Number of Iterations....: 37
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 2.5144326140449122e-14 2.5144326140449122e-14
Constraint violation....: 1.8067145002298446e-15 1.8067145002298446e-15
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5059035597474844e-09 2.5059035597474844e-09
Overall NLP error.......: 2.5059035597474844e-09 2.5059035597474844e-09
Number of objective function evaluations = 57
Number of objective gradient evaluations = 8
Number of equality constraint evaluations = 59
Number of inequality constraint evaluations = 59
Number of equality constraint Jacobian evaluations = 41
Number of inequality constraint Jacobian evaluations = 41
Number of Lagrangian Hessian evaluations = 37
Total seconds in POUNCE = 0.106
EXIT: Optimal Solution Found.
POUNCE 0.10.0: Optimal Solution Found.
Status: Solve_Succeeded
<pyomo.contrib.solver.common.results.Results at 0x25b5d4f0c30>
plot_states(m2);
plot_controls(m2);
Both directions, and the per-point control option#
Perturbed the other way the reactor runs to the hot, high-conversion branch, so
the declared setpoint is a saddle between the two. The last row also exercises
the sequence form of the controls option: one value per free point the
applied profile leaves, here fifty piecewise-constant coolant moves stepping up
at sample ten.
step = [V1_SS] * 10 + [V1_SS * 1.15] * 40 # more coolant from sample 10 on
runs = (("on the setpoint ", 0.6416, 0.5387, None),
("below it ", 0.625, 0.525, None),
("above it ", 0.655, 0.552, None),
("below it, coolant step", 0.625, 0.525, step))
for label, zc0, zt0, v1 in runs:
m = simulate(zc0, zt0, v1=v1)
SolverFactory("pounce").solve(m)
t = m.t.last()
print(f"{label} start ({zc0:.4f}, {zt0:.4f})"
f" end ({pyo.value(m.zc[t]):.5f}, {pyo.value(m.zt[t]):.5f})")
on the setpoint start (0.6416, 0.5387) end (0.64125, 0.53896)
below it start (0.6250, 0.5250) end (0.93045, 0.39719)
above it start (0.6550, 0.5520) end (0.10193, 0.76676)
below it, coolant step start (0.6250, 0.5250) end (0.93386, 0.39436)