Dynamic Systems
This notebook is part of the SECQUOIA Research Group Pyomo Tutorial. It adapts the dynamic-system exercises from the public Pyomo workshop tutorials. The discretized optimization examples use IPOPT; see Setup and Solvers for solver notes.
Learning objectives¶
By the end of this chapter, you should be able to:
represent differential variables with
ContinuousSetandDerivativeVar;connect differential equations, initial conditions, and algebraic constraints in Pyomo;
discretize dynamic models with finite differences and collocation;
use simulation as a check before optimization when SciPy is available;
formulate a small dynamic parameter-estimation problem.
import os
import shutil
import subprocess
import sys
import pyomo.environ as pyo
from pyomo.dae import ContinuousSet, DerivativeVar, Simulator
if "google.colab" in sys.modules and shutil.which("ipopt") is None:
subprocess.run([sys.executable, "-m", "pip", "install", "-q", "idaes-pse"], check=True)
subprocess.run(["idaes", "get-extensions"], check=True)
os.environ["PATH"] = "/root/.idaes/bin:" + os.environ["PATH"]
def solve_with_ipopt(model):
solver = pyo.SolverFactory("ipopt")
if not solver.available(False):
print("IPOPT is not available; the model was built but not solved.")
return None
results = solver.solve(model)
print(f"termination: {results.solver.termination_condition}")
return results
def analytical_z(t):
return (4 * t - 3) / (4 * t + 1)
1. Differential equations by hand¶
The first exercise uses the scalar initial-value problem dz/dt = z**2 - 2*z + 1, with z(0) = -3. The analytical solution is z(t) = (4*t - 3)/(4*t + 1). We will first discretize the derivative manually with a backward finite difference.
Exercise 1.1. Manual finite difference¶
Starting from small_findiff_soln.py, identify the algebraic variables created at each time point and compare the numerical values against the analytical solution.
numpoints = 5
findiff = pyo.ConcreteModel()
findiff.points = pyo.RangeSet(0, numpoints - 1)
findiff.h = pyo.Param(initialize=1.0 / (numpoints - 1))
findiff.z = pyo.Var(findiff.points)
findiff.dzdt = pyo.Var(findiff.points)
findiff.obj = pyo.Objective(expr=1)
@findiff.Constraint(findiff.points)
def zdot(model, i):
return model.dzdt[i] == model.z[i] ** 2 - 2 * model.z[i] + 1
@findiff.Constraint(findiff.points)
def backward_difference(model, i):
if i == 0:
return pyo.Constraint.Skip
return model.dzdt[i] == (model.z[i] - model.z[i - 1]) / model.h
findiff.initial_condition = pyo.Constraint(expr=findiff.z[0] == -3)
result = solve_with_ipopt(findiff)
if result is not None:
print(" i t z_model z_exact")
for i in findiff.points:
t = pyo.value(findiff.h) * i
print(f"{i:2d} {t:5.2f} {pyo.value(findiff.z[i]):8.4f} {analytical_z(t):8.4f}")
termination: optimal
i t z_model z_exact
0 0.00 -3.0000 -3.0000
1 0.25 -1.4721 -1.0000
2 0.50 -0.7267 -0.3333
3 0.75 -0.3026 0.0000
4 1.00 -0.0348 0.2000
2. Pyomo DAE discretization¶
Manual discretizations are useful for learning, but pyomo.dae can generate the finite-difference or collocation equations for you. The continuous model stays close to the differential equation, and the transformation creates the algebraic formulation.
Exercise 2.1. Collocation with pyomo.dae¶
Starting from small_dae_soln.py, change nfe and ncp and observe how the number of discretization points changes.
dae_model = pyo.ConcreteModel()
dae_model.t = ContinuousSet(bounds=(0, 1))
dae_model.z = pyo.Var(dae_model.t)
dae_model.dzdt = DerivativeVar(dae_model.z)
dae_model.obj = pyo.Objective(expr=1)
@dae_model.Constraint(dae_model.t)
def dae_zdot(model, t):
return model.dzdt[t] == model.z[t] ** 2 - 2 * model.z[t] + 1
dae_model.initial_condition = pyo.Constraint(expr=dae_model.z[0] == -3)
pyo.TransformationFactory("dae.collocation").apply_to(
dae_model, nfe=4, ncp=3, scheme="LAGRANGE-RADAU"
)
result = solve_with_ipopt(dae_model)
if result is not None:
time_points = list(dae_model.t)
print(f"discretization points: {len(time_points)}")
for t in time_points[:6]:
print(f"t={float(t):6.4f}, z={pyo.value(dae_model.z[t]):8.4f}, exact={analytical_z(float(t)):8.4f}")
print("...")
t = time_points[-1]
print(f"t={float(t):6.4f}, z={pyo.value(dae_model.z[t]):8.4f}, exact={analytical_z(float(t)):8.4f}")
termination: optimal
discretization points: 13
t=0.0000, z= -3.0000, exact= -3.0000
t=0.0388, z= -2.4773, exact= -2.4630
t=0.1612, z= -1.4212, exact= -1.4317
t=0.2500, z= -1.0000, exact= -1.0000
t=0.2888, z= -0.8571, exact= -0.8561
t=0.4112, z= -0.5115, exact= -0.5123
...
t=1.0000, z= 0.2000, exact= 0.2000
3. Simulation before optimization¶
When a model has fixed inputs and no optimization decisions, simulation is often the fastest way to check signs, units, and initial conditions. The workshop simulation exercise models a consecutive reaction from A to B.
Exercise 3.1. Simulate a reaction model¶
The Simulator interface uses SciPy for this exercise. If SciPy is missing in a local environment, the model-building pattern is still the same and the simulation cell reports the missing dependency.
try:
import scipy # noqa: F401
except Exception as err:
print(f"SciPy is not available, so the simulation step is skipped: {err}")
else:
sim_model = pyo.ConcreteModel()
sim_model.t = ContinuousSet(bounds=(0, 1))
sim_model.a = pyo.Var(sim_model.t)
sim_model.b = pyo.Var(sim_model.t)
sim_model.k1 = pyo.Param(initialize=5)
sim_model.k2 = pyo.Param(initialize=1)
sim_model.dadt = DerivativeVar(sim_model.a)
sim_model.dbdt = DerivativeVar(sim_model.b)
sim_model.a[0].fix(1)
sim_model.b[0].fix(0)
@sim_model.Constraint(sim_model.t)
def a_balance(model, t):
return model.dadt[t] == -model.k1 * model.a[t]
@sim_model.Constraint(sim_model.t)
def b_balance(model, t):
return model.dbdt[t] == model.k1 * model.a[t] - model.k2 * model.b[t]
simulator = Simulator(sim_model, package="scipy")
tsim, profiles = simulator.simulate(integrator="vode", numpoints=25)
print("final simulated concentrations")
for index, var in enumerate(simulator.get_variable_order()):
print(f" {var}: {profiles[-1, index]:.4f}")
final simulated concentrations
a[{t}]: 0.0067
b[{t}]: 0.4514
4. Dynamic parameter estimation¶
Parameter estimation embeds the dynamic model in an optimization problem. The objective penalizes mismatch between model predictions and measurements, while the differential equations constrain the trajectory.
Exercise 4.1. Estimate reaction-rate constants¶
Starting from start_param_est2.py, complete the variables, differential equations, objective, discretization, and solve. A full solution is available in param_est2_soln.py.
a_conc = {0.1: 0.606, 0.2: 0.368, 0.3: 0.223, 0.4: 0.135, 0.5: 0.082,
0.6: 0.050, 0.7: 0.030, 0.8: 0.018, 0.9: 0.011, 1.0: 0.007}
b_conc = {0.1: 0.373, 0.2: 0.564, 0.3: 0.647, 0.4: 0.669, 0.5: 0.656,
0.6: 0.624, 0.7: 0.583, 0.8: 0.539, 0.9: 0.494, 1.0: 0.451}
est = pyo.ConcreteModel()
est.meas_time = pyo.Set(initialize=sorted(a_conc), ordered=True)
est.ameas = pyo.Param(est.meas_time, initialize=a_conc)
est.bmeas = pyo.Param(est.meas_time, initialize=b_conc)
est.time = ContinuousSet(initialize=est.meas_time, bounds=(0, 1))
est.a = pyo.Var(est.time, bounds=(0, 1))
est.b = pyo.Var(est.time, bounds=(0, 1))
est.dadt = DerivativeVar(est.a)
est.dbdt = DerivativeVar(est.b)
est.k1 = pyo.Var(bounds=(0, 10), initialize=4)
est.k2 = pyo.Var(bounds=(0, 10), initialize=1)
@est.Constraint(est.time)
def a_diffeq(model, t):
return model.dadt[t] == -model.k1 * model.a[t]
@est.Constraint(est.time)
def b_diffeq(model, t):
return model.dbdt[t] == model.k1 * model.a[t] - model.k2 * model.b[t]
est.ainit = pyo.Constraint(expr=est.a[0] == 1)
est.binit = pyo.Constraint(expr=est.b[0] == 0)
est.obj = pyo.Objective(
expr=sum((est.a[t] - est.ameas[t]) ** 2 + (est.b[t] - est.bmeas[t]) ** 2 for t in est.meas_time)
)
pyo.TransformationFactory("dae.collocation").apply_to(est, nfe=10, ncp=3, scheme="LAGRANGE-RADAU")
result = solve_with_ipopt(est)
if result is not None:
print(f"k1 = {pyo.value(est.k1):.4f}")
print(f"k2 = {pyo.value(est.k2):.4f}")
print(f"sum of squared errors = {pyo.value(est.obj):.6f}")
termination: optimal
k1 = 5.0035
k2 = 1.0000
sum of squared errors = 0.000001
Exercises to continue¶
Replace the collocation discretization in Exercise 4.1 with finite differences and compare the estimated parameters.
Add bounds to
k1andk2that reflect prior physical knowledge, then solve again.Use the solution trajectory to report residuals at each measured time point.
References¶
Dynamic optimization and DAE discretization patterns are part of Pyomo’s modeling extensions Bynum et al., 2021.
- Bynum, M. L., Hackebeil, G. A., Hart, W. E., Laird, C. D., Nicholson, B. L., Siirola, J. D., Watson, J.-P., & Woodruff, D. L. (2021). Pyomo–optimization modeling in python (Third, Vol. 67). Springer Science & Business Media.