The Klatt-Engell CSTR: the infinite-horizon problem#
The NMPC problem of Lueken, Brandner & Lucia (2023, doi:10.1016/j.ifacol.2023.10.545) with the terminal cost replaced by drto’s terminal segment: one solve from an off-setpoint state, and a training point read off the converged factorization.
models/klatt_engell.py declares the reactor:
four states, two controls, and a tracking stage cost over all six with
Param weights. The steady solve below runs with the paper’s weights
(the two concentrations at one, the rest at zero); after it writes the
equilibrium into the targets, every weight is set to one, so the
horizon problem tracks the computed equilibrium in every state and
control and the tail integrand vanishes there.
from functools import partial
import pyomo.environ as pyo
from pyomo.contrib.solver.common.factory import SolverFactory
import pyomo_pounce # registers the pounce solver with Pyomo
import drto
from models.klatt_engell import klatt_engell
m = klatt_engell(N=1, h=0.005, move_penalty=False)
drto.info(m)
drto registry
transformations (none)
| states: 4, controls: 2 | |
| horizon | t (ContinuousSet, 2 points) |
| states | c_a (free), c_b (free), T_R (free), T_K (free) |
| dynamics | dc_a[t] == F[t]*(c_in - c_a[t]) - k0_1*exp(- e1/(T_R[t] + 273.15))*c_a[t] - k0_3*exp(- e3/(T_R[t] + 273.15))*c_a[t]**2 + w_ca[t] for t in t |
| dynamics | dc_b[t] == - F[t]*c_b[t] + k0_1*exp(- e1/(T_R[t] + 273.15))*c_a[t] - k0_2*exp(- e2/(T_R[t] + 273.15))*c_b[t] + w_cb[t] for t in t |
| dynamics | dT_R[t] == F[t]*(t_in - T_R[t]) + k_w*a_r/(rho*cp_r*v_r)*(T_K[t] - T_R[t]) - (k0_1*exp(- e1/(T_R[t] + 273.15))*c_a[t]*dh1 + k0_2*exp(- e2/(T_R[t] + 273.15))*c_b[t]*dh2 + k0_3*exp(- e3/(T_R[t] + 273.15))*c_a[t]**2*dh3)/(rho*cp_r) + w_tr[t] for t in t |
| dynamics | dT_K[t] == (Q_dot[t] + k_w*a_r*(T_R[t] - T_K[t]))/(m_k*cp_k) + w_tk[t] for t in t |
| controls | F (piecewise_constant, free), Q_dot (piecewise_constant, free) |
| tracking stage cost | cost[t] == q_ca*(c_a[t] - ca_sp)**2 + q_cb*(c_b[t] - cb_sp)**2 + q_tr*(T_R[t] - tr_sp)**2 + q_tk*(T_K[t] - tk_sp)**2 + r_f*(F[t] - f_sp)**2 + r_q*(Q_dot[t] - q_sp)**2 for t in sorted(t)[:-1] |
| terminal cost | term == q_ca*(c_a[0.005] - ca_sp)**2 + q_cb*(c_b[0.005] - cb_sp)**2 + q_tr*(T_R[0.005] - tr_sp)**2 + q_tk*(T_K[0.005] - tk_sp)**2 |
| initial conditions | c_a[0.0] == ca_hat |
| initial conditions | c_b[0.0] == cb_hat |
| initial conditions | T_R[0.0] == tr_hat |
| initial conditions | T_K[0.0] == tk_hat |
| steady-state targets | ca_sp (of c_a), cb_sp (of c_b), tr_sp (of T_R), tk_sp (of T_K) |
| steady-state control targets | f_sp (of F), q_sp (of Q_dot) |
| disturbances | w_ca (piecewise_constant), w_cb (piecewise_constant), w_tr (piecewise_constant), w_tk (piecewise_constant) |
The setpoint#
The paper’s setpoint pair is not an equilibrium of this reactor: k_1 and k_2 share their pre-exponential and activation energy, so holding c_A = 0.7 and c_B = 0.6 at steady state forces a temperature far above the 150 C bound, and a solve with both concentrations fixed converges to local infeasibility. The steady point is instead chosen by the same tracking cost: a clone is reduced, the controls left free, and the reduced stage cost minimized over the steady states within bounds.
Every steady-state target takes the solution’s values, the concentration setpoints included: the tracking cost then vanishes at an attainable equilibrium, which is the point the terminal segment pins. The paper’s 0.7 and 0.6 remain in the problem as the criterion that selected it.
sm = pyo.TransformationFactory("drto.steady_state_optimization").create_using(m)
print(drto.initialize_steady_state(sm))
res = SolverFactory("pounce").solve(sm, tee=True)
print("steady solve:", res.termination_condition.name)
print("steady cost:", pyo.value(sm.drto_objective))
print()
# each steady_state pairing names a component and its target Param; the
# reduced model's records line up by declaration order. The values are
# kept by target name, since the horizon model is built fresh below
reg, reg_s = drto.info(m), drto.info(sm)
setpoint = {}
for kind in ("steady_state", "steady_state_control"):
for here, there in zip(reg.declarations(kind), reg_s.declarations(kind)):
v = there["of"]
value = pyo.value(next(iter(v.values())) if v.is_indexed() else v)
setpoint[here["component"].local_name] = value
print(f"{here['component'].name} = {value:.6g}")
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 : 5
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\tmp23z_7wuu\unknowndws_irsb.nl...
Parsed 7 vars, 5 cons, jac_nnz=18, h_nnz=8 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...: 18
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 8
Total number of variables............................: 7
variables with only lower bounds: 0
variables with lower and upper bounds: 6
variables with only upper bounds: 0
Total number of equality constraints.................: 5
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 1.5937952e+00 3.41e-13 5.40e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.5223143e+00 1.45e-01 1.68e-02 -1.0 1.32e+00 - 9.65e-01 1.00e+00f 1
2 6.4083488e-01 2.30e+01 3.24e-01 -1.0 1.63e+01 - 4.76e-01 1.00e+00f 1
3 7.6759535e-02 1.93e+01 1.49e+00 -1.0 1.06e+02 - 9.01e-01 1.00e+00h 1
4 4.7624828e-01 8.44e+00 1.97e+00 -1.0 8.95e+00 -4.0 7.89e-01 1.00e+00h 1
5 1.2500415e-02 1.49e+01 7.82e-01 -1.0 5.54e+01 - 1.00e+00 1.00e+00h 1
6 -6.7244642e-02 6.44e+00 1.99e-01 -1.0 5.33e+01 - 1.00e+00 1.00e+00h 1
7 4.6376099e-03 2.30e+00 3.34e-02 -1.0 7.80e+02 - 1.00e+00 1.00e+00h 1
8 -6.4454072e-04 1.31e-01 1.09e-02 -1.7 8.85e+01 - 9.83e-01 1.00e+00h 1
9 -6.1648133e-04 5.47e-02 1.85e-04 -2.5 1.15e+02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 -3.0274952e-05 3.45e-03 1.13e-05 -3.8 1.75e+01 - 1.00e+00 1.00e+00h 1
11 8.3692460e-06 8.27e-05 3.35e-07 -5.7 4.06e+00 - 1.00e+00 1.00e+00h 1
12 8.1002096e-06 4.56e-04 2.93e-05 -5.7 3.24e+02 - 1.00e+00 1.00e+00H 1
13 8.1291877e-06 3.25e-05 3.83e-07 -5.7 3.06e+01 - 1.00e+00 1.00e+00H 1
14 8.0461194e-06 5.23e-02 1.86e-06 -8.6 8.64e+01 - 8.81e-01 1.00e+00h 1
15 7.6575166e-06 1.18e+00 3.61e-05 -8.6 4.29e+02 - 8.72e-01 1.00e+00h 1
16 7.6014466e-06 2.04e-01 5.25e-06 -8.6 1.49e+02 - 1.00e+00 8.69e-01h 1
17 7.6068348e-06 3.05e-05 8.33e-09 -8.6 1.98e+00 - 1.00e+00 1.00e+00h 1
18 7.6067523e-06 7.03e-08 1.54e-12 -8.6 1.25e-01 - 1.00e+00 1.00e+00h 1
19 7.6051513e-06 2.48e-05 6.25e-10 -9.0 2.33e+00 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 7.6051526e-06 8.28e-12 4.83e-16 -9.0 7.03e-04 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 20
(scaled) (unscaled)
Objective...............: 7.6051526281730388e-06 7.6051526281730388e-06
Dual infeasibility......: 4.8299345172901641e-16 4.8299345172901641e-16
Constraint violation....: 2.0321522242738865e-12 8.2849282989627682e-12
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 9.0909092452500373e-10 9.0909092452500373e-10
Overall NLP error.......: 9.0909092452500373e-10 9.0909092452500373e-10
Number of objective function evaluations = 25
Number of objective gradient evaluations = 1
Number of equality constraint evaluations = 25
Number of inequality constraint evaluations = 25
Number of equality constraint Jacobian evaluations = 21
Number of inequality constraint Jacobian evaluations = 21
Number of Lagrangian Hessian evaluations = 21
Total seconds in POUNCE = 0.007
EXIT: Optimal Solution Found.
POUNCE 0.11.0: Optimal Solution Found.
Status: Solve_Succeeded
pounce: wrote C:\Users\Devin\AppData\Local\Temp\tmp23z_7wuu\unknowndws_irsb.sol
steady solve: convergenceCriteriaSatisfied
steady cost: 7.605152628173039e-06
ca_sp = 0.701575
cb_sp = 0.597736
tr_sp = 140.82
tk_sp = 139.995
f_sp = 12.9428
q_sp = -715.563
Assemble the horizon#
drto.dynamic_optimization takes the model statement itself. It calls
klatt_engell, discretizes the declared grid at one finite element per
step, appends the terminal segment, and assembles the optimization, so
the four calls this workflow used to spell out are one. The equilibrium
found above and the do-mpc example’s starting state are written into the
model it returns, every tracking term carries unit weight, and the cold
start runs on it. drto.dynamic_optimization declares the
initial-condition Params as pounce sensitivity parameters.
m = drto.dynamic_optimization(
partial(klatt_engell, move_penalty=False), # the page's own option
N=1,
h=0.005,
ncp=2, # do-mpc's mesh
infinite_horizon={"nfe": 3, "ncp": 5},
)
for name, value in setpoint.items(): # the equilibrium from above
m.component(name).set_value(value)
for w in (m.q_ca, m.q_cb, m.q_tr, m.q_tk, m.r_f, m.r_q):
w.set_value(1.0)
# the initial state the horizon starts from (the do-mpc example's)
m.ca_hat.set_value(0.75) # mol/L
m.cb_hat.set_value(0.75) # mol/L
m.tr_hat.set_value(120) # C
m.tk_hat.set_value(120) # C
drto.cold_start_dynamic(m, profile="exponential", time_constant=.3)
drto.info(m)
drto registry
transformations
| states: 4, controls: 2 | |
| horizon | t (ContinuousSet, 3 points) |
| states | c_a (free), c_b (free), T_R (free), T_K (free) |
| dynamics | dc_a[t] == F[t]*(c_in - c_a[t]) - k0_1*exp(- e1/(T_R[t] + 273.15))*c_a[t] - k0_3*exp(- e3/(T_R[t] + 273.15))*c_a[t]**2 + w_ca[t] for t in t |
| dynamics | dc_b[t] == - F[t]*c_b[t] + k0_1*exp(- e1/(T_R[t] + 273.15))*c_a[t] - k0_2*exp(- e2/(T_R[t] + 273.15))*c_b[t] + w_cb[t] for t in t |
| dynamics | dT_R[t] == F[t]*(t_in - T_R[t]) + k_w*a_r/(rho*cp_r*v_r)*(T_K[t] - T_R[t]) - (k0_1*exp(- e1/(T_R[t] + 273.15))*c_a[t]*dh1 + k0_2*exp(- e2/(T_R[t] + 273.15))*c_b[t]*dh2 + k0_3*exp(- e3/(T_R[t] + 273.15))*c_a[t]**2*dh3)/(rho*cp_r) + w_tr[t] for t in t |
| dynamics | dT_K[t] == (Q_dot[t] + k_w*a_r*(T_R[t] - T_K[t]))/(m_k*cp_k) + w_tk[t] for t in t |
| controls | F (piecewise_constant, free), Q_dot (piecewise_constant, free) |
| tracking stage cost | cost[t] == q_ca*(c_a[t] - ca_sp)**2 + q_cb*(c_b[t] - cb_sp)**2 + q_tr*(T_R[t] - tr_sp)**2 + q_tk*(T_K[t] - tk_sp)**2 + r_f*(F[t] - f_sp)**2 + r_q*(Q_dot[t] - q_sp)**2 for t in sorted(t)[:-1] |
| terminal cost | term == q_ca*(c_a[0.005] - ca_sp)**2 + q_cb*(c_b[0.005] - cb_sp)**2 + q_tr*(T_R[0.005] - tr_sp)**2 + q_tk*(T_K[0.005] - tk_sp)**2 |
| initial conditions | c_a[0.0] == ca_hat |
| initial conditions | c_b[0.0] == cb_hat |
| initial conditions | T_R[0.0] == tr_hat |
| initial conditions | T_K[0.0] == tk_hat |
| steady-state targets | ca_sp (of c_a), cb_sp (of c_b), tr_sp (of T_R), tk_sp (of T_K) |
| steady-state control targets | f_sp (of F), q_sp (of Q_dot) |
| disturbances | w_ca (piecewise_constant), w_cb (piecewise_constant), w_tr (piecewise_constant), w_tk (piecewise_constant) |
drto.infinite_horizon: segment=3 elements x 5 Legendre points, beta=1.2, gamma=3.12765494, profile=collocation, horizon=kept, infinite tail appended, disturbances=w_ca fixed at 0.0, w_cb fixed at 0.0, w_tr fixed at 0.0, w_tk fixed at 0.0, terminal_cost=terminal deactivated (the tail owns it), terminal=soft pin z(tau=1)=z_s on 4 states, mu=1000.0drto.parameterize: profiles=F (piecewise_constant), Q_dot (piecewise_constant), w_ca (piecewise_constant), w_cb (piecewise_constant), w_tr (piecewise_constant), w_tk (piecewise_constant)drto.build_objective: objective=sum of 32 weighted cost termsdrto.dynamic_optimization: horizon=kept, tracking_weight=(one stage cost declared), disturbances=w_ca=0, w_cb=0, w_tr=0, w_tk=0, sensitivity=4 initial-condition Params declared
Solve#
From the do-mpc example’s initial state, written into the _hat
Params at assembly, toward the paper’s setpoints c_A = 0.7,
c_B = 0.6 mol/L.
res = SolverFactory("pounce").solve(m, tee=True)
print(res.termination_condition.name)
print("objective", pyo.value(m.drto_objective))
********************************************************************************
#### ### / # /# #/ #### #####
# # # # /# #/ ## / # #
#### # #/ # / # #/# # ####
# # / # /# # /## # #
# ##/ #/# #/ # #### #####
********************************************************************************
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...: 909
Number of nonzeros in inequality constraint Jacobian.: 40
Number of nonzeros in Lagrangian Hessian.............: 198
Total number of variables............................: 201
variables with only lower bounds: 8
variables with lower and upper bounds: 120
variables with only upper bounds: 0
Total number of equality constraints.................: 165
Total number of inequality constraints...............: 8
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 8
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 9.3523418e+02 4.14e+03 1.01e+02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 8.7398172e+02 4.00e+03 9.85e+01 -1.0 4.26e+03 - 1.49e-02 3.21e-02h 1
2 1.1778223e+04 3.67e+02 1.10e+02 -1.0 4.00e+03 - 4.90e-02 1.00e+00h 1
3 1.0886344e+04 3.12e+02 7.87e+01 -1.0 3.73e+02 - 5.53e-01 1.48e-01f 1
4 9.9405280e+03 2.01e+02 4.89e+01 -1.0 3.17e+02 - 6.08e-01 3.56e-01f 1
5 8.8452141e+03 4.13e+00 1.78e+01 -1.0 2.15e+02 - 9.43e-01 1.00e+00f 1
6 8.8909464e+03 3.23e-01 1.98e-01 -1.0 2.46e+01 - 9.98e-01 1.00e+00h 1
7 8.8859825e+03 5.43e-04 3.10e-02 -1.7 1.30e+00 - 1.00e+00 1.00e+00h 1
8 8.8846366e+03 9.49e-05 9.65e-03 -2.5 1.91e-01 - 1.00e+00 1.00e+00f 1
9 8.8844070e+03 1.38e-05 1.02e-03 -3.8 3.11e-02 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 8.8843950e+03 2.26e-07 1.33e-05 -5.7 3.90e-03 - 1.00e+00 1.00e+00h 1
11 8.8843949e+03 5.88e-11 3.31e-09 -8.6 6.29e-05 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 11
(scaled) (unscaled)
Objective...............: 8.8843948948013679e+02 8.8843948948013676e+03
Dual infeasibility......: 3.3109783520705633e-09 3.3109783520705633e-08
Constraint violation....: 2.0176488827219406e-11 5.8804516811505891e-11
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 5.5600286731909703e-09 5.5600286731909702e-08
Overall NLP error.......: 5.5600286731909703e-09 5.5600286731909702e-08
Number of objective function evaluations = 12
Number of objective gradient evaluations = 12
Number of equality constraint evaluations = 12
Number of inequality constraint evaluations = 12
Number of equality constraint Jacobian evaluations = 12
Number of inequality constraint Jacobian evaluations = 12
Number of Lagrangian Hessian evaluations = 12
Total seconds in POUNCE = 0.015
EXIT: Optimal Solution Found.
POUNCE 0.11.0: Optimal Solution Found.
convergenceCriteriaSatisfied
objective 8884.394894801368
One training point#
The label is the first control action, and the objective is the
cost-to-go from the _hat state; the gradient half is the action’s
derivative with respect to the initial state, read from the solve’s
converged factorization through
drto.advanced_step_controller(m, gradient=True).
t0 = m.t.first()
grads = drto.advanced_step_controller(m, gradient=True)
for u in (m.F, m.Q_dot):
print(f"{u.name}[0] = {pyo.value(u[t0]):.6g}")
for p, g in grads[u].items():
print(f" d{u.name}[0]/d{p.name} = {g[u[t0], p]:.6g}")
F[0] = 33.9561
dF[0]/dca_hat = -0.772549
dF[0]/dcb_hat = -1.80919
dF[0]/dtr_hat = -1.12145
dF[0]/dtk_hat = -0.200229
Q_dot[0] = -715.505
dQ_dot[0]/dca_hat = -0.00161395
dQ_dot[0]/dcb_hat = -0.000469408
dQ_dot[0]/dtr_hat = -0.000949226
dQ_dot[0]/dtk_hat = -0.000751104
axes = drto.plot_states(m)
axes = drto.plot_controls(m)
axes = drto.plot_stage_cost(m)
