import sys
import subprocess
from importlib.metadata import version, PackageNotFoundError

REQUIRED_PULP = "3.3.2"
try:
    installed_pulp = version("pulp")
except PackageNotFoundError:
    installed_pulp = None

if installed_pulp != REQUIRED_PULP:
    subprocess.check_call([
        sys.executable, "-m", "pip", "install", "--quiet",
        "--disable-pip-version-check", "--no-cache-dir", "pulp==3.3.2"
    ])
if version("pulp") != REQUIRED_PULP:
    raise RuntimeError("PuLP 3.3.2 se nenainstalovalo. Zkontrolujte výpis instalace výše.")

loaded_pulp = sys.modules.get("pulp")
if loaded_pulp is not None and getattr(loaded_pulp, "__version__", None) != REQUIRED_PULP:
    raise RuntimeError(
        "PuLP 3.3.2 je nainstalované, ale v paměti zůstává jiná verze. "
        "V Colabu zvolte Runtime > Restart session (případně Restart runtime). "
        "V Jupyteru zvolte Kernel > Restart Kernel. Potom spusťte všechny buňky shora."
    )
print("Připraveno: PuLP 3.3.2 je nainstalované. Pokračujte dalšími buňkami.")

Připraveno: PuLP 3.3.2 je nainstalované. Pokračujte dalšími buňkami.

import pulp
import numpy as np
import matplotlib.pyplot as plt

if pulp.__version__ != "3.3.2":
    raise RuntimeError("Spusťte sekci 0 a postupujte podle pokynu k restartu.")
solver = pulp.PULP_CBC_CMD(msg=False)
tol = 1e-5

capacity = {("S", "A"): 500, ("A", "C"): 500, ("C", "T"): 650,
            ("S", "B"): 600, ("B", "D"): 600, ("D", "T"): 350,
            ("C", "D"): 120}
pos = {"S": (0, 0), "A": (1, 1), "B": (1, -1), "C": (2, 1),
       "D": (2, -1), "T": (3, 0)}
fig, ax = plt.subplots(figsize=(7, 3.5))
for (u, v), cap in capacity.items():
    x, y = pos[u]
    xx, yy = pos[v]
    ax.annotate("", pos[v], pos[u], arrowprops=dict(arrowstyle="->", shrinkA=13, shrinkB=13))
    ax.text((x+xx)/2 + (0.13 if x == xx else 0), (y+yy)/2+0.05, str(cap),
            ha="center", backgroundcolor="white")
for node, (x, y) in pos.items():
    ax.text(x, y, node, ha="center", va="center",
            bbox=dict(boxstyle="circle", facecolor="#daf0ed"))
ax.set(xlim=(-0.3, 3.3), ylim=(-1.4, 1.4), title="Kapacity spojů [zásilky/h]")
ax.axis("off")
plt.show()

def route(capacities, omit=None):
    model = pulp.LpProblem("parcel_flow", pulp.LpMaximize)
    flow = {(u, v): pulp.LpVariable(f"x_{u}_{v}", lowBound=0, upBound=limit)
            for (u, v), limit in capacities.items()}
    model += pulp.lpSum(flow[u, v] for u, v in flow if v == "T")
    nodes = {node for edge in flow for node in edge}
    for node in sorted(nodes - {"S", "T", omit}):
        incoming = pulp.lpSum(flow[u, v] for u, v in flow if v == node)
        outgoing = pulp.lpSum(flow[u, v] for u, v in flow if u == node)
        model += incoming == outgoing, f"balance_{node}"
    status = model.solve(solver)
    print(pulp.LpStatus[status])
    if status != pulp.LpStatusOptimal:
        raise RuntimeError(pulp.LpStatus[status])
    return {edge: x.value() for edge, x in flow.items()}

plan = route(capacity)
for edge, amount in plan.items():
    print(" → ".join(edge), amount, "/", capacity[edge])

