Advanced step on the IDAES CSTR#

Advanced-step NMPC runs the expensive solve between samples, at a predicted state, and corrects it the instant the measurement arrives: a backsolve on the solve’s kept factorization, not a re-solve. On a real flowsheet the solve is the cost that matters, so this is the example where the idea earns its keep. The model is the canonical IDAES saponification CSTR.

The advanced-step NMPC controller is Zavala & Biegler (2009), doi:10.1016/j.automatica.2008.06.011: solve at a prediction between samples, correct at the measurement through the NLP’s sensitivity system.

import time

import pyomo.environ as pyo
from pyomo.contrib.solver.common.factory import SolverFactory

import drto
from models.idaes_cstr import DC_START, F_IN, VOLUME, build

# the horizon, the terminal segment, and the assembly in one call. The
# builder holds the feed at the declared sample points and this fixes the
# members collocation adds to it
m = drto.dynamic_optimization(build, N=10, infinite_horizon=True)

# the setpoint comes off a steady simulation, which reduces a declared
# model, so it runs on a second build rather than on the assembled one
declared = build()
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})
drto.initialize_steady_state(ss)
drto.scaled_solve(ss)

cvs = ss.fs.cstr.control_volume
for j, ssp in (("NaOH", m.ss_naoh), ("EthylAcetate", m.ss_ea),
               ("SodiumAcetate", m.ss_sa), ("Ethanol", m.ss_etoh)):
    ssp.set_value(pyo.value(cvs.material_holdup["Liq", j]))
for j, sgn in (("NaOH", 1), ("EthylAcetate", 1),
               ("SodiumAcetate", -1), ("Ethanol", -1)):
    m.mat0[j] = pyo.value(cvs.material_holdup["Liq", j]) + sgn * DC_START * VOLUME
for k in m.eng_ss:
    m.eng_ss[k] = pyo.value(cvs.energy_holdup[k])
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")

setpoint  T = 304.039 K   NaOH = 24.286 mol/m3

The controller#

The terminal segment joins at definition time, the model is cold-started onto the targets, and the assembly declares the five initial-condition Params as pounce sensitivity parameters. The solves run under the solver’s own scaling, which converges on this model, and the correction inherits whatever scaling the prediction solve ran under, since it is a backsolve on that solve’s kept factorization.

drto.cold_start_dynamic(m, profile="exponential", time_constant=3.0)
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_volfs.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_molfs.props.cp_molfs.cstr._flow_vol_inlet_ref[0.0](fs.cstr.control_volume.properties_in[t].temperature - fs.props.temperature_ref) - fs.props.dens_molfs.props.cp_molfs.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

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
drto.parameterize: profiles=fs.cstr.control_volume.heat (piecewise_constant), fs.cstr._flow_vol_inlet_ref (piecewise_constant)
drto.build_objective: objective=sum of 45 weighted cost terms
drto.dynamic_optimization: horizon=kept, tracking_weight=(one stage cost declared), sensitivity=5 initial-condition Params declared

Solve at the prediction#

Directly with pounce: the factorization must live with the model itself, since the correction is a backsolve on it.

tic = time.perf_counter()
SolverFactory("pounce").solve(m)
t_solve = time.perf_counter() - tic

cv = m.fs.cstr.control_volume
fin = m.fs.cstr._flow_vol_inlet_ref
i_heat, i_fin = sorted(cv.heat)[0], sorted(fin)[0]
print(f"solve at the prediction: {t_solve:.2f} s")
print(f"first moves at the prediction: duty = {pyo.value(cv.heat[i_heat]):12.1f} W,"
      f"  flow = {pyo.value(fin[i_fin]):.5f} m3/s")

solve at the prediction: 0.58 s
first moves at the prediction: duty =   -3943647.3 W,  flow = 1.79732 m3/s

The measurement arrives#

The measured holdups differ from the prediction: one more mole of NaOH and half a megajoule of energy. Write them into the initial-condition Params and ask for the corrected solution; the model itself is untouched. pounce flags one caveat, summarized below: the tail’s pin slacks sit exactly at their zero bounds, so any linear step pushes half of each pair out and gets clamped. That flag is about the slacks, not the moves.

