Beispiel: Knapsack Problem

• Wir wollen einen Rucksack mit möglichst wertvollen Gegenständen füllen
• Aber: Wir können nur eine bestimmte Menge an Gewicht tragen.

• Wir haben $n$ Gegenstände.
• Jeder Gegenstand $i$ hat einen Wert $c_i$.
• Jeder Gegenstand $i$ hat ein Gewicht $w_i$.
• Wir können nur eine bestimmte Menge an Gewicht $W$ tragen.

Vollständige ILP Formulierung:

Beispiel

W = 10

          import mip

          model = mip.Model()

          x_A = model.add_var(var_type=mip.BINARY)
          x_B = model.add_var(var_type=mip.BINARY)
          x_C = model.add_var(var_type=mip.BINARY)

          model.objective = mip.maximize(20 * x_A + 15 * x_B + 7 * x_C)
          model.add_constr(8 * x_A + 7 * x_B + 3 * x_C <= 10)

          model.optimize()

          print(x_A.x)
          print(x_B.x)
          print(x_C.x)
        

Anders als bei kürzesten Pfaden ist die Ganzzahligkeit essentiell.

Fraktionale Lösung:

• Sortiere Gegenstände nach $\frac{\text{Wert}}{\text{Gewicht}}$.
• Fülle den Rucksack absteigend auf. Der letzte Gegenstand wird teilweise genommen.

• Wähle A zu 100%: Restkapazität: 2
• Wähle C zu 66%: Restkapazität: 0

Gesamtwert: 24.666

Ein gängiges Pattern beim Modellieren: Binäre Variablen, die eine Zuordnung modellieren.

Das Knapsack Problem ist nur "schwach" NP-schwer.

Wenn alle Gewichte ganzzahlig (mit nicht zu großen Zahlen) sind, können wir das per Dynamic-Programming effizient lösen.

Dynamic-Programming Lösung:

Wir erstellen eine Tabelle: Was ist die beste Lösung für alle Gewichte bis W, mit allen Gegenständen bis zum i-ten?

Wir nehmen entweder:

• Den Eintrag aus der letzten Zeile bei gleichem Gewicht
• Den Eintrag aus der letzten Zeile minus das aktuelle Gewicht + aktueller Wert

\[ \text{opt}(i, w) = \max(\text{opt}(i-1, w), \text{opt}(i-1, w - w_i) + c_i) \]

Dynamic-Programming = Rekursion + Cache

          values = [20, 15, 7]
          weights = [8, 7, 3]

          import numpy as np

          def solve_knapsack(values, weights, limit):
              assert len(values) == len(weights)
              solution = np.zeros(shape=(len(values)+1, limit+1))

              for item_idx, (value, weight) in enumerate(zip(values, weights)):
                  for w in range(limit+1):
                      if w >= weight:
                          solution[item_idx+1, w] = max(
                              solution[item_idx, w],
                              solution[item_idx, w - weight] + value
                          )
                      else:
                          solution[item_idx+1, w] = solution[item_idx, w]
              return solution[len(values), limit]

          solve_knapsack(values, weights, 10)
        

Diese Lösung funktioniert nicht gut, wenn wir zu große Gewichte haben.

In der Regel ist das kein Problem.

https://xkcd.com/287/

          import mip

          prices = {
              "Mixed Fruit": 2.15,
              "French Fries": 2.75,
              "Side Salad": 3.35,
              "Hot Wings": 3.55,
              "Mozzerella Sticks": 4.20,
              "Sampler Plate": 5.80
          }

          m = mip.Model("restaurant order")

          variables = {
              item: m.add_var(var_type=mip.INTEGER) for item in prices
          }

          # m.objective = mip.maximize(sum(variables.values()))

          m.add_constr( sum(variables[item] * prices[item] for item in prices) == 15.05 )

          m.optimize()

          for item, var in variables.items():
              if var.x > 0.5:
                  print(f"{item}: {var.x}")
        

Erweiterung: Bin Packing

