Energy System Benchmark Problems

Optimized grid design and operation is needed as modern society depends on vast quantities of electrical energy. The complexity of the electric grid presents difficult control problems that require powerful solvers and efficient formulations for tractable solutions.

GitHub |

Google Colab

The Gekko Optimization Suite is a machine learning and optimization package in Python for mixed-integer and differential algebraic equations and is capable of solving complex grid design and control problems.

A series of non-dimensional benchmark cases are proposed for grid energy production. These include:

I: load following

II: cogeneration

III: tri-generation

IV: energy storage with constant production

V: energy storage with load following

VI: energy storage with cogeneration

Individual case studies include ramp rate constraints, power production, and energy storage operation as design variables. Variables used in the benchmark problems are defined below.

Symbol	Description
$J$ $J$	objective function
d	demand
g	generation
r	ramp rate
e	storage inventory
e_in, e_out	energy stored, recovered
R	renewable source
s_in, s_out	slack variables for storage switching
$η$ $\eta$	storage efficiency
n	number of generating units
i	subscript indicates product i

Benchmark I: Load Following

The first benchmark problem represents load following, a common scenario in grid systems. The optimizer seeks to match demand and supply with fluctuating demand dynamics. A single generator with ramping constraints attempts to respond to a single load with perfect foresight. The generation and demand match initially, but the generator must ramp in order to ensure this throughout the horizon while minimizing overproduction.

$\begin{align} \min_{r} \ & J = \int_{t=0}^{1} \left [ 1000 \max(0, d - g) + \max(0, g - d) \right ]dt \\ \mathrm{s.t.}\ & \frac{dg}{dt} = r \\ & d=\cos\left({2\pi t}\right) + 3 \label{eqn:benchmark1c}\\ & -1 \leq r \leq 1 \label{eqn:benchmark1d}\\ & g(0)=d(0)=4 \end{align}$

Benchmark II: Cogeneration

In the second benchmark problem, one producer seeks to meet two objectives that are constraining at different times. Benchmark II enhances Benchmark I by replacing the generator with a cogeneration system n=2 that produces (1) electricity and (2) heat in response to electricity demand and a new heat demand profile, again with perfect foresight.

$\begin{align} \min_{r} \ \ & J = \sum_{i=1}^n \int_{t=0}^{1} \, \left [ 1000 \; \max(0, \, d_i - g_i) \vphantom{\int} + \max(0, \, g_i - d_i) \vphantom{\int} \right ]dt \\ \mathrm{s.t.}\ \ & \frac{dg_1}{dt} = r, \enspace g_2 = 2 g_1 \label{eqn:benchmark2b}\\ & d_1=\cos\left({2\pi t}\right) + 3 \\ & d_2=1.5\sin\left({2\pi t}\right) + 7 \\ & -1 \leq r \leq 1 \\ & g_1(0)=d_1(0)=4 \\ & g_2(0)=8, \; d_2(0)=7 \end{align}$

Benchmark III: Tri-generation

The third benchmark problem enhances the previous problems further still, creating a tri-generation system n=3 with two producers, three products, and three demand profiles. The primary producer is the same as the prior benchmark and is ramp-rate constrained to produce the two primary products (e.g., electricity and heat). An additional producer (e.g., a solid oxide electrolysis cell) uses these first two products to make a third product (e.g., hydrogen), thereby utilizing any excess system capacity and maximizing its production while avoiding supply shortages for products one and two. As before, the optimizer has perfect foresight.

$\begin{align} \min_{r_1, r_3} \ \ & J = \sum_{i=1}^n \int_{t=0}^{1} \left [ \vphantom{\int} 1000 \; \max(0, \, d_i - g_i) + \max(0, \, g_i - d_i) \vphantom{\int} \right ]dt - 0.1 \int_{t=0}^{1} g_3 dt \\ \mathrm{s.t.}\ \ & \frac{dg_1}{dt} = r_1, \; \frac{dg_3}{dt} = r_3 \\ & g_2 = 2 g_1 \\ & d_1=\cos\left({2\pi t}\right) + 3 + 2 g_3 \\ & d_2=1.5\sin\left({2\pi t}\right) + 7 + 3 g_3 \\ & d_3=\max\left[0,-0.2\sin\left(2\pi t\right) \right] \\ & -1 \leq r_1 \leq 1, \enspace -1 \leq r_3 \leq 1 \\ & g_1(0)=d_1(0)=4 \\ & g_2(0)=8, \; d_2(0)=7 \\ & g_3(0)=d_3(0)=0 \end{align}$

Benchmark IV: Constant Production with Storage

