Wir suchen eine optimale Lösung für ein allgemeines Lineares Programm:

• Die Nebenbedingungen bilden immer ein Polyeder
• Die beste Lösung ist immer eine der Ecken

Können wir einfach alle Ecken ausprobieren?

Wieviele Ecken hat dieses Polyeder?

\[ \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} x \leq \begin{bmatrix} 1 \\ 1 \end{bmatrix}\\ x \geq 0 \] 4 Ecken

Und dieses?

8 Ecken

Anzahl Ecken unseres Würfels verdoppelt sich mit jeder zusätzlichen Dimension.

Für größere Anzahl an Variablen können wir in keinem Fall alle Ecken ausprobieren.

Iterativer Ansatz: Wir starten in einer Ecke und suchen nach der nächsten, mit der wir uns verbessern können.

Unser Beispielproblem:

Wo starten wir?

Triviale Startlösung: $x = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$

In welche Richtung wandern wir?
Erhöhen wir $x_1$ oder $x_2$?

Zielfunktion war: $\max 0.6 x_1 + 0.5 x_2$.

Ein Schritt Richtung $x_1$ bringt erstmal mehr als ein Schritt Richtung $x_2$.

Wie groß können wir $x_1$ maximal machen?

• $0.8 x_1 + 0.5 x_2 \leq 1500$
• $0.2 x_1 + 0.5 x_2 \leq 900$

$x_1 = \min (\frac{1500}{0.8}, \frac{900}{ 0.2})$
$ = \min (1875, 4500) = 1875$

Neue Lösung: $x = \begin{bmatrix} 1875 \\ 0 \end{bmatrix}$

Zielfunktionswert: $0.6 \cdot 1875 = 1125$

Können wir die Lösung verbessern, indem wir nun $x_2$ erhöhen?

Nicht ganz so einfach, wir dürfen ja nicht aus dem Polyeder rauslaufen.

Wenn wir $x_2$ erhöhen, müssen wir $x_1$ verringern. Nur um wieviel?

$0.8 x_1 + 0.5 x_2 = 1500$

Wenn wir $x_2$ um $1$ erhöhen, verringert sich $x_1$ um $\frac{0.5}{0.8} = 0.625$.

Müssen wir nicht auch an die Nebenbedingung $0.2 x_1 + 0.5 x_2 \leq 900$ denken?

$0.2 \cdot 1875 = 375 < 900$
Die Nebenbedingung ist aktuell noch nicht "aktiv".

Lohnt sich dieser Trade-off?

Zielfunktion war: $\max 0.6 x_1 + 0.5 x_2$.

$x_2$ um $1$ erhöhen, verändert den Zielfunktionswert um
$0.6 \cdot (- 0.625) + 0.5 \cdot 1 = 0.125$

$x_2$ erhöhen macht also Sinn.

Wie weit können wir nun $x_2$ erhöhen?

So lange, bis die nächste Nebenbedingung greift.

$0.2 x_1 + 0.5 x_2 \leq 900$ $x_1 \geq 0$

Wenn wir $x_2$ erhöhen, bis $x_1 = 0$ greift, dann haben wir: \[ 0.8 \cdot 0 + 0.5 \cdot x_2 = 1500 \]

$x_2 =\frac{1500}{0.5} = 3000$.

Die Lösung verletzt
$0.2 \cdot 0 + 0.5 \cdot 3000 = 1500 \not\leq 900$

Was ist die Lösung, wenn wir $x_2$ erhöhen, bis $0.2 x_1 + 0.5 x_2 \leq 900$ greift?

Dann gilt:
$ 0.8 \cdot x_1 + 0.5 \cdot x_2 = 1500\\ 0.2 \cdot x_1 + 0.5 \cdot x_2 = 900 $

bzw:
$ \begin{bmatrix} 0.8 & 0.5 \\ 0.2 & 0.5 \end{bmatrix} x = \begin{bmatrix} 1500 \\ 900 \end{bmatrix} $

Um so ein Gleichungssystem zu lösen gibt es mehrere Möglichkeiten. Zum Beispiel:

• Gauss-Verfahren
• Inverse Matrix berechnen
• Matrix Dekomposition
• Iterative Approximationen

Wenn wir die Matrix $A$ invertieren, bekommen wir die Lösung:

$A x = b \quad \Leftrightarrow\quad x = A^{-1} b$

Die inverse Matrix können wir uns zum Beispiel mit dem Python Numpy Package berechnen lassen:

          import numpy as np

          A = np.array([[0.8, 0.5], [0.2, 0.5]])
          A_inv = np.linalg.inv(A)

          print(A_inv)
        

