Ideal NMPC on the IDAES CSTR

The closed loop the declarations describe, in one call: at each step drto.ideal_nmpc solves the horizon at the measured state, implements the first move on the process, and simulates the process one sample forward. Ideal means the solve is treated as instantaneous, the baseline the delayed and advanced-step loops are judged against. The call builds both sides itself: the input model becomes the controller, a clone becomes the plant, and the cold start, the scaling, and the warm start are the earlier features run in their documented order.

The declared model and its setpoint

The saponification CSTR of models/idaes_cstr.py, as IDAES wrote it. The setpoint comes from the steady reduction at zero duty and nominal feed; the start is hot (318 K) with the conversion knocked back. The model is built once, additive noise terms declared (zero nominally, so the nominal run is untouched by them), and cloned before the first loop consumes the original: the clone is the noisy run’s copy. The setup wraps in a function so the noisy variant at the end can rebuild it. The infinite-horizon transform is part of defining the model; the loop takes it from there.

from functools import partial

import pyomo.environ as pyo

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


def build_cstr_at_setpoint(N=10, h=1, disturbance=False):
    """The declared model with its targets read off a steady simulation.

    A builder in feature 006's sense: it returns the declared,
    undiscretized model, and its first two parameters are N and h with
    defaults. The loop builds both sides from it, so the setpoint is
    written here rather than onto a model the caller holds.
    """
    m = build(N=N, h=h, disturbance=disturbance)
    ss = pyo.TransformationFactory("drto.steady_state_simulation").create_using(
        m, controls={m.fs.cstr.control_volume.heat.name: 0.0,
                     m.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])
    return m


build_noisy_cstr = partial(build_cstr_at_setpoint, disturbance=True)

The loop, one call

Ten samples of closed loop. Inside the call: the process clone in simulation mode at the declared targets, cut to one sample, the dynamic optimization assembly, an exponential cold start (time constant 0.3 s) for the first solves on both sides, and every re-solve warm-started from the previous solution shifted one sample, the warm start options applied at the call site (a warm_start mapping retunes them). scale hands drto.scale the flowsheet’s magnitudes per unit, written once at entry, so both sides carry the factors through their cold starts and every solve; the history reads back in physical units.

import time

tic = time.time()
history = drto.ideal_nmpc(
    build_noisy_cstr,
    steps=10,
    dynamic_optimization={"infinite_horizon": True},
    initialize={"profile": "exponential", "time_constant": 0.3},
    scale={"J": 1e7, "W": 1e6},
)
print(history)
print(f"{time.time() - tic:.0f} s for 10 steps, a controller and a plant solve each")

drto nmpc history: 10 steps, states material_holdup_Liq_NaOH, material_holdup_Liq_EthylAcetate, material_holdup_Liq_SodiumAcetate, material_holdup_Liq_Ethanol, energy_holdup[Liq], moves heat, _flow_vol_inlet_ref
6 s for 10 steps, a controller and a plant solve each

The actual trajectories

With a perfect model and no noise the ideal loop rides its own plan: the holdups run from the knocked-back start onto their setpoints, the energy comes down from the hot start, and the moves settle onto the declared targets, each held over its sample.

axes = drto.plot_states(history)

axes = drto.plot_controls(history)

The same loop under process noise

Every state’s balance carries its own additive noise term: mol/s in each species holdup equation, W in the energy equation, zero nominally, written into the IDAES-generated balances through the control volume’s custom terms. The loop draws each step’s realization normally, zero mean with the given standard deviation per term, reproducible under seed. The plant feels the draws, the controller plans without them, and the holdups settle into a band around their setpoints instead of onto them, the moves working against each draw. The loop runs on the clone taken above, so nothing rebuilds.

noisy = drto.ideal_nmpc(
    build_noisy_cstr,
    steps=10,
    dynamic_optimization={"infinite_horizon": True},
    disturbances={"w_NaOH": 0.5, "w_EthylAcetate": 0.5,
                  "w_SodiumAcetate": 0.5, "w_Ethanol": 0.5,
                  "w_energy": 1e6},
    seed=7,
    initialize={"profile": "exponential", "time_constant": 0.3},
    scale={"J": 1e7, "W": 1e6},
)
print(noisy)
print("w_NaOH draws:", " ".join(f"{w:+.2f}" for w in noisy.realizations["w_NaOH"]))

drto nmpc history: 10 steps, states material_holdup_Liq_NaOH, material_holdup_Liq_EthylAcetate, material_holdup_Liq_SodiumAcetate, material_holdup_Liq_Ethanol, energy_holdup[Liq], moves heat, _flow_vol_inlet_ref
w_NaOH draws: -0.13 -0.11 +0.20 +0.25 +0.15 +0.20 -0.31 +0.12 -0.40 +0.65

axes = drto.plot_states(noisy)

axes = drto.plot_controls(noisy)