$\newcommand{\ket}[1]{\left|{#1}\right\rangle}$
$\newcommand{\bra}[1]{\left\langle{#1}\right|}$
# Challenge 03: Minimum Makespan with Quantum Annealing

## Introduction

The core of this challenge is about translating a real-world problem into a mathematical representation that fits a quantum computer. In this challenge we encourage you to use two representations which even though equivalent, are used in different contexts. 

The first is Quadratic Unconstrained Binary Optimization (QUBO) which is used primarily for the quantum annealers (D-Wave) and the other is Ising Hamiltonian, which is often used with an algorithm called Quantum Approximate Optimization Algorithm (QAOA).

This challenge might require spending some time going through the references (especially 1-3), as we are not aware of any "quick introductions" worth mentioning. If you have a hard time understanding something, feel free to reach out for help on [QOSF Slack](https://join.slack.com/t/qosf/shared_invite/zt-bw59w8b9-WJ~k0~FAMHukTZov4AnLfA).

Unless you have a prior background with these kinds of problems, we encourage you to start from QUBO and then perhaps trying Ising Hamiltonians as well, as there are more materials about it.

It is also worth mentioning that D-Wave allows you to use their machines for free for a limited amount of time, so you can try running your QUBO on a real machine after signing up [here](https://cloud.dwavesys.com/leap/signup/).

## Problem

Niranjan owns a catering service and receives orders one day in advance. Each of the chefs he employs is paid by the hour. In order to optimize the cost of running his business, every night, he looks at all of the dishes he needs to deliver for the next day and asks his employees to work the specific hours he requires them for.

Given an order queue ($Q$) of length $N$ with each order $i$ taking $L_i$ time to prepare, if there are $m$ chefs, design a QUBO/Ising Hamiltonian to help Niranjan find the shortest time ($T$) it would take for the $m$ chefs to prepare all $N$ dishes. Assume each chef can only work on one dish at a time and each dish is worked on by only one chef. We can also safely assume that at any given point, all chefs can be actively working on a dish. 

Since Niranjan wants to ensure that he can run his algorithms on the near term devices, keep track and try to minimize the number of:
1. Gates & Qubits (for Ising Hamiltonian)
2. Variables (for QUBO. For embedding on a real machine also keep track of the lenghts of chains)

**Optional:** 
For the advanced reader, you can also attempt to build the Cost Hamiltonian or Mixer Hamiltonian as described in [7,8].

## Examples

### Case #1: 

*If,* $\quad N = 5 , \quad m = 1, \quad \text{queue} = \{2, 3, 7, 2, 1\} \quad$ *then,* $\quad T=15$

*since* there is only one chef, the shortest time it would take to complete all the orders is simply the sum of the times it would take for each individual order.

### Case #2:

*If,* $\quad N = 5 , \quad m = 2, \quad \text{queue} = \{2, 3, 7, 2, 1\} \quad$ *then,* $\quad T=8$

*since* there are two chefs $m_1$ and $m_2$, the distribution that yields the shortest time is:

$m_1 = \{2, 3, 2\}, \qquad m_2 = \{ 7, 1 \} $

### Case #3:

*If,* $\quad N = 5 , \quad m = 3, \quad \text{queue} = \{2, 3, 2, 2, 1\} \quad$ *then,* $\quad T=4$

*since* there are three chefs $m_1$, $m_2$ and $m_3$, the distribution that yields the shortest time is:

$m_1 = \{3 \}, \qquad m_2 = \{2,2 \}, \qquad m_3 = \{2,1 \} $ 

Note that there could be another possible distribution of tasks, but the shortest time does not change:

$m_1 = \{3,1 \} , \qquad m_2 = \{2,2 \} , \qquad m_3 = \{2 \} $

## Solution by peguerosdc

The third example in the previous section is what gave the hint and motivated this approach.

Defining $m_i$ as the queue of dishes assigned to chef i and defining the total cost of any queue as $cost(m_i)$, we first note that when a chef is assigned with the **optimal queue** $m_{optimal}$ (that is, the queue that satisfies $cost(m_{optimal})=T$), the rest of the chefs can be assigned any combination of the remaining dishes and they are all still valid solutions as long as they fall within the constraint that $cost(m_i) \leq T$ (because of the definition of $m_{optimal}$).

In both of the solutions of the previous example, the **optimal chef** (the one that is assigned with $m_{optimal}$) is always $m_2$, but that doesn't have to be the case as all the chefs are indistinguishable, so it really doesn't matter! That means we can pick any chef as the **optimal chef** (in the following, we will pick $m_1$) and try to minimize its cost subject the constraints that make it the true **optimal chef**:

1. Every dish can be prepared once and only once. This ensures that the solution found is a valid solution.
2. The time it takes to $m_1$ to cook all his dishes is enough time for the other chefs to also cook their dishes: 

$$cost(m_1) = T \geq cost(m_i) \quad \text{for any} \quad i\neq 1$$


To encode this problem into a QUBO/Ising optimisation problem, we use unary encoding to represent the dishes $d_i$ present in a given order Q of N dishes:

$$ \ket{d_0} = \ket{1000...0_N}$$

$$ \ket{d_1} = \ket{0100...0_N}$$

And so on. So, if the order is given by $Q=\{2,3,2,2,1\}$ and chef $m_1$ is preparing the first and the third dishes, his queue/state will be given by:

$$ \ket{m_1} = \ket{10100}$$

And the state $\ket{\Psi}$ of a system of $m$ chefs will be given by the tensor product (the $\otimes$ symbol is omitted for the sake of readability) of $m$ of these states:

$$ \ket{\Psi} = \ket{m_1}\ket{m_2}...\ket{m_m} $$

To find the most optimal solution to $\ket{\Psi}$, we need to build an Ising/QUBO Hamiltonian where the energy represents the cost $m_1$ (which is the quantity we want to minimize) subject to the 2 constraints mentioned above.

### Objective function and constraints

The QUBO Hamiltonian is coded using D-Wave's library and it is described in terms of the binary variables $\{ m_{ij} \}$ which will be 1 if the qubit j of the chef i is in the state $\ket{1}$ and 0 if it is in the state $\ket{0}$.

It is important to mention that even though this library allows us to define the algorithm already in terms of $\{ m_{ij} \}$, we are describing each state in the basis of the $Z$ gate (that is, $\ket{1}$ represents a spin up state and $\ket{0}$ a spin down state), so if we wanted to run the algorithm on a real quantum machine (please someone correct me if I am wrong) or in some other simulators (like pyquil+Grove), we need to replace the variables $m_i$ with the following operators in terms of the $Z$ gate as their eigenvalues translate to the binary description we desire:

$$ m_i \rightarrow \hat{m_i} = \frac{1 - \hat{Z_i}}{2} \quad \text{because} \quad \hat{m_i}\ket{0_i} = 0, \quad \hat{m_i}\ket{1_i} = \ket{1_i}   $$

Following with the $\{ m_i \}$ description (again for the sake of readability), the cost of the optimal chef $m_1$ is encoded as:

$$ H_0 = \vec{Q} \cdot \vec{m_1} $$

Where $\vec{Q}$ is a vector of length $N$ where the entry $i$ corresponds to the cost of the dish $d_i$. Now, it is left to map the two constraints as energy penalties to ensure that the solution is going to be valid and optimal.

#### Constraint 1: every dish can be prepared once and only once

In terms of $\ket{\Psi}$, this constraint means that for any given dish j, only one $m_{ij}$ must be set for all posible values of i. This can be translated to an equality constraint with Lagrange multipliers:

$$ H_1 = \lambda_1 \sum_i^{N} \left( \sum_j^m m_{ji} -1 \right)^2 $$

#### Constraint 2: chef m1 must be the optimal chef

To consider this inequality constraint $cost(m_1) \geq cost(m_i)$ in our optimization problem, we need to translate it to an equality constraint with the aid of one slack variable:

$$ H_2 = \lambda_2 \sum_{i>1}^{m} \left( \vec{Q}\cdot (\vec{m_1} - \vec{m_i}) - s \right)^2 $$

This slack variable is bounded by the worst case scenario in which this constraint is violated, which is given by the case where no dishes are prepared by the optimal chef and all of them are prepared by any other chef $i$, so the value of $s$ is always gonna be capped by $cost(Q)$. Remembering that the algorithm is a **binary** algorithm, we need to encode the slack variable in the $log2( max(s) )$ binary variables that will be represented by ancilla qubits.

#### Hamiltonian

Putting everything together, the objective function is given by $H=H_0 + H_1 + H_2$:

$$ H = \vec{Q} \cdot \vec{m_1} + \lambda_1 \sum_i^{N} \left( \sum_j^m m_{ji} -1 \right)^2 + \lambda_2 \sum_{i>1}^{m} \left( \vec{Q}\cdot (\vec{m_1} - \vec{m_i}) - s \right)^2 $$

Which is trivially quantized to:

$$\hat{H} = \hat{\vec{Q}} \cdot \hat{\vec{m_1}} + \lambda_1 \sum_i^{N} \left( \sum_j^m \hat{m_{ji}} -1 \right)^2 + \lambda_2 \sum_{i>1}^{m} \left( \hat{\vec{Q}}\cdot (\hat{\vec{m_1}} - \hat{\vec{m_i}}) - \hat{s} \right)^2 $$

and its energy eigenvalues are given by:

$$E = cost(m_1) + \text{how much the constraints are violated}$$ 

Which means that the annealing process will (most likely) drive the system to the state where no constraints are violated and $cost(m_1)$ is at its minimum.

Weights $\lambda$ have to be big enough so that even in the worst case scenario (that is, the state for which the energy $E_0$ of $\hat{H_0}$ is maximum), the constraints of $H_1 + H_2$ are still enforced:

$$
\lambda > max(E_0) = sum(\text{dishes})
$$

So, we choose $\lambda=sum(\text{dishes})+1$

### Coding the solution


```python
# math libraries
import math
import numpy as np
# library to manipulate symbolic expressions
from sympy import symbols
# dimod "simulator"
from dimod import ExactSolver
# d-wave libraries
from dwave.system import DWaveSampler, EmbeddingComposite

class MinimumMakespanSolver():

    def __init__(self, dishes, chefs):
        self.dishes = dishes
        self.m = chefs
        self.N = len(dishes)
        # calculate the amount of slack qubits required per chef
        self.ancilla_per_chef = math.ceil(math.log2(sum(dishes)))
        # prepare variables to build the objective function
        amount_qubits = self.N * self.m + (self.m-1)*self.ancilla_per_chef
        self.syms = symbols(f"x:{amount_qubits}")
        # variable to store the QUBO
        self.qubo = None

    def get_index(self, dish, chef):
        """Get the index of the qubit corresponding to a chef's dish"""
        return self.N*chef + dish

    def get_bit(self, dish, chef):
        """Get the symbol of the qubit corresponding to a chef's dish"""
        return self.get_raw_bit(self.get_index(dish, chef))

    def get_raw_bit(self, index):
        """Get the symbol corresponding to qubit 'index' """
        return self.syms[index]

    def get_indexes_of_ancilla(self, chef):
        """Get the indexes of the slack/ancilla qubits used by a chef"""
        return [i for i in range(self.N*self.m + (chef-1)*self.ancilla_per_chef, self.N*self.m + chef*self.ancilla_per_chef) ]

    def cost_of_unprepared_dishes(self, chef):
        """Get the expression that maps to the cost of the state of a certain chef"""
        return sum([cost*self.get_bit(dish, chef) for dish,cost in enumerate(self.dishes)])

    def get_objective(self, weight=1):
        """Builds the objective function enforcing three constraints:
           
           1) Minimize the amount of dishes prepared by chef0 = chef_max

           2) No dishes can be cooked more than once

           3) chef_max should take the highest amount of time
        """
        obj = 0
        # 1 constraint: Minimize the amount of dishes prepared by chef0 = chef_max
        obj += self.cost_of_unprepared_dishes(0)
        # 2 constraint: no dishes can be cooked more than once
        cost_operators = 0
        for dish in range(self.N):
            op = -1
            for chef in range(self.m):
                op += self.get_bit(dish, chef)
            cost_operators += op**2
        obj += weight*cost_operators
        # 3 constraint: chef_max should take the highest amount of time
        cost_operators = 0
        for chef in range(1,self.m):
            # get slack binary variables
            slack = sum( [ (2**i)*self.get_raw_bit(ancilla) for i,ancilla in enumerate(self.get_indexes_of_ancilla(chef))] )
            # get total hamiltonian
            cost_operators += (self.cost_of_unprepared_dishes(0) - self.cost_of_unprepared_dishes(chef) - slack)**2
        obj += weight*cost_operators
        # return expression with all the constraints
        return obj

    def build(self, weight=None):
        """Builds the QUBO for the given problem. Returns the linear and quadratic coefficients.
        If no weight is passed, it is calculated with the worst case scenario
        """
        # get best weight
        weight = weight if weight else sum(self.dishes)+1
        # get the raw function to minimize
        obj = self.get_objective(weight)
        print(f"Raw expression: {obj}")
        # simplify one-variable quadratic terms as linear terms
        for s in self.syms:
            obj = obj.expand().replace(s**2, s)
        # decompose in factors
        terms = obj.expand().as_coefficients_dict()
        # build QUBO
        linear = dict()
        quadratic = dict()
        for term in terms:
            # check how to map each term
            readable = str(term)
            # if the term is in the shape a*xn*xm (that is, with two variables), it is a quadratic term
            if readable.count("x") == 2:
                x0,x1 = readable.split("*")
                quadratic[(x0,x1)] = terms[term]
            # if the term is in the shape a*xn (that is, with one variable), it is a linear term
            elif readable.count("x") == 1:
                linear[(readable, readable)] = terms[term]
        # store the QUBO
        self.qubo = {**linear, **quadratic}
        return self.qubo

    def solve(self, qubo=None, quantum=False):
        """Solves a given QUBO. To use the quantum mode, the D-Wave API keys must be set in the environment.

        Args
        ----
        qubo: coefficients of the problem to solve. If none given, the qubo stored in the class is used.
        quantum: True if meant to be run on a d-wave machine, False if meant to be solved with dimod.ExactSolver
        """
        # get the correct qubo to run
        qubo = qubo if qubo else (self.qubo if self.qubo else self.build())
        # run
        samples = EmbeddingComposite(DWaveSampler(solver={'topology__type': 'chimera'})).sample_qubo(qubo, num_reads=1000) if quantum else ExactSolver().sample_qubo(qubo)
        # get only results with lowest energy as they make up the solution
        results = samples.lowest()
        data = results.record.sample
        # sort the variables correctly
        sorted_data = np.zeros((len(data), self.m, self.N))
        for r,result in enumerate(data):
            for m in range(self.m):
                for j in range(self.N):
                    # get index of this dish
                    i = results.variables.index(f"x{self.get_index(j,m)}")
                    # get value
                    sorted_data[r][m][j] = self.dishes[j] if result[i] else 0.
        return sorted_data
```

And here are some useful functions to show the result:


```python
def visualize_solution(results, queue, m):
    for solution in results:
        for i,chef in enumerate(solution):
            print(f"Chef {i} cooked {chef} with cost {sum(chef)}")
        print("-------")
```

#### Case 1


```python
queue = [2,3,7,2,1]
m = 1
solver = MinimumMakespanSolver(dishes=queue, chefs=m)
solution = solver.solve(quantum=False)
print("\nSOLUTION:")
visualize_solution(solution, queue, m)
```

    Raw expression: 2*x0 + 3*x1 + 7*x2 + 2*x3 + x4 + 16*(x0 - 1)**2 + 16*(x1 - 1)**2 + 16*(x2 - 1)**2 + 16*(x3 - 1)**2 + 16*(x4 - 1)**2
    
    SOLUTION:
    Chef 0 cooked [2. 3. 7. 2. 1.] with cost 15.0
    -------


#### Case 2


```python
queue = [2,3,7,2,1]
m = 2
solver = MinimumMakespanSolver(dishes=queue, chefs=m)
solution = solver.solve(quantum=False)
print("\nSOLUTION:")
visualize_solution(solution, queue, m)
```

    Raw expression: 2*x0 + 3*x1 + 7*x2 + 2*x3 + x4 + 16*(x0 + x5 - 1)**2 + 16*(x1 + x6 - 1)**2 + 16*(x2 + x7 - 1)**2 + 16*(x3 + x8 - 1)**2 + 16*(x4 + x9 - 1)**2 + 16*(2*x0 + 3*x1 - x10 - 2*x11 - 4*x12 - 8*x13 + 7*x2 + 2*x3 + x4 - 2*x5 - 3*x6 - 7*x7 - 2*x8 - x9)**2
    
    SOLUTION:
    Chef 0 cooked [2. 3. 0. 2. 1.] with cost 8.0
    Chef 1 cooked [0. 0. 7. 0. 0.] with cost 7.0
    -------
    Chef 0 cooked [0. 0. 7. 0. 1.] with cost 8.0
    Chef 1 cooked [2. 3. 0. 2. 0.] with cost 7.0
    -------


#### Case 3

Note that this case can't be solved in a classical computer using `ExactSolver`. We need to go full quantum!


```python
queue = [2,3,2,2,1]
m = 3
solver = MinimumMakespanSolver(dishes=queue, chefs=m)
solution = solver.solve(quantum=False)
print("\nSOLUTION:")
visualize_solution(solution, queue, m)
```


```python
queue = [2,3,2,2,1]
m = 3
solver = MinimumMakespanSolver(dishes=queue, chefs=m)
solution = solver.solve(quantum=True)
print("\nSOLUTION:")
visualize_solution(solution, queue, m)
```

    Raw expression: 2*x0 + 3*x1 + 2*x2 + 2*x3 + x4 + 11*(x0 + x10 + x5 - 1)**2 + 11*(x1 + x11 + x6 - 1)**2 + 11*(x12 + x2 + x7 - 1)**2 + 11*(x13 + x3 + x8 - 1)**2 + 11*(x14 + x4 + x9 - 1)**2 + 11*(2*x0 + 3*x1 - 2*x10 - 3*x11 - 2*x12 - 2*x13 - x14 - x19 + 2*x2 - 2*x20 - 4*x21 - 8*x22 + 2*x3 + x4)**2 + 11*(2*x0 + 3*x1 - x15 - 2*x16 - 4*x17 - 8*x18 + 2*x2 + 2*x3 + x4 - 2*x5 - 3*x6 - 2*x7 - 2*x8 - x9)**2
    
    SOLUTION:
    Chef 0 cooked [2. 0. 0. 2. 0.] with cost 4.0
    Chef 1 cooked [0. 3. 0. 0. 0.] with cost 3.0
    Chef 2 cooked [0. 0. 2. 0. 1.] with cost 3.0
    -------


Note that because of the probabilistic nature of the quantum realm, the solution provided by a quantum computer will not always be the most optimal, but still it would be optimal. When I ran my code I was lucky the best solution was provided, but sometimes (just a few, really), the cost can be 5 (which is almost 4, but not 4).

## Resources

\[1\] [D-Wave Problem Solving Handbook](https://docs.dwavesys.com/docs/latest/doc_handbook.html)

\[2\] [Implementation of the Travelling Salesman Problem](https://github.com/mstechly/quantum_tsp_tutorials)

\[3\] [Ising formulations of many NP problems](https://arxiv.org/pdf/1302.5843.pdf)

\[4\] [D-Wave Examples](https://docs.ocean.dwavesys.com/en/stable/getting_started.html#examples)

\[5\] [Quantum Approximate Optimization Algorithm explained](https://www.mustythoughts.com/quantum-approximate-optimization-algorithm-explained)

\[6\] [Job Shop Scheduling Solver based on Quantum Annealing](https://arxiv.org/pdf/1506.08479.pdf)

### References

\[7\] [Quantum Algorithms for Scientific Computing and Approximate Optimization](https://arxiv.org/pdf/1805.03265.pdf)

\[8\] [From the Quantum Approximate Optimization Algorithm to a Quantum Alternating Operator Ansatz](https://arxiv.org/pdf/1709.03489.pdf)

### Other useful resources
\[9\] [D-Wave Getting Started with the System](https://docs.dwavesys.com/docs/latest/doc_getting_started.html)

\[10\] [Qiskit Converters for Quadratic Programs](https://qiskit.org/documentation/tutorials/optimization/2_converters_for_quadratic_programs.html)

\[11\] [D-Wave Support Difference between BQM, Ising, and QUBO problems?](https://support.dwavesys.com/hc/en-us/community/posts/360017439853-Difference-between-BQM-Ising-and-QUBO-problems-)

\[12\] [A QAOA solution to the traveling salesman problem using pyQuil](https://cs269q.stanford.edu/projects2019/radzihovsky_murphy_swofford_Y.pdf)

78eca64b53e2f85092d978b007915e9153125e36
ipynb
Jupyter Notebook
challenge-2021.01-jan/peguerosdc/qosf-monthly-challenge-03-dwave.ipynb
nbaehler/monthly-challenges
59c23b7d0b6e4646d6e558e882dad507f96af77b
Apache-2.0
2020-10-25T03:09:04.000Z
2022-03-29T11:22:58.000Z
challenge-2021.01-jan/peguerosdc/qosf-monthly-challenge-03-dwave.ipynb
nbaehler/monthly-challenges
59c23b7d0b6e4646d6e558e882dad507f96af77b
Apache-2.0
2020-11-04T02:24:40.000Z
2022-01-10T21:12:09.000Z
challenge-2021.01-jan/peguerosdc/qosf-monthly-challenge-03-dwave.ipynb
MaldoAlberto/monthly-challenges
7c669afddf103e93b9bb64f231f5d4baf8d07b96
Apache-2.0
2020-10-25T16:03:00.000Z
2022-03-30T19:11:25.000Z
Qwen/Qwen-72B
1. YES
2. YES