m.mat0["NaOH"] = pyo.value(m.mat0["NaOH"]) + 1.0
m.eng0["Liq"] = pyo.value(m.eng0["Liq"]) + 5.0e5

import warnings

tic = time.perf_counter()
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    est = drto.advanced_step_controller(m)
t_corr = time.perf_counter() - tic
for w in caught:
    assert "clamped" in str(w.message), str(w.message)
if caught:
    print(f"({len(caught)} pounce warning: the tail's pin slacks clamped at "
          f"their zero bounds; the moves are unaffected)")
else:
    print("(no pin-slack warnings this run)")
print(f"correction: {1e3 * t_corr:.1f} ms "
      f"({t_solve / t_corr:.0f}x faster than the solve)")
print(f"corrected first moves:         duty = {est[cv.heat[i_heat]]:12.1f} W,"
      f"  flow = {est[fin[i_fin]]:.5f} m3/s")

(1 pounce warning: the tail's pin slacks clamped at their zero bounds; the moves are unaffected)
correction: 6.8 ms (86x faster than the solve)
corrected first moves:         duty =   -3950674.6 W,  flow = 1.80295 m3/s

Sensitivities#

gradient=True returns each control’s sensitivity to each initial-condition Param: how the first duty and flow moves respond to the measured composition and energy.

g = drto.advanced_step_controller(m, gradient=True)
heat, flow = drto.info(m).components("control")
for hook in (m.mat0["NaOH"], m.eng0["Liq"]):
    Gh, Gf = g[heat][hook], g[flow][hook]
    print(f"d duty[0] / d {hook.name:<12} = {Gh[cv.heat[i_heat], hook]:+11.4g}"
          f"    d flow[0] / d {hook.name:<12} = {Gf[fin[i_fin], hook]:+11.4g}")

d duty[0] / d mat0[NaOH]   =      +55.75    d flow[0] / d mat0[NaOH]   =   +0.000221
d duty[0] / d eng0[Liq]    =    -0.01436    d flow[0] / d eng0[Liq]    =  +1.078e-08

Against the full re-solve#

The honest comparison is the one the closed loop makes: the same model re-solved in place at the measured initial conditions, warm from the prediction’s solution. A freshly initialized model can land on a different local branch of this nonlinear flowsheet; the loop never starts fresh.

tic = time.perf_counter()
SolverFactory("pounce").solve(m)
t_resolve = time.perf_counter() - tic

print(f"re-solve, warm from the prediction: {t_resolve:.2f} s")
print(f"re-solved first moves:         duty = {pyo.value(cv.heat[i_heat]):12.1f} W,"
      f"  flow = {pyo.value(fin[i_fin]):.5f} m3/s")

d_duty = abs(est[cv.heat[i_heat]] - pyo.value(cv.heat[i_heat])) / abs(
    pyo.value(cv.heat[i_heat]))
d_flow = abs(est[fin[i_fin]] - pyo.value(fin[i_fin])) / abs(pyo.value(fin[i_fin]))
print(f"first-move agreement: duty within {100 * d_duty:.2f}%, "
      f"flow within {100 * d_flow:.2f}%")
# the deviation grows along the horizon (first order, and the tail is
# weakly determined); the loop implements only the first move
reg = drto.info(m)
worst = 0.0
for comp in reg.components("control"):
    for vd in comp.values():
        if not vd.fixed:
            worst = max(worst, abs(est[vd] - pyo.value(vd)) / max(1.0, abs(pyo.value(vd))))
print(f"largest deviation over the whole horizon:  {100 * worst:.0f}% "
      f"(late moves, never implemented)")

re-solve, warm from the prediction: 0.40 s
re-solved first moves:         duty =   -3950616.2 W,  flow = 1.80293 m3/s
first-move agreement: duty within 0.00%, flow within 0.00%
largest deviation over the whole horizon:  7% (late moves, never implemented)

Advanced step on the IDAES CSTR

Contents