Damit kommen wir an die Lösung: \[ x = \begin{bmatrix} \frac{5}{3} & -\frac{5}{3} \\ -\frac{2}{3} & \frac{8}{3} \end{bmatrix} \begin{bmatrix} 1500 \\ 900 \end{bmatrix} = \begin{bmatrix} 1000 \\ 1400 \end{bmatrix} \]

Diese erfüllt auch die andere Nebenbedingung:
$x_1 \geq 0$.

Unsere neue Lösung: $x = \begin{bmatrix} 1000 \\ 1400 \end{bmatrix}$

Können wir uns noch weiter verbessern?

Wir könnten $x_2$ noch weiter erhöhen.

Wenn wir $x_2$ weiter erhöhen, müssen wir $x_1$ weiter verringern: \[ 0.2 \cdot x_1 + 0.5 \cdot x_2 = 900 \]

Jeder weitere Schritt Richtung $x_2$ verringert $x_1$ um
$\frac{0.5}{0.2}=2.5$

Wenn wir einen weiteren Schritt Richtung $x_2$ machen, verändert sich also die Zielfunktion um: \[ 0.6 \cdot (-2.5) + 0.5 \cdot 1 = -1.0 \]

Wie können uns also in der Richtung nicht weiter verbessern.

Das Verfahren wollen wir jetzt systematischer machen, um mit beliebigen Problemen umgehen zu können.

Wir gehen von einem beliebigen LP in kanonischer Form aus: \[ \begin{aligned} \max \quad & c^T x \\ \text{subject to} \quad & Ax \leq b \\ & x \geq 0 \end{aligned} \]

Wenn wir ein lineares Gleichungssystem $Sx = b$ haben, mit $S \in \mathbb{R}^{n \times n}$, dann können wir das mit der Inversen lösen: \[ x = S^{-1} b \]

Wichtig sind hier Gleichungen und dass die Matrix quadratisch ist.

Eine Ungleichung können wir mit einer zusätzlichen Variablen in eine Gleichung verwandeln:

\[ a_{i1} x_1 + \dots + a_{in} x_n \leq b_i \Leftrightarrow\\ \quad a_{i1} x_1 + \dots + a_{in} x_n + z_i = b_i, \quad z_i \geq 0 \]

"Slackvariable"

Jede Slackvariable gibt an, wie weit wir von der Grenze der Ungleichung weg sind.

Wir können unsere Nebenbedingungen also als Gleichungen formulieren:

Unsere Nebenbedingungen sind nun:

$S$ ist nicht quadratisch, sondern hat mehr Spalten als Zeilen.

Wenn wir eine Teilmenge $B$ an Spalten auswählen, die so groß ist, wie die Anzahl an Zeilen, dann ist die Teilmatrix $S_B$ quadratisch.

Beispiel

$\begin{bmatrix}x \\ z \end{bmatrix} = \begin{bmatrix} 1875 \\ 0 \\ 0 \\ 525 \end{bmatrix}$

In dieser Lösung produzieren wir also 1875g von Teig A und lassen 525g Schokolade übrig.

Die Auswahl an Spalten $B$ nennen wir Basis.

Die Variablen der Basis sind diejenigen, die einen Wert größer 0 bekommen.

Alle Spalten, die nicht Teil der Basis sind, bezeichnen wir mit $N$.

Alle Variablen, die einer Spalte aus $N$ entsprechen, sind in der Lösung 0.

Jede Basis entspricht einer Ecke vom Polyeder der Nebenbedingungen.

Das geometrische Problem entspricht nun der Frage nach welchen Variablen wir das Gleichungssystem lösen.

Die aktuelle Lösung $\bar b = S_B^{-1} b$ nennen wir Basislösung.

Der Schritt zur nächsten Ecke besteht nun aus folgenden Operationen:

• Entscheide welche neue Spalte zur Basis dazukommt.
• Entscheide welche Spalte aus der Basis entfernt wird.

Diese Entscheidung nennt sich eine "Pivot-Operation"

Neue Basisspalte auswählen:

Für jede Nicht-Basisvariable ermitteln wir, wieviel die Zielfunktion bei einem Schritt in die entsprechende Richtung steigt:

\[ \bar c = c_N - S_N^T {S_B^{-1}}^T c_B \]

          import numpy as np

          Sb = np.array([[0.8, 0], [0.2, 1]])
          Sn = np.array([[0.5, 1], [0.5, 0]])
          cn = np.array([0.5, 0])
          cb = np.array([0.6, 0])

          cn - Sn.T.dot(np.linalg.inv(Sb).T.dot(cb))
        

