diff --git a/examples/quantum_annealing/price_optimization/QUBO-Pricing.ipynb b/examples/quantum_annealing/price_optimization/QUBO-Pricing.ipynb
new file mode 100644
index 000000000..cfd35fd0a
--- /dev/null
+++ b/examples/quantum_annealing/price_optimization/QUBO-Pricing.ipynb
@@ -0,0 +1,2855 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "2c5219f3",
+ "metadata": {},
+ "source": [
+ "# Using quantum annealing on Amazon Braket for price optimization\n",
+ "\n",
+ "Combinatorial Optimization is one of the most important fields in optimization. Practical applications can be found in virtually every industry. Prominent examples include supply chain optimization in transport and logistics, portfolio management in finance, and the optimization of clinical trials in healthcare, among many others. It is also one of the most active research topics in operation research and computer science. However, many practical combinatorial optimization problems are NP-hard and require massive computation costs to find solution of good quality. \n",
+ "\n",
+ "In this notebook, we demonstrate how a quantum annealer on Amazon Braket can be used for price optimization taking into consideration the trade-off between maximizing revenue and minimizing risk. We show how to formulate this problem as a quadratic unconstrained binary optimization problem (QUBO) and use the D-Wave Advantage quantum annealer on Amazon Braket to find close-to-optimal solutions. Overall, this notebook demonstrates that customers can easily leverage quantum computing through Amazon Braket to solve difficult combinatorial optimization challenges in their daily decision-making process."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "10f29dc7",
+ "metadata": {},
+ "source": [
+ "# Table of contents\n",
+ "### 1. Demand model\n",
+ " 1. Create demand dataset\n",
+ " 2. Fit demand model through linear regression\n",
+ "### 2. Price optimization with QUBO\n",
+ " 1. Construct revenue objective\n",
+ " 2. Add penalty for prediction uncertainty\n",
+ " 3. Add equality constraints\n",
+ "### 3. Solve QUBO with Amazon Braket\n",
+ " 1. Evaluate results\n",
+ "### 4. Trade-off between revenue expectation and prediction uncertainty\n",
+ "### 5. Conclusion\n",
+ "### 6. Literature review\n",
+ "-----------------"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "aed8d388",
+ "metadata": {},
+ "source": [
+ "We start by importing Amazon braket SDK and Dwave SDK to be used for quantum annealing in this use case. In addition, we import sklearn to fit the demand estimation model."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "1f31ea86",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Requirement already satisfied: sklearn in /home/ec2-user/anaconda3/envs/Braket/lib/python3.7/site-packages (0.0)\n",
+ "Requirement already satisfied: scikit-learn in /home/ec2-user/anaconda3/envs/Braket/lib/python3.7/site-packages (from sklearn) (1.0.1)\n",
+ "Requirement already satisfied: scipy>=1.1.0 in /home/ec2-user/anaconda3/envs/Braket/lib/python3.7/site-packages (from scikit-learn->sklearn) (1.5.2)\n",
+ "Requirement already satisfied: numpy>=1.14.6 in /home/ec2-user/anaconda3/envs/Braket/lib/python3.7/site-packages (from scikit-learn->sklearn) (1.19.2)\n",
+ "Requirement already satisfied: threadpoolctl>=2.0.0 in /home/ec2-user/anaconda3/envs/Braket/lib/python3.7/site-packages (from scikit-learn->sklearn) (3.0.0)\n",
+ "Requirement already satisfied: joblib>=0.11 in /home/ec2-user/anaconda3/envs/Braket/lib/python3.7/site-packages (from scikit-learn->sklearn) (1.1.0)\n"
+ ]
+ }
+ ],
+ "source": [
+ "!pip install sklearn"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "6718da7a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from braket.aws import AwsDevice\n",
+ "from braket.ocean_plugin import BraketSampler, BraketDWaveSampler\n",
+ "from pyqubo import Binary\n",
+ "# magic word for producing visualizations in notebook\n",
+ "%matplotlib inline\n",
+ "import time\n",
+ "import reprlib\n",
+ "from collections import defaultdict\n",
+ "from itertools import combinations\n",
+ "import math\n",
+ "import networkx as nx\n",
+ "import dwave_networkx as dnx\n",
+ "import minorminer\n",
+ "import dimod\n",
+ "from dimod.binary_quadratic_model import BinaryQuadraticModel\n",
+ "from dwave.system.composites import EmbeddingComposite\n",
+ "import numpy as np\n",
+ "from sklearn.linear_model import LinearRegression\n",
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "from qubo_dynamic_pricing import optimize\n",
+ "\n",
+ "np.random.seed(0)\n",
+ "repr_compact = reprlib.Repr()\n",
+ "repr_compact.maxother=200"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "61be0cc9",
+ "metadata": {},
+ "source": [
+ "__NOTE__: Enter your S3 bucket and key below. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "ce15ec79",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Enter the S3 bucket you created during onboarding in the code below\n",
+ "my_bucket = \"amazon-braket-Your-Bucket-Name\" # the name of the bucket\n",
+ "my_prefix = \"Your-Folder-Name\" # the name of the folder in the bucket\n",
+ "s3_folder = (my_bucket, my_prefix)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "c2fc5e7e",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Device: Device('name': Advantage_system4.1, 'arn': arn:aws:braket:::device/qpu/d-wave/Advantage_system4)\n"
+ ]
+ }
+ ],
+ "source": [
+ "# session and device\n",
+ "device = AwsDevice(\"arn:aws:braket:::device/qpu/d-wave/Advantage_system4\")\n",
+ "print('Device:', device)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "173be395",
+ "metadata": {},
+ "source": [
+ "# 1. Demand model"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5c05da56",
+ "metadata": {},
+ "source": [
+ "The usual goal of price optimization is to maximize the revenue in the next certain period with respect to the corresponding prices, where the revenue can usually be represented by a function of demand and price: \n",
+ "\n",
+ "$$\\max_{{\\mathbf{p}}} R(\\mathbf{p}) = \\sum_{t=T}^{T+n-1}d_tp_t$$\n",
+ "\n",
+ "where $p_t$ and $d_t$ are the price and demand at day $t$.\n",
+ "\n",
+ "In order to do price optimization, we need to create a demand model to estimate demand by a function of price: $d_t=f(\\mathbf{p}_t)$, where $\\mathbf{p}$ is a vector of latest prices. In such a way, the revenue objective can be transformed as a function solely depending on price to do optimization $R=\\sum_{t=T}^{T+n-1}f(\\mathbf{p_t})p_t$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e7de211f",
+ "metadata": {},
+ "source": [
+ "## 1.1 Create training dataset"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "21fe8119",
+ "metadata": {},
+ "source": [
+ "For simplicity, we create a dummy dataset on historical demand and price observations so that a demand estimation model can be trained. \n",
+ "\n",
+ "Here we assume the true demand at day $d_t$ follows a linear model of the latest $n$ days' prices $p_i, i\\in [t-n+1, t-n+2, \\ldots, t]$ plus a noise following normal distribution $\\varepsilon\\sim N(0,\\sigma) $.\n",
+ "\n",
+ "Thus the demand at day $t$ would be\n",
+ "$$d_t = \\sum_{i=t-n+1}^{t}a_{i+n-1-t}p_{i} + b + \\varepsilon $$\n",
+ "where $a_j, j\\in [0, 1, \\ldots, n-1]$ is the elasiticity between the price at day $j+t-n+1$ and demand at day $t$, and\n",
+ "$b$ is a constant.\n",
+ "\n",
+ "We use these assumptions to create our dummy training dataset."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "438e6d85",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "dataset x first 10 samples:\n",
+ "[[ 8 10 10 8 8 10 8]\n",
+ " [10 10 8 8 10 8 13]\n",
+ " [10 8 8 10 8 13 19]\n",
+ " [ 8 8 10 8 13 19 8]\n",
+ " [ 8 10 8 13 19 8 10]\n",
+ " [10 8 13 19 8 10 8]\n",
+ " [ 8 13 19 8 10 8 8]\n",
+ " [13 19 8 10 8 8 16]\n",
+ " [19 8 10 8 8 16 5]\n",
+ " [ 8 10 8 8 16 5 5]]\n",
+ "dataset y first 10 samples:\n",
+ "[176.60983362005817, 166.3135122384355, 122.40886555130055, 139.82753849065216, 136.5494478345374, 139.08672451958816, 171.1141344886614, 149.02456906069798, 153.67576697932017, 185.92865845385478]\n"
+ ]
+ }
+ ],
+ "source": [
+ "a0=[-0.3, -0.5, -1.0, -2.0, -3.0, -3.3, -3.5] # elasticities in linear demand model\n",
+ "b0=300 # constants in the linear demand model \n",
+ "sigma=10. # the standard deviation of the noise\n",
+ "price_levels=[5, 8, 10, 12, 13, 16, 19] # predefined possible prices for each day\n",
+ "probabilities=[0.3, 0.3, 0.2, 0.05, 0.05, 0.05, 0.05] # the probabilities of taking a price choice at a day\n",
+ "n_samples=1000 # the number of sample we want to create\n",
+ "\n",
+ "def create_data_point(p, a, b, sigma):\n",
+ " \"\"\"\n",
+ " estimate the demand\n",
+ " :param p: np.array, (T,)\n",
+ " :param a: np.array, (T,)\n",
+ " :param b: float, the constants\n",
+ " :return: v\n",
+ " \"\"\"\n",
+ " v = np.dot(p,a) + b + np.random.normal(loc=0.0, scale=sigma)\n",
+ " return v\n",
+ "\n",
+ "def create_dataset(a, b, N, price_levels, probabilities, sigma):\n",
+ " \"\"\"\n",
+ " create a dataset for training the demand model\n",
+ " :param a: np.array, (T,)\n",
+ " :param b: float\n",
+ " :param N: int number of samples\n",
+ " :param price_levels: list, price levels\n",
+ " :param probabilities: list, probabilities distribution of the prices\n",
+ " :return:\n",
+ " \"\"\"\n",
+ " t = len(a)\n",
+ " prices = np.random.choice(price_levels, N+t-1, p=probabilities)\n",
+ " data_x = []\n",
+ " data_y = []\n",
+ " for i in range(N):\n",
+ " p = prices[i:i+t]\n",
+ " v = create_data_point(p, a, b, sigma)\n",
+ " data_x.append(\n",
+ " np.expand_dims(p, axis=0)\n",
+ " )\n",
+ " data_y.append(v)\n",
+ "\n",
+ " data_x = np.concatenate(data_x,axis=0)\n",
+ " return data_x, data_y\n",
+ "\n",
+ "data_x, data_y = create_dataset(a0, b0, n_samples, price_levels, probabilities, sigma)\n",
+ "\n",
+ "print(\"dataset x first 10 samples:\")\n",
+ "print(data_x[:10])\n",
+ "print(\"dataset y first 10 samples:\")\n",
+ "print(data_y[:10])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8f68105b",
+ "metadata": {},
+ "source": [
+ "## 1.2. Fit the demand model"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "31302cb4",
+ "metadata": {},
+ "source": [
+ "We use `sklearn` to fit a linear demand model to the training set. This fitted linear demand model will be used for the following price optimization problem."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "08cdb86d",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "fitted elasticities: [-0.3779074420638822, -0.5723441150118157, -0.8832936349398142, -2.044653528191526, -2.9861902495686956, -3.2337479564086022, -3.6095730899274248]\n",
+ "fitted constant: 301.2623906306577\n"
+ ]
+ }
+ ],
+ "source": [
+ "def linear_regression(data_x, data_y):\n",
+ " reg = LinearRegression().fit(data_x, data_y)\n",
+ " a = reg.coef_\n",
+ " b = reg.intercept_\n",
+ " return a, b\n",
+ "\n",
+ "\n",
+ "a, b = linear_regression(data_x, data_y)\n",
+ "a = [i for i in a]\n",
+ "b = b\n",
+ "\n",
+ "print(f'fitted elasticities: {a}')\n",
+ "print(f'fitted constant: {b}')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "694b5278",
+ "metadata": {},
+ "source": [
+ "# 2. Price optimization with QUBO"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c67da567",
+ "metadata": {},
+ "source": [
+ "In this notebook, we leverage a general-purpose mathematical framework: Quadratic Unconstrained Binary Optimization (QUBO), which has been recently proposed as an effective framework that is able to generalize across a large variety of combinatorial optimization use cases in industries [2, 3]. We take this price optimization as an example and showcase how to formulate it into QUBO framework considering prediction uncertainties as well as equality constraints."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "2f6135bd",
+ "metadata": {},
+ "source": [
+ "The formal definition of QUBO is expressed as an optimization problem, Minimize:\n",
+ "$$H=x^TQx$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "277ac4cc",
+ "metadata": {},
+ "source": [
+ "where $x$ is a vector of binary variables representing the decision variables of the combinatorial optimization. $Q$ is a squared matrix of constants inherently encoding the context and constraints of the optimization problems. $Q$ is commonly assumed to be symmetric or upper triangular without loss of generality. The goal is to find the optimal values of the binary decision variables $x$ that minimize the objective $H$.\n",
+ "\n",
+ "As shown above, the format of QUBO is strikingly simple, providing a general-purpose framework for a large class of combinatorial optimization problems. Below we show how to transform the objective function together with the constraints into this QUBO format. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b4ce1392",
+ "metadata": {},
+ "source": [
+ "## 2.1 Construct objective function"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4e5a61db",
+ "metadata": {},
+ "source": [
+ "Assuming tomorrow is day $T$, then the goal of our price optimization is to find the optimal prices $p_t, t\\in[T, T+1, \\ldots, T+n-1]$ in the next $n$ days in order to maximize the revenue in next $n$ days. Based on the demand model we showed above where the demand at a day is correlated with the latest $n$ days' price, we can construct objective function as below:\n",
+ "\n",
+ "$$\\max_{{\\mathbf{p}}}R(\\mathbf{p})=\\sum_{t=T}^{T+n-1}d_tp_t$$\n",
+ "\n",
+ "where $d_t$ is represented by our fitted model above: $d_t = \\sum_{i=t-n+1}^{t}a_{i+n-1-t}p_{i} + b $\n",
+ "\n",
+ "Therefore, we need at least the most recent n days historical prices. Here we randomly choose a historical range of $n$ days' prices as the most recent price data. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "8ab0bc6b",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "[10, 8, 8, 5, 5, 8, 8]"
+ ]
+ },
+ "execution_count": 7,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "p_data = data_x[np.random.randint(0,n_samples)].tolist()\n",
+ "\n",
+ "p_data"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8ced72f4",
+ "metadata": {},
+ "source": [
+ "Now, let's construct the objective function, i.e. the total revenue $R$, for the next $n$ days.\n",
+ "\n",
+ "First, we need to define the price variable $p_t$. According to certain business requirements, the possible prices are usually defined to be a set of discrete price values (e.g. \\\\$5.99, \\\\$9.99, \\\\$14.99). This means the price at each day only has a set of fixed available $m$ options $c_k, k \\in [0, 1, \\ldots, m-1]$, we can represent the price $p_t$ at day $t$ as: $$p_t=\\sum_{k=0}^{m-1}c_{k}x_{t,k}, \\textrm{ subject to } \\sum_{k=0}^{m-1}x_{t,k}=1$$ \n",
+ "\n",
+ "where $x_{t,k}$ is a binary variable with $1$ meaning the price at day $t$ takes option $c_k$, and $0$ meaning otherwise. Therefore, we can use binary variables $x_{t,k}$ to represent the price variables $p_t$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "a668a23a",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "[10,\n",
+ " 8,\n",
+ " 8,\n",
+ " 5,\n",
+ " 5,\n",
+ " 8,\n",
+ " 8,\n",
+ " (Binary(X_006)*Num(19.000000)+Binary(X_005)*Num(16.000000)+Binary(X_004)*Num(13.000000)+Binary(X_003)*Num(12.000000)+Binary(X_002)*Num(10.000000)+Binary(X_000)*Num(5.000000)+Binary(X_001)*Num(8.000000)),\n",
+ " (Binary(X_013)*Num(19.000000)+Binary(X_012)*Num(16.000000)+Binary(X_011)*Num(13.000000)+Binary(X_010)*Num(12.000000)+Binary(X_009)*Num(10.000000)+Binary(X_007)*Num(5.000000)+Binary(X_008)*Num(8.000000)),\n",
+ " (Binary(X_020)*Num(19.000000)+Binary(X_019)*Num(16.000000)+Binary(X_018)*Num(13.000000)+Binary(X_017)*Num(12.000000)+Binary(X_016)*Num(10.000000)+Binary(X_014)*Num(5.000000)+Binary(X_015)*Num(8.000000)),\n",
+ " (Binary(X_027)*Num(19.000000)+Binary(X_026)*Num(16.000000)+Binary(X_025)*Num(13.000000)+Binary(X_024)*Num(12.000000)+Binary(X_023)*Num(10.000000)+Binary(X_021)*Num(5.000000)+Binary(X_022)*Num(8.000000)),\n",
+ " (Binary(X_034)*Num(19.000000)+Binary(X_033)*Num(16.000000)+Binary(X_032)*Num(13.000000)+Binary(X_031)*Num(12.000000)+Binary(X_030)*Num(10.000000)+Binary(X_028)*Num(5.000000)+Binary(X_029)*Num(8.000000)),\n",
+ " (Binary(X_041)*Num(19.000000)+Binary(X_040)*Num(16.000000)+Binary(X_039)*Num(13.000000)+Binary(X_038)*Num(12.000000)+Binary(X_037)*Num(10.000000)+Binary(X_035)*Num(5.000000)+Binary(X_036)*Num(8.000000)),\n",
+ " (Binary(X_048)*Num(19.000000)+Binary(X_047)*Num(16.000000)+Binary(X_046)*Num(13.000000)+Binary(X_045)*Num(12.000000)+Binary(X_044)*Num(10.000000)+Binary(X_042)*Num(5.000000)+Binary(X_043)*Num(8.000000))]"
+ ]
+ },
+ "execution_count": 8,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# here, we use a list to represent the sequence of prices from day T-n to day T+n-1.\n",
+ "\n",
+ "t = len(a) # next number of days to optimize\n",
+ "n_level = len(price_levels) # number of price options\n",
+ "\n",
+ "x = []\n",
+ "p = []\n",
+ "\n",
+ "# get p\n",
+ "for i in range(t):\n",
+ " p_i = 0\n",
+ " for j in range(n_level):\n",
+ " x_ij = Binary(f\"X_{i*n_level+j:03d}\")\n",
+ " x.append(x_ij)\n",
+ " p_i += x_ij*price_levels[j]\n",
+ " p.append(p_i)\n",
+ " \n",
+ "# plus historical prices\n",
+ "all_p = p_data + p\n",
+ "\n",
+ "all_p"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "437d5378",
+ "metadata": {},
+ "source": [
+ "Next, we can input the price variables into our fitted model to estimate the demand, and then use the demand together with the price to represent the revenue $R$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "71dde625",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Revenue:\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "'(((Binary(X_048)*Num(19.000000)+Binary(X_047)*Num(16.000000)+Binary(X_046)*Num(13.000000)+Binary(X...m(12.000000)+Binary(X_009)*Num(10.000000)+Binary(X_007)*Num(5.000000)+Binary(X_008)*Num(8.000000)))'"
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "def get_demand(coeff, b, prices):\n",
+ " assert len(coeff) == len(prices)\n",
+ " d = b\n",
+ " for i in range(len(coeff)):\n",
+ " d += coeff[i]*prices[i]\n",
+ "\n",
+ " return d\n",
+ "\n",
+ "# get d, rev\n",
+ "def get_demands_rev(a,b,hist_p, p):\n",
+ " \"\"\" represent the next n days demands and total revenue, based on:\n",
+ " 1. the fitted coefficients a,b\n",
+ " 2. the historical price hist_p\n",
+ " 3. and the future price decisions represented by binary variable x\n",
+ " \"\"\"\n",
+ " all_p = hist_p + p\n",
+ " t = len(a)\n",
+ " d = [get_demand(\n",
+ " coeff=a,\n",
+ " b=b,\n",
+ " prices=all_p[i+1:i+1+t]\n",
+ " ) for i in range(t)]\n",
+ " rev = np.dot(d, p)\n",
+ " return d, rev\n",
+ "\n",
+ "d, rev = get_demands_rev(a,b,p_data, p)\n",
+ "\n",
+ "print('Revenue:')\n",
+ "repr_compact.repr(rev)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "88ec1a39",
+ "metadata": {},
+ "source": [
+ "## 2.2 Add penalty for prediction uncertainty"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "05742984",
+ "metadata": {},
+ "source": [
+ "In optimization, we usually desire to reduce the risk while maximizing the reward, which means we want to reduce the estimated uncertainty/variance of the predicted demands. Therefore, we add a penalty term in the objective function to represent the uncertainty of the predicted demand. The uncertainty of the predicted demand $d_t$ can be represented by its estimated variance,\n",
+ "\n",
+ "$$H=Revenue-\\beta\\sum_{t=T}^{T+n-1}var(d_t)$$\n",
+ "\n",
+ "where $\\beta$ is the regularization parameter to control the effect of risk estimation. The variance of the demand predictions can be estimated as [10]: $$var(d_t)=\\sigma^2(1+\\vec{p}_t'(\\vec{X}'\\vec{X})^{-1}\\vec{p}_t)$$\n",
+ "\n",
+ "where $\\vec{p}_t'=[1, p_{t-n+1},p_{t-n+2}, \\ldots,p_{t}]$ is the price vector to estimate demand $d_t$. $\\vec{X}$ are the observations in the training dataset used to fit the demand model, where each row is an observation in the training set.\n",
+ "\n",
+ "$$\\vec{X} = \\left[\n",
+ " \\begin{matrix}\n",
+ " 1 & p_{0,1} & p_{1,1} & \\dots & p_{n-1,1}\\\\\n",
+ " 1 & p_{0,2} & p_{1,2} & \\dots & p_{n-1,2}\\\\\n",
+ " \\vdots& \\vdots& \\vdots&& \\vdots\\\\\n",
+ " 1 & p_{0,N}& p_{1,N}& \\dots& p_{n-1,N}\n",
+ " \\end{matrix}\\right].$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "6c3d500e",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Objective with expected revenue and prediction variance:\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "'((((Binary(X_048)*Num(19.000000)+Binary(X_047)*Num(16.000000)+Binary(X_046)*Num(13.000000)+Binary(...X_001)*Num(8.000000))*Num(-0.000535)+Num(0.016381))+Num(1.000000))*Num(100.000000))*Num(-1.000000))'"
+ ]
+ },
+ "execution_count": 10,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "beta=1.\n",
+ "\n",
+ "def get_variance(data_x, p, sigma):\n",
+ " \"\"\"\n",
+ " :param data_x (np.array): [n_samples, n_days]\n",
+ " :param p (list): [n_days]\n",
+ " :return: variance\n",
+ " \"\"\"\n",
+ " n_samples, t = data_x.shape\n",
+ " ones = np.ones((n_samples, 1), dtype=np.float)\n",
+ " x_mat = np.concatenate([ones, data_x], axis=1) # [n_samples, n_days+1]\n",
+ " x_mat = np.linalg.inv(\n",
+ " np.dot(x_mat.T, x_mat)\n",
+ " )\n",
+ " p = np.array([1.]+p)\n",
+ " variance = (sigma**2) * (1. + p.dot(x_mat).dot(p))\n",
+ " return variance\n",
+ "\n",
+ "def get_overall_variance(data_x, hist_p, p, sigma):\n",
+ " all_p = hist_p + p\n",
+ " t = len(p)\n",
+ " var = 0\n",
+ " for i in range(t):\n",
+ " var += get_variance(data_x, all_p[i+1:i+1+t], sigma)\n",
+ "\n",
+ " return var\n",
+ "\n",
+ "objective = rev - beta * get_overall_variance(data_x, p_data, p, sigma)\n",
+ "\n",
+ "print('Objective with expected revenue and prediction variance:')\n",
+ "repr_compact.repr(objective)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "52a8d37f",
+ "metadata": {},
+ "source": [
+ "## 2.3 Add penalty for equality constraints"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0a92e9af",
+ "metadata": {},
+ "source": [
+ "Since the price can only take one value per day, exactly one of the binary variables in a day must be $1$, and the others must be $0$. Formallly, we have equality constraints:\n",
+ "$$\\sum_{k=0}^{m-1}x_{t,k}=1, \\;t\\in[T, T+1, \\cdots, T+n-1]$$\n",
+ "\n",
+ "We can easily incorporate this equality constraint into the objective function by substracting a penalty term $H_p$ with large enough coefficient $L_p$. If the solution satisfies these constraints, the penalty term will be 0, otherwise a large penalty will be imposed.\n",
+ "\n",
+ "$$H_p=L_p\\sum_{t=T}^{t=T+n-1}[(\\sum_{k=0}^{m-1}x_{t,k})-1]^{2}$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "id": "96c8c9be",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Objective with additional equality constraints:\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "'(((((((((((Binary(X_048)*Num(19.000000)+Binary(X_047)*Num(16.000000)+Binary(X_046)*Num(13.000000)+...nary(X_044)+Binary(X_042)+Binary(X_043))+Num(-1.000000))*Num(10000000000000.000000)*Num(-1.000000))'"
+ ]
+ },
+ "execution_count": 11,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# add equalty constraints\n",
+ "Lp=1e13\n",
+ "\n",
+ "for i in range(t):\n",
+ " penalty_i = ((sum(x[i*n_level:(i+1)*n_level]) - 1)**2)*Lp\n",
+ " objective -= penalty_i\n",
+ " \n",
+ "print('Objective with additional equality constraints:')\n",
+ "repr_compact.repr(objective)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6fbf2237",
+ "metadata": {},
+ "source": [
+ "# 3. Solve QUBO with Braket"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9edaa67f",
+ "metadata": {},
+ "source": [
+ "Now we have transformed to the total objective function which:\n",
+ "- maximizes the revenue \n",
+ "- minimizes the prediction uncertainty \n",
+ "- incorporates equality constraints, $\\sum_{k=0}^{m-1}x_{t,k}=1, \\;t\\in[T, T+1, \\cdots, T+n-1]$ \n",
+ "\n",
+ "$$H=Revenue-\\beta\\sum_{t=T}^{T+n-1}var(d_t) - H_p $$\n",
+ "where $H_p=L_p\\sum_{t=T}^{t=T+n-1}[(\\sum_{k=0}^{m-1}x_{t,k})-1]^{2} $ is the penalty term for price equality constraints, and $L_p$ is a constant.\n",
+ "\n",
+ "We can transform this objective function to binary quadratic model, and solve it using Amazon Braket."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "2547ec98",
+ "metadata": {},
+ "source": [
+ "Amazon Braket provides access to quantum annealers from D-Wave Systems Inc. which are specifically designed to solve QUBO problems. A quantum annealer is a hardware implementation of Quantum Annealing (QA) which is a generic algorithm to solve combinatorial optimization problems by taking advantage of quantum fluctuations to escape local minima in complex energy landscapes by effect of tunnelling through barriers separating local minima [5]. Through a simple Amazon Braket API call, the formatted QUBO problem can be embedded to the hardwired topology of the machine and submitted to machine. In turn, the quantum annealer performs an annealing process and obtains a result. However, because of the sparsity of the underlying chip, one typically faces a certain overhead in embedding the original QUBO problem to the D-Wave hardware, because (typically) one logical binary variable maps onto several physical qubits, in order to effectively distill the connectivity of the original QUBO problem. For more details, we refer to the Amazon Braket tutorial notebook available here (https://github.com/aws/amazon-braket-examples/blob/main/examples/quantum_annealing/Dwave_Anatomy.ipynb)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "id": "fdfe3dce",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "\"BinaryQuadraticModel({X_005: -10000000002744.719, X_030: -10000000002576.855, X_036: -100000000021...', 'X_026'): 572.3665512211251, ('X_026', 'X_013'): 906.2470394334479}, 70000000000715.1, 'BINARY')\""
+ ]
+ },
+ "execution_count": 12,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "model = (-objective).compile().to_bqm()\n",
+ "\n",
+ "repr_compact.repr(model)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "id": "54a70c81",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "49"
+ ]
+ },
+ "execution_count": 13,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "len(model)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f149c413",
+ "metadata": {},
+ "source": [
+ "We embed and submit the QUBO task to be solved by a quantum annealer at D-wave backend. Similar to classical simulated annealing, we take several independent shots in order to increase our chances for success in finding a high-quality solution to the QUBO problem. The response contains a number of rows. Each row represents the optimization result of a single shot, including the binary decision variables $[x_T,\\cdots,x_{T+n-1}]$, and corresponding energy which is the negative value of objective function, i.e. $H$. Below we show how to submit a QUBO problem to D-Wave and provide a sample response. Here, the key figure of merit is listed in the column labelled with energy, representing the value for the objective function . We can then easily read off the best bit string found.\n",
+ "\n",
+ "Each row represents an optimized solution. The columns starting with `X_` show the value of the binary decision variables $x_{t,k}$ and the `energy` column shows the negative value of the objective function, i.e. $-H$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "id": "b762f1e9",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "
\n",
+ "Caution: Running the following cell will run a substantial amount of tasks for various senarios taking about 30 minutes. Only uncomment it if you are happy to wait.\n",
+ "