Wir wollen Gegenstände unterschiedlicher Größe ($s_i$) einpacken. Eine Bin hat Größe $B$.
Wie viele Bins brauchen wir mindestens um alle Gegenstände zu verpacken und wie packen wir diese?

Wie modellieren wir das?

Idee: Variable $x_{i,j}$ ob Gegenstand $i$ in Bin $j$ eingepackt wird.

Wir wissen aber gar nicht, wie viele Bins wir am Ende brauchen.

Idee: Wir nehmen mehr Bins als wir brauchen. Zum Beispiel eine pro Gegenstand.

Wir nehmen noch eine Variable $y_j$, die angibt, ob wir die Bin verwenden.

Zielfunktion: Minimiere Anzahl verwendeter Bins

Constraints: Kapazität von verwendeten Bins darf nicht überschritten werden.

Hier haben wir ein häufiges Pattern: Die binäre Variable $y_j$ schaltet die Kapazität der Constraint frei.

Constraints: Jeder Gegenstand muss eingepackt werden.

\[ \begin{align*} \min \quad & \sum_{j=1}^n y_j \\ \text{s.t.} \quad & \sum_{i=1}^n s_i x_{i,j} \leq B y_j \quad & \forall j \\ & \sum_{j=1}^n x_{i,j} = 1 \quad & \forall i \\ & x_{i,j} \in \{0, 1\} \quad & \forall i, \forall j \\ & y_i \in \{0, 1\} \quad & \forall i \end{align*} \]

Nächste Aufgabe: Wir wollen einen Stundenplan erstellen

• Wir haben eine Menge von Kursen $K$
• Wir haben eine Menge von Räumen $R$
  • Jeder Kurs $k \in K$ hat eine Menge zulässige Räume $R(k)$
• Wir haben eine Menge von Zeitfenstern $T$
  • Jeder Kurs $k \in K$ hat eine Menge an zulässigen Zeitfenstern $T(k)$
• Wir haben paarweise Konflikte $C$. $(k_1, k_2) \in C$ sagt, dass $k_1$ und $k_2$ nicht gleichzeitig stattfinden dürfen.

Ziel: Minimiere die Anzahl an Konflikten.

Variablen

• Eine Variable $x_{k, t, r} \in \{0, 1\}$ für jeden Kurs und zulässigen Raum und Zeitslot.
• Eine Variable $y_{k_1, k_2} \in \{0, 1\}$ für jeden Konflikt.

Gängiges Pattern: Eine binäre Variable mit vielen Indizes, die kodiert, welche Ressource einer anderen Ressource zugewiesen wird.

Modellierung mit vielen binären Variablen ist oft besser, als mit wenigen ganzzahligen.

Zielfunktion:

Wir könnten hier auch erweitern, dass manche Konflikte teurer sind als andere.

Constraints: Jeder Kurs braucht einen Zeitslot und einen Raum.

Constraints: Kein Raum darf zu einem Zeitpunkt doppelt gebucht werden.

Constraints: Wir dürfen einen Konflikt nur machen, wenn wir die entsprechende Variable setzen.

\[ \sum_{r \in R(k_1)} x_{k_1, t, r} + \sum_{r \in R(k_2)} x_{k_2, t, r} \leq 1 + y_{k_1, k_2} \quad \forall (k_1, k_2) \in C, \forall t \in T \]

Wenn wir $y$ weglassen bekommen wir einen "harten" Konflikt.

Pattern: Wir haben wieder eine Constraint, die wir mit einer Variable an, oder aus schalten können.

Dadurch wird die Bedingung optional und wir können ihre Verletzung in die Zielfunktion übernehmen.

Beachte: Hier können wir die Nebenbedingung maximal um 1 überschreiten.

Nächstes Beispiel: Machine Scheduling

• Wir haben eine Maschine und $n$ Aufgaben.
• Jede Aufgabe $i$ hat eine Dauer $p_i$.
• Jede Aufgabe $i$ hat eine Deadline $d_i$.

Wir wollen die Anzahl an Aufgaben minimieren, die ihre Deadline reissen.

