Pyomo Homework 4 - Optimization for Decision Science

What this assignment will cover¶

Pyomo Homework 4 is due Monday, September 28, 2026 and is the last homework before Midterm 1 (Wednesday, September 30). It draws on the lectures immediately preceding it:

Lecture	Date	Topic
8	Wednesday, September 16	Parameter Estimation and Optimal Experimental Design
9	Monday, September 21	Stochastic Programming: Introduction (video lecture)
10	Wednesday, September 23	Stochastic Programming: Information and Risk (video lecture)

The course website pages that work these models in full are the best preparation:

Course policies and AI category¶

Read the Artificial Intelligence Policy and Collaboration Policy and Honor Code before starting. Assignment-specific directions control. The categories are No AI, AI permitted after independent work, and AI required.

Unless a problem says otherwise, its category is AI permitted after independent work. Spend about 30 minutes on each top-level problem without AI, solution pages, or another person’s help, stopping early if complete. You may consult lecture notes, textbooks, and nonsolution pages of the course website; bias toward those course sources. Afterward, AI and genuine collaboration, including coding together, are permitted. Everyone must contribute intellectually, understand the work, and verify it.

At the end of each top-level problem, add a concise AI and independent-work report: approximately how long the independent attempt took, how far you got, where you became stuck, any AI or collaborative help used afterward, and how you verified it. If you used no AI, say so. Do not submit prompts or transcripts. Time estimates help the instructor improve the assignment and are not a speed test.

# Import the libraries you need here for the assignment

import sys

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper

    helper.easy_install()
else:
    sys.path.insert(0, "../")
    import helper
helper.set_plotting_style()

# Add your solution here

Problems¶

Problem 1 is below and is final. The remaining problems --- parameter estimation, experimental design, and a stochastic program --- will be added here before the assignment opens.

Problem 1. Pyomo.DAE: reaction kinetics¶

AI category: AI permitted after independent work. Work for about 30 minutes without AI, solution pages, or help from another person, or stop early if you complete the problem.

Moved here from Pyomo Homework 3 for Fall 2026: this system is the one the parameter-estimation and experimental-design material builds on, so it now sits alongside it.

Consider the chemical reaction

A \Leftrightarrow B \Leftrightarrow C

(1)

which is modeled 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*}

(2)

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.

Index analysis¶

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.

Model reformulation¶

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)$ .

Implement index 1 model in Pyomo¶

We will be building a library of functions.

Create model and set initial conditions¶

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.

    m = pyo.ConcreteModel()

    m.t = ContinuousSet(bounds=(0.0, 1))

    # Add your solution here

    return m

Simulate, discretize, and initialize collocation model¶

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 element
    3. If initialize is true, call 'sim.initialize_model()'. This will use the
        Simulator solution to initialize the discretized Pyomo model. Really cool!
    """

    # Add your solution here

    return sim, tsim, profiles

Plot results¶

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
    """

    # Add 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

    m.initcon = pyo.ConstraintList(rule=_init)

    # Solve collocation formulation (no objective, we are just simulating)
    solver = pyo.SolverFactory("ipopt")
    results = solver.solve(m, tee=True)
    assert pyo.check_optimal_termination(results), (
        f"Solve failed: status={results.solver.status}, "
        f"termination={results.solver.termination_condition}"
    )

Simulate and solve Pyomo model with initialization¶

# 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)

Simulate and solve Pyomo model without initialization¶

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.)

# Add your solution here

Discussion: Does initialization matter?¶

Write 1 to 3 sentences for each of the following questions:

What happens if we disable initialization? Does the number of iterations Ipopt needs change?

Why is this specific problem robust to poor initialization? What about this specific DAE system makes it easy to solve?

Degree of Freedom Analysis¶

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 similar number of finite elements. Using the collocation equations from class, show that the discretized Pyomo model has the correct number of variables and algebraic 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.

Note: You might get a strange answer. That is okay.

Problem 1 AI use report¶

At the end of Problem 1, report for the problem as a whole:

approximately how long you spent on the single 30-minute independent attempt for Problem 1;
how far you got during that independent attempt;
where you got stuck, if anywhere;
any AI or collaborative help you used afterward; and
how you verified that help.

If you used no AI, say so explicitly. A few concise bullets are sufficient; do not submit prompts or transcripts. The time estimate gives the instructor useful data for improving the assignment and is not a speed test. Apply the 30-minute guideline once to the entire problem, not separately to each subpart.