# Optimization Sample: Multi-ship loading

In this example, we will take our learnings from the ship-loading sample and generalize to load-balancing between any number of ships. In addition, we'll see how we can make use of the parameter free-solvers to guide us on a selection of parameters which we can use with the parametrized solvers. 

## Pre-requisites

1. [Create an Azure Quantum Workspace](https://github.com/MicrosoftDocs/quantum-docs-private/wiki/Create-quantum-workspaces-with-the-Azure-portal)
2. [Install the `azure-quantum` Python module](https://github.com/MicrosoftDocs/quantum-docs-private/wiki/Use-the-Python-SDK-for-Quantum-Inspired-Optimization).
3. [Complete the ship-loading sample](https://github.com/microsoft/qio-samples/tree/main/samples/ship-loading)


## The Problem

To tackle this problem we will use a PUBO format. As a reminder, these are [cost functions](https://docs.microsoft.com/en-us/azure/quantum/optimization-concepts-cost-functions) where the variables take the values of either 0 or 1 (rather than -1 or 1 for an Ising cost function).
 
In order to balance containers between multiple ships, one option is to define a cost function that:
1. Penalizes variance from a theoretical equal distribution (where an equal distribution is the total weight of the containers divided by the number of ships), and,
2. Penalizes the assignment of the same container on multiple ships

We will create two sub cost-functions $H1$ and $H2$ that we will then sum to evaluate the total cost of a solution. Let's begin with the first cost function, $H1$.

## Penalize variance from equal distribution between ships

Suppose we had 3 containers with respective weights $W0$, $W1$, $W2$, and we define an equal distribution of the container weights to be: 

$$EqDistrib = (W0 + W1 + W2) / 3$$ 

A way to penalize a large variance from the equal distribution for a given ship is to express it in the following way:

$$(W0 + W1 + W2 - EqDistrib)^2$$

Let's take the following example:

|                 |               |
|-----------------|---------------|
|Container weights| 1, 5, 9, 7, 3 |
|Total weight     |      25       |
|Ships            |     A, B, C   |
|EqualDistrib     | 25 / 3 = 8.33 |

Suppose we were to assign those containers to the ships listed, we can calculate the variance from the equal distribution for each of the given ships, shown in the rightmost column:

| Ships\Containers | 1 | 5 | 9 | 7 | 3 |            |        |
|------------------|---|---|---|---|---|------------|--------|
| A                | 0 | 0 | 9 | 0 | 0 | (9-8.33)^2 |= 0.4489|
| B                | 0 | 5 | 0 | 0 | 3 |(5+3-8.33)^2|= 0.1089|
| C                | 1 | 0 | 0 | 7 | 0 |(1+7-8.33)^2|= 0.1089|

As we need to represent our problem in a binary format we need to "encode" the presence ($x_i=1$) or absence ($x_i=0$) of a given container on a ship. To do this, we need to have a label for the weight of each container on each ship. The table below shows how we assign this continuous index by repeating the list of container weights for each ship and assigning a single list of weight labels across all three ships:

|| Ship A |     |     |     |     |  Ship B |     |     |     |     | Ship C   |        |        |        |        |
|---|--------|-----|-----|-----|-----|--------|-----|-----|-----|-----|----------|--------|--------|--------|--------|
|Container weight| 1 |  5  |  9  |  7  |  3  | 1 |  5  |  9  |  7  |  3  | 1 | 5 | 9 | 7 | 3  |
|Weight label|*w<sub>0</sub>*|*w<sub>1</sub>*|*w<sub>2</sub>*|*w<sub>3</sub>*|*w<sub>4</sub>*| *w<sub>5</sub>*|*w<sub>6</sub>*|*w<sub>7</sub>*|*w<sub>8</sub>*|*w<sub>9</sub>*|*w<sub>10</sub>*|*w<sub>11</sub>*|*w<sub>12</sub>*|*w<sub>13</sub>*|*w<sub>14</sub>*|

The cost function $H1$ becomes:

$$ H1 = H_{A} + H_{B} + H_{C} $$

where:

$$ H_{A} = (w_0 x_0 + w_1 x_1 + w_2 x_2 + w_3 x_3 + w_4 x_4 - EqDistrib)^2 $$

$$ H_{B} = (w_5 x_5 + w_6 x_6 + w_7 x_7 + w_8 x_8 + w_9 x_9 - EqDistrib)^2 $$

and 

$$ H_{C} = (w_{10} x_{10} + w_{11} x_{11} + w_{12} x_{12} + w_{13} x_{13} + w_{14} x_{14} - EqDistrib)^2 $$

We can expand the above and group the common terms, for example if we expand $H_{A}$, we get:

$$
\begin{align}
H_{A} &= (\sum_i(w_i x_i) - EqDistrib)^2\\ 
&= (w_0 x_0 + w_1 x_1 + w_2 x_2 + w_3 x_3 + w_4 x_4 - EqDistrib)^2\\
\end{align}
$$ 

To simplify things for the expansion, let's rename the variables as follows:

$$
\begin{align}
w_0 x_0 &= a \\
w_1 x_1 &= b \\
w_2 x_2 &= c \\
w_3 x_3 &= d \\
w_4 x_4 &= e \\
EqDistrib &= f \\
\end{align}
$$ 

So now we have:

$$
\begin{align}
H_{A} &= (\sum_i(w_i x_i) - EqDistrib)^2\\ 
&= (w_0 x_0 + w_1 x_1 + w_2 x_2 + w_3 x_3 + w_4 x_4 - EqDistrib)^2\\
&= (a + b + c + d + e - f)^2\\
&= a^2 + b^2 + c^2 + d^2 + e^2 + f^2 + 2(ab + ac + ad + ae + bc + bd + be + cd + ce + de) - 2(af + bf + cf + df + ef)
\end{align}
$$ 

Substituting our original values back in, this gives us the following:

$$
\begin{align}
H_{A} &= (\sum_i(w_i x_i) - EqDistrib)^2\\ 
&= w_0^2 x_0^2 + w_1^2 x_1^2 + w_2^2 x_2^2 + w_3^2 x_3^2 + w_4^2 x_4^2 + EqDistrib ^2 + 2(w_0 x_0 \cdot w_1 x_1 + w_0 x_0 \cdot w_2 x_2 + w_0 x_0 \cdot w_3 x_3 + w_0 x_0 \cdot w_4 x_4 + w_1 x_1 \cdot w_2 x_2 + w_1 x_1 \cdot w_3 x_3 + w_1 x_1 \cdot w_4 x_4 + w_2 x_2 \cdot w_3 x_3 + w_2 x_2 \cdot w_4 x_4 + w_3 x_3 \cdot w_4 x_4) - 2(w_0 x_0 \cdot EqDistrib +  w_1 x_1 \cdot EqDistrib + w_2 x_2 \cdot EqDistrib + w_3 x_3 \cdot EqDistrib + w_4 x_4 \cdot EqDistrib)
\end{align}
$$ 

We can do the same for $H_{B}$, and $H_{C}$.

## Penalize the assignment of the same container on multiple ships

Using the containers weight encoding above, we can devise a cost function such as this one for the first container:

$$ H_{D} = (w_0 x_0 + w_5 x_5 + w_{10} x_{10} - w_0)^2 $$

As $w_0$, $w_5$ and $w_{10}$ are actually the same value (it is the same container represented across multiple ships) we have: 
$$ H_{D} = (w_0 x_0 + w_0 x_5 + w_0 x_{10} - w_0)^2 $$

If we expand and group the common terms, we get the following:

$$ H_{D} = {w_0}^2 {x_0}^2 + {w_0}^2 {x_5}^2 + {w_0}^2 {x_{10}}^2 + {w_0}^2 +
2 ({w_0}^2 x_0 x_5 + {w_0}^2 x_0 x_{10} + {w_0}^2 x_5 x_{10})
- 2({w_0}^2 x_0 + {w_0}^2 x_5 + {w_0}^2 x_{10}) $$

We can then repeat the above for each container across all ships:

So in addition to:

$$ H_{D} = (w_0 x_0 + w_0 x_5 + w_0 x_{10} - w_0)^2 $$

we also have:

$$ H_{E} = (w_1 x_1 + w_1 x_6 + w_1 x_{11} - w_1)^2 $$ 
$$ H_{F} = (w_2 x_2 + w_2 x_7 + w_2 x_{12} - w_2)^2 $$
$$ H_{G} = (w_3 x_3 + w_3 x_8 + w_3 x_{13} - w_3)^2 $$
$$ H_{H} = (w_4 x_4 + w_4 x_9 + w_4 x_{14} - w_4)^2 $$

Grouping these together into a single cost function, $H2$, we get:

$$ H2 = H_{D} + H_{E} + H_{F}+ H_{G}+ H_{H} $$

which we can expand and group the terms as we did with $H_D$ above.


## Combining our cost functions

The final part of the problem definition is to combine our cost functions $H1$ and $H2$:

$$H = H1 + H2 $$

You will notice that $H1$ and $H2$ have common indices $[i,i]/[m,m]$ and $[i]/[m]$. We will need to be careful to not duplicate them, but sum them, in our final list of terms describing the cost function.


## Solving the Problem in Python

First, we must instantiate a `Workspace` object which allows you to connect to the Workspace you've previously deployed in Azure. Be sure to fill in the settings below which can be retrieved by running `az quantum workspace show`.


```python
# This allows you to connect to the Workspace you've previously deployed in Azure.
# Be sure to fill in the settings below which can be retrieved by running 'az quantum workspace show' in the terminal.
from azure.quantum import Workspace

# Copy the settings for your workspace below
workspace = Workspace (
    subscription_id = "",
    resource_group = "",
    name = "",
    location = ""
)
```

First, let's define a function to add terms according to our definition of the H1 cost function, where we penalized the variance from an equal distribution between the ships.


```python
def AddTermsWeightVarianceCost(start, end, containers, EqDistrib):
    terms: List[Term] = []
    for i,w in enumerate(containers[start:end+1], start):
        # -2*Wi*EqDistrib.xi -2Wi^2.xi (weight variance cost + duplicate container cost)
        terms.append(Term(w=-2*w*EqDistrib - 2*w*w, indices=[i]))
        # Wi^2.xi^2 + Wi^2.xi^2 (weight variance cost + duplicate container cost)
        terms.append(Term(w=2*w*w, indices=[i,i]))

    for c in combinations(range(start, end+1), 2):
        w0 = containers[c[0]]
        w1 = containers[c[1]]
        # 2*Wi*Wj (weight variance cost)
        terms.append(Term(w=2*w0*w1, indices=[c[0],c[1]]))

    return terms
```

Next, let's incorporate our definition of the second part of our cost function, H2, where we wished to penalize the assignment of the same container on multiple ships.


```python
def AddTermsDuplicateContainerCost(start, end, containers):
    terms: List[Term] = []

    # The following is integrated into AddTermsWeightVarianceCost to reduce the number of Terms and speed-up Terms generation
    # for c in combinations(range(start, end+1), 1):
    #     w = containers[c[0]][0]
    #     i1 = containers[c[0]][1]
    #     terms.append(Term(w=w*w, indices=[i1,i1]))              # Wi^2

    # 2.w^2.x_i.x_j terms
    for c in combinations(range(start, end+1), 2):
        w = containers[c[0]][0]
        i1 = containers[c[0]][1]
        i2 = containers[c[1]][1]
        terms.append(Term(w=2*w*w, indices=[i1,i2]))            # Term(w=2*Wm^2, [m,n])

    # The following is integrated into AddTermsWeightVarianceCost to reduce the number of Terms and speed-up Terms generation
    # # for c in combinations(range(start, end+1), 1):
    #     w = containers[c[0]][0]
    #     i1 = containers[c[0]][1]
    #     terms.append(Term(w=-2*w*w, indices=[i1]))              # -2*Wi^2

    # w^2 term
    terms.append(Term(w=containers[start][0]*containers[start][0], indices=[]))

    return terms
```

Combining these together, we can create our cost function definition of the multi-ship problem.


```python
from typing import List
from azure.quantum.optimization import Term
import numpy as np
from itertools import combinations

def createProblemForContainerWeights(containerWeights: List[int], Ships) -> List[Term]:

    terms: List[Term] = []
    containersWithinShip: List[int] = []
    containersAcrossShips: List[int, int] = []
    totalWeight = 0
    EqDistrib = 0

    for c in range (len(containerWeights)):
        totalWeight = totalWeight + containerWeights[c]
    EqDistrib = totalWeight / len(Ships)
    print(Ships)
    print(containerWeights)
    print("Total Weight:", totalWeight)
    print("Equal weight distribution:", EqDistrib)

    # Create container weights in this format:
    # 1  5  9  7  3  - 1  5  9  7  3  - 1   5   9   7   3
    # W0 W1 W2 W3 W4   W5 W6 W7 W8 W9   W10 W11 W12 W13 W14 
    containersWithinShip = containerWeights*len(Ships)

    # Create container weights in this format:
    # 1  1  1  5  5  5  9  9  9  7  7  7  3  3  3
    for i in range(len(containerWeights)):
        for j in range(len(Ships)):
            k = i + j*len(containerWeights)
            containersAcrossShips.append([containersWithinShip[i], k])

    for split in np.array_split(range(len(containersWithinShip)), len(Ships)):
        terms = terms + AddTermsWeightVarianceCost(split[0], split[-1], containersWithinShip, EqDistrib)

    for split in np.array_split(range(len(containersAcrossShips)), len(containerWeights)):
        terms = terms + AddTermsDuplicateContainerCost(split[0], split[-1], containersAcrossShips)

    return terms
```

Before we solve our problem, let's create a function to help us visualize the results of the solver:


```python
def visualize_result(result, containers, ships, target):
    print("\rResult received from: ", target)
    nb_ships = len(ships)
    try:
        config = result['configuration']
        config = list(config.values())
        for ship, sub_config in enumerate(np.array_split(config, nb_ships)):
            shipWeight = 0
            for c,b in enumerate(sub_config):
                shipWeight = shipWeight + b*containers[c]
            print(f'Ship {ships[ship]}: \t' + ''.join(f'{b*containers[c]}' for c,b in enumerate(sub_config)) + ' - ' + str(shipWeight))
    except:
        print('No Configuration')
    try:
        print('Cost: {}'.format(result['cost']))
    except:
        print('No Cost')
    try:
        print('Parameters: {}'.format(result['parameters']))
    except:
        print('No Parameter')
```

Lastly, let's define a function that will submit our problem to Azure Quantum and visualize the results.


```python
def SolveMyProblem(problem, s):
    try:
        # Optimize the problem
        print("Optimizing with:", s.target)
        Job = s.submit(problem)
        Job.wait_until_completed()
        duration = Job.details.end_execution_time - Job.details.begin_execution_time
        if (Job.details.status == "Succeeded"):
            visualize_result(Job.get_results(), containerWeights*len(Ships), Ships, s.target)
            print("Execution duration: ", duration)
        else:
            print("\rJob ID", Job.id, "failed")
        return Job.id
    except BaseException as e:
        print(e)
```

## Submitting our problem to Azure Quantum

Now, we can put it all together. In the following steps we'll:
1. Define the list of containers and their weights
2. Instantiate the problem, creating our list of terms
3. Submit these terms to a parameter-free solver
4. Use the parameters returned to submit the problem to a parametrized solver


```python
# This array contains a list of the weights of the containers:
containerWeights = [3, 8, 3, 4, 1, 5, 2, 2, 7, 9, 5, 4, 8, 9, 4, 6, 8, 7, 6, 2, 2, 9, 4, 6, 3, 8, 5, 7, 2, 4, 9, 4]
Ships = ["A", "B", "C", "D", "E"]
```


```python
# Create the Terms for this list of containers:
terms = createProblemForContainerWeights(containerWeights,Ships)
```


```python
from azure.quantum.optimization import Problem, ProblemType

# Create the Problem to submit to the solver:
nbTerms = len(terms)
problemName = f'Balancing {str(len(containerWeights))} containers between {str(len(Ships))} Ships ({nbTerms:,} terms)'
print(problemName)
problem = Problem(name=problemName, problem_type=ProblemType.pubo, terms=terms)
```


```python
from azure.quantum.optimization import SimulatedAnnealing

# Try to call a solver with different timeout value and see if it affects the results
jobid = SolveMyProblem(problem, SimulatedAnnealing(workspace, timeout=10))
# jobid = SolveMyProblem(problem, SimulatedAnnealing(workspace, timeout=20))
# jobid = SolveMyProblem(problem, SimulatedAnnealing(workspace, timeout=30))
# jobid = SolveMyProblem(problem, SimulatedAnnealing(workspace))
```


```python
# Try using the parameters returned by the parameter free versions and observe the significant performance improvement

# First use the job id to view the parameters selected by the parameter free solver
job = workspace.get_job(jobid)
results = job.get_results()

print(results)
```


```python
# From the results, let's extract the beta_start, beta_stop, restarts, and sweeps parameters selected

beta_start = results["parameters"]["beta_start"]
beta_stop = results["parameters"]["beta_stop"]
restarts = results["parameters"]["restarts"]
sweeps = results["parameters"]["sweeps"]
```


```python
# Now let's call the solver again, this time the parametrized version, using these parameters

jobid = SolveMyProblem(problem, SimulatedAnnealing(workspace, timeout=5, beta_start=beta_start, beta_stop=beta_stop, restarts=restarts, sweeps=sweeps))

```


```python
from azure.quantum.optimization import ParallelTempering, Tabu, HardwarePlatform, QuantumMonteCarlo

# Here's how we could experiment with different solvers from Microsoft

jobid = SolveMyProblem(problem, SimulatedAnnealing(workspace, platform=HardwarePlatform.FPGA, timeout=5))
# jobid = SolveMyProblem(problem, Tabu(workspace, timeout=5))
# jobid = SolveMyProblem(problem, ParallelTempering(workspace, timeout=60))
# jobid = SolveMyProblem(problem, QuantumMonteCarlo(workspace))
```


```python
from azure.quantum.optimization.oneqbit import PathRelinkingSolver

# And how we can submit the same jobs to solvers by third-party providers, such as 1QBit
# Note: PathRelinkingSolver is only available if the 1QBit provider is enabled in your quantum workspace

# jobid = SolveMyProblem(problem, PathRelinkingSolver(workspace))
```


```python

```

d427344b494d31432ef3283a53ddea7c3f3b3eac
ipynb
Jupyter Notebook
samples/multiship-loading-sample/multiship-loading.ipynb
anraman/qio-samples
7576a78749f169b510f16b13e0787c97932e2b78
MIT
samples/multiship-loading-sample/multiship-loading.ipynb
anraman/qio-samples
7576a78749f169b510f16b13e0787c97932e2b78
MIT
2021-09-24T20:38:26.000Z
2021-09-24T20:38:26.000Z
samples/multiship-loading-sample/multiship-loading.ipynb
anraman/qio-samples
7576a78749f169b510f16b13e0787c97932e2b78
MIT
Qwen/Qwen-72B
1. YES
2. YES