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 nejprve sekci 0 a restartujte jádro podle pokynů.")
print(pulp.__version__)

3.3.2

model = pulp.LpProblem("mixtures", pulp.LpMaximize)

x = pulp.LpVariable("x", lowBound=0, cat="Continuous")
y = pulp.LpVariable("y", lowBound=0, cat="Continuous")

model += 40*x + 30*y, "profit"

model += 2*x + y <= 6, "machine"

model += x + y <= 4.5, "material"

print(model)
print("Dolní mez x:", x.lowBound)
print("Dolní mez y:", y.lowBound)

mixtures:
MAXIMIZE
40*x + 30*y + 0
SUBJECT TO
machine: 2 x + y <= 6

material: x + y <= 4.5

VARIABLES
x Continuous
y Continuous

Dolní mez x: 0
Dolní mez y: 0

solver = pulp.PULP_CBC_CMD(msg=False)
status = model.solve(solver)
print("Kód stavu:", status)
print("Stav:", pulp.LpStatus[status])

Kód stavu: 1
Stav: Optimal

if status == pulp.LpStatusOptimal:
    print("Směs A [l]:", x.value())
    print("Směs B [l]:", y.value())
    print("Zisk [Kč]:", pulp.value(model.objective))
else:
    print("Nebylo potvrzeno optimum. Stav:", pulp.LpStatus[status])

Směs A [l]: 1.5
Směs B [l]: 3.0
Zisk [Kč]: 150.0

if status == pulp.LpStatusOptimal:
    print("Čas [min]:", 2*x.value() + y.value(), "z 6")
    print("Surovina [kg]:", x.value() + y.value(), "ze 4.5")
    print("Rezerva času [min]:", 6 - (2*x.value() + y.value()))
    print("Rezerva suroviny [kg]:", 4.5 - (x.value() + y.value()))
    print("Zisk přepočtený ze zadání [Kč]:", 40*x.value() + 30*y.value())

Čas [min]: 6.0 z 6
Surovina [kg]: 4.5 ze 4.5
Rezerva času [min]: 0.0
Rezerva suroviny [kg]: 0.0
Zisk přepočtený ze zadání [Kč]: 150.0

integer_model = pulp.LpProblem("whole_containers", pulp.LpMaximize)
x_i = pulp.LpVariable("x_i", lowBound=0, cat="Integer")
y_i = pulp.LpVariable("y_i", lowBound=0, cat="Integer")
integer_model += 40*x_i + 30*y_i, "profit"
integer_model += 2*x_i + y_i <= 6, "machine"
integer_model += x_i + y_i <= 4.5, "material"
status_i = integer_model.solve(pulp.PULP_CBC_CMD(msg=False))
print(pulp.LpStatus[status_i])
if status_i == pulp.LpStatusOptimal:
    print("x =", x_i.value(), "y =", y_i.value())
    print("Zisk [Kč]:", pulp.value(integer_model.objective))

Optimal
x = 2.0 y = 2.0
Zisk [Kč]: 140.0

mixed_model = pulp.LpProblem("mixed_containers", pulp.LpMaximize)
x_m = pulp.LpVariable("x_m", lowBound=0, cat="Integer")
y_m = pulp.LpVariable("y_m", lowBound=0, cat="Continuous")
mixed_model += 40*x_m + 30*y_m, "profit"
mixed_model += 2*x_m + y_m <= 6, "machine"
mixed_model += x_m + y_m <= 4.5, "material"
status_m = mixed_model.solve(pulp.PULP_CBC_CMD(msg=False))
print(pulp.LpStatus[status_m])
if status_m == pulp.LpStatusOptimal:
    print("x =", x_m.value(), "y =", y_m.value())
    print("Zisk [Kč]:", pulp.value(mixed_model.objective))

Optimal
x = 1.0 y = 3.5
Zisk [Kč]: 145.0

volume_model = pulp.LpProblem("total_volume", pulp.LpMaximize)
x_v = pulp.LpVariable("x_v", lowBound=0, cat="Continuous")
y_v = pulp.LpVariable("y_v", lowBound=0, cat="Continuous")
volume_model += x_v + y_v, "volume"
volume_model += 2*x_v + y_v <= 6, "machine"
volume_model += x_v + y_v <= 4.5, "material"
status_v = volume_model.solve(pulp.PULP_CBC_CMD(msg=False))
print(pulp.LpStatus[status_v])
if status_v == pulp.LpStatusOptimal:
    print("x =", x_v.value(), "y =", y_v.value())
    print("Objem [l]:", pulp.value(volume_model.objective))

