Domácí úkol ke cvičení 3: Malý výrobní projekt v PuLP¶
Řešení.
Tento domácí úkol procvičuje samostatné sestavení, ověření a porovnání výrobních plánů. Otevřete notebook v Jupyteru nebo Google Colabu a spouštějte buňky shora dolů pomocí Shift+Enter. Používáme Python 3.12 nebo novější a PuLP 3.3.2. Nejprve spusťte instalační buňku ze C2 v sekci 0.
V zadání doplňujte buňky označené Doplňte a komentáře TODO. Samotné spuštění nevyplněné šablony model nevyřeší. Před opětovným sestavením modelu spusťte znovu buňku, která jej vytváří; před odevzdáním restartujte jádro a spusťte celý doplněný notebook shora. Ve verzi s řešením jsou všechny kroky vyplněné. Data jsou umělá a úplně uvedená níže.
0. Příprava prostředí — spusťte jako první¶
Tato buňka nainstaluje PuLP 3.3.2 do prostředí aktuálního notebooku, pokud správná verze chybí. Stačí ji spustit; příkazy není potřeba upravovat.
Pokud už máte v paměti PuLP 4, buňka po instalaci výpočet zastaví a vyzve
k restartu relace/jádra. V Colabu otevřete nabídku Runtime (Běhové prostředí)
a zvolte restart relace; v Jupyteru Kernel → Restart Kernel. Potom spusťte
notebook znovu shora. Samotné opakování import pulp starou verzi z paměti neodstraní.
Po hlášení „Připraveno“ můžete pokračovat. Instalace vyžaduje připojení k internetu.
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
if pulp.__version__ != "3.3.2":
raise RuntimeError("Spusťte sekci 0 a postupujte podle pokynu k restartu.")
print("PuLP:", pulp.__version__)
TOL = 1e-7
PuLP: 3.3.2
1. Zadání: co máme rozhodnout?¶
Malá výrobna připravuje tři kapalné směsi A, B a C. Rozhodujeme, kolik litrů každé směsi vyrobit v jedné dávce. K dispozici je dohromady 18 minut stroje a 16 kg společné suroviny.
| Údaj | A | B | C |
|---|---|---|---|
| Příspěvek k zisku [Kč/l] | 40 | 30 | 50 |
| Čas stroje [min/l] | 2 | 1 | 3 |
| Spotřeba suroviny [kg/l] | 1 | 2 | 2 |
| Nejvyšší odbyt [l] | 4 | 6 | 5 |
Směsi jsou dělitelné: lze vyrobit i část litru. Vše vyrobené se prodá, nejvýše však do uvedené meze odbytu. Nevyrábíme do zásoby. Spotřeby i příspěvky jsou přímo úměrné množství. Příspěvek k zisku již zohledňuje náklady na daný litr; další náklady neuvažujeme. Data jsou uměle vytvořená pro toto cvičení.
Doplňte: Definujte proměnné včetně jednotek a oborů. Na papír nebo do následující buňky zapište účelovou funkci a všechna omezení. Proč nejvyšší odbyt není povinné výrobní množství?
Označme $x_A,x_B,x_C$ vyrobené množství jednotlivých směsí v litrech. Úloha je $$\begin{aligned} \max\quad &40x_A+30x_B+50x_C\\ \text{za podmínek}\quad &2x_A+x_B+3x_C\le18 &&\text{(minuty stroje)},\\ &x_A+2x_B+2x_C\le16 &&\text{(kg suroviny)},\\ &0\le x_A\le4,\quad 0\le x_B\le6,\quad 0\le x_C\le5. \end{aligned}$$ Proměnné jsou spojité. Hodnota účelové funkce je v Kč za dávku. Meze odbytu jsou horní meze: méně vyrobit smíme, více ne.
2. Data oddělená od modelu¶
Seznam products obsahuje názvy výrobků. Slovník přiřazuje klíči hodnotu: například profit["A"] je příspěvek k zisku z litru směsi A. Ve vnořeném slovníku use nejprve vybereme zdroj, potom výrobek: use["machine"]["A"] je počet minut stroje na litr A.
Doplňte ústně: V jakých jednotkách jsou capacity["material"], use["material"]["B"] a demand["B"]? Která čísla jsou společná pro celou dávku?
products = ["A", "B", "C"]
resources = ["machine", "material"]
profit = {"A": 40, "B": 30, "C": 50} # Kč/l
use = {
"machine": {"A": 2, "B": 1, "C": 3}, # min/l
"material": {"A": 1, "B": 2, "C": 2}, # kg/l
}
capacity = {"machine": 18, "material": 16}
units = {"machine": "min", "material": "kg"}
demand = {"A": 4, "B": 6, "C": 5} # l
print(profit["A"])
print(use["machine"]["A"])
40 2
3. Proměnné a účelová funkce¶
Vytvoříme nový model a prázdný slovník x. Cyklus pro každý výrobek vytvoří jednu proměnnou. K x["A"] pak přistupujeme stejně jako k údaji profit["A"]; v x jsou ale objekty proměnných, jejichž hodnoty teprve hledáme.
model = pulp.LpProblem("production", pulp.LpMaximize)
x = {}
for p in products:
x[p] = pulp.LpVariable(f"x_{p}", lowBound=0, cat="Continuous")
print(x)
{'A': x_A, 'B': x_B, 'C': x_C}
Výraz profit[p] * x[p] je příspěvek výrobku p. Zápis
pulp.lpSum(profit[p] * x[p] for p in products) tyto členy sečte přes celý seznam výrobků.
Doplňte: Přidejte do modelu účelovou funkci. Výraz má odpovídat součtu $40x_A+30x_B+50x_C$. Použijte lpSum, aby zápis fungoval i při změně seznamu výrobků.
model += pulp.lpSum(profit[p] * x[p] for p in products), "profit"
4. Od jednoho omezení k celé skupině¶
Nejprve sestavíme pouze výraz pro spotřebu stroje. Následující buňka jej nepřidává jako omezení; jen ukazuje tři členy součtu.
machine_use = 2*x["A"] + x["B"] + 3*x["C"]
print(machine_use)
2*x_A + x_B + 3*x_C
Doplňte: Tentýž výraz sestavte pomocí use["machine"][p] a lpSum a vypište jej. Potom nahraďte klíč "machine" proměnnou r: cyklem přes resources přidejte pro každý zdroj jednu kapacitní podmínku. Pravá strana je capacity[r]. Omezení pojmenujte f"capacity_{r}".
Druhým cyklem přidejte pro každý výrobek podmínku x[p] <= demand[p], pojmenovanou f"sales_{p}". Nezápornost už zajišťuje lowBound=0. V modelu tak bude pět pojmenovaných omezení a tři nezáporné proměnné.
print(pulp.lpSum(use["machine"][p] * x[p] for p in products))
for r in resources:
model += (
pulp.lpSum(use[r][p] * x[p] for p in products) <= capacity[r],
f"capacity_{r}",
)
for p in products:
model += x[p] <= demand[p], f"sales_{p}"
2*x_A + x_B + 3*x_C
5. Řešení a první čtení výsledku¶
Doplňte: Vypište model a porovnejte jej s formulací z první sekce. Potom jej vyřešte pomocí pulp.PULP_CBC_CMD(msg=False) a vypište pulp.LpStatus[model.status].
Pro tuto malou úlohu pokračujte k optimálnímu plánu pouze při model.status == pulp.LpStatusOptimal. Jinak výpočet zastavte pomocí raise RuntimeError(...) a vraťte se k modelu. Po úspěšném řešení uložte číselné hodnoty do slovníku plan a hodnotu účelové funkce do base_profit. Symbolické proměnné x a nalezená čísla plan jsou dvě různé věci.
print(model)
model.solve(pulp.PULP_CBC_CMD(msg=False))
base_status = pulp.LpStatus[model.status]
print("Stav:", base_status)
if model.status != pulp.LpStatusOptimal:
raise RuntimeError("Nejprve zkontrolujte model a stav řešiče.")
plan = {p: x[p].value() for p in products}
base_profit = pulp.value(model.objective)
print("Plán [l]:", plan)
print("Příspěvek k zisku [Kč]:", base_profit)
production:
MAXIMIZE
40*x_A + 30*x_B + 50*x_C + 0.0
SUBJECT TO
capacity_machine: 2 x_A + x_B + 3 x_C <= 18
capacity_material: x_A + 2 x_B + 2 x_C <= 16
sales_A: x_A <= 4
sales_B: x_B <= 6
sales_C: x_C <= 5
VARIABLES
x_A Continuous
x_B Continuous
x_C Continuous
Stav: Optimal
Plán [l]: {'A': 4.0, 'B': 4.0, 'C': 2.0}
Příspěvek k zisku [Kč]: 380.0
6. Kontrola původních podmínek¶
Stav Optimal je informace od řešiče o zadaném modelu. Ještě ověříme, zda nalezená čísla splňují původní podmínky zadání. Kontrolu počítáme ze vstupních dat a hodnot plan, nikoli jen z interních údajů o omezeních v PuLP. Přitom musíme porovnat i slovníky s tabulkou zadání; chybný údaj převzatý do modelu i kontroly by samotný přepočet neodhalil.
Doplňte: Pro každý zdroj sestavte řádek tabulky: spotřeba, kapacita, rezerva, jednotka a splnění podmínky. Rezerva je kapacita minus spotřeba. Následně ověřte meze 0 <= plan[p] <= demand[p] a znovu spočítejte příspěvek k zisku.
Zde už sčítáme čísla, proto stačí vestavěné sum(...). Použijte absolutní toleranci TOL = 1e-7: pro kapacity v minutách, respektive kg, pro množství v litrech a pro shodu hodnot účelové funkce v Kč. Pro jiná měřítka dat je nutné toleranci zvolit znovu. Nezaokrouhlujte hodnoty před kontrolou; dvě desetinná místa použijte jen při výpisu. assert podmínka zastaví buňku, pokud podmínka není splněna.
print("Zdroj Spotřeba Kapacita Rezerva Jednotka Splněno")
for r in resources:
used = sum(use[r][p] * plan[p] for p in products)
reserve = capacity[r] - used
ok = used <= capacity[r] + TOL
print(f"{r:9} {used:9.2f} {capacity[r]:12.2f} {reserve:12.2f} {units[r]:8} {ok}")
assert ok, f"Překročena kapacita: {r}"
print("\nVýrobek Množství [l] Mez odbytu [l] Rezerva [l] Splněno")
for p in products:
ok = -TOL <= plan[p] <= demand[p] + TOL
print(f"{p:9} {plan[p]:12.2f} {demand[p]:16.2f} {demand[p]-plan[p]:13.2f} {ok}")
assert ok, f"Nesplněny meze proměnné: {p}"
checked_profit = sum(profit[p] * plan[p] for p in products)
assert abs(checked_profit - pulp.value(model.objective)) <= TOL
print("\nPřepočtený příspěvek k zisku [Kč]:", checked_profit)
print("Shoda s hodnotou účelové funkce z řešiče: ano")
Zdroj Spotřeba Kapacita Rezerva Jednotka Splněno machine 18.00 18.00 0.00 min True material 16.00 16.00 0.00 kg True Výrobek Množství [l] Mez odbytu [l] Rezerva [l] Splněno A 4.00 4.00 0.00 True B 4.00 6.00 2.00 True C 2.00 5.00 3.00 True Přepočtený příspěvek k zisku [Kč]: 380.0 Shoda s hodnotou účelové funkce z řešiče: ano
Interpretace: Které zdroje jsou plně využité? U kterých výrobků zbývá prostor do meze odbytu? Zaručuje samotné splnění kontrol optimalitu?
Plán je 4 l A, 4 l B a 2 l C s příspěvkem 380 Kč. Stroj spotřebuje 18 minut a surovina 16 kg, takže oba zdroje jsou plně využité. Rezervy odbytu jsou 0 l u A, 2 l u B a 3 l u C. Nezápornost i všechny horní meze jsou splněny.
Kontroly dokládají přípustnost nalezeného plánu a správný výpočet hodnoty účelové funkce v tomto řešení. Samy o sobě neprokazují optimalitu; tu zde přebíráme ze stavu řešiče.
7. Samostatná změna: více času stroje¶
Výrobna může pro stejnou dávku využít 22 minut stroje místo 18. Ostatní údaje zůstávají stejné. Zatím neuvažujeme cenu za dodatečný čas.
Nejprve odhadněte: Může se nejlepší dosažitelný příspěvek k zisku zhoršit? Musí se zvýšit výroba všech směsí? Zdůvodněte svou odpověď ještě před výpočtem.
Doplňte: Vytvořte scenario_capacity = capacity.copy() a změňte pouze položku "machine". Kopie zachová data původní varianty. Sestavte nový model scenario_model s novými proměnnými scenario_x stejným postupem jako výše, vyřešte jej a ověřte stav. Uložte scenario_plan a scenario_profit. Původní model, plan a capacity zachovejte pro porovnání.
Původní plán zůstává přípustný, takže nejlepší dosažitelný příspěvek nemůže klesnout. Není však předem zaručeno, že vzroste nebo že vzroste množství každého výrobku.
scenario_capacity = capacity.copy()
scenario_capacity["machine"] = 22
scenario_model = pulp.LpProblem("production_more_time", pulp.LpMaximize)
scenario_x = {}
for p in products:
scenario_x[p] = pulp.LpVariable(f"scenario_x_{p}", lowBound=0, cat="Continuous")
scenario_model += pulp.lpSum(profit[p] * scenario_x[p] for p in products), "profit"
for r in resources:
scenario_model += (
pulp.lpSum(use[r][p] * scenario_x[p] for p in products) <= scenario_capacity[r],
f"capacity_{r}",
)
for p in products:
scenario_model += scenario_x[p] <= demand[p], f"sales_{p}"
scenario_model.solve(pulp.PULP_CBC_CMD(msg=False))
scenario_status = pulp.LpStatus[scenario_model.status]
print("Stav:", scenario_status)
if scenario_model.status != pulp.LpStatusOptimal:
raise RuntimeError("Nejprve zkontrolujte změněný model a stav řešiče.")
scenario_plan = {p: scenario_x[p].value() for p in products}
scenario_profit = pulp.value(scenario_model.objective)
Stav: Optimal
Doplňte kontrolu nové varianty: Zopakujte tabulku ze sekce 6 s scenario_plan, scenario_capacity a scenario_model. Použijete-li omylem původní kapacitu, budete kontrolovat jinou úlohu. Potom porovnejte obě varianty v jedné tabulce.
print("Zdroj Spotřeba Kapacita Rezerva Jednotka Splněno")
for r in resources:
used = sum(use[r][p] * scenario_plan[p] for p in products)
reserve = scenario_capacity[r] - used
ok = used <= scenario_capacity[r] + TOL
print(f"{r:9} {used:9.2f} {scenario_capacity[r]:12.2f} {reserve:12.2f} {units[r]:8} {ok}")
assert ok, f"Překročena kapacita: {r}"
print("\nVýrobek Množství [l] Mez odbytu [l] Rezerva [l] Splněno")
for p in products:
ok = -TOL <= scenario_plan[p] <= demand[p] + TOL
print(f"{p:9} {scenario_plan[p]:12.2f} {demand[p]:16.2f} {demand[p]-scenario_plan[p]:13.2f} {ok}")
assert ok, f"Nesplněny meze proměnné: {p}"
checked_profit = sum(profit[p] * scenario_plan[p] for p in products)
assert abs(checked_profit - pulp.value(scenario_model.objective)) <= TOL
print("\nPřepočtený příspěvek k zisku [Kč]:", checked_profit)
print("Shoda s hodnotou účelové funkce z řešiče: ano")
Zdroj Spotřeba Kapacita Rezerva Jednotka Splněno machine 22.00 22.00 0.00 min True material 16.00 16.00 0.00 kg True Výrobek Množství [l] Mez odbytu [l] Rezerva [l] Splněno A 4.00 4.00 0.00 True B 2.00 6.00 4.00 True C 4.00 5.00 1.00 True Přepočtený příspěvek k zisku [Kč]: 420.0 Shoda s hodnotou účelové funkce z řešiče: ano
print("Výrobek Původně [l] Nově [l] Rozdíl [l]")
for p in products:
print(f"{p:9} {plan[p]:11.2f} {scenario_plan[p]:10.2f} {scenario_plan[p]-plan[p]:12.2f}")
print("Příspěvek k zisku [Kč]:", base_profit, "→", scenario_profit)
print("Změna příspěvku [Kč]:", scenario_profit - base_profit)
Výrobek Původně [l] Nově [l] Rozdíl [l] A 4.00 4.00 0.00 B 4.00 2.00 -2.00 C 2.00 4.00 2.00 Příspěvek k zisku [Kč]: 380.0 → 420.0 Změna příspěvku [Kč]: 40.0
Komentář: Jak se změnila skladba výroby? Proč neporovnáváme jen hodnotu účelové funkce? Co ještě potřebujeme vědět, abychom rozhodli o pořízení dodatečného času stroje?
Nový plán je 4 l A, 2 l B a 4 l C, příspěvek 420 Kč. Výroba B klesla o 2 l a výroba C vzrostla o 2 l; A zůstává na mezi odbytu. Nahrazení 2 l B dvěma litry C spotřebuje o 4 minuty stroje více, ale stejné množství suroviny. Příspěvek vzroste o 40 Kč. Oba zdroje jsou opět plně využité, rezervy odbytu jsou 0, 4 a 1 l.
Dodatečné čtyři minuty přinesou v tomto modelu 40 Kč před započtením ceny tohoto času. Pořízení celého balíčku času je při ostatních předpokladech výhodné, stojí-li méně než 40 Kč; při ceně 40 Kč je čistý přírůstek nulový. Tento výpočet neříká, že každé další čtyři minuty přinesou dalších 40 Kč.
8. Krátký závěr projektu¶
Doplňte vlastními slovy krátký závěr, který obsahuje doporučený plán, jeho kontrolu a srovnání s novou variantou. Připojte předpoklady, které omezují použití výsledku, zdroj dat a návod k opakování výpočtu. Neopisujte jen výstup Optimal.
Při přípravě vlastního zápočtového projektu lze použít stejnou stavbu: zadání a model → výpočet → kontrola a interpretace → jedna změna zadání. Samostatný projekt musí obsahovat vaše vlastní zpracování a uvedení použitých zdrojů a pomoci.
Příklad stručného zápisu¶
Pro původní dávku doporučujeme vyrobit 4 l A, 4 l B a 2 l C; příspěvek k zisku je 380 Kč. CBC oznámil stav Optimal. Dosazením do původních podmínek jsme ověřili obě kapacity, nezápornost, všechny meze odbytu i hodnotu účelové funkce s absolutní tolerancí $10^{-7}$ v příslušných jednotkách.
Při zvýšení času stroje z 18 na 22 minut doporučujeme 4 l A, 2 l B a 4 l C. I tento plán prošel kontrolami a CBC oznámil Optimal; příspěvek vzroste na 420 Kč. Rozdíl 40 Kč ještě nezahrnuje případnou cenu dodatečného času. Výsledek závisí na dělitelnosti směsí, stálých spotřebách a příspěvcích na litr a možnosti prodat vše do uvedených mezí. Například prodej výhradně v celých lahvích by vyžadoval nové zdůvodnění oborů proměnných.
Data: umělá tabulka z domácího úkolu ke C3 kurzu 01LIP; všechny vstupy jsou v sekci 2, žádný externí soubor není potřeba. Postup a kód: tento řešený domácí úkol ke C3. Při převzetí do vlastní práce uveďte tuto předlohu i další využité zdroje a pomoc.
Opakování: Python 3.12 nebo novější, PuLP 3.3.2 a CBC dodaný s touto instalací. Spusťte instalační buňku v sekci 0, případně podle jejího hlášení restartujte jádro/relaci a spusťte všechny buňky shora. Používáme PULP_CBC_CMD(msg=False) s výchozími nastaveními; msg=False pouze skryje protokol. Konkrétní verzi CBC použitého při výpočtu zobrazí následující buňka.
# Verze prostředí a řešiče pro záznam výpočtu.
print("Python:", sys.version.split()[0])
print("PuLP:", pulp.__version__)
cbc_info = subprocess.run(
[pulp.PULP_CBC_CMD().path, "-stop"],
capture_output=True, text=True, check=True,
)
for line in cbc_info.stdout.splitlines():
if "Version:" in line:
print("CBC", line.strip())
Python: 3.12.14 PuLP: 3.3.2 CBC Version: 2.10.3