This notebook contains material from ND-Pyomo-Cookbook; content is available on Github. The text is released under the CC-BY-NC-ND-4.0 license, and code is released under the MIT license.

4.6 Unit Commitment

Keywords: semi-continuous variables, cbc usage, gdp, disjunctive programming

4.6.1 Imports

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import display, HTML

import shutil
import sys
import os.path

if not shutil.which("pyomo"):
    !pip install -q pyomo
    assert(shutil.which("pyomo"))

if not (shutil.which("cbc") or os.path.isfile("cbc")):
    if "google.colab" in sys.modules:
        !apt-get install -y -qq coinor-cbc
    else:
        try:
            !conda install -c conda-forge coincbc 
        except:
            pass

assert(shutil.which("cbc") or os.path.isfile("cbc"))
import pyomo.environ as pyo
import pyomo.gdp as gdp

4.6.2 Problem statement

A set of $N$ electrical generating units are available to meet a required demand $d_t$ for time period $t \in 1, 2, \ldots, T$. The power generated by unit $n$ for time period $t$ is denoted $x_{n,t}$. Each generating unit is either off, $x_{n,t} = 0$ or else operating in a range $[p_n^{min}, p_n^{max}]$. The incremental cost of operating the generator during period $t$ is $a_n x_{n,t} + b_n$. A binary variable variable $u_{n,t}$ indicates the operational state of a generating unit.

The unit commmitment problem is then

\begin{align*} \min \sum_{n\in N} \sum_{t\in T} a_n x_{n,t} + b_n u_{n,t} \end{align*}

subject to

\begin{align*} \sum_{n\in N} x_{n,t} & = d_t \qquad \forall t \in T \\ p_{n}^{min}u_{n,t} & \leq x_{n,t} \qquad \forall n \in N, \ \forall t \in T \\ p_{n}^{max}u_{n,t} & \geq x_{n,t} \qquad \forall n \in N, \ \forall t \in T \\ \end{align*}

where we use the short-cut notation $T = [1, 2, \ldots T]$ and $N = [1, 2, \ldots, N]$.

This is a minimal model. A realistic model would include additional constraints corresponding to minimum up and down times for generating units, limits on the rate at which power levels can change, maintenance periods, and so forth.

4.6.3 Model

4.6.3.1 Demand

In [2]:
# demand
T = 20
T = np.array([t for t in range(0, T)])
d = np.array([100 + 100*np.random.uniform() for t in T])

fig, ax = plt.subplots(1,1)
ax.bar(T+1, d)
ax.set_xlabel('Time Period')
ax.set_title('Demand')
Out[2]:
Text(0.5, 1.0, 'Demand')

4.6.3.2 Generating Units

In [3]:
# generating units
N = 5
pmax = 2*max(d)/N
pmin = 0.6*pmax

N = np.array([n for n in range(0, N)])
a = np.array([0.5 + 0.2*np.random.randn() for n in N])
b = np.array([10*np.random.uniform() for n in N])

p = np.linspace(pmin, pmax)

fig, ax = plt.subplots(1,1)
for n in N:
    ax.plot(p, a[n]*p + b[n])
ax.set_xlim(0, pmax)
ax.set_ylim(0, max(a*pmax + b))
ax.set_xlabel('Unit Production')
ax.set_ylabel('Unit Operating Cost')
ax.grid()

4.6.3.3 Pyomo model 1: Conventional implementation for emi-continuous variables

In [4]:
def unit_commitment():
    m = pyo.ConcreteModel()

    m.N = pyo.Set(initialize=N)
    m.T = pyo.Set(initialize=T)

    m.x = pyo.Var(m.N, m.T, bounds = (0, pmax))
    m.u = pyo.Var(m.N, m.T, domain=pyo.Binary)
    
    # objective
    m.cost = pyo.Objective(expr = sum(m.x[n,t]*a[n] + m.u[n,t]*b[n] for t in m.T for n in m.N), sense=pyo.minimize)
    
    # demand
    m.demand = pyo.Constraint(m.T, rule=lambda m, t: sum(m.x[n,t] for n in N) == d[t])
    
    # semi-continuous
    m.lb = pyo.Constraint(m.N, m.T, rule=lambda m, n, t: pmin*m.u[n,t] <= m.x[n,t])
    m.ub = pyo.Constraint(m.N, m.T, rule=lambda m, n, t: pmax*m.u[n,t] >= m.x[n,t])
    return m
   
m = unit_commitment()
pyo.SolverFactory('cbc').solve(m).write()