Optimal
x = 0.0 y = 4.5
Objem [l]: 4.5

infeasible_model = pulp.LpProblem("impossible_order", pulp.LpMaximize)
x_f = pulp.LpVariable("x_f", lowBound=0, cat="Continuous")
y_f = pulp.LpVariable("y_f", lowBound=0, cat="Continuous")
infeasible_model += 40*x_f + 30*y_f, "profit"
infeasible_model += 2*x_f + y_f <= 6, "machine"
infeasible_model += x_f + y_f <= 4.5, "material"
infeasible_model += x_f + y_f >= 5, "order"
status_f = infeasible_model.solve(pulp.PULP_CBC_CMD(msg=False))
print(pulp.LpStatus[status_f])

Infeasible

unbounded_model = pulp.LpProblem("no_upper_bound", pulp.LpMaximize)
u = pulp.LpVariable("u", lowBound=0, cat="Continuous")
unbounded_model += u, "objective"
status_u = unbounded_model.solve(pulp.PULP_CBC_CMD(msg=False))
print(pulp.LpStatus[status_u])

Unbounded

practice_model = pulp.LpProblem("four_kg", pulp.LpMaximize)
x_p = pulp.LpVariable("x_p", lowBound=0, cat="Continuous")
y_p = pulp.LpVariable("y_p", lowBound=0, cat="Continuous")
practice_model += 40*x_p + 30*y_p, "profit"
practice_model += 2*x_p + y_p <= 6, "machine"
practice_model += x_p + y_p <= 4, "material"
status_p = practice_model.solve(pulp.PULP_CBC_CMD(msg=False))
print(pulp.LpStatus[status_p])
if status_p == pulp.LpStatusOptimal:
    print("x =", x_p.value(), "y =", y_p.value())
    print("Zisk [Kč]:", pulp.value(practice_model.objective))
    print("Čas [min]:", 2*x_p.value() + y_p.value(), "z 6")
    print("Surovina [kg]:", x_p.value() + y_p.value(), "ze 4")

Optimal
x = 2.0 y = 2.0
Zisk [Kč]: 140.0
Čas [min]: 6.0 z 6
Surovina [kg]: 4.0 ze 4

import random
import hashlib
import json

SEED = 20261016
N_PRODUCTS = 6000
N_RESOURCES = 600
RESOURCES_PER_PRODUCT = 12
rng = random.Random(SEED)

profit = [rng.randint(10, 100) for j in range(N_PRODUCTS)]
rows = [[] for i in range(N_RESOURCES)]
for j in range(N_PRODUCTS):
    for i in rng.sample(range(N_RESOURCES), RESOURCES_PER_PRODUCT):
        rows[i].append((j, rng.randint(1, 10)))
capacity = [0.5 * sum(a for j, a in row) for row in rows]

fingerprint = hashlib.sha256(
    json.dumps([profit, rows, capacity]).encode()
).hexdigest()[:16]
print("n =", N_PRODUCTS, "m =", N_RESOURCES)
print("nnz =", sum(len(row) for row in rows))
print("seed =", SEED, "SHA-256 (16) =", fingerprint)

n = 6000 m = 600
nnz = 72000
seed = 20261016 SHA-256 (16) = 9aafeb8801e6d560

%pip install pulp==3.3.2 pyomo==6.10.1 scipy==1.18.1 highspy==1.15.1 matplotlib

import os
import platform
import subprocess
import warnings
import numpy as np
import scipy
from scipy.optimize import linprog, OptimizeWarning
from scipy.sparse import csr_matrix
import pyomo.environ as pyo
import highspy
from importlib.metadata import version

HIGHS_OPTIONS = dict(threads=1, solver="simplex", parallel="off",
                     simplex_strategy=1, presolve="on", random_seed=0)
print("Python:", platform.python_version())
for package in ("pulp", "pyomo", "scipy", "highspy", "numpy"):
    print(package + ":", version(package))
print("HiGHS přes highspy:", highspy.Highs().version())
# SciPy obsahuje vlastní HiGHS; neveřejný import slouží jen k výpisu verze.
try:
    from scipy.optimize._highspy import _core
    scipy_highs = ".".join(str(getattr(_core, "HIGHS_VERSION_" + part))
                           for part in ("MAJOR", "MINOR", "PATCH"))
except (ImportError, AttributeError):
    scipy_highs = "verze nezjištěna"