Variablen:

Was wir suchen ist vor allem die Reihenfolge in der wir die Aufgaben bearbeiten. Für jede der $n$ Aufgaben gibt es $n$ mögliche Positionen, an der sie ausgeführt werden kann.

• Eine Variable $x_{i, j} \in \{0, 1\}$: Aufgabe $i$ wird an $j$-ter Stelle ausgeführt.
• Eine Variable $y_j \in \{0, 1\}$: Der Job, der an Position $j$ ausgeführt wird, hat seine Deadline verpasst.
• Eine Variable $t_j \geq 0$: Der Zeitpunkt, an dem die Aufgabe an der $j$-ten Stelle begonnen wird.

Wir haben hier ein echtes MILP, da wir ganzzahlige mit fraktionalen Variablen mischen.

Zielfunktion:

Constraints: Jede Aufgabe wird genau einmal ausgeführt.

Constraints: An jeder Position wird genau eine Aufgabe ausgeführt.

Constraints: Jede Aufgabe kann erst nach der vorherigen starten.

Mit $\sum_{i=1}^n p_i x_{i,j}$ nehmen wir nur genau die Verarbeitungszeit, die an der Position $j$ relevant ist.

Constraints: Eine Aufgabe darf ihre Deadline nur reissen, wenn die entsprechende $y$ Variable gesetzt ist:

\[ t_j + \sum_{i=1}^n p_i x_{i,j} \leq \sum_{i=1}^n d_i x_{i,j} + M y_j \quad \forall j \in \{1, \dots, n\} \]

Die gleiche Variable kommt hier zweimal vor. Das macht die Nebenbedingung verständlicher.

Das kann auch so formuliert werden:

\[ t_j + \sum_{i=1}^n (p_i - d_i) x_{i,j} \leq M y_j \quad \forall j \in \{1, \dots, n\} \]

Oder, sogar so (da wir wissen dass sich die $x$ zu 1 addieren):

\[ \sum_{i=1}^n (t_j + p_i - d_i) x_{i,j} \leq M y_j \quad \forall j \in \{1, \dots, n\} \]

Das zweite Pattern ist der Term $M y_j$. Hier schalten wir die Nebenbedingung aus, wenn $y_j = 1$.

$M$ ist eine "große" Zahl, so dass die Bedingung einfach zu erfüllen wird.

Mit diesen "Big-M" Constraints können wir gut ein
"if-then" modellieren.

Sie machen das Modell aber auch schwerer lösbar.

Beispiel: $M=10000$. Ein Job wird so geplant, dass er $10$ Zeiteinheiten zu spät ist.

Mit $y=0.001$ wird die Nebenbedingung in der LP-Relaxierung ausgeschaltet, kostet in der Zielfunktion aber nicht viel.

Wenn wir solche Constraints brauchen, sollten wir immer versuchen, das $M$ so klein wie möglich zu wählen.

Hier: Was ist die Zeit, die ein Job im schlimmsten Fall zu spät sein kann?

Der Job mit der frühesten Deadline wird als allerletzter bearbeitet.

Noch besser: Für jede Position $j$ überlegen wir uns, wie viel der Job an der Position maximal zu spät sein kann:

$P_j$: Sei die Bearbeitungszeit der $j$ längsten Jobs

Wir können mit diesem Pattern auch ein "Entweder-Oder" modellieren:

Eine der Bedingungen ist aus und die andere eingeschaltet.

Beispiel: Ich möchte, dass die Variable $x$ entweder kleiner als 10 oder größer als 100 ist:

Wenn wir die Zeit in ein Raster diskretisieren können, erhalten wir hier potentiell ein besser lösbares Modell.
(ohne Big-M Constraints)

In der Praxis gibt es diverse Varianten dieses Problems:

• Minimiere die Menge an Minuten, welche die Aufgaben zu spät sind.
• Aufgaben können erst ab einem bestimmten Release Zeitpunkt bearbeitet werden.
• Wir haben mehrere Maschinen.
• Aufgaben können unterbrochen und "geparkt" werden.
• ...

