Unit commitment problem¶
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 this notebook we show how to formulate the optimization problem as a QUBO that can lateron solved with quantum algorithms, quantum annealing, quantum-inspired methods or classical heuristics
In [34]:
# 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
import timeit
In [ ]:
In [ ]:
In [35]:
## XXXS
cost=[2,1]
demand=[1]
minup=[1,1]
mindown=[0,0]
maxup=[0,0]
maxdown=[0,0]
maxgen=[1,1]
startcost=[1,1]
start_time = timeit.default_timer()
xxxs,mapsxxxs=createSATqubo(cost, maxgen, demand, minup, mindown, maxup, maxdown, startcost)
end_time = timeit.default_timer()
print('Duration in seconds: {}'.format(end_time - start_time))
print("Size of Qubo: ",len(xxxs),"x",len(xxxs))
print(xxxs)
print(mapsxxxs)
bestenergy,res=dwaveGreedySolver(xxxs,1000)
print(bestenergy,res)
SAT formulation for constraint qubo: [[-2, 1], [-1, 2], [-4, 3], [-3, 4]] Number of formulas: 4 maxVariable: 4 Constraint QUBO: [[ 2. -4. 0. 0.] [ 0. 2. 0. 0.] [ 0. 0. 2. -4.] [ 0. 0. 0. 2.]] Cost QUBO: [[ 4. -4. 0. 0.] [ 0. 2. 0. 0.] [ 0. 0. 3. -4.] [ 0. 0. 0. 2.]] penaltyDemand: 100 Demand QUBO: [[-96. -4. 200. 0.] [ 0. 2. 0. 0.] [ 0. 0. -97. -4.] [ 0. 0. 0. 2.]] Duration in seconds: 0.0022689440002068295 Size of Qubo: 4 x 4 [[-96. -4. 200. 0.] [ 0. 2. 0. 0.] [ 0. 0. -97. -4.] [ 0. 0. 0. 2.]] [[0], [2]] #nonzeros/#quboEntries 7 / 16 sparsity 0.4375 -99.0 [0, 0, 1, 1]
In [36]:
## XXS
cost=[2,1]
demand=[1]
minup=[1,1]
mindown=[0,0]
maxup=[0,0]
maxdown=[0,0]
maxgen=[1,1]
startcost=[1,1]
start_time = timeit.default_timer()
xxs,mapsxxs=createSATqubo(cost, maxgen, demand, minup, mindown, maxup, maxdown, startcost)
#xxs,mapsxxs=createSATquboPenalty(cost, maxgen, demand, minup, mindown, maxup, maxdown, startcost,100,1,1,1) #100,50,10,10
end_time = timeit.default_timer()
print('Duration in seconds: {}'.format(end_time - start_time))
print("Size of Qubo: ",len(xxs),"x",len(xxs))
print(xxs)
print(mapsxxs)
bestenergyxxs,resxxs=dwaveGreedySolver(xxs,1000)
print(bestenergyxxs,resxxs)
SAT formulation for constraint qubo: [[-2, 1], [-1, 2], [-4, 3], [-3, 4]] Number of formulas: 4 maxVariable: 4 Constraint QUBO: [[ 2. -4. 0. 0.] [ 0. 2. 0. 0.] [ 0. 0. 2. -4.] [ 0. 0. 0. 2.]] Cost QUBO: [[ 4. -4. 0. 0.] [ 0. 2. 0. 0.] [ 0. 0. 3. -4.] [ 0. 0. 0. 2.]] penaltyDemand: 100 Demand QUBO: [[-96. -4. 200. 0.] [ 0. 2. 0. 0.] [ 0. 0. -97. -4.] [ 0. 0. 0. 2.]] Duration in seconds: 0.0017852540004241746 Size of Qubo: 4 x 4 [[-96. -4. 200. 0.] [ 0. 2. 0. 0.] [ 0. 0. -97. -4.] [ 0. 0. 0. 2.]] [[0], [2]] #nonzeros/#quboEntries 7 / 16 sparsity 0.4375 -99.0 [0, 0, 1, 1]
In [37]:
## XS
cost=[2,1]
demand=[1,2,1]
#demand=[1,2,1,2,1,2,1,2,1,2,1,2,1,2,1]
#minup=[2,3]
minup=[2,3]
mindown=[0,0]
maxup=[0,0]
maxdown=[0,0]
maxgen=[1,1]
startcost=[1,1]
start_time = timeit.default_timer()
xs,mapsxs=createSATqubo(cost, maxgen, demand, minup, mindown, maxup, maxdown, startcost)
#xs,mapsxs=createSATquboPenalty(cost, maxgen, demand, minup, mindown, maxup, maxdown, startcost,100,1,1,1) #100,50,10,10
end_time = timeit.default_timer()
print('Duration in seconds: {}'.format(end_time - start_time))
print("Size of Qubo: ",len(xs),"x",len(xs))
print(xs)
print(mapsxs)
bestenergyxs,resxs=dwaveGreedySolver(xs,1000)
print(bestenergyxs,resxs)
SAT formulation for constraint qubo: [[-4, 1], [-4, 2], [-5, 2], [-5, 3], [-1, 4], [-2, 4, 5], [-3, 5], [-9, 6], [-9, 7], [-9, 8], [-6, 9], [-7, 9], [-8, 9]] Number of formulas: 13 maxVariable: 9 Constraint QUBO: [[ 2. 0. 0. -4. 0. 0. 0. 0. 0. 0.] [ 0. -1. 0. -4. -2. 0. 0. 0. 0. 3.] [ 0. 0. 2. 0. -4. 0. 0. 0. 0. 0.] [ 0. 0. 0. 7. 0. 0. 0. 0. 0. -3.] [ 0. 0. 0. 0. 2. 0. 0. 0. 0. 2.] [ 0. 0. 0. 0. 0. 2. 0. 0. -4. 0.] [ 0. 0. 0. 0. 0. 0. 2. 0. -4. 0.] [ 0. 0. 0. 0. 0. 0. 0. 2. -4. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 6. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -3.]] Cost QUBO: [[ 4. -1. 0. -4. 0. 0. 0. 0. 0. 0.] [ 0. 2. -1. -4. -2. 0. 0. 0. 0. 3.] [ 0. 0. 5. 0. -4. 0. 0. 0. 0. 0.] [ 0. 0. 0. 7. 0. 0. 0. 0. 0. -3.] [ 0. 0. 0. 0. 2. 0. 0. 0. 0. 2.] [ 0. 0. 0. 0. 0. 3. -1. 0. -4. 0.] [ 0. 0. 0. 0. 0. 0. 4. -1. -4. 0.] [ 0. 0. 0. 0. 0. 0. 0. 4. -4. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 6. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -3.]] penaltyDemand: 100 Demand QUBO: [[ -96. -1. 0. -4. 0. 200. 0. 0. 0. 0.] [ 0. -298. -1. -4. -2. 0. 200. 0. 0. 3.] [ 0. 0. -95. 0. -4. 0. 0. 200. 0. 0.] [ 0. 0. 0. 7. 0. 0. 0. 0. 0. -3.] [ 0. 0. 0. 0. 2. 0. 0. 0. 0. 2.] [ 0. 0. 0. 0. 0. -97. -1. 0. -4. 0.] [ 0. 0. 0. 0. 0. 0. -296. -1. -4. 0.] [ 0. 0. 0. 0. 0. 0. 0. -96. -4. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 6. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -3.]] Duration in seconds: 0.002686586999516294 Size of Qubo: 10 x 10 [[ -96. -1. 0. -4. 0. 200. 0. 0. 0. 0.] [ 0. -298. -1. -4. -2. 0. 200. 0. 0. 3.] [ 0. 0. -95. 0. -4. 0. 0. 200. 0. 0.] [ 0. 0. 0. 7. 0. 0. 0. 0. 0. -3.] [ 0. 0. 0. 0. 2. 0. 0. 0. 0. 2.] [ 0. 0. 0. 0. 0. -97. -1. 0. -4. 0.] [ 0. 0. 0. 0. 0. 0. -296. -1. -4. 0.] [ 0. 0. 0. 0. 0. 0. 0. -96. -4. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 6. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -3.]] [[0, 1, 2], [5, 6, 7]] #nonzeros/#quboEntries 27 / 100 sparsity 0.27 -595.0 [0, 1, 0, 0, 1, 1, 1, 1, 1, 0]
In [38]:
dffull, solutionVector, bufferSupply, solutionCostTotal, TotalSupply=printSolution(resxs, mapsxs, cost, minup, demand, maxgen)
total_violate, ratio_violate, num_constraints, unmatched_demand, unmatched_demand_ratio = checkMinupMinDownMaxupMaxdownChecks(minup,mindown,maxup,maxdown,solutionVector,TotalSupply,demand)
Total costs of production: 5 ['010', '111'] unit: 0 --> Minup VIOLATION -1 unit: 1 --> Minup Correct no mindown constraint no mindown constraint no maxup constraint no maxup constraint no maxdown constraint no maxdown constraint demand = supply --> constraint matched unmatched demand ratio = 1.0 Demand = 1 | Supply = 1 demand = supply --> constraint matched unmatched demand ratio = 1.0 Demand = 2 | Supply = 2 demand = supply --> constraint matched unmatched demand ratio = 1.0 Demand = 1 | Supply = 1 Number of demand constraint breaks = 0
In [39]:
## S
cost=[2,1]
demand=[1,2,1,2,1]
minup=[2,3]
mindown=[0,0]
maxup=[0,0]
maxdown=[0,0]
maxgen=[1,1]
startcost=[1,1]
start_time = timeit.default_timer()
#s,mapss=createQUBOCostMaxgenDemandMinupMindownMaxupMaxDownStartcostIterativeSort(cost, maxgen, demand, minup, mindown, maxup, maxdown, startcost)
## createSATqubo(cost,maxgen,demand,minup,mindown,maxup,maxdown,startcost,PDemand,PCost,PStart,PConstr)
s,mapss=createSATqubo(cost, maxgen, demand, minup, mindown, maxup, maxdown, startcost) #100,50,10,10
end_time = timeit.default_timer()
print('Duration in seconds: {}'.format(end_time - start_time))
print("Size of Qubo: ",len(s),"x",len(s))
print(mapss)
bestenergys,ress=dwaveGreedySolver(s,1000)
print(bestenergys,ress)
dffull, solutionVector, bufferSupply, solutionCostTotal, TotalSupply=printSolution(ress, mapss, cost, minup, demand, maxgen)
total_violate, ratio_violate, num_constraints, unmatched_demand, unmatched_demand_ratio = checkMinupMinDownMaxupMaxdownChecks(minup,mindown,maxup,maxdown,solutionVector,TotalSupply,demand)
SAT formulation for constraint qubo: [[-6, 1], [-6, 2], [-7, 2], [-7, 3], [-8, 3], [-8, 4], [-9, 4], [-9, 5], [-1, 6], [-2, 6, 7], [-3, 7, 8], [-4, 8, 9], [-5, 9], [-15, 10], [-15, 11], [-15, 12], [-16, 11], [-16, 12], [-16, 13], [-17, 12], [-17, 13], [-17, 14], [-10, 15], [-11, 15, 16], [-12, 15, 16, 17], [-13, 16, 17], [-14, 17]]
Number of formulas: 27
maxVariable: 17
Constraint QUBO:
[[ 2. 0. 0. 0. 0. -4. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.]
[ 0. -1. 0. 0. 0. -4. -2. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 3.
0. 0. 0. 0. 0. 0.]
[ 0. 0. -1. 0. 0. 0. -4. -2. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
3. 0. 0. 0. 0. 0.]
[ 0. 0. 0. -1. 0. 0. 0. -4. -2. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 3. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 2. 0. 0. 0. -4. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 7. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -3.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 5. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 2.
-3. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 5. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
2. -3. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 2. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 2. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 2. 0. 0. 0. 0. -4. 0. 0. 0.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. -4. -2. 0. 0.
0. 0. 3. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. -4. -2. -2. 0.
0. 0. 0. 3. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. -4. -2. 0.
0. 0. 0. 0. 0. 3.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 2. 0. 0. -4. 0.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 12. 0. 0. 0.
0. 0. -3. -3. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 8. 0. 0.
0. 0. 2. 2. -3. -3.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 2. 0.
0. 0. 0. 0. 2. 2.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -3.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
-3. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. -3. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. -3. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. -3. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. -3.]]
Cost QUBO:
[[ 4. -1. 0. 0. 0. -4. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.]
[ 0. 2. -1. 0. 0. -4. -2. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 3.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 2. -1. 0. 0. -4. -2. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
3. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 2. -1. 0. 0. -4. -2. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 3. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 5. 0. 0. 0. -4. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 7. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -3.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 5. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 2.
-3. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 5. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
2. -3. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 2. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 2. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 3. -1. 0. 0. 0. -4. 0. 0. 0.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -1. 0. 0. -4. -2. 0. 0.
0. 0. 3. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -1. 0. -4. -2. -2. 0.
0. 0. 0. 3. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. -1. 0. -4. -2. 0.
0. 0. 0. 0. 0. 3.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 4. 0. 0. -4. 0.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 12. 0. 0. 0.
0. 0. -3. -3. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 8. 0. 0.
0. 0. 2. 2. -3. -3.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 2. 0.
0. 0. 0. 0. 2. 2.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -3.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
-3. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. -3. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. -3. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. -3. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. -3.]]
penaltyDemand: 100
Demand QUBO:
[[ -96. -1. 0. 0. 0. -4. 0. 0. 0. 200. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. -298. -1. 0. 0. -4. -2. 0. 0. 0. 200. 0.
0. 0. 0. 0. 0. 3. 0. 0. 0. 0. 0. 0.]
[ 0. 0. -98. -1. 0. 0. -4. -2. 0. 0. 0. 200.
0. 0. 0. 0. 0. 0. 3. 0. 0. 0. 0. 0.]
[ 0. 0. 0. -298. -1. 0. 0. -4. -2. 0. 0. 0.
200. 0. 0. 0. 0. 0. 0. 3. 0. 0. 0. 0.]
[ 0. 0. 0. 0. -95. 0. 0. 0. -4. 0. 0. 0.
0. 200. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 7. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. -3. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 5. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 2. -3. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 5. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 2. -3. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 2. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 2. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. -97. -1. 0.
0. 0. -4. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -299. -1.
0. 0. -4. -2. 0. 0. 0. 0. 3. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -99.
-1. 0. -4. -2. -2. 0. 0. 0. 0. 3. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
-299. -1. 0. -4. -2. 0. 0. 0. 0. 0. 0. 3.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. -96. 0. 0. -4. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 12. 0. 0. 0. 0. 0. -3. -3. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 8. 0. 0. 0. 0. 2. 2. -3. -3.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 2. 0. 0. 0. 0. 0. 2. 2.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. -3. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. -3. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. -3. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. -3. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -3. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -3.]]
Duration in seconds: 0.009583523000401328
Size of Qubo: 24 x 24
[[0, 1, 2, 3, 4], [9, 10, 11, 12, 13]]
#nonzeros/#quboEntries 73 / 576
sparsity 0.1267361111111111
-1105.0 [1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1]
Total costs of production: 11
['11011', '01110']
unit: 0 --> Minup Correct
unit: 1 --> Minup Correct
no mindown constraint
no mindown constraint
no maxup constraint
no maxup constraint
no maxdown constraint
no maxdown constraint
demand = supply --> constraint matched
unmatched demand ratio = 1.0
Demand = 1 | Supply = 1
demand = supply --> constraint matched
unmatched demand ratio = 1.0
Demand = 2 | Supply = 2
demand = supply --> constraint matched
unmatched demand ratio = 1.0
Demand = 1 | Supply = 1
demand = supply --> constraint matched
unmatched demand ratio = 1.0
Demand = 2 | Supply = 2
demand = supply --> constraint matched
unmatched demand ratio = 1.0
Demand = 1 | Supply = 1
Number of demand constraint breaks = 0
new we solve the QUBO with the inbuild QAOA solver¶
In [40]:
import pygrnd
from pygrnd.optimize.qaoa import *
from pygrnd.optimize.sat_ucp import *
from qiskit.visualization import plot_histogram
import dimod
import qiskit
import numpy as np
In [41]:
## show solutions landscape
qaoaLandscape(xxs,10,1000)
0.0 0.0 1000 -96.0 0.0 0.3141592653589793 0011 -99.0
In [42]:
# solve example set xxs
# xs is a 4x4 qubo that requires 4 qubits
vec, counts, obj, prob, qc, res1, res2, bestBetas, bestGammas = QAOAoptimizeMaxCount(xxs,5,1000)
Selected device: qasm_simulator with 1000 shots
Trying 5 layer
Generating inital random parameters beta and gamma
Starting with betas: [3.0316061425719845, -0.1519250911700416, -2.060103847622278, 1.6365896789072494, -5.345778860338202]
Starting with gammas: [4.000753853832176, 5.4180642003749675, 1.9802957639121512, -6.1007873218954725, -5.06725191258352]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [ 3.03160614 -0.15192509 -2.06010385 1.63658968 -5.34577886] Best Gamma [ 4.00075385 5.4180642 1.98029576 -6.40582669 -5.06725191]
Now run the QAOA with the found parameters
----------------------------------------------
Optimum = -98.0 with probability = 0.274
----------------------------------------------
Depth: 67
Gate counts: OrderedDict([('rz', 35), ('cx', 30), ('ry', 20), ('barrier', 11), ('h', 4), ('measure', 4)])
In [43]:
plot_histogram(counts)
Out[43]:
In [44]:
## do not approx expectation value by average over all solutions but only use the max count bitstring
Nmax=[]
probamax=[]
Optmax=[]
for i in range(1,10):
vec, counts, obj, prob, qc, res1, res2, bestBetas, bestGammas = QAOAoptimizeMaxCount(xxs,i,1000)
Nmax.append(i)
probamax.append(prob)
Optmax.append(obj)
Selected device: qasm_simulator with 1000 shots
Trying 1 layer
Generating inital random parameters beta and gamma
Starting with betas: [-0.7667554870727287]
Starting with gammas: [4.5194823454833575]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [-0.84534792] Best Gamma [4.50818364]
Now run the QAOA with the found parameters
----------------------------------------------
Optimum = 5.0 with probability = 0.161
----------------------------------------------
Depth: 15
Gate counts: OrderedDict([('rz', 7), ('cx', 6), ('h', 4), ('ry', 4), ('measure', 4), ('barrier', 3)])
Selected device: qasm_simulator with 1000 shots
Trying 2 layer
Generating inital random parameters beta and gamma
Starting with betas: [0.037155049742974455, 1.2405130243676892]
Starting with gammas: [-1.680187571835848, 4.85981031127122]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [0.03715505 1.24051302] Best Gamma [-1.76419695 4.85981031]
Now run the QAOA with the found parameters
----------------------------------------------
Optimum = 0.0 with probability = 0.465
----------------------------------------------
Depth: 28
Gate counts: OrderedDict([('rz', 14), ('cx', 12), ('ry', 8), ('barrier', 5), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 3 layer
Generating inital random parameters beta and gamma
Starting with betas: [-3.577795290290057, -1.5984094739294914, -2.1195462834196173]
Starting with gammas: [-3.914556116068772, 4.400492652220558, 4.94262406307908]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [-3.61754857 -1.61616958 -2.07244525] Best Gamma [-4.05591509 4.36382188 4.99754211]
Now run the QAOA with the found parameters
----------------------------------------------
Optimum = 0.0 with probability = 0.214
----------------------------------------------
Depth: 41
Gate counts: OrderedDict([('rz', 21), ('cx', 18), ('ry', 12), ('barrier', 7), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 4 layer
Generating inital random parameters beta and gamma
Starting with betas: [-5.437406140851462, 5.062048853057593, -1.5534827142613077, 3.271463467953126]
Starting with gammas: [-0.9971621017678665, 0.624268335250231, -6.10650184737705, -3.0482833068603226]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [-5.43740614 5.06204885 -1.55348271 3.27146347] Best Gamma [-0.9971621 0.62426834 -6.10650185 -3.04828331]
Now run the QAOA with the found parameters
----------------------------------------------
Optimum = 5.0 with probability = 0.323
----------------------------------------------
Depth: 54
Gate counts: OrderedDict([('rz', 28), ('cx', 24), ('ry', 16), ('barrier', 9), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 5 layer
Generating inital random parameters beta and gamma
Starting with betas: [4.566441888593328, 2.826431188375512, 2.8708372286428485, 5.715785858239107, 0.383098022160711]
Starting with gammas: [-1.0646764825398582, -6.02892611048706, -3.5653629096297172, -0.41168512015416336, -0.39524598184121906]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [4.56644189 2.82643119 2.87083723 5.71578586 0.38309802] Best Gamma [-1.09129339 -6.17964926 -3.56536291 -0.41168512 -0.39524598]
Now run the QAOA with the found parameters
----------------------------------------------
Optimum = 0.0 with probability = 0.155
----------------------------------------------
Depth: 67
Gate counts: OrderedDict([('rz', 35), ('cx', 30), ('ry', 20), ('barrier', 11), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 6 layer
Generating inital random parameters beta and gamma
Starting with betas: [-0.7501843620847204, -2.9499382215046586, 1.605129775372836, 2.6188301925297743, 0.27599249728486086, -3.883161968725844]
Starting with gammas: [5.033181403624903, -5.733415756279973, 5.184477665297733, 5.016778885420768, -2.1952952017001914, 3.693598592306362]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [-0.75018436 -2.94993822 1.60512978 2.61883019 0.28979212 -3.88316197] Best Gamma [ 5.0331814 -5.73341576 5.18447767 5.01677889 -2.25017758 3.78593856]
Now run the QAOA with the found parameters
----------------------------------------------
Optimum = 5.0 with probability = 0.15
----------------------------------------------
Depth: 80
Gate counts: OrderedDict([('rz', 42), ('cx', 36), ('ry', 24), ('barrier', 13), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 7 layer
Generating inital random parameters beta and gamma
Starting with betas: [-2.0917459230119944, -3.1734689500386053, 1.5016518745427865, -6.093750496286678, 4.894818319447703, -2.264229731331918, 1.0325039615140916]
Starting with gammas: [-0.9398240006580894, 4.191913784443285, -1.8538154238820734, 3.6824390737045487, -2.554012641717735, 6.26632401689767, -1.4365238980449018]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [-2.09174592 -3.17346895 1.50165187 -6.0937505 4.89481832 -2.26422973
1.03250396] Best Gamma [-0.9868152 4.19191378 -1.85381542 3.68243907 -2.55401264 6.26632402
-1.4365239 ]
Now run the QAOA with the found parameters
----------------------------------------------
Optimum = -99.0 with probability = 0.338
----------------------------------------------
Depth: 93
Gate counts: OrderedDict([('rz', 49), ('cx', 42), ('ry', 28), ('barrier', 15), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 8 layer
Generating inital random parameters beta and gamma
Starting with betas: [5.683139696815974, -5.74380033350016, 3.1156574360578304, 5.557578255495642, -4.271808793470049, -3.5908471150410737, -1.7190955836369533, 5.955729312522351]
Starting with gammas: [1.227473073329218, 2.6575062121551625, -0.30693752834871546, -3.5508491363649632, 6.039749238259079, 0.1672965832420692, 2.194498352012694, 4.274055218889998]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [ 5.6831397 -5.74380033 3.11565744 5.55757826 -4.27180879 -3.59084712
-1.71909558 5.95572931] Best Gamma [ 1.22747307 2.79038152 -0.30693753 -3.55084914 6.03974924 0.16729658
2.19449835 4.27405522]
Now run the QAOA with the found parameters
----------------------------------------------
Optimum = -99.0 with probability = 0.175
----------------------------------------------
Depth: 106
Gate counts: OrderedDict([('rz', 56), ('cx', 48), ('ry', 32), ('barrier', 17), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 9 layer
Generating inital random parameters beta and gamma
Starting with betas: [2.490943023336392, -5.514960516966169, -2.122216994538338, 2.0406473519879764, -3.376822891061745, 3.761596442840677, -4.901810758707899, -2.9723422174701084, 4.520494429584939]
Starting with gammas: [6.265816969617411, -5.7199617242016245, -3.8828176972673414, 5.33838872696537, 0.5069595192707865, -3.3879414364496485, 5.5485789277729225, 3.0433555434063386, -4.708802862571192]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [ 2.48116846 -5.50436362 -2.115895 2.03824622 -3.39548599 3.78238611
-4.89932895 -2.98876982 4.49800886] Best Gamma [ 6.30044702 -5.75157493 -3.88530266 5.36789305 0.50976139 -3.3882701
5.49832815 3.0467905 -4.85254764]
Now run the QAOA with the found parameters
----------------------------------------------
Optimum = -99.0 with probability = 0.217
----------------------------------------------
Depth: 119
Gate counts: OrderedDict([('rz', 63), ('cx', 54), ('ry', 36), ('barrier', 19), ('h', 4), ('measure', 4)])
In [45]:
bestenergyxxs
Out[45]:
-99.0
In [46]:
plt.axhline(y=bestenergyxxs, color='g', linestyle=':')
plt.ylabel("Cost")
plt.xlabel("Iterations N")
plt.title("Example set XXS with 4 qubits: Convergence vs iterations")
plt.plot(Nmax,Optmax,'r+')
Out[46]:
[<matplotlib.lines.Line2D at 0x7efd6d98c8e0>]
In [47]:
## approximate the expectation value with all measured bitstrings and take the average
#vec, counts, expectationValue, prob, qc, res1, res2, bestBetas, bestGammas, optimum = QAOAoptimizeExpectation(xxs,5,1000)
Nexp=[]
probaexp=[]
Optexp=[]
expVexp=[]
for i in range(1,10):
vec, counts, expectationValue, prob, qc, res1, res2, bestBetas, bestGammas, optimum = QAOAoptimizeExpectation(xxs,i,1000)
Nexp.append(i)
probaexp.append(prob)
expVexp.append(expectationValue)
Optexp.append(optimum)
Selected device: qasm_simulator with 1000 shots
Trying 1 layer
Generating inital random parameters beta and gamma
Starting with betas: [-2.2949524076404204]
Starting with gammas: [1.0057312668756282]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [-1.91825229] Best Gamma [1.3494441]
Now run the QAOA with the found parameters
----------------------------------------------
Expectation value = -45.916
----------------------------------------------
Optimum = 2.0 with probability = 0.154
----------------------------------------------
Depth: 15
Gate counts: OrderedDict([('rz', 7), ('cx', 6), ('h', 4), ('ry', 4), ('measure', 4), ('barrier', 3)])
Selected device: qasm_simulator with 1000 shots
Trying 2 layer
Generating inital random parameters beta and gamma
Starting with betas: [-3.4171331931550886, -5.968806991022075]
Starting with gammas: [5.563187289718794, 3.685419923888709]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [-3.20442738 -6.31699821] Best Gamma [6.08671872 3.76833311]
Now run the QAOA with the found parameters
----------------------------------------------
Expectation value = -47.0
----------------------------------------------
Optimum = 5.0 with probability = 0.094
----------------------------------------------
Depth: 28
Gate counts: OrderedDict([('rz', 14), ('cx', 12), ('ry', 8), ('barrier', 5), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 3 layer
Generating inital random parameters beta and gamma
Starting with betas: [5.587253081714104, 3.219839915022419, 4.794235418644313]
Starting with gammas: [0.055296185026889155, -1.0334799131639203, -0.12671593314843665]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [5.24531162 3.00579502 4.61636948] Best Gamma [ 0.05850503 -1.12092744 -0.14006672]
Now run the QAOA with the found parameters
----------------------------------------------
Expectation value = -47.798
----------------------------------------------
Optimum = -96.0 with probability = 0.106
----------------------------------------------
Depth: 41
Gate counts: OrderedDict([('rz', 21), ('cx', 18), ('ry', 12), ('barrier', 7), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 4 layer
Generating inital random parameters beta and gamma
Starting with betas: [2.33362566581156, -2.3149122870140415, -5.219349854025362, -0.5786531916704627]
Starting with gammas: [2.719067639388051, -1.4596955077863045, -4.772902309092033, -5.49686887617206]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [ 2.37347287 -2.2960179 -5.24222439 -0.59036493] Best Gamma [ 2.85685637 -1.52673015 -4.73190274 -5.26543614]
Now run the QAOA with the found parameters
----------------------------------------------
Expectation value = -50.239
----------------------------------------------
Optimum = 7.0 with probability = 0.156
----------------------------------------------
Depth: 54
Gate counts: OrderedDict([('rz', 28), ('cx', 24), ('ry', 16), ('barrier', 9), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 5 layer
Generating inital random parameters beta and gamma
Starting with betas: [4.688903196933991, 1.1118646992746122, -5.022928549221403, 4.05737167118977, -4.726261978440847]
Starting with gammas: [-5.598428498044197, 5.90780887552846, 3.115335165959136, 1.8256050343483086, -5.696589986720251]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [ 4.65460811 1.12068786 -5.17185383 4.08943002 -4.97337726] Best Gamma [-5.59566451 5.88350684 3.12934639 1.83013783 -5.69913852]
Now run the QAOA with the found parameters
----------------------------------------------
Expectation value = -53.619
----------------------------------------------
Optimum = -94.0 with probability = 0.224
----------------------------------------------
Depth: 67
Gate counts: OrderedDict([('rz', 35), ('cx', 30), ('ry', 20), ('barrier', 11), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 6 layer
Generating inital random parameters beta and gamma
Starting with betas: [4.672661823327026, 0.9793597768297149, 3.0187377941040463, 1.7873965930050133, -1.034638305383841, -5.956697842263989]
Starting with gammas: [2.1683174575731705, -0.07717023832717018, 0.7112019636403772, -1.032289144610913, 0.5600720055374264, -3.6751740194434603]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [ 4.71405564 0.99519478 3.05428607 1.8035227 -1.05137607 -5.99151041] Best Gamma [ 2.17427149 -0.08036919 0.68705965 -1.04199339 0.56292494 -3.67191739]
Now run the QAOA with the found parameters
----------------------------------------------
Expectation value = -59.69
----------------------------------------------
Optimum = 3.0 with probability = 0.178
----------------------------------------------
Depth: 80
Gate counts: OrderedDict([('rz', 42), ('cx', 36), ('ry', 24), ('barrier', 13), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 7 layer
Generating inital random parameters beta and gamma
Starting with betas: [5.188414879137735, -3.2869483131469304, -0.8289447212057253, 5.356146779893553, 2.782210346027684, -0.36642694991945035, -1.3321774046908175]
Starting with gammas: [-0.5538597890519679, 0.4974173040528207, 0.9818832548812928, 2.8128096425450053, 4.804032149684469, -1.9080563557298138, -2.5268024246564087]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [ 5.34533673 -3.29489068 -0.85608114 5.28372598 2.82084766 -0.36742846
-1.3535569 ] Best Gamma [-0.56131464 0.4947271 0.99623635 2.70302766 4.74996978 -1.93708471
-2.56154841]
Now run the QAOA with the found parameters
----------------------------------------------
Expectation value = -50.308
----------------------------------------------
Optimum = -99.0 with probability = 0.207
----------------------------------------------
Depth: 93
Gate counts: OrderedDict([('rz', 49), ('cx', 42), ('ry', 28), ('barrier', 15), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 8 layer
Generating inital random parameters beta and gamma
Starting with betas: [4.01500124960393, 1.16415940466998, 4.089047823659499, -5.93343103052138, -1.8469380605781351, -2.0159076802647062, 2.718697019314506, -3.2314371367614747]
Starting with gammas: [2.7597924543644616, -1.4243533701860596, 1.0458248459461164, 2.6980406293786636, 5.745201978092473, 0.2433993154765055, 0.23854202718096396, -0.1270785452009715]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [ 3.84682168 1.17644705 4.06336394 -5.96055566 -1.88109157 -2.04031013
2.7411552 -3.24323639] Best Gamma [ 2.76551159 -1.43443028 1.01700237 2.6766355 5.95220755 0.25639667
0.24108732 -0.1270501 ]
Now run the QAOA with the found parameters
----------------------------------------------
Expectation value = -52.167
----------------------------------------------
Optimum = -96.0 with probability = 0.251
----------------------------------------------
Depth: 106
Gate counts: OrderedDict([('rz', 56), ('cx', 48), ('ry', 32), ('barrier', 17), ('h', 4), ('measure', 4)])
Selected device: qasm_simulator with 1000 shots
Trying 9 layer
Generating inital random parameters beta and gamma
Starting with betas: [0.1288389842307769, -1.1192571688090043, 0.46822355100712354, 3.530310658430258, 3.9888689531656176, 5.634427873008759, -4.369177762993817, 1.3117427927014749, -1.0866371079569976]
Starting with gammas: [-0.5034447667448658, 3.1235259440930285, -5.09156561420877, -0.3938870372962331, 5.893412831059214, -5.336840667694324, 3.853060250975613, 3.4682243920611846, 4.247078231594186]
Optimize FIRST round with random initialisation
Optimize SECOND round with the found initialization
Best Beta [ 0.12893305 -1.12701925 0.47559974 3.54636413 4.00736317 5.6526842
-4.48477679 1.31487418 -1.09109445] Best Gamma [-0.50433483 3.11916511 -5.06507953 -0.39463995 5.87406004 -5.4550044
3.85055358 3.55086832 4.24754685]
Now run the QAOA with the found parameters
----------------------------------------------
Expectation value = -68.636
----------------------------------------------
Optimum = -97.0 with probability = 0.244
----------------------------------------------
Depth: 119
Gate counts: OrderedDict([('rz', 63), ('cx', 54), ('ry', 36), ('barrier', 19), ('h', 4), ('measure', 4)])
In [49]:
plt.axhline(y=bestenergyxxs, color='g', linestyle=':')
plt.ylabel("Cost")
plt.xlabel("Iterations N")
plt.title("Example set XXS with 4 qubits: Convergence vs iterations")
plt.plot(Nexp,expVexp,'r+')
plt.plot(Nexp,Optexp,'b*')
Out[49]:
[<matplotlib.lines.Line2D at 0x7efd6d7292e0>]
In [58]:
plt.axhline(y=bestenergyxxs, color='g', linestyle=':')
plt.plot(Nexp,expVexp,'r+',label="Expectation approx")
plt.plot(Nexp,Optexp,'g+',label="Expectation approx - optimum")
plt.plot(Nmax,Optmax,'b+',label="Max count - optimum")
plt.legend()
plt.ylabel("Cost")
plt.xlabel("Iterations N")
plt.title("Example set XXS with 4 qubits: Convergence vs iterations")
Out[58]:
Text(0.5, 1.0, 'Example set XXS with 4 qubits: Convergence vs iterations')
In [ ]:
In [ ]:
In [ ]: