Pyomo Homework 1

Course policies and submitting this assignment¶

Read the Artificial Intelligence Policy and Collaboration Policy and Honor Code before starting. Then please see the assignment-specific directions below.

Submit two files to Canvas:

1. One scanned PDF containing all handwritten work for Problems 1, 2, and 4-A. Put the work in problem order and label each part. A clear phone photo assembled into a single PDF is fine.

2. A copy of this notebook, with your code cells run.

This assignment is graded on completion. Published answers, if provided, are so you can check your own work. Attempt each problem before you look.

How this assignment is organized¶

This homework has four problems. The first two are practice exam problems: pencil and paper, no solver, no code cell, written in the exact format of the in-person midterm, so the exam is not the first time you see a question like this. Problems 3 and 4 combine independent attempts with disclosed AI use.

Problem 1. Student diet¶

Pencil and paper, practice exam problem. Attempt every part before you open the linked answer page.

Background. You want to save money eating while remaining healthy. A healthy diet requires at least $P = 6$ units of protein, $C = 15$ units of carbohydrates, $F = 5$ units of fats and $V = 7$ units of vitamins per day. Due to compounding factors (e.g., blizzard during Lent), campus only has these options:

Instruction. Consider the following Python and Pyomo code, which should NOT contain any mistakes.

# Data. food_options is the table above, read from ../data/student_diet.csv
nutrients = food_options.columns.to_list()[0:4]     # ['P', 'C', 'F', 'V']
foods = food_options.index.to_list()                # ['takeaway', 'vegetables', 'bread']
food_info = food_options[nutrients].stack().to_dict()   # {('takeaway', 'P'): 3, ...}
price = food_options["price"].to_dict()                 # {'takeaway': 5, ...}
nutrient_requirements = {"P": 6, "C": 15, "F": 5, "V": 7}

m = pyo.ConcreteModel()

## Define sets
m.FOOD = pyo.Set(initialize=foods)
m.NUTRIENTS = pyo.Set(initialize=nutrients)

## Define parameters
m.needs = pyo.Param(m.NUTRIENTS, initialize=nutrient_requirements, units=u.nutrient)
m.food_info = pyo.Param(
    m.FOOD, m.NUTRIENTS, initialize=food_info, units=u.nutrient / u.serving
)
m.price = pyo.Param(m.FOOD, initialize=price, units=u.USD / u.serving)

## Define variables
m.food_eaten = pyo.Var(
    m.FOOD, initialize=1.0, domain=pyo.NonNegativeReals, units=u.serving
)


## Define constraints
# C1
@m.Constraint(m.NUTRIENTS)
def diet_min(b, n):
    return sum(b.food_info[f, n] * b.food_eaten[f] for f in b.FOOD) >= b.needs[n]


## Define objective
@m.Objective(sense=pyo.minimize)
def cost(b):
    return sum(b.food_eaten[f] * b.price[f] for f in b.FOOD)

Solving this model with Ipopt gives:

Units of takeaway eaten = 1.0
Units of vegetables eaten = 1.5
Units of bread eaten = 3.0
Total cost = 12.50 USD/day

1-A. Translate the Pyomo code into a mathematical model (symbols and equations) using set notation similar to our in-class examples¶

Next to each set, parameter, variable, and equation, write a few-word description.

Sets.

Parameters (Data). Hint: state the units.

Variables. Hint: state the units.

Objective. (write mathematical equations)

Constraints. Write mathematical equations. Write a brief description for C1. You are filling in the missing Python comments here.

1-B. Degree of freedom analysis¶

Count the continuous variables, the equality constraints, and the inequality constraints separately, and report the degrees of freedom. Then answer: why is it acceptable for this model to have more inequality constraints than variables?

1-C. Problem size¶

Let $N_f$ represent the number of food options and $N_n$ represent the number of nutrients. Using these symbols, determine the variables and constraints in your model. After each header below, give a number with a brief justification.

Number of continuous variables:

Number of integer/discrete variables:

Number of equality constraints:

Number of inequality constraints:

Number of variable bounds: As in Lecture 2’s degree of freedom analysis, report these separately from the inequality constraints above, split into lower-only, upper-only, and both.

Number of parameters (data values) you must supply:

1-D. Classify the problem¶

Is this an LP, QP, NLP, MILP or MINLP? Justify your classification by referring to the objective, the constraints and the variable domains. Is the problem convex? What does that buy you?

1-E. Now check your work in Pyomo, and refresh your Pyomo fundamentals¶

Only after you have attempted 1-A through 1-D on paper, open Continuous Optimization: Linear Programming --- Continuous Optimization: Linear Programming (class website). It is the same diet problem, worked in full, with the answers to 1-A through 1-D.