Ein großer Vorteil von MILP Modellierungen ist, dass wir oft neue Anforderungen in existierende Modelle integrieren können.

Das ist ein großer Vorteil im agilen Entwicklungsprozess.

Modellierungs-Tipp:

Das Problem so diskretisieren, dass wir binäre Variablen verwenden können.

Modellierungs-Tipp:

Nebenbedingungen verwenden, wo alle Variablen nur aufaddiert werden:

Modellierungs-Tipp:

Nebenbedingungen verwenden, wo alle konstanten Terme eine 1 sind:

Modellierungs-Tipp:

Am besten beide Tipps von gerade zusammen:

Was machen wir, wenn wir eine nicht lineare Funktion abbilden wollen?

Beispiel:

• Wir produzieren zwei Güter A und B.
• Wir können maximal 300 Einheiten insgesamt produzieren.
• Produkt B bringt 9 Euro je Einheit.
• Produkt A bringt 10 Euro für die ersten 100 Einheiten, 8 Euro für die nächsten 100 Einheiten und anschließend 6 Euro je Einheit.

Mit konstanten Einnahmen je Einheit A wäre es einfach:

          import mip

          model = mip.Model()

          prod_A = model.add_var(lb=0)
          prod_B = model.add_var(lb=0)

          model.add_constr(prod_A + prod_B <= 300)
          model.objective = mip.maximize(10 * prod_A + 9 * prod_B)

          model.optimize()
        

Eigentlich wollen wir folgendes lösen:

\[ \begin{align*} \max \quad & f(x_A) + 9 x_B \\ \text{s.t.} \quad & x_A + x_B \leq 300 \\ & x_A, x_B \geq 0 \end{align*} \]

mit

\[ f(x) = \begin{cases} 10 x & \text{if } x \leq 100 \\ 1000 + 8 (x - 100) & \text{if } 100 < x \leq 200 \\ 1800 + 6 (x - 200) & \text{if } 200 < x \end{cases} \]

Die Einnahmen $f$ für A sind eine
stückweise lineare Funktion
(piecewise linear function)

Eine stückweise lineare Funktion hat:

• $n$ lineare Segmente
• Break points $b_1, \ldots, b_{n+1}$
• An break point $i$ den Wert $f(b_i)$

$x$ muss sich auf einem Segment zwischen zwei break points befinden:

Wenn wir das Segment kennen, können wir $x$ als lineare Kombination der break points darstellen:

Beispiel:

$x = 140$

Break points $b_2 = 100$ und $b_3 = 200$.

$x = 0.6 b_2 + 0.4 b_3 = 60 + 80$

Innerhalb des Segments ist $f$ linear:

\[ = z_k f(b_k) + z_{k+1} f(b_{k+1}) \]

Beispiel:

$x = 140$

$f(140) = f(60 + 80) = f(0.6 b_2 + 0.4 b_3)$

$= 0.6 f(b_2) + 0.4 f(b_3)$

$= 0.6 \cdot 1000 + 0.4 \cdot 1800 = 1320$

Wir können jetzt $f$ mithilfe der break points modellieren:

\[ \begin{align*} \max \quad & {\color{red} 0 z_1 + 1000 z_2 + 1800 z_3 + 2400 z_4} + 9 x_B \\ \text{s.t.} \quad & x_A + x_B \leq 300 \\ & {\color{red} x_A = 0 z_1 + 100 z_2 + 200 z_3 + 300 z_4} \\ & {\color{red} z_1 + z_2 + z_3 + z_4 = 1} \\ & x_A, x_B \geq 0 \\ & {\color{red} z_1, z_2, z_3, z_4 \geq 0 } \end{align*} \]

Aktuell kann das Modell die Segmente noch beliebig vermischen, z.B.:

$z_1 = 0.2 \quad z_2 = 0.3 \quad z_3 = 0 \quad z_4 = 0.5$

Wir müssen sicherstellen, dass nur genau ein Segment ausgewählt wird.