The fourth benchmark problem models a hybrid system with a single generator with constant production constraints coupled with energy storage that together must meet an oscillating electricity demand. The goal of the problem is to minimize the required power production and use energy storage to capture excess generation serve the oscillating energy demand while keeping the base-load generator production constant. As before, the model has perfect foresight. In order to prevent the energy storage from charging and discharging simultaneously without requiring mixed-integer variables, slack variables are used to control when the storage charges and discharges, allowing it to switch modes in a way that is both continuous and differentiable. This allows the modeling language to use automatic differentiation to provide exact gradients to the solver.

$\begin{align} \min_g \ \ & g \\ \mathrm{s.t.} & \frac{de}{dt} = e_{\text{in}} - e_{\text{out}} \cdot \eta \\ & e_{\text{in}} = g - d + s_{\text{in}} \\ & e_{\text{out}} = d - g + s_{\text{out}} \\ & g - d = s_{\text{out}} - s_{\text{in}} \\ & s_{\text{out}}, \, s_{\text{in}} \geq 0, \enspace e_{\text{out}} \cdot e_{\text{in}} \leq 0 \\ & g + e_{\text{out}}/\eta - e_{\text{in}} \geq d \\ & e \geq 0, \enspace \eta = 0.7 \\ & d = 10 - 2 \sin(2 \, \pi \, t) \\ & e(0) = e(1) = 0 \end{align}$

Benchmark V: Load Following with Storage

The fifth benchmark combines energy storage with a load-following problem similar to Benchmark I. The first-half of the time horizon is nearly identical to Benchmark I, but now the excess energy can be stored. This allows the system to meet a higher demand in the second half of the time horizon without needing extremes in generation. The solver minimizes the ramping needs and operates more flexibly by storing and then recovering the overproduction caused by the ramping constraints. Energy storage allows this generator to meet the load without requiring significant overproduction.

$\begin{align} \min_{r} \ \ & J = \sum_{i=1}^n \int_{t=0}^{1} \, [ 1000 \, \max(0, \, \Psi) \vphantom{\int} + \, \max(0, \, -\Psi) \vphantom{\int} ]dt \\ & \text{where} \ \Psi = d - g - R + e_{\text{out}}/\eta - e_{\text{in}} \\ \mathrm{s.t.}\ \ & \frac{de}{dt} = e_{\text{in}} - e_{\text{out}} \cdot \eta \\ & e_{\text{in}} = g + R - d + s_{\text{in}} \\ & e_{\text{out}} = d - g - R + s_{\text{out}} \\ & g + R - d = s_{\text{out}} - s_{\text{in}} \\ & s_{\text{out}}, \, s_{\text{in}} \geq 0, \enspace e_{\text{out}} \cdot e_{\text{in}} \leq 0\\ & g + R + e_{\text{out}}/\eta - e_{\text{in}} \geq d \\ & e \geq 0, \enspace \eta = 0.85 \\ & d = 7 - 2 \sin(2 \, \pi \, t) \\ & R = \begin{cases} 3 + 3 \cos(4 \, \pi \, t) & \frac{1}{4} \leq t \leq \frac{3}{4} \\ 0 & \mbox{otherwise} \\ \end{cases} \\ & \frac{dg}{dt} = r, \enspace -1\leq r \leq 1 \\ & e(0) = e(1) = 0 \end{align}$

Benchmark VI: Cogeneration with Dual Storage

The sixth benchmark problem is a combination of Benchmarks II and V where the ramp rate of the generator is the manipulated variable but now must meet both electrical (1) and heat (2) demand with the use of both electrical and thermal storage. A renewable generation source (such as solar PV) is added to the system as an auxiliary electrical energy source that cannot be controlled. The objective is the same as in Benchmark IV, to minimize power production, but this time while meeting both the heat and power demands. The model again has perfect foresight. Slack variables are used in the same way as in Benchmarks IV and V.

$\begin{align} \min_{r_1} \ \ & g_1 \\ \mathrm{s.t.}\ \ & \frac{de_i}{dt} = e_{\text{in},i} - e_{\text{out},i}\cdot \eta_i \\ & e_{\text{in},i} = g_i + R_i - d_i + s_{\text{in},i} \\ & e_{\text{out},i} = d_i - g_i - R_i + s_{\text{out},i} \\ & g_i + R_i - d_i = s_{\text{out},i} - s_{\text{in},i} \\ & s_{\text{out},i}, \, s_{\text{in},i} \geq 0, \enspace e_{\text{out},i} \cdot e_{\text{in},i} \leq 0 \\ & g_i + R_i + e_{\text{out},i}/\eta_i - e_{\text{in},i} \geq d_i \\ & e_i \geq 0, \enspace \eta_1 = 0.7, \enspace \eta_2 = 0.8 \\ & d_1 = 10 - 2 \sin(2 \, \pi \, t) \\ & d_2 = 15 + 1.5 \cos(2 \, \pi \, t) \\ & R_1 = \begin{cases} 2 + 2 \cos(4 \, \pi \, t) & \frac{1}{4} \leq t \leq \frac{3}{4} \\ 0 & \mbox{otherwise} \\ \end{cases} \\ & \frac{dg_1}{dt} = r_1, \enspace -3\leq r_1 \leq 3 \\ & R_2 = 0, \enspace g_2 = 1.5 \cdot g_1 \\ & e_1(0) = 0, \enspace e_2(0) = e_2(1) = 0.5 \end{align}$

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