fig, ax = plt.subplots(max(N)+1, 1, figsize=(8, 1.5*max(N)+1))
for n in N:
    ax[n].bar(T+1, [m.x[n,t]() for t in T])
    ax[n].set_xlim(0, max(T)+2)
    ax[n].set_ylim(0, 1.1*pmax)
    ax[n].plot(ax[n].get_xlim(), np.array([pmax, pmax]), 'r--')
    ax[n].plot(ax[n].get_xlim(), np.array([pmin, pmin]), 'r--')
    ax[n].set_title('Unit ' + str(n+1))
fig.tight_layout()
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 1018.71533244
  Upper bound: 1018.71533244
  Number of objectives: 1
  Number of constraints: 200
  Number of variables: 180
  Number of binary variables: 100
  Number of integer variables: 100
  Number of nonzeros: 180
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.06
  Wallclock time: 0.06
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
    Black box: 
      Number of iterations: 0
  Error rc: 0
  Time: 0.07869720458984375
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

4.6.3.4 Pyomo model 2: GDP implementation

In [5]:
def unit_commitment_gdp():
    m = pyo.ConcreteModel()

    m.N = pyo.Set(initialize=N)
    m.T = pyo.Set(initialize=T)

    m.x = pyo.Var(m.N, m.T, bounds = (0, pmax))
    
    # demand
    m.demand = pyo.Constraint(m.T, rule=lambda m, t: sum(m.x[n,t] for n in N) == d[t])
    
    # representing the semicontinous variables as disjuctions
    m.sc1 = gdp.Disjunct(m.N, m.T, rule=lambda d, n, t: d.model().x[n,t] == 0)
    m.sc2 = gdp.Disjunct(m.N, m.T, rule=lambda d, n, t: d.model().x[n,t] >= pmin)
    m.sc = gdp.Disjunction(m.N, m.T, rule=lambda m, n, t: [m.sc1[n,t], m.sc2[n,t]])
    
    # objective. Note use of the disjunct indicator variable
    m.cost = pyo.Objective(expr = sum(m.x[n,t]*a[n] + m.sc2[n,t].indicator_var*b[n] for t in m.T for n in m.N), sense=pyo.minimize)

    # alternative formulation. But how to access the indicator variable?
    #m.semicontinuous = gdp.Disjunction(m.N, m.T, rule=lambda m, n, t: [m.x[n,t]==0, m.x[n,t] >= pmin])
    pyo.TransformationFactory('gdp.chull').apply_to(m)
    return m
   
m_gdp = unit_commitment_gdp()
pyo.SolverFactory('cbc').solve(m_gdp).write()
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 863.60019688
  Upper bound: 863.60019688
  Number of objectives: 1
  Number of constraints: 20
  Number of variables: 100
  Number of binary variables: 200
  Number of integer variables: 200
  Number of nonzeros: 100
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.02
  Wallclock time: 0.03
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
    Black box: 
      Number of iterations: 0
  Error rc: 0
  Time: 0.04323315620422363
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

4.6.3.5 There is a problem here!

Why are the results different? Somehow it appears values of the indicator variables are being ignored.

In [16]:
for n in N:
    for t in T:
        print("n = {0:2d}  t = {1:2d}   {2} {3}  {4:5.2f}".format(n, t, m_gdp.sc1[n,t].indicator_var(), m_gdp.sc2[n,t].indicator_var(), m.x[n,t]()))