Füge eine binäre Variable $y_1, \ldots, y_n$ je Segment ein.

• Wähle nur ein Segment:
  $y_1 + y_2 + y_3 = 1$
• Verwende $z_1$ nur, wenn Segment $1$ aktiv ist:
  $z_1 \leq y_1$
• Verwende $z_2$ nur, wenn Segment $1$ oder $2$ aktiv ist:
  $z_2 \leq y_1 + y_2$
• Verwende $z_3$ nur, wenn Segment $2$ oder $3$ aktiv ist:
  $z_3 \leq y_2 + y_3$
• Verwende $z_4$ nur, wenn Segment $3$ aktiv ist:
  $z_4 \leq y_3$

\[ \begin{align*} \max \quad & 0 z_1 + 1000 z_2 + 1800 z_3 + 2400 z_4 + 9 x_B \\ \text{s.t.} \quad & x_A + x_B \leq 300 \\ & x_A = 0 z_1 + 100 z_2 + 200 z_3 + 300 z_4 \\ & z_1 + z_2 + z_3 + z_4 = 1 \\ & {\color{red} y_1 + y_2 + y_3 = 1 }\\ & {\color{red} z_1 \leq y_1 }\\ & {\color{red} z_2 \leq y_1 + y_2 }\\ & {\color{red} z_3 \leq y_2 + y_3 }\\ & {\color{red} z_4 \leq y_3 }\\ & x_A, x_B \geq 0 \\ & z_1, z_2, z_3, z_4 \geq 0 \\ & {\color{red} y_1, y_2, y_3 \in \{0, 1\}} \end{align*} \]

          import mip

          model = mip.Model(sense=mip.MAXIMIZE)

          prod_A = model.add_var(lb=0)
          prod_B = model.add_var(lb=0)

          model.add_constr(prod_A + prod_B <= 300)

          z_A1 = model.add_var(lb=0, obj= 0)
          z_A2 = model.add_var(lb=0, obj= 1000)
          z_A3 = model.add_var(lb=0, obj= 1800)
          z_A4 = model.add_var(lb=0, obj= 2400)
          y_A1 = model.add_var(var_type=mip.BINARY)
          y_A2 = model.add_var(var_type=mip.BINARY)
          y_A3 = model.add_var(var_type=mip.BINARY)

          model.add_constr(z_A1 + z_A2 + z_A3 + z_A4 == 1)
          model.add_constr(y_A1 + y_A2 + y_A3 == 1)
          model.add_constr(z_A1 <= y_A1)
          model.add_constr(z_A2 <= y_A1 + y_A2)
          model.add_constr(z_A3 <= y_A2 + y_A3)
          model.add_constr(z_A4 <= y_A3)
          model.add_constr(prod_A == z_A1 * 0 + z_A2 * 100 + z_A3 * 200 + z_A4 * 300)

          model.objective += 9 * prod_B

          model.optimize()

          print(prod_A.x)
          print(prod_B.x)
        

Die Abbildung der stückweise linearen Funktion ist zwar etwas umständlich, aber sehr modular.

Kein anderer Teil des Modells braucht Zugriff auf die Hilfsvariablen $z$ und $y$.

Was tun wir, wenn eine Funktion nicht stückweise linear ist?

Wir können die Funktion stückweise linear approximieren:

Je nach benötigter Genauigkeit erhalten wir mehr Segmente:

Wie gehen wir mit Unsicherheiten um?

Beispiel: Ein Energieunternehmen muss jedes Jahr entscheiden, wie viel Gas es zum aktuellen Preis kaufen möchte, ohne den zukünftigen Preis und Bedarf zu kenne.

• Dieses Jahr kostet eine Einheit Gas €5,-
• Der Bedarf dieses Jahr sind 100 Einheiten
• Wenn wir Gas lagern wollen müssen wir pauschal €100,- bezahlen
• Für jede gelagerte Einheit müssen wir €0.5,- bezahlen.