m = GEKKO(remote=False)
t = np.linspace(0, 1, 101)
m.time = t

m.options.SOLVER = 3
m.options.IMODE = 6
m.options.NODES = 2
m.options.CV_TYPE = 1
m.options.MAX_ITER = 1000

p = m.SV(10) # production (constant)
s = m.Var(0.1, lb=0) # storage inventory
stored = m.SV() # store energy rate
recovery = m.SV() # recover energy rate
vx = m.SV(lb=0) # recover slack variable
vy = m.SV(lb=0) # store slack variable

eps = 0.85 # Storage efficiency

d = m.MV((-20*np.sin(np.pi*t/12*24)+100)/10)
d_h = m.MV((15*np.cos(np.pi*t/12*24)+150)/10)

p_h_initial = m.Intermediate(p*1.5)

p_h = m.SV(p_h_initial)
s_h = m.Var(0.5,lb=0)
stored_h = m.SV()
recovery_h = m.SV()

#renewable energy source
renewable = (20*np.cos(np.pi*t/6*24)+20)/10
center = np.ones(len(t))
num = len(t)
center[0:int(num/4)] = 0
center[-int(num/4):] = 0
renewable *= center
r = m.Param(renewable)

r1 = m.MV(ub=3,lb=-3)
r1.STATUS=1

m.periodic(s_h)

zx = m.SV(lb=0)
zy = m.SV(lb=0)

eps_h = 0.8 # heat storage efficiency

net = m.Intermediate(d-r)
m.Equations([p + r + recovery/eps - stored >= d,
p + r - d == vx - vy,
stored == p + r - d + vy,
recovery == d - p - r + vx,
s.dt() == stored - recovery/eps,
p.dt() == r1,
stored * recovery <= 0,
p_h + recovery_h/eps_h - stored_h >= d_h,
p_h - d_h == zx - zy,
stored_h == p_h - d_h + zy,
recovery_h == d_h - p_h + zx,
s_h.dt() == stored_h - recovery_h/eps_h,
stored_h * recovery_h <= 0,
p_h == 1.5 * p])
m.Minimize(p)
m.solve()

# Plot solution
fig, axes = plt.subplots(5, 1, figsize=(5, 5.1), sharex=True)
axes = axes.ravel()

ax = axes[0]
ax.plot(t, d, 'r-', label='Demand 1 ($d_1$)')
ax.plot(t, p,'b:', label='Production 1 ($g_1$)',lw=2)
ax.plot(t, net, 'k--', label='Net ($d_1-R_1$)')

ax = axes[1]
ax.plot(t,r, 'b-',label='Source 1 ($R_1$)')
ax.plot(t,r1, 'k--', label='Ramp Rate ($r$)')

ax = axes[2]
ax.plot(t,s, 'k-', label='Storage 1 ($e_1$)')
ax.plot(t,stored,'g--',label='Stored ($e_{\text{in},1}$)')
ax.plot(t,recovery,'b:',label='Recovered ($e_{\text{out},1}$)',lw=2)

ax = axes[3]
ax.plot(t,d_h, 'r-', label='Demand 2 ($d_2$)')
ax.plot(t[1:], p_h.value[1:],'b:',\
label='Production 2 ($g_2$)',lw=2)

ax = axes[4]
ax.plot(t,s_h, 'k-', label='Storage 2 ($e_2$)')
ax.plot(t,stored_h,'g--',label='Stored ($e_{\text{in},2}$)')
ax.plot(t[1:],recovery_h.value[1:],'b:',\
label='Recovered ($e_{\text{out},2}$)',lw=2)
ax.set_xlabel('Time')

for ax in axes:
ax.legend(loc='center left',\
bbox_to_anchor=(1,0.5),frameon=False)
ax.grid()
ax.set_xlim(0, 1)
plt.tight_layout()
plt.savefig('grid_energy6.png', dpi=600,\
bbox_inches = 'tight')
plt.show()

[$[Get Code]]

Reference

Gates, N.S., Hill, D.C., Billings, B.W., Powell, K.M., Hedengren, J.D., Benchmarks for Grid Energy Management with Python Gekko, 60th IEEE Conference on Decision and Control (CDC), Tutorial Session: Open Source Software for Control Applications, Austin, TX, Dec 13-15, 2021. Preprint

Acknowledgement

Financial support from the Department of Energy Nuclear Energy University Program Grant Number DE-NE0008866.