# Import the libraries you need here for the assignment
# YOUR SOLUTION HERE
Consider the chemical reaction $$A \Leftrightarrow B \Leftrightarrow C$$
which is modeling with the following differential algebraic equations:
$$\begin{align*} \frac{dz_A}{dt} &= -p_1 z_A(t) + p_2 z_B(t),\quad z_A(0)=1 \\ \frac{dz_B}{dt} &= p_1 z_A(t) - (p_2 + p_3) z_B(t) + p_4 z_C(t), \quad z_B(0)=0 \\ 1 &= z_A(t) + z_B(t) + z_C(t) \end{align*}$$where $p_1=4$, $p_2=2$, $p_3=40$, and $p_4=20$ are parameters with the appropriate units. $z_A(t)$, $z_B(t)$, and $z_C(t)$ are time varying concentrations of species $A$, $B$, and $C$ respectively.
Determine the index of the above differential algebraic equation (DAE) system above.
Tip: do this on paper. On Gradescope, there will be a separate assignment for you to turn in your handwritten work.
Apply the index reduction algorithm from class as needed. Ultimately identify two versions of the model: one that is index 1 and another that is index 0. Find a consistent initial condition $z_C(0)$.
We will be building a library of functions.
def create_model():
''' Create index 1 model and set initial conditions
Return:
m: Pyomo model
'''
# Tip: Set time to go from 0 to 1 when creating the model.
# YOUR SOLUTION HERE
return m
def simulate_discretize_model(m,NFE,initialize):
''' Simulation, discretize, and initialize the Pyomo model
Arguments:
m: Pyomo model
NFE: number of finite elements to consider (integer)
initialize: if True, initialize the discretized model with the
integrator solution (boolean)
Returns:
sim: Simulator object from Pyomo.DAE
tsim: Timesteps returned from simulator
profiles: Results returned from simulator
Overall Steps:
1. Create Pyomo.DAE simulator and integrate with casadi/idas
2. Transform model using 'dae.collocation' strategy. Use 3 collocation points
per finite elemebt
3. If initialize is true, call 'sim.initialize_model()'. This will use the
Simulator solution to initialize the discretized Pyomo model. Really cool!
'''
# YOUR SOLUTION HERE
return sim, tsim, profiles
def plot_result(m, sim, tsim, profiles, include_model_values):
""" Plot the results from the simulator (and optionally Pyomo model)
Arguments:
m: Pyomo model
sim: Pyomo.DAE simulator
tsim: timesteps from simulator
profiles: results from simulation
include_model_values: if True, also plot the values from the Pyomo model m
Returns:
nothing
Actions/Steps/Tips:
1. Plot the results stored in tsim and profiles as solid lines. Recycle code from class.
2. If 'include_model_values' is true, plot za, zb, and zc values stored in Pyomo model 'm'.
Use a solid symbol.
3. Add a legend and axes labels
"""
# YOUR SOLUTION HERE
# Tip: Do not forget to include `plt.show()` (assuming you imported matplotlib.pyplot as plt)
def solve_model(m):
""" Solve discretized model with Ipopt
Arguments:
m: Pyomo model
Returns:
nothing
"""
# Specify initial conditions
def _init(m):
yield m.za[0] == 1
yield m.zb[0] == 0
model.initcon = ConstraintList(rule=_init)
# Solve collocation formulation (no objective, we are just simulating)
solver = SolverFactory('ipopt')
solver.solve(model,tee=True)
# Create Pyomo model
model = create_model()
# Initialize discretized model with simulation result?
init = True
# Number of finite elements
NFE = 6
# Simulate model
sim, tsim, profiles = simulate_discretize_model(model,NFE,init)
# Plot simulation results
plot_result(model, sim, tsim, profiles,True)
# Solve collocation formulation with Ipopt
solve_model(model)
# Plot results
plot_result(model, sim, tsim, profiles, True)
Repeat the analysis from above, but do NOT initialize the discretized model with the simulation result. (Tip: you only need to change two small things.)
# YOUR SOLUTION HERE
Write 1 to 3 sentences for each of the following questions:
What happens if we disable initialization? Does the number of interactions Ipopt need change?
Why is this specific problem robust to poor initialization? What about this specific DAE system makes it easy to solve?
Please answer on paper and turn in via Gradescope.
Rerun the code above for two different numbers of finite elements. Record the total number of variables and equality constraints.
How many degrees of freedom are in the simulation problem? (1 sentence is fine.)
Choose $N=3$ or a simular number of finite elements. Using the collocation equations from class, show that the discretized Pyomo model has the correct number of variables and algebriac equations. This will be a little tedious, but good to practice once on a simple model. You might need to do this when debugging a DAE model for research.