print("HiGHS uvnitř SciPy:", scipy_highs)
print("OS:", platform.system(), platform.release(), platform.machine())
cpu = platform.processor() or platform.machine()
if platform.system() == "Darwin":
    query = subprocess.run(["sysctl", "-n", "machdep.cpu.brand_string"],
                           capture_output=True, text=True)
    if query.returncode == 0:
        cpu = query.stdout.strip()
print("CPU:", cpu, "logických CPU:", os.cpu_count())

Python: 3.12.14
pulp: 3.3.2
pyomo: 6.10.1
scipy: 1.18.1
highspy: 1.15.1
numpy: 2.5.3
HiGHS přes highspy: 1.15.1
HiGHS uvnitř SciPy: 1.12.0
OS: Darwin 25.5.0 arm64
CPU: Apple M5 Pro logických CPU: 18

def prepare_pulp():
    model = pulp.LpProblem("large_production", pulp.LpMinimize)
    x = [pulp.LpVariable(f"q_{j:05}", lowBound=0, upBound=1)
         for j in range(N_PRODUCTS)]
    model += pulp.lpSum(-profit[j] * x[j] for j in range(N_PRODUCTS))
    for i, row in enumerate(rows):
        model += pulp.lpSum(a * x[j] for j, a in row) <= capacity[i]
    solver = pulp.HiGHS(msg=False, **HIGHS_OPTIONS)

    def solve():
        status = model.solve(solver)
        ok = status == pulp.LpStatusOptimal and model.sol_status == pulp.LpSolutionOptimal
        return ok, [q.value() for q in x], -pulp.value(model.objective)
    return solve

def prepare_pyomo():
    model = pyo.ConcreteModel()
    model.x = pyo.Var(range(N_PRODUCTS), bounds=(0, 1))
    model.objective = pyo.Objective(
        expr=pyo.quicksum(-profit[j] * model.x[j] for j in range(N_PRODUCTS)),
        sense=pyo.minimize)
    model.resources = pyo.Constraint(range(N_RESOURCES), rule=lambda m, i:
        pyo.quicksum(a * m.x[j] for j, a in rows[i]) <= capacity[i])
    solver = pyo.SolverFactory("appsi_highs")
    solver.options.update(HIGHS_OPTIONS)

    def solve():
        result = solver.solve(model)
        ok = result.solver.termination_condition == pyo.TerminationCondition.optimal
        return ok, [pyo.value(model.x[j]) for j in range(N_PRODUCTS)], -pyo.value(model.objective)
    return solve

def make_arrays():
    coefficients, columns, starts = [], [], [0]
    for row in rows:
        for j, a in row:
            columns.append(j)
            coefficients.append(a)
        starts.append(len(coefficients))
    A = csr_matrix((np.array(coefficients, dtype=float),
                    np.array(columns, dtype=np.int32),
                    np.array(starts, dtype=np.int32)),
                   shape=(N_RESOURCES, N_PRODUCTS))
    c = -np.array(profit, dtype=float)
    b = np.array(capacity, dtype=float)
    return c, A, b

def prepare_scipy():
    c, A, b = make_arrays()

    def solve():
        # SciPy předá threads do HiGHS, ale upozorní na dodatečnou volbu.
        with warnings.catch_warnings():
            warnings.filterwarnings("ignore", category=OptimizeWarning,
                                    message=r"Unrecognized options detected: \{'threads': 1\}\..*")
            result = linprog(c, A_ub=A, b_ub=b, bounds=(0, 1),
                             method="highs-ds", options={"threads": 1})
        return result.success and result.status == 0, result.x, -result.fun if result.success else None
    return solve

def check_highs(status):
    if status != highspy.HighsStatus.kOk:
        raise RuntimeError(f"HiGHS: {status}")

def prepare_highspy():
    c, A, b = make_arrays()
    lp = highspy.HighsLp()
    lp.num_col_, lp.num_row_ = N_PRODUCTS, N_RESOURCES
    lp.col_cost_ = c
    lp.col_lower_, lp.col_upper_ = np.zeros(N_PRODUCTS), np.ones(N_PRODUCTS)
    lp.row_lower_, lp.row_upper_ = np.full(N_RESOURCES, -highspy.kHighsInf), b
    lp.a_matrix_.format_ = highspy.MatrixFormat.kRowwise
    lp.a_matrix_.start_, lp.a_matrix_.index_, lp.a_matrix_.value_ = A.indptr, A.indices, A.data
    solver = highspy.Highs()
    check_highs(solver.setOptionValue("output_flag", False))
    for option, value in HIGHS_OPTIONS.items():
        check_highs(solver.setOptionValue(option, value))

    def solve():
        check_highs(solver.passModel(lp))
        check_highs(solver.run())
        ok = solver.getModelStatus() == highspy.HighsModelStatus.kOptimal
        return ok, solver.getSolution().col_value, -solver.getObjectiveValue()
    return solve

