Branch and Bound¶
Algorytm branch and bound polega na systematycznym przeszukiwaniu przestrzeni rozwiązań problemu optymalizacji kombinatorycznej przy jednoczesnym ograniczaniu liczby badanych przypadków. W ogólności, składa się z 3 rodzajów operacji:
Rozgałęzianie (branching): Problem jest dzielony na mniejsze podproblemy, które odpowiadają różnym podzbiorom przestrzeni rozwiązań.
Wyznaczanie ograniczeń (bounding): Dla każdego podproblemu oblicza się ograniczenie (np. dolne lub górne), które informuje o tym, jaką najlepszą wartość funkcji celu można osiągnąć w obrębie danego poddrzewa.
Odrzucanie niekorzystnych gałęzi (pruning): Jeśli obliczone ograniczenie dla danego podproblemu nie pozwala na uzyskanie lepszego wyniku niż już znalezione rozwiązanie, cały podproblem jest pomijany, co znacząco zmniejsza przestrzeń poszukiwań.
Branch and Bound dla modelu Isinga.¶
Pokażemy heurystyczną wersję tego algorytmu dla modelu Isinga. Nie daje nam gwarancji znalezienia stanu podstawowego. Zeletą jest łatwośc implementacji. Jest to punk wyjścia do budowania innych algorytmów o podobnej strukturze (np. branch and bound w przestrzeni probabilistycznej.)
Pseudokod wygląda w następująco:
Poniżej jest graficzna reprezentacja jak takie przeszukiwanie drzewa może wyglądać.
# Ustawienie ścierzek dla plików pomocniczych
from setup import setup_paths
setup_paths()
# Implementacja prostego Branch and Bound
import numpy as np
from typing import NamedTuple
from math import inf
from tqdm import tqdm
class State(NamedTuple):
configuration: list
energy: float
# Funkcja pomocniczna, warto zauważyć że są prowadzone pewne modyfikacje
def calculate_energy(J: np.ndarray, h: np.ndarray, state: list):
n = len(h)
k = len(state)
state = state + [0] * (n - k)
state = np.array(state)
return state @ J @ state.T + state @ h
def branch_and_bound(J, h, num_states, num_returned_elements = 1):
n = len(h)
states = [State([], inf)]
for _ in tqdm(range(n), desc="Wykonywanie branch and bound"):
# Branching
temp = []
for state in states:
plus = state.configuration + [1]
plus_energy = calculate_energy(J, h, plus)
temp.append(State(plus, plus_energy))
minus = state.configuration + [-1]
minus_energy = calculate_energy(J, h, minus)
temp.append(State(minus, minus_energy))
# Bounding
temp.sort(key=lambda x: x.energy)
if len(temp) > num_states:
temp = temp[:num_states]
states = temp
return states[:num_returned_elements]
# test dla bardzo małej instancji (24 spiny)
from funkcje_pomocnicze import read_instance, small_pegasus
J, h = read_instance(small_pegasus.path, convention="dwave")
solutions = branch_and_bound(J, h, 100)
solution = solutions[0]
print(f"Otrzymana energia: {solution.energy}")
print(f"Stan podstawowy: {small_pegasus.best_energy}")
Wykonywanie branch and bound: 100%|██████████| 24/24 [00:00<00:00, 915.06it/s]
Otrzymana energia: -39.0 Stan podstawowy: -39.0
Testy algorytmu¶
W następnych komurkach przeprowadzimy testy jak ten algorytm się zachowuje dla różnych rozmiarów instancji oraz różnych wartości
# wykresy zachowania algorytmu
import os
import time
import matplotlib.pyplot as plt
import numpy as np
from Setup import setup
setup()
from funkcje_pomocnicze import P2, P4, P8, read_instance
from IPython.utils.io import capture_output
from tqdm import tqdm
gaps = {}
times = {}
is_automated_run = "PYTEST_CURRENT_TEST" in os.environ
instances = [P2, P4] if is_automated_run else [P2, P4, P8]
beam_widths = [10, 100] if is_automated_run else [10, 100, 500]
with tqdm(total=len(instances) * len(beam_widths), desc="zbieranie danych") as pbar:
for instance in instances:
for m in beam_widths:
J, h = read_instance(instance.path, convention="dwave")
with capture_output() as captured:
start = time.time()
solutions = branch_and_bound(J, h, m)
solution = solutions[0]
end = time.time()
gap = (instance.best_energy - solution.energy)/instance.best_energy
gaps[(instance.name, m)] = gap * 100
times[(instance.name, m)] = end - start
pbar.update(1)
instancje = [instance.name for instance in instances]
x = np.arange(len(instancje))
width = min(0.8 / max(len(beam_widths), 1), 0.25)
offsets = np.linspace(-width * (len(beam_widths) - 1) / 2, width * (len(beam_widths) - 1) / 2, len(beam_widths))
fig, ax = plt.subplots(figsize=(12, 6))
bars_by_width = []
for offset, m in zip(offsets, beam_widths):
values = [gaps[(instance_name, m)] for instance_name in instancje]
bars = ax.bar(x + offset, values, width, label=f"{m} trzymanych stanów")
bars_by_width.append(bars)
ax.set_ylim((-3, 30))
ax.set_xlabel("Rozmiar instancji")
ax.set_xticks(x)
ax.set_xticklabels(instancje, rotation=45)
ax.set_ylabel("Luka energetyczna [%]")
ax.set_title("Luka energetyczna dla różnych rozmiarów instancji")
ax.axhline(0, color='black', linestyle='--')
ax.legend()
for bars in bars_by_width:
for bar in bars:
if bar.get_height() == 0:
center = bar.get_x() + bar.get_width()/2
ax.plot(center, 0, marker='*', markersize=10, color=bar.get_facecolor())
plt.grid(axis="y")
plt.show()
fig2, ax2 = plt.subplots(figsize=(12, 6))
for offset, m in zip(offsets, beam_widths):
values = [times[(instance_name, m)] for instance_name in instancje]
ax2.bar(x + offset, values, width, label=f"{m} trzymanych stanów")
ax2.set_xlabel("Rozmiar instancji")
ax2.set_xticks(x)
ax2.set_xticklabels(instancje, rotation=45)
ax2.set_ylabel("Czas [s]")
ax2.set_title("Czas do znalezienia rozwiązania dla różnych rozmiarów instancji")
ax2.legend()
plt.grid(axis="y")
plt.show()
zbieranie danych: 100%|██████████| 9/9 [06:50<00:00, 45.64s/it]