Optimal
S → A 500.0 / 500
A → C 500.0 / 500
C → T 500.0 / 650
S → B 350.0 / 600
B → D 350.0 / 600
D → T 350.0 / 350
C → D 0.0 / 120

def balances(plan):
    net = {}
    for (u, v), amount in plan.items():
        net[u] = net.get(u, 0) + amount
        net[v] = net.get(v, 0) - amount
    return net

def check_plan(plan, capacities):
    net = balances(plan)
    assert all(-tol <= plan[e] <= limit + tol for e, limit in capacities.items())
    assert all(abs(value) <= tol for node, value in net.items() if node not in ("S", "T"))
    assert abs(net["S"] + net["T"]) <= tol
    return -net["T"]

print(balances(plan))
print("Tok [zásilky/h]:", check_plan(plan, capacity))

{'S': 850.0, 'A': 0.0, 'C': 0.0, 'T': -850.0, 'B': 0.0, 'D': 0.0}
Tok [zásilky/h]: 850.0

suspect = route(capacity, omit="C")
print(balances(suspect))

Optimal
{'S': 350.0, 'A': 0.0, 'C': 650.0, 'T': -1000.0, 'B': 0.0, 'D': 0.0}

for edge in [("A", "C"), ("D", "T")]:
    changed = capacity.copy()
    changed[edge] += 100
    result = check_plan(route(changed), changed)
    print(edge, result, "Δ =", result - check_plan(plan, capacity))

Optimal
('A', 'C') 850.0 Δ = 0.0
Optimal
('D', 'T') 950.0 Δ = 100.0

rng = np.random.default_rng(41)
width = 48
collection, hubs, delivery = [
    [f"{prefix}{j:02d}" for j in range(width)] for prefix in ["A", "B", "C"]
]
large = {("S", node): int(rng.integers(850, 1350)) for node in collection}
for left, right in [(collection, hubs), (hubs, delivery)]:
    for j, u in enumerate(left):
        for offset in range(8):
            large[u, right[(j + offset) % width]] = int(rng.integers(90, 215))
for node in delivery:
    large[node, "T"] = int(rng.integers(850, 1350))

large_plan = route(large)
baseline = check_plan(large_plan, large)
print("Spoje:", len(large), "| Tok [zásilky/h]:", baseline)
print("Největší odchylka bilance:", max(abs(n) for v, n in balances(large_plan).items()
                                        if v not in ("S", "T")))

Optimal
Spoje: 864 | Tok [zásilky/h]: 52651.0
Největší odchylka bilance: 0.0

hub_load = {h: sum(x for (u, v), x in large_plan.items() if v == h) for h in hubs}
closed = max(hubs, key=hub_load.get)
after = {edge: 0 if closed in edge else limit for edge, limit in large.items()}
revised = check_plan(route(after), after)
print("Centrum:", closed, "| Původní průtok centra:", hub_load[closed])
print("Tok po výpadku:", revised, "| Ztráta [zásilky/h]:", baseline - revised)

Optimal
Centrum: B40 | Původní průtok centra: 1319.0
Tok po výpadku: 51779.0 | Ztráta [zásilky/h]: 872.0

rng = np.random.default_rng(42)
length = rng.integers(10, 91, 150)
mass = rng.uniform(1, 24, 150).round(1)
# První zásilka je krabice. V další variantě k ní přidáme tubus.
length[0], mass[0] = 40, 4
manual = mass + 0.2 * length > 24

def parcel_plot():
    for mask, marker, label in [(~manual, "o", "Automatická linka"),
                                (manual, "^", "Ruční zpracování")]:
        plt.scatter(length[mask], mass[mask], marker=marker, label=label, alpha=0.8)
    plt.xlabel("Nejdelší strana [cm]")
    plt.ylabel("Hmotnost [kg]")
    plt.legend()

parcel_plot()
plt.show()

