Pyomo Homework 4

How this assignment is organized¶

This assignment connects four ideas needed later in the course:

formulate, index-reduce, discretize, and initialize a DAE model;
estimate parameters from noisy dynamic data;
use parameter covariance to design the next experiment; and
formulate two-stage stochastic and risk-averse models with scenario Blocks.

Problems 2 and 3 use one deliberately small biochemical oxygen demand model so that the estimation--design loop remains visible.

#	Problem	Format
1	Pyomo.DAE: reaction kinetics	mixed --- index analysis, model reformulation, and the degree-of-freedom analysis are pencil and paper; everything else is notebook/code
2	Parameter estimation with `parmest`	notebook/code
3	Design the next measurement with Pyomo.DoE	notebook/code
4	Two-stage newsvendor under uncertainty	notebook/code

Course policies and submitting this assignment¶

Read the Artificial Intelligence Policy and Collaboration Policy and Honor Code before starting.

Submit two files to Canvas:

One scanned PDF containing all handwritten work from Problem 1: the index analysis, the model reformulation (including the consistent initial condition $z_C(0)$ ), and the degree-of-freedom analysis. Put the work in order and label each part. A clear phone photo assembled into a single PDF is fine.
A copy of this notebook, with your code cells run.

The three AI categories are:

No AI.
AI permitted after independent work --- spend the stated time on your own first, then AI and genuine collaboration (including coding together) are permitted. Everyone must still contribute intellectually, understand the work, and verify it.
AI required --- independent work first, then using an AI tool is part of the problem.

Each problem below states its own category and independent-work window --- see the note at the start of each problem rather than a single rule for the whole assignment.

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

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; it is part of the handwritten PDF you submit to Canvas.

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

Tip: do this on paper, together with the index analysis above --- both go in the same handwritten PDF.

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; it is part of the handwritten PDF you submit to Canvas.

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.

Problem 2. Parameter estimation with `parmest`¶

AI category: AI permitted after independent work. Spend about 20 minutes planning the model and API calls before using AI.

The Rooney--Biegler biochemical oxygen demand model is

y(t;\theta)=\theta_1\left(1-\exp(-\theta_2t)\right).

(3)

Use six observations with a measurement standard deviation of 1.0:

$t$ (days)	1	2	3	4	5	7
$y$	8.3	10.3	19.0	16.0	15.6	19.8

The two parameters are the asymptotic demand $\theta_1$ and rate constant $\theta_2$ . We use six observations because covariance estimation requires more experiments than unknown parameters.

2-A. Fit the parameters and estimate covariance¶

Use the current Pyomo 6.10 Experiment/Estimator workflow:

create one labeled experiment per observation;
use obj_function="SSE_weighted";
call theta_est(); and
call cov_est().

Report the estimates, covariance matrix, parameter standard deviations, and correlation coefficient. Explain in one or two sentences why two observations could fit two parameters but cannot provide the covariance requested here.

# Add your solution here

2-B. Check the fit¶

Plot the observations and fitted curve on the same axes. Then answer:

Which time region is most informative about the rate constant?
Which time region is most informative about the asymptote?
Does the fitted curve make the scatter in the measurements disappear? Why should it not?

# Add your solution here

Problem 2 AI use report¶

In concise bullets, report your independent attempt, any AI or collaborative help used afterward, and how you verified that help. If you used no AI, say so. Do not submit prompts or a transcript.

Problem 3. Design the next measurement with Pyomo.DoE¶

AI category: AI permitted after independent work. Spend about 20 minutes deriving the sensitivities and predicting a useful measurement time before using AI.

For the BOD model,

\frac{\partial y}{\partial\theta_1}=1-\exp(-\theta_2t), \qquad \frac{\partial y}{\partial\theta_2}=\theta_1t\exp(-\theta_2t).

(4)

Use the covariance from Problem 2 as the prior covariance. Choose one additional measurement time in $0.5\le t\le 5$ days, again with measurement standard deviation 1.0.

3-A. Predict before optimizing¶

Plot both sensitivities over the design interval. On paper, predict where a new measurement should be placed if the goal is to maximize the determinant of the updated Fisher information matrix. Explain the tradeoff in two or three bullets.

# Add your solution here

3-B. Optimize the design¶

Use DesignOfExperiments with:

the fitted parameters as nominal values;
$\mathbf{F}_{\mathrm{prior}}=\mathbf{V}_{\theta}^{-1}$ ;
one experiment whose hour is a design variable on $[0.5,5]$ ; and
ObjectiveLib.determinant.

Report the selected time and the prior and updated log determinants. Confirm the result by evaluating

\mathbf{F}_{\mathrm{new}}(t)=\mathbf{V}_{\theta}^{-1}+\mathbf{q}(t)\mathbf{q}(t)^\mathsf{T}

(5)

on a dense grid, where $\mathbf{q}(t)$ is the sensitivity vector above.

# Add your solution here

3-C. Interpret the design¶

Answer in concise bullets:

Why does the optimal time not occur at the latest possible measurement?
How did the prior covariance affect the answer?
What practical factors are absent from this design objective?
What would change if measurement variance depended on time?

Problem 3 AI use report¶

In concise bullets, report your independent attempt, any AI or collaborative help used afterward, and how you verified that help. If you used no AI, say so. Do not submit prompts or a transcript.

Problem 4. Two-stage newsvendor under uncertainty¶

AI category: AI permitted after independent work. Spend about 25 minutes deriving the deterministic-equivalent model before using AI.

A vendor orders newspapers before demand is known:

order cost: $c=\$10$ per paper;
selling price: $q=\$25$ per paper;
salvage value: $r=\$5$ per unsold paper; and
demand: $D\sim\mathrm{Uniform}[50,150]$ .

The order quantity $x$ is the first-stage decision. Sales and leftovers are scenario-specific recourse decisions.

4-A. Formulate the deterministic equivalent¶

Write the optimization model on paper. Then explain why:

$x$ belongs on the root model;
each scenario’s sales and leftovers belong in a scenario Block; and
no explicit nonanticipativity constraints are needed with this structure.

Derive the risk-neutral continuous solution from the critical fractile before coding.

4-B. Solve a sample-average approximation¶

Draw 200 demand scenarios from the uniform distribution using NumPy seed 60499. Build and solve the deterministic equivalent with root-level $x$ and indexed scenario Blocks. Report the SAA order quantity and compare it with the analytical value.

# Add your solution here

4-C. Quantify the value of stochastic information¶

For the continuous uniform distribution, calculate:

the expected-value solution and its expected result (EEV);
the recourse problem value (RP);
the wait-and-see value (WS);
the value of the stochastic solution (VSS); and
the expected value of perfect information (EVPI).

Use the profit convention, so $\mathrm{VSS}=\mathrm{RP}-\mathrm{EEV}$ and $\mathrm{EVPI}=\mathrm{WS}-\mathrm{RP}$ .

# Add your solution here

4-D. Add CVaR¶

Extend the block model to minimize a 50--50 blend of expected loss and 90% CVaR of loss. Report the risk-averse order quantity and compare its empirical profit distribution with the risk-neutral solution.

In concise bullets, explain why a model can have a lower expected profit but be preferable to a risk-averse decision maker.

# Add your solution here

Problem 4 AI use report¶

In concise bullets, report your independent attempt, any AI or collaborative help used afterward, and how you verified that help. If you used no AI, say so. Do not submit prompts or a transcript.

How this assignment is organized¶

Course policies and submitting this assignment¶

Problem 1. Pyomo.DAE: reaction kinetics¶

Index analysis¶

Model reformulation¶

Implement index 1 model in Pyomo¶

Create model and set initial conditions¶

Simulate, discretize, and initialize collocation model¶

Plot results¶

Simulate and solve Pyomo model with initialization¶

Simulate and solve Pyomo model without initialization¶

Discussion: Does initialization matter?¶

Degree of Freedom Analysis¶

Problem 1 AI use report¶

Problem 2. Parameter estimation with parmest¶

2-A. Fit the parameters and estimate covariance¶

2-B. Check the fit¶

Problem 2 AI use report¶

Problem 3. Design the next measurement with Pyomo.DoE¶

3-A. Predict before optimizing¶

3-B. Optimize the design¶

3-C. Interpret the design¶

Problem 3 AI use report¶

Problem 4. Two-stage newsvendor under uncertainty¶

4-A. Formulate the deterministic equivalent¶

4-B. Solve a sample-average approximation¶

4-C. Quantify the value of stochastic information¶

4-D. Add CVaR¶

Problem 4 AI use report¶

Problem 2. Parameter estimation with `parmest`¶