n =  0  t =  0   1.0 0.0  76.13
n =  0  t =  1   1.0 0.0  45.86
n =  0  t =  2   1.0 0.0  45.86
n =  0  t =  3   1.0 0.0  75.96
n =  0  t =  4   1.0 0.0  45.86
n =  0  t =  5   1.0 0.0  45.86
n =  0  t =  6   1.0 0.0  45.86
n =  0  t =  7   1.0 0.0  73.80
n =  0  t =  8   1.0 0.0  68.79
n =  0  t =  9   1.0 0.0  61.04
n =  0  t = 10   1.0 0.0  47.89
n =  0  t = 11   1.0 0.0  56.14
n =  0  t = 12   1.0 0.0  45.86
n =  0  t = 13   1.0 0.0  45.86
n =  0  t = 14   1.0 0.0  45.86
n =  0  t = 15   1.0 0.0  47.08
n =  0  t = 16   1.0 0.0  45.86
n =  0  t = 17   1.0 0.0  53.14
n =  0  t = 18   1.0 0.0  45.86
n =  0  t = 19   1.0 0.0  47.56
n =  1  t =  0   1.0 0.0   0.00
n =  1  t =  1   1.0 0.0   0.00
n =  1  t =  2   1.0 0.0   0.00
n =  1  t =  3   1.0 0.0   0.00
n =  1  t =  4   1.0 0.0   0.00
n =  1  t =  5   1.0 0.0   0.00
n =  1  t =  6   1.0 0.0   0.00
n =  1  t =  7   1.0 0.0   0.00
n =  1  t =  8   1.0 0.0   0.00
n =  1  t =  9   1.0 0.0   0.00
n =  1  t = 10   1.0 0.0   0.00
n =  1  t = 11   1.0 0.0   0.00
n =  1  t = 12   1.0 0.0   0.00
n =  1  t = 13   1.0 0.0   0.00
n =  1  t = 14   1.0 0.0   0.00
n =  1  t = 15   1.0 0.0   0.00
n =  1  t = 16   1.0 0.0   0.00
n =  1  t = 17   1.0 0.0   0.00
n =  1  t = 18   1.0 0.0   0.00
n =  1  t = 19   1.0 0.0   0.00
n =  2  t =  0   1.0 0.0   0.00
n =  2  t =  1   1.0 0.0   0.00
n =  2  t =  2   1.0 0.0   0.00
n =  2  t =  3   1.0 0.0   0.00
n =  2  t =  4   1.0 0.0   0.00
n =  2  t =  5   1.0 0.0   0.00
n =  2  t =  6   1.0 0.0   0.00
n =  2  t =  7   1.0 0.0   0.00
n =  2  t =  8   1.0 0.0   0.00
n =  2  t =  9   1.0 0.0   0.00
n =  2  t = 10   1.0 0.0   0.00
n =  2  t = 11   1.0 0.0   0.00
n =  2  t = 12   1.0 0.0   0.00
n =  2  t = 13   1.0 0.0   0.00
n =  2  t = 14   1.0 0.0   0.00
n =  2  t = 15   1.0 0.0   0.00
n =  2  t = 16   1.0 0.0   0.00
n =  2  t = 17   1.0 0.0   0.00
n =  2  t = 18   1.0 0.0   0.00
n =  2  t = 19   1.0 0.0   0.00
n =  3  t =  0   1.0 0.0  76.43
n =  3  t =  1   1.0 0.0  61.93
n =  3  t =  2   1.0 0.0  64.26
n =  3  t =  3   1.0 0.0  76.43
n =  3  t =  4   1.0 0.0  57.86
n =  3  t =  5   1.0 0.0  67.42
n =  3  t =  6   1.0 0.0  65.12
n =  3  t =  7   1.0 0.0  76.43
n =  3  t =  8   1.0 0.0  76.43
n =  3  t =  9   1.0 0.0  76.43
n =  3  t = 10   1.0 0.0  76.43
n =  3  t = 11   1.0 0.0  76.43
n =  3  t = 12   1.0 0.0  72.97
n =  3  t = 13   1.0 0.0  72.25
n =  3  t = 14   1.0 0.0  72.77
n =  3  t = 15   1.0 0.0  76.43
n =  3  t = 16   1.0 0.0  70.91
n =  3  t = 17   1.0 0.0  76.43
n =  3  t = 18   1.0 0.0  72.74
n =  3  t = 19   1.0 0.0  76.43
n =  4  t =  0   1.0 0.0   0.00
n =  4  t =  1   1.0 0.0  45.86
n =  4  t =  2   1.0 0.0  45.86
n =  4  t =  3   1.0 0.0   0.00
n =  4  t =  4   1.0 0.0   0.00
n =  4  t =  5   1.0 0.0  45.86
n =  4  t =  6   1.0 0.0  45.86
n =  4  t =  7   1.0 0.0   0.00
n =  4  t =  8   1.0 0.0  45.86
n =  4  t =  9   1.0 0.0   0.00
n =  4  t = 10   1.0 0.0   0.00
n =  4  t = 11   1.0 0.0   0.00
n =  4  t = 12   1.0 0.0  45.86
n =  4  t = 13   1.0 0.0   0.00
n =  4  t = 14   1.0 0.0   0.00
n =  4  t = 15   1.0 0.0  45.86
n =  4  t = 16   1.0 0.0   0.00
n =  4  t = 17   1.0 0.0   0.00
n =  4  t = 18   1.0 0.0  45.86
n =  4  t = 19   1.0 0.0  45.86