Als Ergebnis erhalten wir \[ \begin{bmatrix} 0.125 \\ -0.75 \end{bmatrix} \]

$x_2$ um 1 steigern, steigert den Zielfunktionswert um $0.125$.
$z_1$ steigern verschlechtert dagegen die Zielfunktion.

Diese Werte $\bar c$ werden reduzierte Kosten genannt.

Wenn alle reduzierten Kosten negativ sind, haben wir eine optimale Lösung gefunden, da keine Verbesserung mehr möglich ist.

Wir wählen die Variable $i \in N$ mit den größten reduzierten Kosten aus. (also $x_2$)

Das ist nicht notwendigerweise die beste Wahl. Es gibt ausgefeiltere Strategien für die Wahl der neuen Basisspalte.

Nun müssen wir uns eine Variable $j \in B$ aussuchen, welche die Basis verlässt.

Voraussetzung: Keine Nebenbedingung darf verletzt werden.

Also darf keine Variable negativ werden.

Sei $\bar S = S_B^{-1} S_N$.

Wenn wir Spalte $j$ aus der Basis entfernen, dann bekommt die neue Variable den Wert: \[ \frac{\bar b_j}{\bar S_{j i}} \]

Die aktuelle Basislösung $\bar b$ ist positiv. Damit der neue Wert nicht negativ wird, muss also $\bar S_{j i} > 0$ sein.

Wir dürfen mit der neuen Variable nur minimal weit laufen, um nicht über die nächste Ecke hinauszuschießen.

Wähle $j$ mit kleinstem nicht negativen neuen Variablenwert: \[ j = \operatornamewithlimits{arg\, min}_{j \in B, \bar S_{j i} > 0} \frac{\bar b_j}{\bar S_{j i}} \]

          import numpy as np

          Sb = np.array([[0.8, 0], [0.2, 1]])
          Sn = np.array([[0.5, 1], [0.5, 0]])
          b = np.array([1500, 900])

          Sb_inv = np.linalg.inv(Sb)
          basic_solution = Sb_inv.dot(b)
          S_bar = Sb_inv.dot(Sn)

          new_values = basic_solution / S_bar[:, 0] # Beware of division by zero
          j = np.argmin(new_values[new_values > 0])
        

In dem Beispiel, kann $x_2$ als neuen Wert entweder $3000$ oder $1400$ annehmen.

Aktuelle Basis war $B = [0, 3]$. Also muss $z_2$ rausfliegen.

Neue Basis ist $B = [0, 1]$.

$\bar b = S_B^{-1} b = \begin{bmatrix} 1000 \\ 1400 \end{bmatrix}$

Um zu Prüfen, ob es noch weitergeht, berechnen wir wieder die reduzierten Kosten:

\[ \bar c = c_N - S_N^T {S_B^{-1}}^T c_B = \begin{bmatrix} -\frac{2}{3} \\ -\frac{1}{3} \end{bmatrix} \]

Wir haben also eine optimale Lösung gefunden.

Die reduzierten Kosten sagen uns aber noch mehr:

• Wenn wir den Abstand zur ersten Ungleichung um 1 erhöhen, sinkt die Zielfunktion um $\frac{2}{3}$.
• Umgekehrt: Könnten wir über die Ungleichung um 1 hinaus laufen, würde der Zielfunktionswert um $\frac{2}{3}$ steigen.

          import numpy as np

          Sb = np.array([[0.8, 0.5], [0.2, 0.5]])
          cb = np.array([0.6, 0.5])
          b_bar = np.linalg.inv(Sb).dot(np.array([1501, 900]))
          cb.dot(b_bar)

          # => 1300.6666666666665
        

Ein Extra Gramm Mehl bringt also $0.667$ Cent mehr, ein Extra Gramm Schokolade nur $0.333$ Cent.

Der Mehrwert steigt linear, bis die Lösung ungültig wird und wir eine neue Ecke suchen müssen.

Wir haben nun schon die wichtigsten Schritte für einen funktionierenden Simplex-Algorithmus.

Ist dieser denn besser, als alle Ecken auszuprobieren?

Schlechte Nachricht: Im schlimmsten Fall muss der Algorithmus alle Ecken erkunden. (Klee-Minty cube)

In der Praxis passiert das aber nicht. Daher funktioniert der Algorithmus fasst immer gut.

Allgemein müssen wir uns noch um ein paar Sonderfälle kümmern, damit beliebige LPs gerechnet werden können.