import gc
import math
from time import perf_counter
from statistics import median

PREPARE = {"PuLP": prepare_pulp, "Pyomo": prepare_pyomo,
           "SciPy linprog": prepare_scipy, "highspy": prepare_highspy}
REPEATS = 3

def run_once(name):
    gc.collect()
    start = perf_counter()
    solve = PREPARE[name]()
    built = perf_counter()
    optimal, values, reported_profit = solve()
    finished = perf_counter()

    if not optimal:
        raise RuntimeError(f"{name}: optimalita nebyla potvrzena")
    if not all(v is not None and math.isfinite(v) and -1e-7 <= v <= 1+1e-7 for v in values):
        raise RuntimeError(f"{name}: porušeny meze proměnných")
    violation = max(0.0, max(
        (sum(a * values[j] for j, a in row) - capacity[i]) / max(1.0, capacity[i])
        for i, row in enumerate(rows)))
    if violation > 1e-7:
        raise RuntimeError(f"{name}: porušena kapacitní omezení")
    objective = sum(p * v for p, v in zip(profit, values))
    if not math.isclose(objective, reported_profit, rel_tol=1e-7, abs_tol=1e-6):
        raise RuntimeError(f"{name}: nesouhlasí přepočet zisku")
    if not 0.5 * sum(profit) - 1e-6 <= objective <= sum(profit) + 1e-6:
        raise RuntimeError(f"{name}: zisk je mimo známé meze")
    return dict(library=name, build_s=built-start, solve_s=finished-built,
                total_s=finished-start, objective=objective, status="Optimal")

# Jeden zahřívací běh pro každou knihovnu; nepatří do souhrnu.
for name in PREPARE:
    run_once(name)

runs = []
names = list(PREPARE)
for repeat in range(REPEATS):
    order = names[repeat:] + names[:repeat]
    for name in order:
        result = run_once(name)
        runs.append(result)
        print(f"{repeat+1}: {name:13s} {result['total_s']:.3f} s  {result['status']}")
reference = runs[0]["objective"]
for run in runs:
    if not math.isclose(run["objective"], reference, rel_tol=1e-7, abs_tol=1e-6):
        raise RuntimeError("Knihovny vrátily rozdílné optimální hodnoty")

summary = []
for name in names:
    selected = [r for r in runs if r["library"] == name]
    totals = [r["total_s"] for r in selected]
    summary.append(dict(library=name, count=len(selected), status="Optimal",
                        total_s=median(totals), min_s=min(totals), max_s=max(totals),
                        build_s=median(r["build_s"] for r in selected),
                        solve_s=median(r["solve_s"] for r in selected),
                        objective=selected[-1]["objective"]))
summary.sort(key=lambda r: r["total_s"])

1: PuLP          2.351 s  Optimal

1: Pyomo         2.623 s  Optimal

1: SciPy linprog 2.173 s  Optimal

1: highspy       2.176 s  Optimal

2: Pyomo         2.620 s  Optimal

2: SciPy linprog 2.172 s  Optimal

2: highspy       2.169 s  Optimal

2: PuLP          2.337 s  Optimal

3: SciPy linprog 2.169 s  Optimal

3: highspy       2.176 s  Optimal

from IPython.display import Markdown, display

header = "| Knihovna | Sestavení [s] | Řešení + načtení [s] | Celkem [s] | Min–max celkem [s] | Zisk [Kč] |"
lines = [header, "|---|---:|---:|---:|---:|---:|"]
for row in summary:
    lines.append(f"| {row['library']} | {row['build_s']:.3f} | {row['solve_s']:.3f} | "
                 f"{row['total_s']:.3f} | {row['min_s']:.3f}–{row['max_s']:.3f} | "
                 f"{row['objective']:.3f} |")
display(Markdown("\n".join(lines)))

import matplotlib.pyplot as plt

medians = [r["total_s"] for r in summary]
errors = [[r["total_s"]-r["min_s"] for r in summary],
          [r["max_s"]-r["total_s"] for r in summary]]