def separator(records):
    model = pulp.LpProblem("sorting")
    a, b, c = [pulp.LpVariable(name) for name in ("a", "b", "c")]
    model += 0
    for side, weight, by_hand in records:
        score = a * int(side) + b * float(weight) + c
        model += (score >= 1) if by_hand else (score <= -1)
    status = model.solve(solver)
    print(pulp.LpStatus[status])
    if status == pulp.LpStatusOptimal:
        return a.value(), b.value(), c.value()
    return None

records = list(zip(length, mass, manual))
a, b, c = separator(records)
scores = a * length + b * mass + c
assert min(scores[manual]) >= 1 - tol and max(scores[~manual]) <= -1 + tol
parcel_plot()
xs = np.array([length.min(), length.max()])
if abs(b) > tol:
    plt.plot(xs, -(a * xs + c) / b, color="black")
else:
    plt.axvline(-c/a, color="black")
plt.show()

Optimal

separator(records + [(40, 4, True)])

Infeasible

rng = np.random.default_rng(43)
load = np.linspace(5, 100, 300)
power = 85 + 1.65 * load + rng.normal(0, 3, len(load))
changed_power = power.copy()
changed_power[150] += 90

def fit(readings, criterion):
    model = pulp.LpProblem("server_power", pulp.LpMinimize)
    a, b = pulp.LpVariable("slope"), pulp.LpVariable("intercept")
    if criterion == "sum":
        bounds = [pulp.LpVariable(f"e_{i}", lowBound=0) for i in range(len(load))]
        model += pulp.lpSum(bounds)
    else:
        t = pulp.LpVariable("t", lowBound=0)
        bounds = [t] * len(load)
        model += t
    for u, p, bound in zip(load, readings, bounds):
        residual = a * float(u) + b - float(p)
        model += bound >= residual
        model += bound >= -residual
    status = model.solve(solver)
    if status != pulp.LpStatusOptimal:
        raise RuntimeError(pulp.LpStatus[status])
    predicted = a.value() * load + b.value()
    errors = abs(predicted - readings)
    checked = errors.sum() if criterion == "sum" else errors.max()
    assert np.isclose(checked, pulp.value(model.objective), rtol=1e-5)
    return predicted, errors.mean(), errors.max()

fig, axes = plt.subplots(1, 2, figsize=(10, 4), sharex=True, sharey=True)
for ax, readings, title in zip(axes, [power, changed_power],
                               ["Původní měření", "Jedno zvýšené měření"]):
    ax.scatter(load, readings, s=7, color="0.65")
    print(title, "| průměrná / největší absolutní odchylka [W]")
    for criterion, label, style in [("sum", "Součet", "-"), ("max", "Maximum", "--")]:
        predicted, mean_error, max_error = fit(readings, criterion)
        ax.plot(load, predicted, style, label=label)
        print(label, round(mean_error, 2), "/", round(max_error, 2))
    ax.set(title=title, xlabel="Vytížení [%]")
    ax.legend()
axes[0].set_ylabel("Příkon [W]")
fig.tight_layout()
plt.show()

Původní měření | průměrná / největší absolutní odchylka [W]
Součet 2.4 / 8.69
Maximum 2.53 / 8.13

Jedno zvýšené měření | průměrná / největší absolutní odchylka [W]
Součet 2.68 / 86.2
Maximum 39.42 / 46.94

Cvičení 4: Tok v síti, třídění zásilek a proložení přímky¶

0. Příprava prostředí — spusťte jako první¶

Co si dnes vyzkoušíme¶

1. Kolik zásilek projde sítí dep?¶

2. Od jednoho překladiště k celé síti¶

Matematický model maximálního toku¶

3. Kontrola výsledku a hledání chyby¶

Co změří kontrola bilance?¶

4. Který spoj posílit? Doplňte dvě varianty¶

5. Noční přeprava v síti se 146 uzly¶

Výpadek třídicího centra¶

6. Které zásilky prošly automatickou linkou?¶

LP pro oddělení dvou skupin¶

7. Jak souvisí příkon serveru s jeho vytížením?¶

Dvě účelové funkce pro stejná měření¶