Classes with completion bonus and "if and only if" reduction
Class: BIP · Links: if and only if (two), CNF · Script: python/fam07_6_classesbonus.py
Difficulty: ★★★★☆ · Time: 45–60 min
Problem 7.6
A company has \(n\) jobs executable on a machine with availability \(a\). For each job \(j\), \(t_j\) is the time and \(r_j\) the revenue. The jobs are partitioned into \(q \ge 2\) classes. The availability is reduced by \(u > 0\) minutes if and only if the machine executes jobs of at least two different classes. For each class \(c\) an extra revenue \(v_c > 0\) is obtained if and only if all the jobs of the class are executed. Maximise the total revenue.
The problem in words. Two "if and only if"s: one rewards (completion), one penalises (mixing). For each it suffices to impose one direction: the objective imposes the other.
Model
Variables. \(n + q + 1\) binary: \(x_j\) (job executed), \(y_c\) (class complete), \(z\) (jobs of at least two classes).
\[
\begin{aligned}
\max ~~ \sum_{j=1}^{n} r_j\, x_j + \sum_{c=1}^{q} v_c\, y_c & & \\
\text{subject to} \quad x_j - y_c &\ge 0, & \forall c \in \{1, 2, \dots, q\},\ \forall j \in \mathscr{J}_c, \\
x_j + x_i - z &\le 1, & \forall c, g \in \{1, 2, \dots, q\},\ c < g,\ \forall j \in \mathscr{J}_c,\ \forall i \in \mathscr{J}_g, \\
\sum_{j=1}^{n} t_j\, x_j + u\, z &\le a, & \\
x_j &\in \{0, 1\}, & \forall j \in \{1, 2, \dots, n\}, \\
y_c &\in \{0, 1\}, & \forall c \in \{1, 2, \dots, q\}, \\
z &\in \{0, 1\}. &
\end{aligned}
\]
- the objective maximises revenues of the jobs plus bonuses of the complete classes;
- the all constraints: if a class is declared complete, all its jobs are executed (\(n\) constraints);
- the mixed constraints: if two jobs of different classes are executed, \(z = 1\) (\(\sum_{c<g} |\mathscr{J}_c|\,|\mathscr{J}_g|\) constraints);
- the availability constraint with the reduction \(u z\) (\(1\) constraint);
- the domain constraints.
Link between the variables: four implications
\(y_c\), from the constraint. \(y_c \Rightarrow (x_j \,\mathtt{AND}\, x_i \,\mathtt{AND}\, \dots)\): the expression \((x_j \,\mathtt{AND}\, \dots) \,\mathtt{OR}\, \mathtt{NOT}\,y_c\) becomes by distributivity the CNF \((x_j \,\mathtt{OR}\, \mathtt{NOT}\,y_c) \,\mathtt{AND}\, \dots\), i.e.\ \(x_j \ge y_c\). \(y_c\), from the optimum. If all jobs are executed then \(y_c = 1\): setting \(y_c = 1\) stays feasible and increases the objective by \(v_c > 0\).
\(z\), from the constraint. \(x_j \,\mathtt{AND}\, x_i \Rightarrow z\) for every mixed pair: De Morgan gives \(\mathtt{NOT}\,x_j \,\mathtt{OR}\, \mathtt{NOT}\,x_i \,\mathtt{OR}\, z\), i.e.\ \(x_j + x_i - z \le 1\). \(z\), from the optimum. If the executed jobs lie in a single class, \(z = 0\): setting \(z = 0\) stays feasible, frees \(u\) minutes and does not change the objective (where \(z\) does not appear) — "there exists an optimum", not "in every optimum".
The model in gurobipy
pairs = [(j, i) for c in range(q) for g in range(c + 1, q)
for j in J[c] for i in J[g]]
m = gp.Model("classes_bonus"); m.Params.OutputFlag = 0
x = m.addVars(n, vtype=GRB.BINARY, name="x")
y = m.addVars(q, vtype=GRB.BINARY, name="y")
z = m.addVar(vtype=GRB.BINARY, name="z")
m.setObjective(gp.quicksum(r[j] * x[j] for j in range(n))
+ gp.quicksum(v[c] * y[c] for c in range(q)), GRB.MAXIMIZE)
m.addConstrs((x[j] - y[c] >= 0 for c in range(q) for j in J[c]), name="all")
m.addConstrs((x[j] + x[i] - z <= 1 for (j, i) in pairs), name="mixed")
m.addConstr(gp.quicksum(t[j] * x[j] for j in range(n)) + u * z <= a, name="availability")
m.optimize()
The instance
\(n = 6\), \(q = 3\): \(\mathscr{J}_1 = \{1, 2\}\), \(\mathscr{J}_2 = \{3, 4\}\), \(\mathscr{J}_3 = \{5, 6\}\), \(a = 50\), \(u = 10\).
| \(j=1\) | \(j=2\) | \(j=3\) | \(j=4\) | \(j=5\) | \(j=6\) | |
|---|---|---|---|---|---|---|
| \(r_j\) | 10 | 5 | 20 | 12 | 10 | 22 |
| \(t_j\) | 5 | 15 | 25 | 15 | 10 | 38 |
| \(c=1\) | \(c=2\) | \(c=3\) | |
|---|---|---|---|
| \(v_c\) | 5 | 4 | 10 |
The model written on the data of the instance:
\[
\begin{array}{rrrrrrrrrrr c l}
\max & 10x_1 & +5x_2 & +20x_3 & +12x_4 & +10x_5 & +22x_6 & +5y_1 & +4y_2 & +10y_3 & & & \\
\text{subject to} & x_1 & & & & & & -y_1 & & & & \ge & 0\\
& & x_2 & & & & & -y_1 & & & & \ge & 0\\
& & & x_3 & & & & & -y_2 & & & \ge & 0\\
& & & & x_4 & & & & -y_2 & & & \ge & 0\\
& & & & & x_5 & & & & -y_3 & & \ge & 0\\
& & & & & & x_6 & & & -y_3 & & \ge & 0\\
& x_1 & & +x_3 & & & & & & & -z & \le & 1\\
& x_1 & & & +x_4 & & & & & & -z & \le & 1\\
& & x_2 & +x_3 & & & & & & & -z & \le & 1\\
& & x_2 & & +x_4 & & & & & & -z & \le & 1\\
& x_1 & & & & +x_5 & & & & & -z & \le & 1\\
& x_1 & & & & & +x_6 & & & & -z & \le & 1\\
& & x_2 & & & +x_5 & & & & & -z & \le & 1\\
& & x_2 & & & & +x_6 & & & & -z & \le & 1\\
& & & x_3 & & +x_5 & & & & & -z & \le & 1\\
& & & x_3 & & & +x_6 & & & & -z & \le & 1\\
& & & & x_4 & +x_5 & & & & & -z & \le & 1\\
& & & & x_4 & & +x_6 & & & & -z & \le & 1\\
& 5x_1 & +15x_2 & +25x_3 & +15x_4 & +10x_5 & +38x_6 & & & & +10z & \le & 50\\
& x_1, & x_2, & x_3, & x_4, & x_5, & x_6 & & & & & \in & \{0, 1\}\\
& & & & & & & y_1, & y_2, & y_3 & & \in & \{0, 1\}\\
& & & & & & & & & & z & \in & \{0, 1\}
\end{array}
\]
Constructive heuristic: the primal bound
Class by class; from the second class on, the first executed job also pays \(u\).
- Steps 1–2. Class 1: \(x[1] = x[2] = 1\), \(ra = 30\); class complete, \(y[1] = 1\).
- Step 3. Class 2: \(t_3 + u = 35 > 30\), skipped. Step 4. \(t_4 + u = 25 \le 30\): \(x[4] = 1\), \(z = 1\), \(ra = 5\).
- Steps 5–6. Class 3: \(t_5, t_6 > 5\), skipped.
Revenue \(10 + 5 + 12 + 5 = 32\): \(z(\mathit{MILP}) \ge 32\).
LP relaxation and dual: the dual bound
With \(\pi_j \le 0\) (all), \(\lambda_{ji} \ge 0\) (mixed), \(\mu \ge 0\) (availability):
\[
\begin{aligned}
\min ~~ \sum_{c<g}\ \sum_{j \in \mathscr{J}_c}\ \sum_{i \in \mathscr{J}_g} \lambda_{ji} + a\, \mu & & \\
\text{subject to} \quad \pi_j + \sum_{i \notin \mathscr{J}_c} \lambda_{ji} + t_j\, \mu &\ge r_j, & \forall c \in \{1, 2, \dots, q\},\ \forall j \in \mathscr{J}_c, \\
-\sum_{j \in \mathscr{J}_c} \pi_j &\ge v_c, & \forall c \in \{1, 2, \dots, q\}, \\
-\sum_{c<g}\ \sum_{j \in \mathscr{J}_c}\ \sum_{i \in \mathscr{J}_g} \lambda_{ji} + u\, \mu &\ge 0, & \\
\pi_j &\le 0, & \forall j \in \{1, 2, \dots, n\}, \\
\lambda_{ji} &\ge 0, & \forall c, g \in \{1, 2, \dots, q\},\ c < g,\ \forall j \in \mathscr{J}_c,\ \forall i \in \mathscr{J}_g, \\
\mu &\ge 0. &
\end{aligned}
\]
The same dual, written on the data of the instance:
\[
\begin{array}{rrrrrrrrrrrrrrrrrrrr c l}
\min & & & & & & & \lambda_{13} & +\lambda_{14} & +\lambda_{23} & +\lambda_{24} & +\lambda_{15} & +\lambda_{16} & +\lambda_{25} & +\lambda_{26} & +\lambda_{35} & +\lambda_{36} & +\lambda_{45} & +\lambda_{46} & +50\mu & & \\
\text{subject to} & \pi_1 & & & & & & +\lambda_{13} & +\lambda_{14} & & & +\lambda_{15} & +\lambda_{16} & & & & & & & +5\mu & \ge & 10\\
& & \pi_2 & & & & & & & +\lambda_{23} & +\lambda_{24} & & & +\lambda_{25} & +\lambda_{26} & & & & & +15\mu & \ge & 5\\
& & & \pi_3 & & & & +\lambda_{13} & & +\lambda_{23} & & & & & & +\lambda_{35} & +\lambda_{36} & & & +25\mu & \ge & 20\\
& & & & \pi_4 & & & & +\lambda_{14} & & +\lambda_{24} & & & & & & & +\lambda_{45} & +\lambda_{46} & +15\mu & \ge & 12\\
& & & & & \pi_5 & & & & & & +\lambda_{15} & & +\lambda_{25} & & +\lambda_{35} & & +\lambda_{45} & & +10\mu & \ge & 10\\
& & & & & & \pi_6 & & & & & & +\lambda_{16} & & +\lambda_{26} & & +\lambda_{36} & & +\lambda_{46} & +38\mu & \ge & 22\\
& -\pi_1 & -\pi_2 & & & & & & & & & & & & & & & & & & \ge & 5\\
& & & -\pi_3 & -\pi_4 & & & & & & & & & & & & & & & & \ge & 4\\
& & & & & -\pi_5 & -\pi_6 & & & & & & & & & & & & & & \ge & 10\\
& & & & & & & -\lambda_{13} & -\lambda_{14} & -\lambda_{23} & -\lambda_{24} & -\lambda_{15} & -\lambda_{16} & -\lambda_{25} & -\lambda_{26} & -\lambda_{35} & -\lambda_{36} & -\lambda_{45} & -\lambda_{46} & +10\mu & \ge & 0\\
& \pi_1, & \pi_2, & \pi_3, & \pi_4, & \pi_5, & \pi_6 & & & & & & & & & & & & & & \le & 0\\
& & & & & & & \lambda_{13}, & \lambda_{14}, & \lambda_{23}, & \lambda_{24}, & \lambda_{15}, & \lambda_{16}, & \lambda_{25}, & \lambda_{26}, & \lambda_{35}, & \lambda_{36}, & \lambda_{45}, & \lambda_{46} & & \ge & 0\\
& & & & & & & & & & & & & & & & & & & \mu & \ge & 0
\end{array}
\]
A hand-built dual solution. The bonus of every class loaded on one job: \(\bar\pi_1 = -5\), \(\bar\pi_3 = -4\), \(\bar\pi_5 = -10\); \(\bar\lambda = 0\); \(\bar\mu = \max_j (r_j - \bar\pi_j)/t_j = \max\{3, \tfrac{1}{3}, \tfrac{24}{25}, \tfrac{4}{5}, 2, \tfrac{11}{19}\} = 3\); value \(150\): \(32 \le z(\mathit{MILP}) \le 150\).
What the solver says. \(z(\mathit{LP}) = 5280/113 = 46.7\). Integer optimum \(42\): class 3 alone, complete, jobs 5 and 6 (\(48 \le 50\), \(z = 0\)), revenue \(10 + 22 + 10\). Heuristic gap \(24\%\).
| \(LB\) | \(UB\) (hand dual) | \(z(\mathit{LP})\) | \(z(\mathit{MILP})\) | heuristic gap |
|---|---|---|---|---|
| 32 | 150 | \(5280/113\) | 42 | \(23.8\%\) |
Additional considerations
- The optimality direction for \(y_c\) is imposed with \(\sum_{j \in \mathscr{J}_c} x_j - y_c \le |\mathscr{J}_c| - 1\) (\(q\) optimality-preserving, not valid, constraints).
- The optimality direction for \(z\): \(z \le \sum_{j \notin \mathscr{J}_c} x_j\) for every \(c\).
- An aggregated form with class variables \(w_c \ge x_j\) and \(\sum_c w_c - 1 \le (q-1) z\) replaces the pairs: same integer set, different relaxation.
Additional modelling questions
7.6.1 — Penalty for a class started and not finished
Starting a class without completing it costs \(w = 3\).
A worked variant: at least one job per class
A union agreement requires executing at least one job of every class.
A covering constraint per class, \(\sum_{j \in \mathscr{J}_c} x_j \ge 1\) for every \(c\) (\(q\) linear constraints). With \(q \ge 2\) classes all touched, for every pair of classes there is a mixed pair of executed jobs, and constraints force \(z = 1\): the reduction \(u\) is certain and it can be substituted in the capacity constraint (\(a - u\) instead of \(a\), eliminating \(z\)). On the instance the optimum drops from \(42\) to \(40\): for instance jobs \(1\), \(3\) and \(5\) (\(5 + 25 + 10 + 10 = 50\)), revenue \(40\), no complete class.
"At least one job per class" adds a family \(\omega_c \le 0\) to the dual, one per class, which enters the objective: being non-positive, it lowers it. Each \(\omega_c\) is pushed down as far as the columns of that class allow. The heuristic starts from the shortest job of each class — so it is feasible by construction — and then fills the remaining time.
| value | what it is | |
|---|---|---|
| \(\mathit{UB}\) | \(107\) | dual certificate built by hand |
| \(\mathit{LB}\) | \(32\) | heuristic solution |
| \(z(\mathit{LP})\) | \(\frac{309}{7}\) | relaxation without the bounds |
| \(z(\mathit{LP}^+)\) | \(\frac{309}{7}\) | relaxation with the bounds |
| \(z(\mathit{MILP})\) | \(40\) | optimum of the MILP |
Code
Full script: python/fam07_6_classesbonus.py;
notebook: notebooks/fam07_6_classesbonus.ipynb.
Show the complete script — python/fam07_6_classesbonus.py (241 lines)
"""Problem 7.6 -- Classes with completion bonus and if-and-only-if reduction.
Two "if and only if"s: each with one direction imposed by the constraint (via
CNF) and the other by the optimum -- the general scheme to model an iff.
"""
import gurobipy as gp
import numpy as np
import pandas as pd
from gurobipy import GRB
from euristiche import best_fit, first_fit, matrice, next_fit
from mip import (ammissibile, dualita_forte, due_rilassamenti, frazione,
nuovo_modello, registra_bound, rilassamenti, risolvi,
stampa_soluzione, valuta)
from stile import CICLO, ROSSO, intestazione, plt, salva_dati, salva_figura
from esteso import salva_modello
R = range
# ---------- 1. MODEL AND INSTANCE ----------
intestazione("6. Bonus if the whole class is executed; reduction u if and only if >= 2 classes")
r6 = [10, 5, 20, 12, 10, 22]
t6 = [5, 15, 25, 15, 10, 38]
J6 = [[0, 1], [2, 3], [4, 5]]
v6 = [5, 4, 10]
a6, u6 = 50, 10
salva_dati(pd.DataFrame({"job": R(1, 7), "r": r6, "t": t6,
"class": [c + 1 for j in R(6) for c in R(3) if j in J6[c]]}), "fam07_6_lavori")
salva_dati(pd.DataFrame({"class": R(1, 4), "v": v6}), "fam07_6_classi")
def coppie(J):
return [(j, i, c, g) for c in R(len(J)) for g in R(c + 1, len(J)) for j in J[c] for i in J[g]]
def modello_6(r, t, J, v, a, u):
n, q = len(r), len(J)
m = nuovo_modello("classi_premio")
x = m.addVars(n, vtype=GRB.BINARY, name="x")
y = m.addVars(q, vtype=GRB.BINARY, name="y")
z = m.addVar(vtype=GRB.BINARY, name="z")
m.setObjective(gp.quicksum(r[j] * x[j] for j in R(n)) + gp.quicksum(v[c] * y[c] for c in R(q)),
GRB.MAXIMIZE)
m.addConstrs((x[j] - y[c] >= 0 for c in R(q) for j in J[c]), name="all")
m.addConstrs((x[j] + x[i] - z <= 1 for (j, i, c, g) in coppie(J)), name="miste")
m.addConstr(gp.quicksum(t[j] * x[j] for j in R(n)) + u * z <= a, name="disponibilita")
return m, x, y, z
def duale_6(r, t, J, v, a, u):
"""min sum lam_ji + a mu; pi_j + sum lam + t_j mu >= r_j; -sum_{J_c} pi_j >= v_c;
-sum lam + u mu >= 0; pi <= 0, lam >= 0, mu >= 0."""
n, q = len(r), len(J)
cp = coppie(J)
d = nuovo_modello("duale_classi_premio")
pi = d.addVars(n, lb=-GRB.INFINITY, ub=0.0, name="pi")
lam = d.addVars([(j, i) for (j, i, _, _) in cp], name="lam")
mu = d.addVar(name="mu")
d.setObjective(lam.sum() + a * mu, GRB.MINIMIZE)
for j in R(n):
d.addConstr(pi[j] + gp.quicksum(lam[jj, ii] for (jj, ii, _, _) in cp if jj == j or ii == j)
+ t[j] * mu >= r[j], name=f"rc_x[{j}]")
d.addConstrs((-gp.quicksum(pi[j] for j in J[c]) >= v[c] for c in R(q)), name="rc_y")
d.addConstr(-lam.sum() + u * mu >= 0, name="rc_z")
return d
def euristica_6(r, t, J, v, a, u):
"""Class by class: from the second class on, the first job also pays the reduction u."""
n, q = len(r), len(J)
x, y, z, ra, passi = [0] * n, [0] * q, 0, a, []
for c in R(q):
cnt = 0
for j in J[c]:
if c == 0 or z == 1:
if t[j] <= ra:
x[j], ra, cnt = 1, ra - t[j], cnt + 1
passi.append(f"Class {c + 1}: t[{j + 1}] = {t[j]} <= ra; x[{j + 1}] = 1, ra = {ra}.")
else:
passi.append(f"Class {c + 1}: t[{j + 1}] = {t[j]} > ra = {ra}; job {j + 1} is skipped.")
else:
if t[j] + u <= ra:
x[j], z, ra, cnt = 1, 1, ra - t[j] - u, cnt + 1
passi.append(f"Class {c + 1}, reduction not applied yet: t[{j + 1}] + u = {t[j] + u} <= ra; "
f"x[{j + 1}] = 1, z = 1, ra = {ra}.")
else:
passi.append(f"Class {c + 1}, reduction not applied yet: t[{j + 1}] + u = {t[j] + u} > ra = {ra}; "
f"job {j + 1} is skipped.")
if cnt == len(J[c]):
y[c] = 1
passi.append(f"All the jobs of class {c + 1} are executed: y[{c + 1}] = 1 (bonus v = {v[c]}).")
return x, y, z, passi
m6, x6, y6, z6 = modello_6(r6, t6, J6, v6, a6, u6)
salva_modello(m6, "fam07_6_primale")
# ---------- 2. THE LP RELAXATION ----------
zlp6, zlp6r, _ = rilassamenti(m6)
# ---------- 3. THE DUAL OF THE RELAXATION (UPPER BOUND: IT IS A MAXIMUM) ----------
d6 = duale_6(r6, t6, J6, v6, a6, u6)
salva_modello(d6, "fam07_6_duale")
pi_mano = {f"pi[{J6[c][0]}]": -v6[c] for c in R(3)} # the first job of every class carries the bonus
mu_mano = max((r6[j] - pi_mano.get(f"pi[{j}]", 0)) / t6[j] for j in R(6))
mano = dict(pi_mano, mu=mu_mano)
ub6, viol = valuta(d6, mano)
assert viol <= 1e-9
print(f"Hand-built dual solution: pi_1 = -5, pi_3 = -4, pi_5 = -10, lam = 0, "
f"mu = max_j (r_j - pi_j)/t_j = {frazione(mu_mano)} -> ub = {frazione(ub6)}")
dualita_forte(d6, zlp6)
# ---------- 4. CONSTRUCTIVE HEURISTIC (LOWER BOUND: IT IS A MAXIMUM) ----------
xe, ye, ze, passi = euristica_6(r6, t6, J6, v6, a6, u6)
print("Class-by-class heuristic:")
for i, s in enumerate(passi, 1):
print(f" Step {i}. {s}")
lb6 = sum(r6[j] * xe[j] for j in R(6)) + sum(v6[c] * ye[c] for c in R(3))
print(f" lb = {lb6} (x = {xe}, y = {ye}, z = {ze})")
# ---------- 5. OPTIMAL SOLUTION OF THE MILP ----------
z6v = risolvi(m6)
print("Optimal solution of the MILP:")
stampa_soluzione(m6, solo_non_nulle=True)
riga = registra_bound("6 classes bonus", ub6, lb6, zlp6, zlp6r, z6v, senso="max")
salva_dati(pd.DataFrame([riga]), "fam07_6_bound")
# ---------- 6. ADDITIONAL MODELLING QUESTIONS ----------
varianti = {}
def variante(nome, m):
z = risolvi(m)
print(f" {nome:70s} z = {frazione(z)}")
return z
# 6a: at least one job per class
m, x, y, z = modello_6(r6, t6, J6, v6, a6, u6)
m.addConstrs((gp.quicksum(x[j] for j in J6[c]) >= 1 for c in R(3)), name="at_least_one")
varianti["6a"] = variante("6a. At least one job per class (hence z = 1)", m)
# 6b: penalty w per class started and not completed
w6 = 3
m, x, y, z = modello_6(r6, t6, J6, v6, a6, u6)
st = m.addVars(3, vtype=GRB.BINARY, name="s")
m.addConstrs((st[c] >= x[j] for c in R(3) for j in J6[c]), name="started")
m.update()
m.setObjective(m.getObjective() - w6 * gp.quicksum(st[c] - y[c] for c in R(3)), GRB.MAXIMIZE)
varianti["6b"] = variante("6b. Penalty 3 per class started and not completed (s_c >= x_j)", m)
salva_dati(pd.DataFrame({"variant": list(varianti), "z": list(varianti.values())}), "fam07_6_varianti")
# ---------- 7. THE SANDWICH ON THE VARIANT 6a ----------
intestazione("6a. The sandwich on the variant: at least one job per class")
def modello_6a(r, t, J, v, a, u):
mm_, xx, yy, zz = modello_6(r, t, J, v, a, u)
mm_.addConstrs((gp.quicksum(xx[j] for j in J[c]) >= 1 for c in R(len(J))), name="almeno_uno")
return mm_, xx, yy, zz
def duale_6a(r, t, J, v, a, u):
"""To the dual of 7.6 one adds omega_c <= 0 for every constraint
sum_{j in J_c} x_j >= 1 (a >= direction in a maximisation). The right-hand
side is 1, so omega enters the objective --- and, being non-positive, lowers it."""
nn, q = len(r), len(J)
cp = coppie(J)
d = nuovo_modello("duale_classi_premio_6a")
pi = d.addVars(nn, lb=-GRB.INFINITY, ub=0.0, name="pi")
lam = d.addVars([(j, i) for (j, i, _, _) in cp], name="lam")
mu = d.addVar(name="mu")
om = d.addVars(q, lb=-GRB.INFINITY, ub=0.0, name="omega")
d.setObjective(lam.sum() + a * mu + om.sum(), GRB.MINIMIZE)
for c in R(q):
for j in J[c]:
d.addConstr(pi[j] + gp.quicksum(lam[jj, ii] for (jj, ii, _, _) in cp
if jj == j or ii == j)
+ t[j] * mu + om[c] >= r[j], name=f"rc_x[{j}]")
d.addConstrs((-gp.quicksum(pi[j] for j in J[c]) >= v[c] for c in R(q)), name="rc_y")
d.addConstr(-lam.sum() + u * mu >= 0, name="rc_z")
return d
m6a, x6a, y6a, z6a = modello_6a(r6, t6, J6, v6, a6, u6)
salva_modello(m6a, "fam07_6a_primale")
# -- feasible heuristic: the base one, repaired by adding the missing job --
print("Constructive heuristic: the new constraint forces every class, so one starts from the")
print("shortest job of each class --- so the solution is feasible by")
print("construction --- and then the others are added by decreasing ratio r_j/t_j.")
def tempo_6a(scelti):
"""The time spent: the jobs plus the reduction u, which triggers with two classes."""
classi = {c for c in R(3) for j in J6[c] if j in scelti}
return sum(t6[j] for j in scelti) + (u6 if len(classi) >= 2 else 0)
scelti_6a = [min(J6[c], key=lambda j: t6[j]) for c in R(3)]
print(f" one job per class: {sorted(j + 1 for j in scelti_6a)}, "
f"time {tempo_6a(scelti_6a)} out of {a6}")
for j in sorted(set(R(6)) - set(scelti_6a), key=lambda j: -r6[j] / t6[j]):
if tempo_6a(scelti_6a + [j]) <= a6:
scelti_6a.append(j)
print(f" job {j + 1} enters as well, time {tempo_6a(scelti_6a)}")
y_6a = [1 if all(j in scelti_6a for j in J6[c]) else 0 for c in R(3)]
z_6a = 1 if len({c for c in R(3) if any(j in scelti_6a for j in J6[c])}) >= 2 else 0
lb6a = sum(r6[j] for j in scelti_6a) + sum(v6[c] * y_6a[c] for c in R(3))
sol_6a = ({f"x[{j}]": 1 for j in scelti_6a}
| {f"y[{c}]": y_6a[c] for c in R(3)} | {"z": z_6a})
assert ammissibile(m6a, sol_6a), "the heuristic solution of the variant must be feasible"
print(f" jobs {sorted(j + 1 for j in scelti_6a)} -> lb = {frazione(lb6a)}")
# -- dual certificate: omega discounts every class --
d6a = duale_6a(r6, t6, J6, v6, a6, u6)
salva_modello(d6a, "fam07_6a_duale")
pi_6a = {f"pi[{J6[c][0]}]": -v6[c] for c in R(3)}
mu_6a = max((r6[j] - pi_6a.get(f"pi[{j}]", 0)) / t6[j] for j in R(6))
om_6a = {c: min(0.0, max(r6[j] - pi_6a.get(f"pi[{j}]", 0) - t6[j] * mu_6a for j in J6[c]))
for c in R(3)}
mano_6a = dict(pi_6a, mu=mu_6a) | {f"omega[{c}]": om_6a[c] for c in R(3)}
ub6a, viol_6a = valuta(d6a, mano_6a)
assert viol_6a <= 1e-9, viol_6a
print("Dual solution by hand: pi and mu as in the base problem; then omega_c is lowered as")
print(" far as the columns of class c allow, that is")
print(" omega_c = max_{j in J_c} (r_j - pi_j - t_j mu), which is non-positive.")
print(f" omega = {[frazione(om_6a[c]) for c in R(3)]}")
print(f" -> ub = {frazione(ub6a)} (with omega = 0 it would be {frazione(a6 * mu_6a)})")
zlp6a, zlp6ar, _ = due_rilassamenti(m6a, d6a)
z6a_val = risolvi(m6a)
riga_6a = registra_bound("6a at least one job per class", ub6a, lb6a, zlp6a, zlp6ar, z6a_val,
senso="max")
salva_dati(pd.DataFrame([riga_6a]), "fam07_6a_bound")
assert lb6a <= z6a_val <= zlp6a + 1e-9 <= ub6a + 1e-9
print("Fine.")
print("Done.")