fig, ax = plt.subplots(figsize=(7, 3.5), layout="constrained")
ax.barh([r["library"] for r in summary], medians, xerr=errors,
        color="#315e85", capsize=4)
ax.invert_yaxis()
ax.set_xlim(left=0)
ax.set_xlabel("Celkový čas [s] — medián ze 3 běhů")
ax.set_title(f"{N_PRODUCTS} proměnných, {N_RESOURCES} kapacitních omezení")
ax.set_axisbelow(True)
ax.grid(axis="x", alpha=0.2)
ax.spines[["top", "right"]].set_visible(False)
plt.show()

Knihovna	Jak se v ní pracuje
PuLP	Vytvoříme pojmenované proměnné a zapisujeme výrazy a omezení podobně jako na papíře. Umí zadat LP i celočíselné a smíšené modely.
SciPy: `scipy.optimize.linprog`	Lineární program zadáváme pomocí vektoru koeficientů účelové funkce, matic omezení a mezí.
Pyomo	Modelovací knihovna pro řadu optimalizačních úloh, včetně lineárních, celočíselných i nelineárních modelů.
OR-Tools	Sada optimalizačních nástrojů; rozhraní `MPSolver` umožňuje sestavovat LP a smíšené celočíselné modely.

Model	Obor $x$	Obor $y$	Nalezený plán $(x,y)$ [l]	Zisk [Kč]
LP	spojitá	spojitá	(1,5; 3)	150
ILP	celočíselná	celočíselná	(2; 2)	140
MILP	celočíselná	spojitá	(1; 3,5)	145

Plán $(x,y)$ [l]	Čas [min]	Surovina [kg]	Celkový objem [l]
(0; 4,5)	4,5	4,5	4,5
(1; 3,5)	5,5	4,5	4,5

Stav	Kód	Co z něj čteme
`Optimal`	1	V našich dokončených výpočtech bylo nalezeno optimum; nemusí být jediné.
`Infeasible`	-1	Omezení nelze splnit současně.
`Unbounded`	-2	Účelovou funkci lze zlepšovat bez konečné meze.
`Not Solved`	0	Model zatím není vyřešen; tento stav nedokazuje nepřípustnost.
`Undefined`	-3	Není k dispozici jednoznačný závěr; je třeba prohlédnout zprávu řešiče.

Knihovna	Sestavení [s]	Řešení + načtení [s]	Celkem [s]	Min–max celkem [s]	Zisk [Kč]
SciPy linprog	0.006	2.166	2.172	2.169–2.173	230781.446
highspy	0.011	2.164	2.176	2.169–2.176	230781.446
PuLP	0.157	2.195	2.351	2.337–2.352	230781.446
Pyomo	0.067	2.553	2.622	2.620–2.623	230781.446

Cvičení 2: První model v PuLP krok za krokem¶

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

Co si dnes vyzkoušíme¶

1. Které knihovny připadají v úvahu?¶

Načtení knihovny¶

2. Zadání: dvě kapalné směsi¶

Model na papíře¶

3. Sestavení modelu po krocích¶

Krok 1: prázdná maximalizační úloha¶

Krok 2: proměnné a jejich obory¶

Krok 3: účelová funkce¶

Krok 4: omezení na čas stroje¶

Krok 5: omezení na surovinu¶

Krok 6: kontrola zapsaného modelu¶

4. Spuštění řešiče a čtení výsledku¶

Nejprve stav¶

Potom hodnoty proměnných a účelové funkce¶

Kontrola dosazením do zadání¶

5. Stejná data, jiné obory proměnných¶

Obě proměnné celočíselné: ILP¶

Jen jedna proměnná celočíselná: MILP¶

6. Více optimálních řešení¶

7. Model bez přípustného řešení¶

8. Krátké srovnání: neomezená úloha¶

Přehled stavů v PuLP 3.3.2¶

9. Krátký samostatný úkol¶

10. Pro zajímavost: PuLP a další knihovny na větší úloze¶

Zadání a data¶

Příprava knihoven¶

Stejný model ve čtyřech rozhraních¶

PuLP: pojmenované proměnné a výrazy¶

Pyomo: proměnné, účelová funkce a omezení jako části modelu¶

SciPy a highspy: maticový zápis¶

Co přesně měříme?¶

Výsledky měření¶

Jak srovnání číst¶

Údaj na jeden litr	Směs A	Směs B
Příspěvek k zisku [Kč/l]	40	30
Čas stroje [min/l]	2	1
Spotřeba společné suroviny [kg/l]	1	1