Merit order¶
The ucp is a mixed- integer combinatorial optimization problem including uncertain supply from renewable energies (e.g. wind, solar), potential machine failure or demand. The objective is to allocate power ressources to match a certain demand at all times producing minimal cost.
In the easiest way the problem is equivalent to the knapsack problem: https://en.wikipedia.org/wiki/Knapsack_problem
More information here: https://ercim-news.ercim.eu/en128/special/energy-economics-fundamental-modelling-with-quantum-algorithms
- Start with an easy example
- introduce resolution
- introduce slack variable to formulate unequalitites
- minimum/maximum up/down unit commitment problem including satisfiablitiy formulation to formulate QUBO
In [4]:
# import pygrnd and other libraries needed
# we build on top of the open source framework qiskit (qiskit.org)
import pygrnd
from pygrnd.qc.helper import *
from pygrnd.qc.brm import brm
from pygrnd.qc.brm_oracle import brmoracle
from pygrnd.qc.QAE import qae
from pygrnd.optimize.sat_ucp import *
from pygrnd.optimize.meritorder import *
from pygrnd.optimize.bruteforce import *
from pygrnd.optimize.MonteCarloSolver import *
from pygrnd.optimize.qaoa import *
from qiskit import execute
from qiskit import Aer
from math import pi
import math
import cmath
import random
from qiskit.visualization import plot_histogram
import matplotlib.pyplot as plt
from scipy.stats import norm
import networkx as nx
from IPython.display import Image
Motivation QUBO (Quadratic Unconstrained Binary Optimization)¶
- quadratic unconstrained binary optimization
- Minimize/maximize $\langle x | Q | x \rangle = \sum_{ij} x_i Q_{ij} x_j$
- Binary variables $x_i \in \{0,1\} \iff x_i = x_i^2$
- entries in symmetric matrix Q are real numbers
- solves combinatorial problems
constraints need to be encoded as penalty terms
Constraints for knapsack encoded in parameters $Q_{ij}$
$Q = Q_{cost} + Q_{constraint}$
the optimal solution is hard to find
- QUBOs can be solved by different Ising solver, annealer, quantum annealer, quantum simulator
QUBO stands equivalent with Ising spin model
Cost are encoded on the main diagonal $c_i$
Knapsack Problem as QUBO¶
- Minimize -$\sum_i x_i v_i$ s.t. $\sum_i x_i w_i = W$
- Binary variables $x_i\in\{0,1\}$, i.e., take an element or do not take it at all
- Maximize sum of selected values $v_i \in \mathbb{R}_{\geq 0}$
- Respect constraint $\sum_i x_i w_i \leq W$ with weights $w_i \in \mathbb{R}_{\geq 0}$ and maximum weight $W\geq 0$
- Introduce slack variable to obtain equality $\sum_i x_i w_i +s = W$ with $s\in \mathbb{R}$
QUBO constraint for demand¶
D = W = demand
$ w_i = weights_i$ = $maxgen_i$ of each unit i
$ \sum_i w_i = D $
$\implies (\sum_i w_i - D)^2 = 0 $
$ \implies \sum_i w_i \sum w_j - 2*D \sum_i w_i + D^2 = 0$
for $i = j: \quad (\sum_i w_i )^2 - 2*D*\sum_i w_i $ # main diagonal matrix elements
for $i \neq j: \quad \sum_i w_i * \sum_i w_j $
we ignore constant ("offset", $D^2$) - needs to be added/subtracted from the solution
Rewrite Constraints To Obtain QUBO Formulation¶
- $\sum_i x_i w_i =W$ can be written as $\left( \sum_i x_i w_i -W \right)^2$
- Solve QUBO -$\sum_i v_i x_i + P \left( \sum_i x_i w_i -W \right)^2$
- Find appropriate penalty factor $P$
- Use $x^2_i = x_i$ for binary variables $x_i$
Knapsack With Resolution¶
- $0/1$ value for $x_i$ should be more fine-grained
- Solution: Split $x_i$ in several parts and represent it as $0...0$ to $1...1$ with $\frac{1}{2^{b}-1}$ per part
- Example: $b=3$ leads to $\frac{0}{7}, \frac{1}{7}, \ldots, \frac{7}{7}$
- Components of Costs and Weights must be weighted by $\frac{1}{7}, \frac{2}{7}, \frac{4}{7}$ for $(x_0,x_1,x_2)$
In [5]:
# Weight constraint violated. P too small.
M=QUBO_knapsack([1,3,4,5],[1,3,4,5],7,1)
print(M)
solver(M)
[[-12. 3. 4. 5.] [ 3. -30. 12. 15.] [ 4. 12. -36. 20.] [ 5. 15. 20. -40.]]
16it [00:00, 20834.79it/s]
Out[5]:
(-42.0, matrix([[0, 1, 1, 0]]))
In [6]:
# Weight constraint satisfied.
M=QUBO_knapsack([1,3,4,5],[1,3,4,5],7,2)
print(M)
solver(M)
## penalty finetuning is art
[[-25. 6. 8. 10.] [ 6. -63. 24. 30.] [ 8. 24. -76. 40.] [ 10. 30. 40. -85.]]
16it [00:00, 68061.73it/s]
Out[6]:
(-91.0, matrix([[0, 1, 1, 0]]))
In [7]:
# Two solutions are possible: 10/01 and 11/10. Both should lead to 5/3 as result.
# Penalty 3 should be necessary, 2 leads to wrong results.
values_split=splitParameters([1,2],2)
weights_split=splitParameters([1,2],2)
M=QUBO_knapsack(values_split,weights_split,5/3,3)
print(M)
res,v=solver(M)
print(v)
res=recombineSolution(values_split, [1,0,0,1], 2)
print(sum(res))
res=recombineSolution(values_split, [1,1,1,0], 2)
print(sum(res))
[[-2.66666667 0.66666667 0.66666667 1.33333333] [ 0.66666667 -4.66666667 1.33333333 2.66666667] [ 0.66666667 1.33333333 -4.66666667 2.66666667] [ 1.33333333 2.66666667 2.66666667 -6.66666667]]
16it [00:00, 48489.06it/s]
[[0 1 1 0]] 1.6666666666666665 1.6666666666666665
Slack Variable for Upper Weight Bound¶
- Transform $\sum_i w_i x_i \leq W$ to $\sum_i w_i x_i + s =W$ with a slack value $s \in \mathbb{R}_{\geq 0}$
- Decompose $s=\frac{s_0}{2^c-1}+ \frac{2s_1}{2^c-1}+\frac{4s_2}{2^c-1}+\ldots$ with sufficient resolution of $c$ bits
- Use the normal fractional method from above and just extend slack variable with weight $W$ and value $0$
In [8]:
QUBO_knapsack_slack_resolution([1,2],[1,1],3,3,10,True)
512it [00:00, 86037.01it/s]
[0, 0, 0, 0, 0, 0, 1, 1, 1] -90.0 recombined fractions: [0.0, 0.0, 1.0] recombined values: [0.0, 0.0, 0.0] recombined weights: [0.0, 0.0, 3.0] total/slack value: 0.0 0.0 real/slack weight: 0.0 3.0 total/demanded/diff weight: 3.0 3 0.0
Out[8]:
0.0
In [9]:
# With penalty 7.1 we start to see good results
for i in range(1,100):
buffer=i*0.1
res=QUBO_knapsack_slack_resolution([1,2],[1,2],5/7,3,buffer)
if abs(res-5/7)<0.1:
print(i,buffer,abs(res-5/7))
break
512it [00:00, 79927.19it/s] 512it [00:00, 97875.38it/s] 512it [00:00, 102627.65it/s] 512it [00:00, 114440.91it/s] 512it [00:00, 97964.68it/s] 512it [00:00, 116495.80it/s] 512it [00:00, 108245.56it/s] 512it [00:00, 113473.38it/s] 512it [00:00, 113731.79it/s] 512it [00:00, 96260.87it/s] 512it [00:00, 118462.25it/s] 512it [00:00, 92703.81it/s] 512it [00:00, 102090.97it/s] 512it [00:00, 103229.52it/s] 512it [00:00, 116660.35it/s] 512it [00:00, 119344.43it/s] 512it [00:00, 103913.85it/s] 512it [00:00, 116174.39it/s] 512it [00:00, 95274.34it/s] 512it [00:00, 121094.15it/s] 512it [00:00, 109025.93it/s] 512it [00:00, 120361.15it/s] 512it [00:00, 108158.33it/s] 512it [00:00, 113449.40it/s] 512it [00:00, 122608.26it/s] 512it [00:00, 117702.58it/s] 512it [00:00, 117137.60it/s] 512it [00:00, 119563.70it/s] 512it [00:00, 114203.55it/s] 512it [00:00, 125981.68it/s] 512it [00:00, 116787.23it/s] 512it [00:00, 119490.52it/s] 512it [00:00, 121244.56it/s] 512it [00:00, 121760.14it/s] 512it [00:00, 125952.12it/s] 512it [00:00, 119463.93it/s] 512it [00:00, 112604.67it/s] 512it [00:00, 125334.64it/s] 512it [00:00, 105678.05it/s] 512it [00:00, 120869.23it/s] 512it [00:00, 122196.63it/s] 512it [00:00, 123135.53it/s] 512it [00:00, 121560.27it/s] 512it [00:00, 121587.80it/s] 512it [00:00, 116869.86it/s] 512it [00:00, 115580.39it/s] 512it [00:00, 118129.91it/s] 512it [00:00, 120957.74it/s] 512it [00:00, 115998.68it/s] 512it [00:00, 121822.31it/s] 512it [00:00, 121032.73it/s] 512it [00:00, 118214.45it/s] 512it [00:00, 125466.44it/s] 512it [00:00, 120916.87it/s] 512it [00:00, 117451.52it/s] 512it [00:00, 126226.04it/s] 512it [00:00, 107724.29it/s] 512it [00:00, 111187.93it/s] 512it [00:00, 115035.55it/s] 512it [00:00, 112876.93it/s] 512it [00:00, 105134.81it/s] 512it [00:00, 117220.72it/s] 512it [00:00, 110008.90it/s] 512it [00:00, 112286.73it/s] 512it [00:00, 106400.62it/s] 512it [00:00, 120638.37it/s] 512it [00:00, 106802.79it/s] 512it [00:00, 117054.60it/s] 512it [00:00, 114930.89it/s] 512it [00:00, 119125.96it/s] 512it [00:00, 113858.42it/s] 512it [00:00, 109192.23it/s] 512it [00:00, 111720.09it/s] 512it [00:00, 106263.73it/s] 512it [00:00, 97118.47it/s] 512it [00:00, 106332.13it/s] 512it [00:00, 114440.91it/s] 512it [00:00, 105439.37it/s] 512it [00:00, 113497.37it/s] 512it [00:00, 116331.73it/s] 512it [00:00, 96386.16it/s] 512it [00:00, 128185.02it/s] 512it [00:00, 104363.30it/s] 512it [00:00, 107767.53it/s] 512it [00:00, 120638.37it/s] 512it [00:00, 110797.84it/s] 512it [00:00, 115208.35it/s] 512it [00:00, 104530.94it/s] 512it [00:00, 119152.40it/s] 512it [00:00, 113371.54it/s] 512it [00:00, 106511.44it/s] 512it [00:00, 121498.37it/s] 512it [00:00, 113079.02it/s] 512it [00:00, 113073.06it/s] 512it [00:00, 117974.16it/s] 512it [00:00, 113858.42it/s] 512it [00:00, 113858.42it/s] 512it [00:00, 110746.41it/s] 512it [00:00, 98063.09it/s]
example 1¶
We only consider 3 units and demand = 20 = D
One time step
No ramp up/down costs, no min/max up/down times
| unit | maxgen (MW) | cost (€) |
|---|---|---|
| 1 | 15 | 2 |
| 2 | 5 | 10 |
| 3 | 5 | 10 |
In [10]:
#small example
demand = 20
penalty = 1
M=QUBO_knapsack([2,10,10],[15,5,5],demand,penalty)
print(M)
res,vec=solver(M)
print("offset = ",demand**2)
print("Solution = ",res+demand**2,"| Vector = ",vec)
[[-373. 75. 75.] [ 75. -165. 25.] [ 75. 25. -165.]]
8it [00:00, 20674.33it/s]
offset = 400 Solution = 12.0 | Vector = [[1 0 1]]
In [11]:
# Here a bigger example
# Modelling of up parameter only - no time dependecies
maxgen=[1000,800,1000,700,350]
cost=[6000,22000,34000,26600,23100]
demand=1000
print(cost)
[6000, 22000, 34000, 26600, 23100]
In [12]:
penalty=2
offset=penalty*((demand)**2)
print("Demand = ",demand)
print("Penalty = ",penalty)
units=["Nuc","Lign","Coal","CC","GT"]
print(units)
M=QUBO_knapsack(cost,maxgen,demand,penalty)
print(M)
print(solver(M))
res,vec=solver(M)
print("offset = ",offset)
print("Solution = ",res+offset,"| Vector = ",vec)
Demand = 1000 Penalty = 2 ['Nuc', 'Lign', 'Coal', 'CC', 'GT'] [[-1994000. 1600000. 2000000. 1400000. 700000.] [ 1600000. -1898000. 1600000. 1120000. 560000.] [ 2000000. 1600000. -1966000. 1400000. 700000.] [ 1400000. 1120000. 1400000. -1793400. 490000.] [ 700000. 560000. 700000. 490000. -1131900.]]
32it [00:00, 55188.21it/s]
(-1994000.0, matrix([[1, 0, 0, 0, 0]]))
32it [00:00, 87211.00it/s]
offset = 2000000 Solution = 6000.0 | Vector = [[1 0 0 0 0]]
In [ ]: