Pyomo Mini-Project: Receding Horizon Stochastic Control#

Assignment Goals#

The primary purpose of this assignment is to give you hands-on-keyboard experience with optimization modeling and analysis. In this assignment, I incorporated many of the best practices I share with my Ph.D. students during their first year using Pyomo for research. Specifically:

Planning is the most important (and often most overlooked step) in computer programming. Two best practices are to (1) write pseudocode and (2) write the optimization model on paper and perform degree of freedom analysis before you use the computer. We practiced tip (2) in class. If you cannot do (1) and (2), no amount of Googling Pyomo/Python syntax will save you.
Functions (or classes) allow us to reuse most of our code. Below I gave you function definitions and doc strings to show you a good way to “modularize” this analysis. Once the functions defined below work, performing the analysis is easy.
Solving your optimization problem a few times is (often) not a sufficient research result. Who cares? Application-driven research is all about the analysis. Optimization is a tool to better answer analysis questions. By giving you >50% of the code for this assignment, I am intentionally emphasizing the analysis (not just coding). You will see we can take a really simple optimization model and ask/answer some interesting analysis questions.

You should start this assignment early. Almost everyone will make at least one Python/Pyomo syntax mistake. Give yourself plenty of time to ask questions on Slack. Please also be generous with helping your classmates. This assignment (and grades in general in this class) is not a competition. You are welcome (encouraged!) to post screenshots of your code in Slack when asking/answering questions. Please refrain from posting your discussion responses. You are welcome to discuss the analysis questions on Slack; everyone should answer the discussion questions in their own words.

# This code cell installs packages on Colab

import sys
if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/CBE60499/main/notebooks/helper.py"
    import helper
    helper.install_idaes()
    helper.install_ipopt()
    helper.download_data(['Prices_DAM_ALTA2G_7_B1.csv'])

import pandas as pd
import pyomo.environ as pyo
import numpy as np
import matplotlib.pyplot as plt

ca_data = pd.read_csv('https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/data/Prices_DAM_ALTA2G_7_B1.csv',names=['price'])

def prepare_data_array(day=2,num_days=8):
    ''' Prepare numpy array of prices
    
    Arguments:
        day: day to start (0 = January 1, 2015)
        num_days: numbers of days to include in the vector
        
    Returns:
        numpy array of prices
    '''

    return ca_data["price"][(day)*24:24*(day+num_days)].to_numpy()

In the previous notebook, we considered model predictive control (MPC) of a battery energy storage system interacting with the electricity market. In essence, we posed an optimization problem:

Given a forecasting for energy prices, what charging and discharging strategy maximizes net market revenue?

In this homework mini-project, you will consider two modifications to make the example more realistic:

Add a constraint to prevent simultanous charging and discharging.
Receding horizon control
Uncertainty

The goal is to help you develop Python and Pyomo programming skills.

Model Refinement and Code Reorganization#

Revised Model#

Propose a new mathematical model that includes a constraint to prevent charging and discharging. Hint: using either integer variables or a disjunction will ensure the objective and constraints remaining linear.

Below is the model description copied from our in-class example. Please edit such that this notebook describes your complete new model.

Sets

Timesteps: $\mathcal{T} = \{0, 1, ..., N\}$

Timesteps without the initial time: $\mathcal{T}' = \{1, ..., N\}$

Variables

$E_t$, energy stored in battery at time $t$, units: MWh
$d_t$, battery discharge power (sold to market) during time interval [t-1, t), units: MW
$c_t$, battery charge power (purchased from the market) during time interval [t-1, t), units: MW

Parameters

$\pi_t$: Energy price during time interval [t-1, t), units: $/MW
$\eta$: Round trip efficiency, units: dimensionless
$c_{max}$ Maximum charge power, units: MW
$d_{max}$ Maximum discharge power, units: MW
$E_{max}$ Maximum storage energy, units: MWh
$E_{0}$ Energy in storage at time $t=0$, units: MWh
$\Delta t = 1$ hour, Timestep for grid decisions and prices (fixed)

Optimization Problem $$