Im nächsten Jahr wissen wir jedoch noc nicht, wie hoch der Bedarf und der Preis ist.

Wie entscheiden wir die Einkaufsmenge dieses Jahr?

Idee: Unsicherheit in Szenarien verpacken.

Vorhersage guter Szenarien ist nicht einfach.

Benötigen gute Modelle / Analysen zur Vorhersage (z.B. mit Machine Learning)

Wie feingranular wollen wir die Szenarien unterteilen?

Beispiel:

Wettermodell sagt Mild, Kalt, Sehr kalt mit W.-keit 0.4, 0.4 und 0.2 vorher.

Bedarfsmodell sagt Durchschnitt basierend auf Wetter vorher.

Beispiel:

Historische Bedarfe

          demand_1 = np.mean(data[data < np.percentile(x, 33)])
          demand_2 = np.mean(data[data < np.percentile(x, 66)])
          demand_3 = np.mean(data[data >= np.percentile(x, 66)])
        

Gegeben die Szenarien, wie vewenden wir diese?

• Ein Teil des Modells bildet die Entscheidung im ersten Jahr ab.
• Im zweiten Jahr bekommen wir ein Teilmodell je Szenario.

Variablen:

• Eingekaufte Einheiten in Jahr 1: $x_\text{buynow} \geq 0$
• Eingelagerte Einheiten in Jahr 1: $x_\text{store} \geq 0$
• Lager verwendet: $x_\text{usestore} \in \{0, 1\}$
• Eingekaufte Einheiten in Jahr 2 je Szenario:
  $x_\text{buylater}^1, x_\text{buylater}^2, x_\text{buylater}^3 \geq 0$

Alles was wir jetzt einkaufen wir gelagert oder verbraucht:

\[ x_\text{buynow} = x_\text{store} + \text{demand now} \]

Wenn wir etwas lagern, müssen wir das Lager verwenden:

\[ x_\text{store} \leq x_\text{usestore} \cdot \text{max demand} \]

(Big-M constraint)

In jedem Szenario können wir entscheiden, wieviel noch gebraucht wird:

• \[ x_\text{store} + x_\text{buylater}^1 \geq \text{scenario 1 demand} \]
• \[ x_\text{store} + x_\text{buylater}^2 \geq \text{scenario 2 demand} \]
• \[ x_\text{store} + x_\text{buylater}^3 \geq \text{scenario 3 demand} \]

In jedem Szenario bezahlen wir nur den Erwartungswert:

\[ \min 5 x_\text{buynow} + 100 x_\text{usestore} + 0.5 x_\text{store} + \\ 0.4 \cdot 5 x_\text{buylater}^1 + 0.4 \cdot 7 x_\text{buylater}^2 + 0.2 \cdot 9 x_\text{buylater}^3 \]

          import mip

          storage_base_cost = 100
          storage_unit_cost = 0.5
          cost_now = 5
          demand_now = 100

          scenarios = [{"cost": 5, "demand": 100, "prob": 0.4},
                       {"cost": 7, "demand": 150, "prob": 0.4},
                       {"cost": 9, "demand": 200, "prob": 0.2}]

          model = mip.Model()

          buy_now = model.add_var(lb = 0)
          store_now = model.add_var(lb = 0)
          use_storage = model.add_var(var_type=mip.BINARY)

          model.add_constr(buy_now == demand_now + store_now)
          model.add_constr(store_now <= use_storage * max(s["demand"] for s in scenarios))

          model.objective = mip.minimize(buy_now * cost_now + storage_base_cost * use_storage + storage_unit_cost * store_now)

          buy_later_vars = []
          for scenario in scenarios:
              buy_later = model.add_var(lb = 0)
              buy_later_vars.append(buy_later)
              model.add_constr(store_now + buy_later >= scenario["demand"])
              model.objective += scenario["prob"] * scenario["cost"] * buy_later

          model.optimize()
        

Stochastische Modelle funktionieren auch für reine LPs
(Ganzzahligkeit ist nicht notwendig für Abbildung der Szenarien)