---
author:
- "[Eugeniy E. Mikhailov](http://physics.wm.edu/~evmik/) and [Greg
  Bentsen](https://physics.wm.edu/~gbentsen/)"
title: "Numerical Optimization IV: Combinatorial optimization"
---

## Logistics and Agenda

Midterm 2

-   Take-home posted Monday Oct 26
    -   report and code are due Sunday Nov 1st at 11:59pm
    -   you can work in groups but you must write **your own report**
-   In-class exam Monday Nov 2

\

Agenda

-   This week:  Combinatorial optimization and linear algebra methods
-   Next week:  midterm and neural networks as base for machine learning

## Today

Optimization in multiple dimensions

-   Revisiting simulated annealing
-   Back pack problem
-   Traveling salesman problem

## Simulated annealing a.k.a. Metropolis--Hastings algorithm

$$
\boxed{P(\mathrm{accept})=\min\left(1,e^{-\Delta E/T}\right)}
$$

with

$$
\Delta E=E_{\mathrm{new}}-E_{\mathrm{old}}
$$

-   Downhill moves $\Delta E < 0$: always accepted
-   Uphill moves $\Delta E >0$: accept with probability
    $e^{-\Delta E/T}$

Uphill moves are especially probable at **high temperature** $T$, but
occur less frequently at low temperature.

**Gradually lower $T$:** explore first, settle later. Good news: we can
find global minimum, bad news it might take infinite amount of time.

\

Question: what is the good starting temperature and final temperature?

. . .

-   $T_{\mathrm{start}}$ should be about the size of the largest "hill"
-   $T_{\mathrm{final}}$ should be about the size of the "measurement"
    noise

## Combinatorial optimization problem statement

We still want to optimize (minimize) our multidimensional energy or cost
or merit function $E$

::: definition
Find $\vec{x}=x_1, x_2, x_3 \cdots$ that minimizes $E(\vec{x})$ subject
to $g(\vec{x})=0$, $h(\vec{x}) \le 0$
:::

The only difference is that values of the $\vec{x}$ are discrete,
i.e. any component of the $\vec{x}$ can take a **countable** set of
different values.

In this case, we cannot run our golden search algorithm or anything else
which assumes continuous space for the $\vec{x}$.

. . .

Well, actually we can but constructing a proper constraining function
will be a nightmare.

. . .

Instead, we have to find a method to search through discrete sets of all
possible input values, i.e. try all possible combinations of $\vec{x}$
components.

Hence, the name `\alert{combinatorial optimization}`{=tex}.

## Example: Backpack problem

::: columns
::: {.column width="50% "}
Suppose you have a backpack with the given size (volume) and a set of
objects with given volumes and monetary (or sentimental) values.

Our job is to find a subset of items that still fits in the backpack and
has the maximum combined value.

For simplicity, we will assume that every item occurs only once.

Our job is to maximize $$\begin{aligned}
    E(\vec{x}) = \sum \text{value}_i x_i = \overrightarrow{\text{values}} \cdot \vec{x}
\end{aligned}$$ Subject to the following constraints $$\begin{aligned}
    \sum \text{volume}_i x_i = \overrightarrow{\text{volumes}} \cdot \vec{x}
 \le \text{BackpackSize}
\end{aligned}$$

Where $x_i=(0$ or $1)$, i.e. it reflects whether we take this object or
not.
:::

::: {.column width="50%"}
![Backpack illustration (generated by
ChatGPT)](./pics/backpack_problem.png){width="100%"}
:::
:::

## Brute force backpack optimization: binary encoding

Recall that the $\vec{x}$ component is either 0 or 1.

So, $\vec{x}$ is a combination of zeros and ones
$\vec{x}=[0,1,0,1,\cdots,1,1,0,1,1]$.

How would we generate all possible combinations of $\vec{x}$ components?

. . .

-   $\vec{x}$ looks like a binary number.
-   let's start with $\vec{x}=[0,0,0,0,\cdots,0,0]$
-   every new component will be generated by adding 1 to the previous
    $x$ according to binary addition rules
    -   for example $x_{next}=[1,0,1,\cdots,1,1,0,1,1]+1 =
                            [1,0,1,\cdots,1,1,1,0,0]$
    -   recall that $1 + 1 = 10_2$
-   for every new $\vec{x}$, we check to see if the items fit into the
    backpack and if the new packed value is larger than the previous
-   once we have tried all $2^N$ combinations of $\vec{x}$, we are done

Notice that if there are $N$ objects then the number of all possible
combinations is $2^N$. Thus the time of the optimization grows
`\alert{exponentially}`{=tex} with the number $N$ of items to chose, but
we will find the `\alert{global}`{=tex} optimum.

## Backpack optimization with binary encoding: test run

-   My computer sorts through 20 items in 4.4 seconds, but
    `\alert{30}`{=tex} items would take 1000 times longer, i.e. about
    `\alert{75}`{=tex} minutes to solve.

``` python
#!/usr/bin/env python3
import numpy as np

def backpack_binary(max_capacity, values, volumes):
    """ Solve backpack problem using binary coding approach."""
    n = len(values)
    max_value = 0
    best_combination = np.zeros(n, dtype=int)

# Try all possible combinations using binary representation

# From 0 to 2^n - 1 (all possible binary numbers with n bits)

for i in range(2**n):
        selected_items = np.zeros(n, dtype=int)
        
        code = bin(i)[2:].zfill(n) # string representing bits
        for i,s in enumerate(code):
            selected_items[i] = 1 if s=="1" else 0
        total_volume = np.sum(volumes * selected_items)
        total_value = np.sum(values * selected_items)

        if total_volume <= max_capacity:
            if total_value > max_value:
                max_value = total_value
                best_combination = selected_items.copy()
    return max_value, best_combination


def print_solution(values, volumes, max_capacity, max_value, selected_items):
    print(f"Backpack Capacity: {max_capacity}")
    print(f"Maximum Value: {max_value}")
    print(f"Selected Items: {selected_items}")
    print(f"{values=}")
    print(f"{volumes=}")
    print(f"Total number of objects: {len(values)}")


if __name__ == "__main__":
    print("====== Example 1 ===========")
    values1 = np.array([12, 10, 20, 15], dtype=int)
    volumes1 = np.array([2, 1, 3, 2], dtype=int)
    max_capacity1 = 5
    
    max_val1, selected1 = backpack_binary(max_capacity1, values1, volumes1)
    print_solution(values1, volumes1, max_capacity1, max_val1, selected1)
    
    print("====== Example 2 ===========")
    N=20
    values2 = (100*np.random.uniform(size=(N))).astype(int)
    volumes2 = (80*np.random.uniform(size=(N))).astype(int)
    max_capacity2 = 50
    
    import time
    t0 = time.time()
    max_val2, selected2 = backpack_binary(max_capacity2, values2, volumes2)
    print_solution(values2, volumes2, max_capacity2, max_val2, selected2)
    print(f"Execution time: {time.time() - t0}")
```

For realization of this algorithm, have a look at the
[backpack_binaray.py](code/backpack_binaray.py)