\[\begin{align*} \max_{\mathbf{E},\mathbf{d},\mathbf{c}} \quad & \psi := \sum_{t \in \mathcal{T}'} \pi_{t} \Delta t (d_{t} - c_{t}) \\ \mathrm{s.t.} \quad & E_{t} = E_{t-1} + \Delta t \left( c_{t} \sqrt{\eta} - \frac{d_{t}}{\sqrt{\eta}} \right), ~~ \forall ~ t \in \mathcal{T}' \\ & E_{0} = E_{N} \\ & 0 \leq c_{t} \leq c_{max}, ~~\forall ~ t \in \mathcal{T}' \\ & 0 \leq d_{t} \leq d_{max}, ~~\forall ~ t \in \mathcal{T}' \\ & 0 \leq E_{t} \leq E_{max}, ~~\forall ~ t \in \mathcal{T}' \end{align*}\]

Degree of Freedom Analysis#

Perform degree of freedom analysis. Identify:

Total number of continous variables
Total number of integer variables (or disjunctions if using)
Total number of equality constraints
Total number of inequality constraints
Degrees of freedom

Tip: perform this on paper then typeset below.

Number of continous variables:

Number of integer variables:

Number of equality constraints:

Number of inequality constraints:

Degrees of freedom. Please discuss in a few sentences. There are a few ways to do this; we are less concerned about you getting the “right” answer and more about giving an answer with a supporting argument.

Python Code#

Programmers organize their code into functions, often following the rule “1 task/action = 1 function”. Here are a few of the reasons you should always do this:

Functions make your code easier to use
More specifically, functions force you to abstract/generalize your code
Writing doc strings (required for this course!) forces you to document your code and think about the structure.
Functions allow your code to best tested piece by piece

To show you how to do this, we have provided skeletons for several functions in this assignment. Based on the doc string, complete the functions.

def build_model(price_data, horizon_length, include_pbc=True, include_disjunction=True):
    '''
    Builds HORIZON pyomo model
    
    Inputs:
    price_data: dictionary of price data
    horizon_length: number of timesteps to evaluate
    include_pbc: Boolean, include periodic boundary condition constraint
    
    Returns:
    m: HORIZON pyomo model
    '''
    
    m = pyo.ConcreteModel()
    
    # Add your solution here
    
    return m

def extract_solution(m):
    '''
    Function to extract the solution from the solved pyomo model
    
    Inputs:
        m: Solved pyomo model
    
    Outputs:
        c_control: numpy array of charging power values, MW
        d_control: numpy array of discharging power values, MW
        E_control: numpy array of energy stored values, MWh
        t: numpy array of time indices from the Pyomo model, hr
    '''
    
    c_control = np.empty(m.N)
    d_control = np.empty(m.N)
    E_control = np.empty(m.N)
    t = np.empty(m.N)

    # Tip: We added a negative sign to d to help with plotting in the code discussed during class.
    # Add your solution here
    
    return c_control, d_control, E_control, t

def plot_solution(c, d, E, t=None):
    '''
    Function to plot the solution of the HORIZON problem
    
    Inputs:
        c_control: numpy array of charging power values, MW
        d_control: numpy array of discharging power values, MW
        E_control: numpy array of energy stored values, MWh
        t: numpy array of time indices from the Pyomo model, hr. (default=None)
    
    Actions. Creates the following plots: 
        State of charge (energy) verus time
        Power to the grid verus time
    '''
    
    # Plot the state of charge (E)
    plt.figure()

    N = len(c)
    assert N > 0, "Input t must have length greate than 0"
    assert len(d) == N, "All inputs should be the same length"
    assert len(E) == N, "All inputs should be the same length"
    
    if t is None:
        # add E0
        t = np.array(range(1,N+1))
        t_ = np.array(range(0,N+1))
    else:
        t_ = np.concatenate(([0],  t))

    # Add your solution here
    
    # Plot charging/discharging plot
    plt.figure()

    # double up first data point to make the step plot
    c_control_ = np.concatenate(([c[0]],  c))
    d_control_ = np.concatenate(([d[0]],  d))

    plt.step(t_,c_control_,'r.-',where='pre')
    plt.step(t_,d_control_,'g.-',where='pre')
    plt.xlabel('Time (hr)')
    plt.ylabel('Power from Grid (MW)')
    plt.grid(True)
    plt.show()
    
    return

Optimal Control Over 3 Day Horizon#

Now let’s use our functions for some analysis. Calculate the optimal battery operation strategy over 3 days starting on January 4, 2015.

three_days = prepare_data_array(4,3)

Original Model#

Hint: After solving your model, plot the solution. Does the solution make sense? If not, debug your model before attempting the rolling horizon optimization.

# Add your solution here

With Constraint Enforcing No Simultaneous Charging or Discharging#

Hint: After solving your model, plot the solution and visually compare to the previous solution. Does the result make sense? Don’t overthink this.

# Add your solution here

Comparison#

Write Python code below to compare the solutions. Print out any timesteps where the charge, discharge, and energy profiles are different by more than 10$^{-4}$ MW (or MWh).

def compare_solutions(solution1,solution2):
    '''
    Function to compare outputs between the original model and model with disjuncutions
    
    Inputs:
    original: list of original charging/discharging data
    disjunction: list of disjunction model results charging/discharging data
    t: list of timesteps
    
    Action:
    prints same/different + timestep
    '''
    
    assert len(solution1) == len(solution2), "solution 1 and solution2 must be the same length"
    
    n = len(solution1)
    
    # Add your solution here
    print("Done.")

# Charge
compare_solutions(c_control_og, c_control_disj)

# Discharge
compare_solutions(d_control_og, d_control_disj)

# Energy
compare_solutions(E_control_og, E_control_disj)

Receding Horizon Control#

We now have functions to build, solve, and analyze the battery problem for a single problem. Now we can consider receeding horizon control.

Please ask questions during class.

Nomenclature: Let $u^*_{i}$ represent the optimal control action from time $t=i$ to $t=i+1$. For the battery system, $u_{i} = [c_i, d_i]^{T}$ (charging and discharging).

Algorithm:

Consider time $t=m$.
Compute the optimal control actions $u^*=[u^*_m, u^*_{m+1}, ..., u^*_{m+N-1}]^T$ from time $t=m$ until $t=N+m$ where $N$ is the length of the planning horizon. In our battery example, we considered $N = 24$.
Simulate the system for one timestep using control action $u^*_{m}$.
(Optional). Update the system model or forecasts using new observations.
Step $t=m+1$ and goto Step 2.

Write Pseudocode#

Pseudocode is an important planning tool for complex computer programs. The idea is to write out the code on paper or a whiteboard using generic programming syntax. As an example, consider calculating the 3rd through 10th elements of the Fibonacci sequence:

# Governing equation where f[i] is the ith element of the sequence
# f[i] = f[i-1] + f[i-2]

# Algorithm:

# Previous previous number in the sequence, f[i-2]
n_2prev = 1

# Previous number in the sequence, f[i-1]
n_1prev = 1

# Set counter
j = 3

while j <= 10 do:
    # Apply Fib. formula
    n = n1_prev + n2_prev
    
    # Print to screen
    print(n)
    
    # Shift history.
    # Previous number becomes previous previous number
    n2prev = n1prev
    
    # Current number becomes previous number
    n1prev = n
    
    # Increment counter
    j = j + 1

Here we sketched out the main algorithm (with comments!) in generic syntax. We did not worry about the correct way to write a while loop in Python. Once we figure out the algorithm we can worry about those details.

Write pseudocode for the receding horizon control example on paper or a whiteboard. Your algorithm should include the following main steps:

Load price data from a text file
Loop over timesteps. Keep track of implemented control actions (charge and discharge), battery SOC, and net revenue.
Make a plot of battery SOC and control actions that were actually implemented.

Scan your pseudocode as a PDF and turn in via Gradescope.

Python Code#

Tip: You should use the same create_model function throughout the assignment. In your model, define price to be a mutable parameter. This means we can just update the Pyomo model in the function below. We do not need to recreate it.

def update_model(model, new_price_data, new_E0):
    """ Update the Pyomo model with new price data
    
    Arguments:
        model: pyomo model
        new_price_data: numpy array or list containing price data
        new_E0: (float) value for initial battery charge
    """
    
    # Add your solution here

def receding_horizon_simulation(price_data, simulation_length, num_hours, include_pbc=True,verbose=True):  
    '''
    Perform receding horizon simulation
    
    Inputs:
        price_data: numpy array price data for entire time period of interest
        simulation_length: number of total hours to evaluate
        num_hours: horizon length for each optimization problem
        include_pbc: boolean, include periodic boundary condition
        verbose: boolean, if true, print out message each hour
    
    Returns:
        c_control: charging decision implemented
        d_control: discharging decision implemented
        E_control: state-of-charge profile
        revenue_control: total revenue over simulation horizon
    
    '''
    
    # Check input data length
    assert len(price_data) >= simulation_length + num_hours, "Input price data is too short"
    
    assert simulation_length > 0, "Must evaluate more than 0 days"
    
    assert num_hours > 0, "Must consider planning horizon longer than 0 days"

    nS = simulation_length
        
    ## Preallocate vectors to store results
    c_control = np.empty(nS)
    d_control = np.empty(nS)
    E_control = np.empty(nS)
    revenue_control = np.empty(nS)

    # Add your solution here

    return c_control, d_control, E_control, revenue_control

Simulate Performance Over 3 Days#

Perform receding horizon control for the first 3 days in January 2015. Make a plot to show the implemented control action.

Tips:

The implemented control action is a concationation of the first timestep for each optimization problem.
How can you resuse the plot function above?

three_days = prepare_data_array(0,4)

# Tip: do not overwrite these results. We will use them later.
# If you reuse this code, change the variable names used to save the results.
c_control_pbc, d_control_pbc, E_control_pbc, revenue_control_pbc = receding_horizon_simulation(three_days, 24*3, 
                                                                               24, include_pbc=True,
                                                                              verbose=True)

plot_solution(c_control_pbc, d_control_pbc, E_control_pbc)

Impact of Horizon Length#

Repeat the sensitivity analysis for horizon lengths of N=2 through N=48 (hours). For each simulation, compute the following:

Revenue from first 30 days
Total energy purchased from the market during the first 30 days
Total energy sold to the market during the first 30 days
Total computation time

# Study this cell to see how to time calculations in Python
import time

time_start = time.perf_counter()
first_month = prepare_data_array(0,34)
time_elapsed = (time.perf_counter() - time_start)

print("Elapsed Time",time_elapsed,"s")

rev = []
charge = []
discharge = []
timer = []

# Tip: This will take a few minutes to finish. While developing your solutions, 
# only simulate the first few values of N. This will be helpful if you need to
# restart your kernel and rerun the cells. Once you are happy with your code
# for the entire assignment, you can set N to the values given here then
# restart and run all.
N = [2, 4, 6, 8, 12, 24, 36, 48]

# Add your solution here

Make four well labeled plots to show these trends.

plt.plot(N,rev,'o-',color="purple",linewidth=3)
plt.xlabel('N',fontsize=18)
plt.ylabel('Revenue, $',fontsize=18)
plt.grid(True)
plt.show()

plt.plot(N,charge,'ro-',linewidth=3)
plt.xlabel('N',fontsize=18)
plt.ylabel('Charging, [MWh]',fontsize=18)
plt.grid(True)
plt.show()

plt.plot(N,discharge,'go-',linewidth=3)
plt.xlabel('N',fontsize=18)
plt.ylabel('Discharging, [MWh]',fontsize=18)
plt.grid(True)
plt.show()

plt.plot(N,timer,'ko-',linewidth=3)
plt.xlabel('N',fontsize=18)
plt.ylabel('Computational Time, [s]',fontsize=18)
plt.grid(True)
plt.show()

Discussion How long of a simulation horizon do you recommend? Discuss the trade-off between revenue and computation time.

Impact of Periodic Boundary Condition#

We will now explore the impact of the periodic boundary constraint.

First, simulate the performance over 3 days without enforcing the periodic boundary constraint. Plot the implemented control action and calculate the total revenue. Use $N=4$.

c_control2, d_control2, E_control2, revenue_control = receding_horizon_simulation(three_days, 24*3, 4,
                                                                                  include_pbc=False,verbose=True)

plot_solution(c_control2, d_control2, E_control2)

Next, simulate with the periodic boundary constraint using $N=4$. Plot the results.

c_control3, d_control3, E_control3, revenue_control = receding_horizon_simulation(three_days, 24*3, 4,
                                                                                  include_pbc=True,verbose=False)
plot_solution(c_control3, d_control3, E_control3)

Now let’s compare.

# Charge
compare_solutions(c_control2, c_control3)

# Discharge
compare_solutions(d_control2, d_control3)

# Energy
compare_solutions(E_control2, E_control3)

Discussion Compare the control actions and revenue with and without the periodic boundary constraint. How are they different (if at all)? Offer an explanation for any differences. (The purpose of this questions is to get you thinking about how to critically analyze the optimization results. This skill will be really important for the semester projects.)

Next, repeat the horizon length sensitivity analysis without the peridoic boundary constraint.

# Add your solution here

Discussion Compare these result to the analysis with periodic boundary constraints. Based on these results, make a recommendation to a battery operator.

Forecast Uncertainty (required for CBE 60499, optional for CBE 40499)#

So far, we have assumed the battery operator can perfectly forecast the market price. This is not realistic. What is the impact of uncertainty?

Generate white noise with mean zero and standard deviation 5 $/MWh using np.random.normal() and add this to the historical price. We will treat this as a price forecast. Simulate the first three days with $N=24$. Calculate how much less revenue you make with price uncertainty compared to the perfect information case. Hint: Each time you run the simulation, you will get a different answer. Repeat the analysis 10 times and record the average.

Are you enrolled in CBE 60499 (are you a graduate student)? Please answer yes/no:

revenue_unc_predicted = []
revenue_unc_actual = []

# Perform simulation with perfect information
c_control, d_control, E_control, revenue_control = receding_horizon_simulation(three_days, 
                                                                                       24*3, 24, include_pbc=True,
                                                                                  verbose=False)
# We can also calculate the revenue after the fact
rev_check = 0
for t in range(len(c_control)):
    # Why is this slightly different? In the instructor/TA solutions for receding_horizon_simulation,
    # we added a negative sign to d_control for plotting.
    rev_check -=  three_days[t]*(d_control[t] + c_control[t])

revenue_perfect_information = np.sum(revenue_control)
print("Revenue with perfect price forecast =",revenue_perfect_information,"=",rev_check,"$")

# Tip: normalize the revenue from the uncertain simulations by `revenue_perfect_information`
# This will make the interpretation much easier!

n = len(three_days)
nsim = 10

for i in range(nsim):

    # Add random noise to data
    three_days_with_noise = three_days + 5*np.random.normal(size=n)

    # Add your solution here
    
print("Average Revenue =",np.mean(revenue_unc_actual)/revenue_perfect_information*100,"%")

Discussion Compare the results with your earlier analysis. How much less revenue do you get with Gaussian noise with standard deviation 5 $/MWh.

Next, make a plot showing the average revenue as a function of the standard deviation of the uncertainty. Hint: Normalize the revenue by dividing by the perfect information revenue (no uncertainty).

# Add your solution here

Discussion What is the impact with the larger uncertainty?

Discussion Why does uncertainty in the price decrease the revenue (on average)?

Nonlinear and Stochastic Optimization

Pyomo Mini-Project: Receding Horizon Stochastic Control

Contents

Pyomo Mini-Project: Receding Horizon Stochastic Control#

Assignment Goals#

Model Refinement and Code Reorganization#

Revised Model#

Degree of Freedom Analysis#

Python Code#

Optimal Control Over 3 Day Horizon#

Original Model#

With Constraint Enforcing No Simultaneous Charging or Discharging#

Comparison#

Receding Horizon Control#

Write Pseudocode#

Python Code#

Simulate Performance Over 3 Days#

Impact of Horizon Length#

Impact of Periodic Boundary Condition#

Forecast Uncertainty (required for CBE 60499, optional for CBE 40499)#