Wir haben angenommen, dass wir zu Beginn eine gültige Lösung haben.

In unserem Fall war das $x=0$, also die Basis nur aus Slack Variablen.

Es kann aber sein, dass $x=0$ nicht gültig ist:

(Wir müssen mindestens 400g von Teig A herstellen)

Wie kommen wir nun an eine initial gültige Lösung?

Wir können eine Hilfvariable $z' \geq 0$ einführen mit $-x_1 - z' = -400$.

Dann suchen wir zuerst eine Lösung, die nur die Werte von $z'$ minimiert.

Die Suche nach der gültigen Startlösung nennt sich Phase 1 des Simplex Algorithmus und ist potentiell genauso aufwendig wie Phase 2.

Es kann sein, dass wir keine gültige Startlösung finden können (also $z' > 0$ bleibt).

Dann ist das LP unlösbar (infeasible).

Beispiel

Es gibt noch einen weiteren Spezialfall:

Die Lösung ist unbeschränkt (unbounded)

Den Fall können wir im Simplex bemerken, wenn wir keine Spalte zum verlassen der Basis finden, weil alle $\bar S_{j i} < 0$ sind.

In der Praxis können "degenerierte" Ecken Probleme machen:

Hier kann es zu Endlosschleifen kommen. Durch eine gute Pivot Regel und numerische Tricks lässt sich das in der Praxis vermeiden.

Der teuerste Schritt beim Simplex ist die Berechnung von $S_B^{-1}$.

Der Trick ist, die Inverse nicht jeden Schritt neu zu berechnen, sondern nur updates zu machen.

Zum Beispiel durch LU Dekomposition und Forrest-Tomlin Updates

Lower Triangular Matrix (L)

Eine (quadratische) Matrix, die nur Einträge auf oder unter der Diagonalen hat:

\[ L = \begin{bmatrix} l_{11} & 0 & 0 & \dots & 0 \\ l_{21} & l_{22} & 0 & \dots & 0 \\ l_{31} & l_{32} & l_{33} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & l_{n3} & \dots & l_{nn} \end{bmatrix} \]

Angenommen wir wollen $L x = b$ lösen.

\[ \begin{align*} &l_{11} x_1 &= b_1 \\ &l_{21} x_1 + l_{22} x_2 &= b_2 \\ &l_{31} x_1 + l_{32} x_2 + l_{33} x_3 &= b_3 \\ &\vdots \\ &l_{n1} x_1 + l_{n2} x_2 + l_{n3} x_3 + \dots + l_{nn} x_n &= b_n \end{align*} \]

Nach $x_1$ lösen ist einfach: \[ x_1 = \frac{b_1}{l_{11}} \]

Wenn wir $x_1$ kennen ist auch nach $x_2$ lösen einfach:

Analog können wir das gesamte Gleichungssystem schnell lösen.

Eine LU-Dekomposition einer Matrix ist die Aufteilung der Matrix in eine Lower Triangular Matrix (L) und eine Upper Triangular Matrix (U)

\[ \begin{bmatrix} a_{11} & a_{12} & a_{13} & \dots & a_{1n} \\ a_{21} & a_{22} & a_{23} & \dots & a_{2n} \\ a_{31} & a_{32} & a_{33} & \dots & a_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & a_{n3} & \dots & a_{nn} \end{bmatrix} = \begin{bmatrix} l_{11} & 0 & 0 & \dots & 0 \\ l_{21} & l_{22} & 0 & \dots & 0 \\ l_{31} & l_{32} & l_{33} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & l_{n3} & \dots & l_{nn} \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13} & \dots & u_{1n} \\ 0 & u_{22} & u_{23} & \dots & u_{2n} \\ 0 & 0 & u_{33} & \dots & u_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & u_{nn} \end{bmatrix} \]

Das Gleichungssystem $A x = b$ können wir nun schnell lösen:

• $Ly = b$ ist einach zu lösen.
• Analog ist $Ux = y$ einfach zu lösen.
• So bekommen wir die Lösung für $LU x = b$
  und damit für $Ax = b$.

Wenn wir erstmal eine LU Dekomposition haben und in der Original Matrix ändert sich nur eine Zeile und eine Spalte, so kann man die Dekomposition mit wenig Aufwand updaten.

Leistungsfähige LP Solver:

• COIN-OR CLP (Open Source)
• CPLEX (kommerziell)
• GLPK (Open Source)
• GLOP (Open Source)
• Gurobi (kommerziell)
• SoPlex (Open Source)