#!/usr/bin/env python3 # erzeuge_dualitaets_bruecke.py """ Erzeugt die Doppeltafel zum starken Dualitaetssatz: bilder_04/kap_lp_dualitaets_bruecke.svg ERGAENZT das vorhandene kap04_primal_dual.svg, ersetzt es nicht: Jenes stellt die beiden Probleme schematisch gegenueber - wer fragt was. Dieses zeigt ihre Geometrie und damit den Satz selbst. Links der primale Mengenraum (x1, x2): ein beschraenktes Polyeder, das Maximum in einer Ecke. Rechts der duale Preisraum (y1, y2): ein nach oben offener Bereich, das Minimum ebenfalls in einer Ecke. Die beiden Bereiche haben nichts miteinander gemein - andere Achsen, andere Einheiten, der eine beschraenkt, der andere nicht. Trotzdem steht an beiden Optima dieselbe Zahl. Das ist der Satz, und das ist der einzige Grund fuer die Doppeltafel. Der duale Raum hat drei Variablen; gezeichnet wird die Ebene y3 = 0. Das ist keine Vereinfachung, sondern folgt aus dem komplementaeren Schlupf: Das Marktlimit x1 <= 8 ist bei x1* = 7,5 nicht bindend, also ist y3* = 0. Der Selbsttest rechnet nach, dass das Minimum in dieser Ebene wirklich das Minimum des vollen Dualproblems ist - sonst zeigte die rechte Tafel einen zu grossen Wert und der optische Beweis waere keiner. Modell und Loesung stammen aus dem Kapitelprogramm Dualitaet_Nachweis.py; hier wird nichts nachgebaut und nichts eingetragen. Aufruf (aus dem Repository-Wurzelverzeichnis): python3 bilder_04/erzeuge_dualitaets_bruecke.py Benoetigt: numpy, scipy, matplotlib """ from __future__ import annotations import itertools import os import sys import numpy as np HIER = os.path.dirname(os.path.abspath(__file__)) BASIS = os.path.dirname(HIER) sys.path.insert(0, HIER) sys.path.insert(0, os.path.join( BASIS, "Operations_Research_mit_Python_Version_04_Programme")) from stil_04 import FARBEN, speichere # noqa: E402 import matplotlib.pyplot as plt # noqa: E402 from Dualitaet_Nachweis import (A_PRIMAL, B_PRIMAL, # noqa: E402 C_PRIMAL, loese) def ecken(A: np.ndarray, b: np.ndarray) -> np.ndarray: """Ecken des Bereichs {x >= 0 : A x <= b}, gegen den Uhrzeigersinn. Alle Paare von Randgeraden schneiden (die Achsen zaehlen mit), die unzulaessigen Schnittpunkte verwerfen, den Rest um den Schwerpunkt sortieren. Fuer zwei Dimensionen ist das vollstaendig.""" zeilen = np.vstack([A, [[-1.0, 0.0], [0.0, -1.0]]]) rechte = np.concatenate([b, [0.0, 0.0]]) punkte = [] for i, j in itertools.combinations(range(len(zeilen)), 2): M = zeilen[[i, j]] if abs(np.linalg.det(M)) < 1e-9: continue p = np.linalg.solve(M, rechte[[i, j]]) if np.all(zeilen @ p <= rechte + 1e-9): punkte.append(p) punkte = np.unique(np.round(punkte, 9), axis=0) mitte = punkte.mean(axis=0) winkel = np.arctan2(punkte[:, 1] - mitte[1], punkte[:, 0] - mitte[0]) return punkte[np.argsort(winkel)] def zeichne(x_stern, z_stern, y_stern, w_stern) -> None: figur, achsen = plt.subplots(1, 2, figsize=(10.2, 4.6)) beschriftung = f"{z_stern:,.0f} €".replace(",", ".") # --- links: der primale Mengenraum ------------------------------------ links = achsen[0] polygon = ecken(A_PRIMAL, B_PRIMAL) links.fill(polygon[:, 0], polygon[:, 1], color=FARBEN["haupt"], alpha=0.18, zorder=2) links.plot(np.append(polygon[:, 0], polygon[0, 0]), np.append(polygon[:, 1], polygon[0, 1]), "-", color=FARBEN["haupt"], linewidth=1.8, zorder=3) x1 = np.linspace(0, 11, 200) # Ein Niveau unter dem Optimum, das Optimum selbst, eines darueber - das # obere liegt sichtbar ausserhalb des Polyeders, und genau das ist die # Aussage: weiter geht es nicht. for niveau, stil in ((1500.0, ":"), (z_stern, "--"), (3000.0, ":")): # 150 x1 + 250 x2 = niveau links.plot(x1, (niveau - C_PRIMAL[0] * x1) / C_PRIMAL[1], stil, color=FARBEN["warnung"], linewidth=1.4, zorder=4) links.annotate("", xy=(4.2, 4.2), xytext=(2.4, 3.1), arrowprops=dict(arrowstyle="-|>", color=FARBEN["warnung"], linewidth=2.0), zorder=5) links.text(4.35, 4.25, r"$\mathbf{c} = (150, 250)^\top$", fontsize=8.5, color=FARBEN["warnung"]) links.plot(*x_stern, "*", color=FARBEN["gut"], markersize=17, zorder=6) links.annotate(f"$\\mathbf{{x}}^* = ({x_stern[0]:.1f}\\,;\\," f"{x_stern[1]:.0f})$\n$Z^* = {beschriftung}$", xy=x_stern, xytext=(4.55, 1.05), textcoords="data", fontsize=9, color=FARBEN["gut"], fontweight="bold", ha="left", zorder=7, arrowprops=dict(arrowstyle="-", color=FARBEN["gut"], linewidth=1.0, alpha=0.7)) links.set_xlim(0, 11) links.set_ylim(0, 9) links.set_xlabel("$x_1$ — Bots Typ A") links.set_ylabel("$x_2$ — Bots Typ B") links.set_title("primal: Mengen — maximieren, beschränkt", fontsize=10.5, color=FARBEN["text"]) # --- rechts: der duale Preisraum, Ebene y3 = 0 ------------------------ rechts = achsen[1] g1 = np.linspace(0, 55, 320) g2 = np.linspace(0, 55, 320) yy1, yy2 = np.meshgrid(g1, g2) # A^T y >= c, hier ohne y3: 2 y1 + 4 y2 >= 150 und 5 y1 + 6 y2 >= 250 zulaessig = ((A_PRIMAL[0, 0] * yy1 + A_PRIMAL[1, 0] * yy2 >= C_PRIMAL[0] - 1e-9) & (A_PRIMAL[0, 1] * yy1 + A_PRIMAL[1, 1] * yy2 >= C_PRIMAL[1] - 1e-9)) rechts.contourf(yy1, yy2, zulaessig.astype(float), levels=[0.5, 1.5], colors=[FARBEN["zweit"]], alpha=0.18, zorder=2) for spalte, farbe in ((0, FARBEN["zweit"]), (1, FARBEN["zweit"])): # A[0, spalte] * y1 + A[1, spalte] * y2 = c[spalte] rechts.plot(g1, (C_PRIMAL[spalte] - A_PRIMAL[0, spalte] * g1) / A_PRIMAL[1, spalte], "-", color=farbe, linewidth=1.8, zorder=3) for niveau, stil in ((1200.0, ":"), (w_stern, "--"), (3600.0, ":")): # 40 y1 + 60 y2 = niveau rechts.plot(g1, (niveau - B_PRIMAL[0] * g1) / B_PRIMAL[1], stil, color=FARBEN["warnung"], linewidth=1.4, zorder=4) rechts.plot(y_stern[0], y_stern[1], "*", color=FARBEN["gut"], markersize=17, zorder=6) rechts.annotate(f"$\\mathbf{{y}}^* = ({y_stern[0]:.2f}\\,;\\," f"{y_stern[1]:.2f})$\n$W^* = " f"{f'{w_stern:,.0f}'.replace(',', '.')}$ €", xy=(y_stern[0], y_stern[1]), xytext=(14, 16), textcoords="offset points", fontsize=9, color=FARBEN["gut"], fontweight="bold", zorder=7) rechts.annotate("zulässig: nach oben offen —\nteurer geht immer", xy=(0.04, 0.97), xycoords="axes fraction", fontsize=8.5, color=FARBEN["text"], va="top") rechts.set_xlim(0, 55) rechts.set_ylim(0, 55) rechts.set_xlabel("$y_1$ — Preis je vCPU-Stunde") rechts.set_ylabel("$y_2$ — Preis je GB RAM") rechts.set_title("dual: Preise — minimieren, unbeschränkt", fontsize=10.5, color=FARBEN["text"]) figur.suptitle(f"Zwei verschiedene Räume, ein gemeinsamer Wert: " f"{beschriftung}", fontsize=11.5, color=FARBEN["text"], y=1.01) speichere(figur, "kap_lp_dualitaets_bruecke") if __name__ == "__main__": primal, dual = loese() x_stern, z_stern = primal.x, -primal.fun y_stern, w_stern = dual.x, dual.fun print(f" primal: x* = ({x_stern[0]:.4f}, {x_stern[1]:.4f}) " f"Z* = {z_stern:.4f}") print(f" dual: y* = ({y_stern[0]:.4f}, {y_stern[1]:.4f}, " f"{y_stern[2]:.4f}) W* = {w_stern:.4f}") if abs(z_stern - w_stern) > 1e-6: raise SystemExit( f"Dualitaets-Bruecke: Z* = {z_stern:.6f}, W* = {w_stern:.6f}. " f"Die Tafel behauptet einen gemeinsamen Wert - er ist keiner.") # Die rechte Tafel zeigt die Ebene y3 = 0. Das ist nur zulaessig, wenn # das Optimum dort liegt. if abs(y_stern[2]) > 1e-9: raise SystemExit( f"Dualitaets-Bruecke: y3* = {y_stern[2]:.6f} ist nicht null. Die " f"gezeichnete Ebene y3 = 0 enthaelt das Optimum dann nicht.") print(f" starker Dualitätssatz bestätigt: " f"{z_stern:,.0f} €".replace(",", ".")) print(f" y3* = 0 — die gezeichnete Ebene enthält das Optimum") zeichne(x_stern, z_stern, y_stern, w_stern)