notebooks/merit_order.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1eb92fd8",
   "metadata": {},
   "source": [
    "# Merit order\n",
    "\n",
    "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.\n",
    "\n",
    "In the easiest way the problem is equivalent to the knapsack problem:\n",
    "https://en.wikipedia.org/wiki/Knapsack_problem\n",
    "\n",
    "More information here: \n",
    "https://ercim-news.ercim.eu/en128/special/energy-economics-fundamental-modelling-with-quantum-algorithms\n",
    "\n",
    "1. Start with an easy example\n",
    "2. introduce resolution\n",
    "3. introduce slack variable to formulate unequalitites\n",
    "4. minimum/maximum up/down unit commitment problem including satisfiablitiy formulation to formulate QUBO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5197507d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import pygrnd and other libraries needed\n",
    "# we build on top of the open source framework qiskit (qiskit.org)\n",
    "import pygrnd\n",
    "\n",
    "from pygrnd.qc.helper import *\n",
    "from pygrnd.qc.brm import brm\n",
    "from pygrnd.qc.brm_oracle import brmoracle\n",
    "from pygrnd.qc.QAE import qae\n",
    "\n",
    "from pygrnd.optimize.sat_ucp import *\n",
    "from pygrnd.optimize.meritorder import *\n",
    "\n",
    "from pygrnd.optimize.bruteforce import *\n",
    "from pygrnd.optimize.MonteCarloSolver import *\n",
    "from pygrnd.optimize.qaoa import *\n",
    "\n",
    "#from qiskit import execute\n",
    "#from qiskit import Aer\n",
    "\n",
    "from math import pi\n",
    "import math\n",
    "import cmath\n",
    "import random\n",
    "from qiskit.visualization import plot_histogram\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm\n",
    "import networkx as nx\n",
    "from IPython.display import Image"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc1dbde7",
   "metadata": {},
   "source": [
    "# Motivation QUBO (Quadratic Unconstrained Binary Optimization)\n",
    "\n",
    "- quadratic unconstrained binary optimization\n",
    "- Minimize/maximize $\\langle x | Q | x \\rangle = \\sum_{ij} x_i Q_{ij} x_j$\n",
    "- Binary variables $x_i \\in \\{0,1\\} \\iff x_i = x_i^2$ \n",
    "- entries in symmetric matrix Q are real numbers\n",
    "- solves combinatorial problems\n",
    "- constraints need to be encoded as penalty terms\n",
    "\n",
    "- Constraints for knapsack encoded in parameters $Q_{ij}$\n",
    "- $Q = Q_{cost} + Q_{constraint}$\n",
    "\n",
    "- the optimal solution is hard to find\n",
    "- QUBOs can be solved by different Ising solver, annealer, quantum annealer, quantum simulator\n",
    "- QUBO stands equivalent with Ising spin model\n",
    "\n",
    "- Cost are encoded on the main diagonal $c_i$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19c8b3a4",
   "metadata": {},
   "source": [
    "# Knapsack Problem as QUBO\n",
    "\n",
    "- Minimize -$\\sum_i x_i v_i$ s.t. $\\sum_i x_i w_i = W$\n",
    "- Binary variables $x_i\\in\\{0,1\\}$, i.e., take an element or do not take it at all\n",
    "- Maximize sum of selected values $v_i \\in \\mathbb{R}_{\\geq 0}$\n",
    "- 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$\n",
    "- Introduce slack variable to obtain equality $\\sum_i x_i w_i +s = W$ with $s\\in \\mathbb{R}$\n",
    "\n",
    "# QUBO constraint for demand\n",
    "\n",
    "- D = W = demand\n",
    "\n",
    "- $ w_i = weights_i$ = $maxgen_i$ of each unit i\n",
    "\n",
    "- $ \\sum_i w_i = D $\n",
    "\n",
    "- $\\implies (\\sum_i w_i - D)^2 = 0 $\n",
    "\n",
    "- $ \\implies \\sum_i w_i \\sum w_j - 2*D \\sum_i w_i + D^2 = 0$\n",
    "\n",
    "   - for $i = j: \\quad (\\sum_i w_i )^2 - 2*D*\\sum_i w_i $ # main diagonal matrix elements\n",
    "\n",
    "   - for $i \\neq j: \\quad \\sum_i w_i * \\sum_i w_j $\n",
    "\n",
    "- we ignore constant (\"offset\", $D^2$) - needs to be added/subtracted from the solution\n",
    "\n",
    "\n",
    "# Rewrite Constraints To Obtain QUBO Formulation\n",
    "\n",
    "- $\\sum_i x_i w_i =W$ can be written as $\\left( \\sum_i x_i w_i -W \\right)^2$\n",
    "- Solve QUBO -$\\sum_i v_i x_i + P \\left( \\sum_i x_i w_i -W \\right)^2$\n",
    "- Find appropriate penalty factor $P$\n",
    "- Use $x^2_i = x_i$ for binary variables $x_i$\n",
    "\n",
    "# Knapsack With Resolution\n",
    "\n",
    "- $0/1$ value for $x_i$ should be more fine-grained\n",
    "- Solution: Split $x_i$ in several parts and represent it as $0...0$ to $1...1$ with $\\frac{1}{2^{b}-1}$ per part\n",
    "- Example: $b=3$ leads to $\\frac{0}{7}, \\frac{1}{7}, \\ldots, \\frac{7}{7}$\n",
    "- 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)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "405f9a6f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-12.   3.   4.   5.]\n",
      " [  3. -30.  12.  15.]\n",
      " [  4.  12. -36.  20.]\n",
      " [  5.  15.  20. -40.]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "0it [00:00, ?it/s]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "16it [00:00, 30545.68it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-42.0, matrix([[0, 1, 1, 0]]))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Weight constraint violated. P too small.\n",
    "M=QUBO_knapsack([1,3,4,5],[1,3,4,5],7,1)\n",
    "print(M)\n",
    "solver(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f5dadefd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-25.   6.   8.  10.]\n",
      " [  6. -63.  24.  30.]\n",
      " [  8.  24. -76.  40.]\n",
      " [ 10.  30.  40. -85.]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "16it [00:00, 24546.04it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-91.0, matrix([[0, 1, 1, 0]]))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Weight constraint satisfied.\n",
    "M=QUBO_knapsack([1,3,4,5],[1,3,4,5],7,2)\n",
    "print(M)\n",
    "solver(M)\n",
    "\n",
    "## penalty finetuning is art"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d6989c17",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-2.66666667  0.66666667  0.66666667  1.33333333]\n",
      " [ 0.66666667 -4.66666667  1.33333333  2.66666667]\n",
      " [ 0.66666667  1.33333333 -4.66666667  2.66666667]\n",
      " [ 1.33333333  2.66666667  2.66666667 -6.66666667]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "16it [00:00, 29615.56it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1 1 0]]\n",
      "1.6666666666666665\n",
      "1.6666666666666665\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# Two solutions are possible: 10/01 and 11/10. Both should lead to 5/3 as result.\n",
    "# Penalty 3 should be necessary, 2 leads to wrong results.\n",
    "\n",
    "values_split=splitParameters([1,2],2)\n",
    "weights_split=splitParameters([1,2],2)\n",
    "\n",
    "M=QUBO_knapsack(values_split,weights_split,5/3,3)\n",
    "print(M)\n",
    "res,v=solver(M)\n",
    "print(v)\n",
    "res=recombineSolution(values_split, [1,0,0,1], 2)\n",
    "print(sum(res))\n",
    "res=recombineSolution(values_split, [1,1,1,0], 2)\n",
    "print(sum(res))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14848470",
   "metadata": {},
   "source": [
    "# Slack Variable for Upper Weight Bound\n",
    "\n",
    "- 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}$\n",
    "- 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\n",
    "- Use the normal fractional method from above and just extend slack variable with weight $W$ and value $0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2a7504e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "512it [00:00, 84699.99it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 0, 0, 0, 0, 0, 1, 1, 1] -90.0\n",
      "recombined fractions: [0.0, 0.0, 1.0]\n",
      "recombined values: [0.0, 0.0, 0.0]\n",
      "recombined weights: [0.0, 0.0, 3.0]\n",
      "total/slack value: 0.0 0.0\n",
      "real/slack weight: 0.0 3.0\n",
      "total/demanded/diff weight: 3.0 3 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "QUBO_knapsack_slack_resolution([1,2],[1,1],3,3,10,True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f9dc650c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "512it [00:00, 77117.24it/s]\n",
      "512it [00:00, 84251.39it/s]\n",
      "512it [00:00, 86795.07it/s]\n",
      "512it [00:00, 95554.14it/s]\n",
      "512it [00:00, 94916.40it/s]\n",
      "512it [00:00, 86826.65it/s]\n",
      "512it [00:00, 87356.45it/s]\n",
      "512it [00:00, 84783.59it/s]\n",
      "512it [00:00, 82671.84it/s]\n",
      "512it [00:00, 76852.29it/s]\n",
      "512it [00:00, 80013.55it/s]\n",
      "512it [00:00, 80150.92it/s]\n",
      "512it [00:00, 89418.87it/s]\n",
      "512it [00:00, 83193.88it/s]\n",
      "512it [00:00, 80478.33it/s]\n",
      "512it [00:00, 103021.52it/s]\n",
      "512it [00:00, 82792.95it/s]\n",
      "512it [00:00, 84473.43it/s]\n",
      "512it [00:00, 90891.08it/s]\n",
      "512it [00:00, 92452.37it/s]\n",
      "512it [00:00, 78167.06it/s]\n",
      "512it [00:00, 94536.17it/s]\n",
      "512it [00:00, 90431.79it/s]\n",
      "512it [00:00, 79627.86it/s]\n",
      "512it [00:00, 89713.98it/s]\n",
      "512it [00:00, 91788.50it/s]\n",
      "512it [00:00, 79787.61it/s]\n",
      "512it [00:00, 73173.08it/s]\n",
      "512it [00:00, 61619.00it/s]\n",
      "512it [00:00, 66640.30it/s]\n",
      "512it [00:00, 80403.00it/s]\n",
      "512it [00:00, 86074.94it/s]\n",
      "512it [00:00, 84699.99it/s]\n",
      "512it [00:00, 72766.46it/s]\n",
      "512it [00:00, 87792.14it/s]\n",
      "512it [00:00, 95164.57it/s]\n",
      "512it [00:00, 70745.63it/s]\n",
      "512it [00:00, 75205.17it/s]\n",
      "512it [00:00, 75382.04it/s]\n",
      "512it [00:00, 86679.46it/s]\n",
      "512it [00:00, 84586.56it/s]\n",
      "512it [00:00, 79456.97it/s]\n",
      "512it [00:00, 80040.39it/s]\n",
      "512it [00:00, 71807.79it/s]\n",
      "512it [00:00, 74963.65it/s]\n",
      "512it [00:00, 80502.46it/s]\n",
      "512it [00:00, 80013.55it/s]\n",
      "512it [00:00, 77917.48it/s]\n",
      "512it [00:00, 93825.74it/s]\n",
      "512it [00:00, 104632.80it/s]\n",
      "512it [00:00, 79009.70it/s]\n",
      "512it [00:00, 78312.44it/s]\n",
      "512it [00:00, 84433.58it/s]\n",
      "512it [00:00, 83036.26it/s]\n",
      "512it [00:00, 79196.18it/s]\n",
      "512it [00:00, 91250.26it/s]\n",
      "512it [00:00, 81208.73it/s]\n",
      "512it [00:00, 88319.29it/s]\n",
      "512it [00:00, 79965.88it/s]\n",
      "512it [00:00, 80526.61it/s]\n",
      "512it [00:00, 81824.49it/s]\n",
      "512it [00:00, 78538.70it/s]\n",
      "512it [00:00, 79924.21it/s]\n",
      "512it [00:00, 77988.22it/s]\n",
      "512it [00:00, 80013.55it/s]\n",
      "512it [00:00, 80988.22it/s]\n",
      "512it [00:00, 81153.49it/s]\n",
      "512it [00:00, 77838.40it/s]\n",
      "512it [00:00, 84446.86it/s]\n",
      "512it [00:00, 82633.66it/s]\n",
      "512it [00:00, 94552.82it/s]\n",
      "512it [00:00, 86164.73it/s]\n",
      "512it [00:00, 79950.99it/s]\n",
      "512it [00:00, 79047.51it/s]\n",
      "512it [00:00, 78395.34it/s]\n",
      "512it [00:00, 81116.70it/s]\n",
      "512it [00:00, 93503.01it/s]\n",
      "512it [00:00, 83306.84it/s]\n",
      "512it [00:00, 78130.09it/s]\n",
      "512it [00:00, 81408.83it/s]\n",
      "512it [00:00, 84393.76it/s]\n",
      "512it [00:00, 81896.26it/s]\n",
      "512it [00:00, 86837.19it/s]\n",
      "512it [00:00, 76981.78it/s]\n",
      "512it [00:00, 81696.86it/s]\n",
      "512it [00:00, 82780.19it/s]\n",
      "512it [00:00, 81257.89it/s]\n",
      "512it [00:00, 78132.93it/s]\n",
      "512it [00:00, 77225.39it/s]\n",
      "512it [00:00, 75751.65it/s]\n",
      "512it [00:00, 79554.11it/s]\n",
      "512it [00:00, 78241.11it/s]\n",
      "512it [00:00, 75358.24it/s]\n",
      "512it [00:00, 77610.54it/s]\n",
      "512it [00:00, 80796.25it/s]\n",
      "512it [00:00, 84390.44it/s]\n",
      "512it [00:00, 78118.72it/s]\n",
      "512it [00:00, 78492.77it/s]\n",
      "512it [00:00, 79799.47it/s]\n"
     ]
    }
   ],
   "source": [
    "# With penalty 7.1 we start to see good results\n",
    "for i in range(1,100):\n",
    "    buffer=i*0.1\n",
    "    res=QUBO_knapsack_slack_resolution([1,2],[1,2],5/7,3,buffer)\n",
    "    if abs(res-5/7)<0.1:\n",
    "        print(i,buffer,abs(res-5/7))\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "139c3d02",
   "metadata": {},
   "source": [
    "# example 1\n",
    "\n",
    "We only consider 3 units and demand = 20 = D\n",
    "\n",
    "One time step\n",
    "\n",
    "No ramp up/down costs, no min/max up/down times\n",
    "\n",
    "| unit | maxgen (MW) | cost (€) |\n",
    "|------|-------------|----------|\n",
    "|  1   |     15      |     2    |\n",
    "|  2   |      5      |    10    |\n",
    "|  3   |      5      |    10    |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "77485122",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-373.   75.   75.]\n",
      " [  75. -165.   25.]\n",
      " [  75.   25. -165.]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "8it [00:00, 42581.77it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "offset =  400\n",
      "Solution =  12.0 | Vector =  [[1 0 1]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "#small example\n",
    "demand = 20\n",
    "penalty = 1\n",
    "\n",
    "M=QUBO_knapsack([2,10,10],[15,5,5],demand,penalty)\n",
    "print(M)\n",
    "res,vec=solver(M)\n",
    "\n",
    "print(\"offset = \",demand**2)\n",
    "print(\"Solution = \",res+demand**2,\"| Vector = \",vec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9c3105f6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6000, 22000, 34000, 26600, 23100]\n"
     ]
    }
   ],
   "source": [
    "# Here a bigger example\n",
    "# Modelling of up parameter only - no time dependecies\n",
    "\n",
    "maxgen=[1000,800,1000,700,350]\n",
    "cost=[6000,22000,34000,26600,23100]\n",
    "demand=1000\n",
    "\n",
    "print(cost)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0e74ab9b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Demand =  1000\n",
      "Penalty =  2\n",
      "['Nuc', 'Lign', 'Coal', 'CC', 'GT']\n",
      "[[-1994000.  1600000.  2000000.  1400000.   700000.]\n",
      " [ 1600000. -1898000.  1600000.  1120000.   560000.]\n",
      " [ 2000000.  1600000. -1966000.  1400000.   700000.]\n",
      " [ 1400000.  1120000.  1400000. -1793400.   490000.]\n",
      " [  700000.   560000.   700000.   490000. -1131900.]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "32it [00:00, 64652.08it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-1994000.0, matrix([[1, 0, 0, 0, 0]]))\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "32it [00:00, 45130.37it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "offset =  2000000\n",
      "Solution =  6000.0 | Vector =  [[1 0 0 0 0]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "penalty=2\n",
    "offset=penalty*((demand)**2)\n",
    "\n",
    "print(\"Demand = \",demand)\n",
    "print(\"Penalty = \",penalty)\n",
    "\n",
    "units=[\"Nuc\",\"Lign\",\"Coal\",\"CC\",\"GT\"]\n",
    "print(units)\n",
    "\n",
    "M=QUBO_knapsack(cost,maxgen,demand,penalty)\n",
    "print(M)\n",
    "print(solver(M))\n",
    "\n",
    "res,vec=solver(M)\n",
    "\n",
    "print(\"offset = \",offset)\n",
    "print(\"Solution = \",res+offset,\"| Vector = \",vec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6b6e516",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "q",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  },
  "vscode": {
   "interpreter": {
    "hash": "59ee0f4e51d5fd22062f9f81a021fbeb1bb3e4bf9568561b7ac3ccf818c074f1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}