Compare your answers against the notebook. Then run the notebook in Google Colab. Experiment with the notebook by changing the model and resolve.

Write one thing that you learned from checking your answer or changing and resolving the model.

1-F. AI use report¶

At the end of Problem 1, report:

• approximately how long the independent attempt took,

• how far you got

• where you got stuck, if anywhere,

• any AI or collaborative help you used afterward

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

Problem 2. Portfolio optimization¶

Pencil and paper, practice exam problem. Attempt every part before you open the linked answer page.

Background. You have a fixed amount of money and $N$ assets to put it in. Each asset’s return fluctuates from day to day. You want to choose what fraction of your money goes into each asset so that the variance of the portfolio’s return is as small as possible, while still achieving at least a required expected return $\rho$. This is the Markowitz mean/variance model, and it is the problem that won Markowitz a Nobel Prize in 1990.

The data are five years of daily adjusted closing prices for $N = 5$ market indices: DJI, GSPC (S&P 500), IXIC (NASDAQ Composite), RUT (Russell 2000) and VIX. From 1258 one-day return rates the code estimates $\bar{r}$, the average one-day return of each index, and $\Sigma_r$, the $5 \times 5$ covariance matrix of those returns.

Instruction. Consider the following Pyomo code, which should NOT contain any mistakes.

def create_portfolio_model(rho, r_avg, cov):
    # Create the Markowitz mean/variance portfolio model in Pyomo
    #
    # Arguments:
    #     rho: required expected return of the portfolio (float)
    #     r_avg: average one-day return of each asset (pandas Series)
    #     cov: covariance matrix of the one-day returns (pandas DataFrame)

    m = pyo.ConcreteModel()

    m.ASSETS = pyo.Set(initialize=list(r_avg.index))

    m.rho = pyo.Param(initialize=rho, mutable=True, units=1 / u.day)
    m.r_avg = pyo.Param(m.ASSETS, initialize=r_avg.to_dict(), units=1 / u.day)
    m.cov = pyo.Param(
        m.ASSETS,
        m.ASSETS,
        initialize={(i, j): cov.loc[i, j] for i in m.ASSETS for j in m.ASSETS},
        units=1 / u.day**2,
    )

    m.x = pyo.Var(
        m.ASSETS, domain=pyo.NonNegativeReals, initialize=0.0, units=u.dimensionless
    )

    @m.Objective(sense=pyo.minimize)
    def OBJ(b):
        return sum(b.x[i] * b.cov[i, j] * b.x[j] for i in b.ASSETS for j in b.ASSETS)

    # C1
    @m.Constraint()
    def required_return(b):
        return sum(b.r_avg[i] * b.x[i] for i in b.ASSETS) >= b.rho

    # C2
    @m.Constraint()
    def budget(b):
        return sum(b.x[i] for i in b.ASSETS) == 1

    assert_units_consistent(m)

    return m

Solving with rho = 0.0008 per day gives:

Standard deviation of the return rate = 0.004580 per day

Optimal allocation of funds:
    DJI   30.62%
  GSPC   18.51%
  IXIC   42.93%
   RUT    0.02%
   VIX    7.93%

Expected return = 0.000800 per day

2-A. Translate the Pyomo code into a mathematical model (symbols and equations) using set notation similar to our in-class examples¶

Next to each set, parameter, variable, and equation, write a few-word description. After C1 and C2, write a brief description. You are filling in the missing Python comments here.

Sets.

Parameters (Data). Hint: state the units of each one, and say which one is mutable and why.

Variables.

Objective. Write it twice: once with explicit sums over the index sets, and once in matrix notation.

Constraints.

2-B. Degree of freedom analysis¶

Count the variables, equality constraints and inequality constraints separately and report the degrees of freedom. Is C1 active at the reported solution? How can you tell from the printed output alone?

2-C. Problem size¶

Let $N$ represent the number of assets. Using this symbol, determine the variables and constraints in your model. After each header below, give a number with a brief justification.

Number of continuous variables:

Number of integer/discrete variables:

Number of equality constraints:

Number of inequality constraints:

Number of parameters (data values) you must supply:

Finally: compare how this model scales with how the diet model in Problem 1 scales. They are not the same, and the difference is the interesting part.

2-D. Classify the problem¶

Is this an LP, QP, NLP, MILP or MINLP? Justify your classification. Where, exactly, is the nonlinearity, and what does it cost you?

2-E. The Hessian of the objective¶

In this part we explore the Hessian $\nabla_{xx} z$ of the portfolio variance. The Hessian is what determines whether the problem is convex, and how hard it is to solve numerically.

We will review the needed linear algebra, including eigenvalues, after the modeling portion of the course. For now, use the definitions and prompts here; this is practice, not an assumption that you already remember every detail.

1. Starting from $z = \sum_{i \in \mathcal{S}} \sum_{j \in \mathcal{S}} x_i \Sigma_{i,j} x_j$, compute $\partial z / \partial x_k$ and then $\partial^2 z / \partial x_k \partial x_l$. Use the symmetry of $\Sigma_r$ to simplify. Show that $\nabla_{xx} z = 2 \Sigma_r$.

2. State the condition on the eigenvalues of $\Sigma_r$ under which the problem is convex, and explain in one or two sentences why a covariance matrix always satisfies it.

3. The eigenvalues of $\Sigma_r$ for the five indices above are $9.59 \times 10^{-7}$, $8.14 \times 10^{-6}$, $1.71 \times 10^{-5}$, $1.02 \times 10^{-4}$ and $7.55 \times 10^{-3}$. Interpret the smallest one. What does it say about the portfolio, and what does it say about the numerics?

2-F. Check your work¶

Only after you have attempted 2-A through 2-E on paper, open Portfolio Optimization (on the class website). It works the same model in full, sweeps $\rho$ to trace the efficient frontier, and closes with a section on how much the estimated $\bar{r}$ can be trusted.

Check your answers against these solutions. Write one thing that you learned from this problem.

2-G. AI use report¶

Please report:

• approximately how long the independent attempt took,

• how far you got

• where you got stuck, if anywhere,

• any AI or collaborative help you used afterward

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

Problem 3. Nonlinear programs: initialization and problem formulation are very important!¶

AI category: AI permitted after independent work. For each of 3-A, 3-B, and 3-C, work for 10 minutes without AI, solution pages, or help from another person, or stop early if you complete the subproblem. Three subparts at 10 minutes is the same 30-minute independent budget the syllabus sets for every top-level problem; this assignment just spends it in three pieces. Record your attempt before using AI. You may then use AI, but you must test and verify its suggestions. Complete the report in 3-D after finishing the entire problem.

3-A. Alternative initialization¶

Effective initialization can be critical for solving nonlinear problems, since they can have several local solutions and numerical difficulties. Solve the Rosenbrock problem using different initial values for the x (and optionally y) variable(s). Write a loop that varies the initial value from 2.0 to 6.0, solves the problem, and prints the solution for each iteration of the loop.

As elaborated here, the Rosenbrock problem is a classic “hard” test case for optimization algorithms. Your results may surprise you (and show the effectiveness of Pyomo and Ipopt!).

3-B. Evaluation errors¶

Consider the following problem with initial values $x$=5, $y$=5.

$\min_{x,y} f(x,y)=(x-1.01)^2+y^2$

s.t. $y=\sqrt{x-1.0}$

1. Formulate this Pyomo model and solve using IPOPT. You should get a list of errors from the solver. Add the IPOPT solver option solver.options[‘halt_on_ampl_error’]=‘yes’ to find the problem. Hint: the error output might be ordered strangely, look this up in the console output. What did you discover? How might you fix this?

Question Answers

Fill in here

2. Add bounds $x \geq 1$ to fix this problem and resolve. Comment on the number of iterations and the quality of solution. Note, the problem still occurs because $x \geq 1$ is not enforced exactly, and small numerical values still cause the error.

Discussion

Fill in here

3. Think about other solutions for this problem and attempt to implement one of these solutions. Hint: $x \geq 1.001$.

3-C. Alternative formulations¶

Consider the following problem with initial values $x$=5, $y$=5.

$\min_{x,y} f(x,y)=(x-1.01)^2+y^2$

s.t. $\frac{x-1}{y}=1$

Note, the solution to this problem is $x$=1.005 and $y$=0.005. There are several ways that the problem above can be reformulated. Some examples are shown below. Which ones do you expect to be better? Why?

1. $\min_{x,y} f(x,y)=(x-1.01)^2+y^2$

s.t. $\frac{x-1}{y}=1$

2. $\min_{x,y} f(x,y)=(x-1.01)^2+y^2$

s.t. $\frac{x}{y+1}=1$

3. $\min_{x,y} f(x,y)=(x-1.01)^2+y^2$

s.t. $y=x-1$

Implement Pyomo models for each formulation and solve with IPOPT.

Formulation 1¶

Formulation 2¶

Formulation 3¶

Note the number of iterations and quality of solutions. What can you learn about the problem formulation from these examples?

Discussion

Fill in here

Bounds and initialization can be very helpful when solving nonlinear optimization problems. Resolve the original problem below, but add bounds, $y \geq 0$. Note the number of iterations and quality of solution, and compare with what you found for Formulation 1.

Discussion

Fill in here

3-D. AI use report¶

For each of 3-A, 3-B, and 3-C, briefly report:

• approximately how long the independent attempt took;

• how far you got during the 10-minute independent period;

• where you got stuck, if anywhere;

• how AI was helpful after the independent period; and

• how you verified any AI-generated explanation or code.

If you did not use AI for a subproblem, say so explicitly. Concise bullets are sufficient; do not submit prompts or transcripts. The time estimates help the instructor improve the assignment.

Problem 4. Reactor design (Bynum et al., 2021; Bequette, 2003)¶

AI category: AI permitted after independent work. Work on each of 4-A, 4-B, and 4-C for 10 minutes without AI, solution pages, or help from another person, or stop early if you complete the subproblem. As in Problem 3, three subparts at 10 minutes is the syllabus’s 30-minute independent budget, spent in three pieces. Record your attempt before using AI. You may then use AI, provided you verify it.

Here we will consider a chemical reactor designed to produce product B from reactant A using a reaction scheme known as the Van de Vusse reaction:

$A \overset{k_1}{\rightarrow} B \overset{k_2}{\rightarrow} C$

$2A \overset{k_3}{\rightarrow} D$

Under appropriate assumptions, $F$ is the volumetric flowrate through the tank. The concentration of component A in the feed is $c_{Af}$, and the concentrations in the reactor are equivalent to the concentrations of each component flowing out of the reactor, given by $c_A$, $c_B$, $c_C$, and $c_D$.

If the reactor is too small, we will not produce sufficient quantity of B, and if the reactor is too large, much of B will be further reacted to form the undesired product C. Therefore, our goal is to solve for the reactor volume that maximizes the outlet concentration for product B.

The steady-state mole balances for each of the four components are given by

$0=\frac{F}{V}c_{Af}-\frac{F}{V}c_A-k_1c_A-2k_3c^2_A$

$0=-\frac{F}{V}c_{B}+k_1c_A-k_2c_B$

$0=-\frac{F}{V}c_{C}+k_2c_B$

$0=-\frac{F}{V}c_{D}+k_3c^2_A$

The known parameters for the system are:

$c_{Af}=10000\frac{\mathrm{gmol}}{\mathrm{m}^3}$

$k_1=\frac{5}{6}\mathrm{min}^{-1}$

$k_2=\frac{5}{3}\mathrm{min}^{-1}$

$k_3=\frac{1}{6000}\frac{\mathrm{m}^3}{\mathrm{gmol}~\mathrm{min}}$

Since the volumetric flowrate $F$ always appears as the numerator over the reactor volume $V$, it is common to consider this ratio as a single variable, called the space-velocity $SV$.

4-A. Degree-of-freedom analysis and classification (pencil and paper)¶

Write the model in symbols, then complete the Lecture 2 counting table. Count variable bounds separately rather than subtracting them. Submit this work in the handwritten PDF.

Classify the problem as LP, QP, NLP, MILP, or MINLP and justify your answer. Then classify the continuous problem as convex or nonconvex, identifying the expressions that determine your answer.

4-B. Implement and solve in Pyomo¶

Implement the model in Pyomo, solve it with Ipopt, check the termination condition before reading variable values, and report the optimal space velocity and outlet concentration of B with units.

4-C. Multistart initialization¶

Choose at least five meaningfully different initial values for the space velocity and concentrations. Solve from each starting point, record termination and the resulting objective, and explain whether the runs support a claim of local or global optimality. A multistart study is evidence, not a proof.

References: Bynum, M. L., Hackebeil, G. A., Hart, W. E., Laird, C. D., Nicholson, B. L., Siirola, J. D., Watson, J.-P., and Woodruff, D. L. Pyomo — Optimization Modeling in Python, Third Edition. Springer Optimization and Its Applications, Vol. 67, 2021. (§7.4.4, p. 106)

B.W. Bequette. Process control: modeling, design, and simulation. Prentice Hall, 2003.

4-D. AI use report¶

For each of 4-A, 4-B, and 4-C, briefly report:

• approximately how long the independent attempt took;

• how far you got during the 10-minute independent period;

• where you got stuck, if anywhere;

• how AI was helpful after the independent period; and

• how you verified any AI-generated mathematical claims, explanations, or code.

If you did not use AI for a subproblem, say so explicitly. Concise bullets are sufficient; do not submit prompts or transcripts. The time estimates help the instructor improve the assignment.