# Variable binaire: installer antenne ou non
    model.build = Var(model.Locations, within=Binary)
    
    # Objectif: minimiser coût total
    model.total_cost = Objective(
        expr=sum(model.cost[loc] * model.build[loc] for loc in model.Locations),
        sense=minimize
    )
    
    # Contrainte: chaque zone doit être couverte
    def coverage_rule(model, zone):
        return sum(model.covers[loc, zone] * model.build[loc] 
                  for loc in model.Locations) >= 1
    model.coverage_constraint = Constraint(model.Zones, rule=coverage_rule)
    
    # Contrainte de redondance (optionnel): double couverture pour zones critiques
    model.CriticalZones = Set(within=model.Zones, initialize=[1, 10, 20, 30])
    
    def redundancy_rule(model, zone):
        return sum(model.covers[loc, zone] * model.build[loc] 
                  for loc in model.Locations) >= 2
    model.redundancy = Constraint(model.CriticalZones, rule=redundancy_rule)
    
    return model


# === Agriculture - Planification de cultures ===

def crop_planning():
    """Planification optimale des cultures agricoles"""
    model = ConcreteModel(name="Crop Planning")
    
    # Cultures et parcelles
    model.Crops = Set(initialize=['Wheat', 'Corn', 'Soybean', 'Rice'])
    model.Fields = RangeSet(1, 10)
    model.Seasons = Set(initialize=['Spring', 'Summer', 'Fall'])
    
    # Surface des parcelles (hectares)
    model.field_area = Param(model.Fields, initialize={
        i: 50 + i*5 for i in model.Fields
    })
    
    # Rendement par culture et saison (kg/ha)
    yield_data = {
        ('Wheat', 'Spring'): 3000, ('Wheat', 'Summer'): 0, ('Wheat', 'Fall'): 2500,
        ('Corn', 'Spring'): 0, ('Corn', 'Summer'): 5000, ('Corn', 'Fall'): 0,
        ('Soybean', 'Spring'): 2000, ('Soybean', 'Summer'): 2500, ('Soybean', 'Fall'): 1800,
        ('Rice', 'Spring'): 4000, ('Rice', 'Summer'): 4500, ('Rice', 'Fall'): 0,
    }
    model.yield_rate = Param(model.Crops, model.Seasons, initialize=yield_data)
    
    # Prix de vente ($/kg)
    model.price = Param(model.Crops, initialize={
        'Wheat': 0.25, 'Corn': 0.20, 'Soybean': 0.40, 'Rice': 0.30
    })
    
    # Coût de production ($/ha)
    model.prod_cost = Param(model.Crops, initialize={
        'Wheat': 500, 'Corn': 600, 'Soybean': 450, 'Rice': 700
    })
    
    # Besoin en eau (m³/ha)
    model.water_need = Param(model.Crops, initialize={
        'Wheat': 450, 'Corn': 550, 'Soybean': 400, 'Rice': 1200
    })
    
    # Disponibilité en eau par saison
    model.water_available = Param(model.Seasons, initialize={
        'Spring': 15000, 'Summer': 20000, 'Fall': 12000
    })
    
    # Variables: surface cultivée
    model.planted = Var(model.Crops, model.Fields, model.Seasons, 
                       within=NonNegativeReals)
    
    # Objectif: maximiser profit
    def profit_rule(model):
        revenue = sum(model.price[c] * model.yield_rate[c,s] * model.planted[c,f,s]
                     for c in model.Crops for f in model.Fields for s in model.Seasons)
        cost = sum(model.prod_cost[c] * model.planted[c,f,s]
                  for c in model.Crops for f in model.Fields for s in model.Seasons)
        return revenue - cost
    
    model.profit = Objective(rule=profit_rule, sense=maximize)
    
    # Contrainte: surface par parcelle
    def field_capacity_rule(model, f, s):
        return sum(model.planted[c,f,s] for c in model.Crops) <= model.field_area[f]
    model.field_capacity = Constraint(model.Fields, model.Seasons, 
                                     rule=field_capacity_rule)
    
    # Contrainte: disponibilité en eau
    def water_constraint_rule(model, s):
        return sum(model.water_need[c] * model.planted[c,f,s]
                  for c in model.Crops for f in model.Fields) <= model.water_available[s]
    model.water_limit = Constraint(model.Seasons, rule=water_constraint_rule)
    
    # Contrainte: rotation des cultures (pas même culture 2 saisons consécutives)
    # Simplifié: entre Spring et Summer
    def rotation_rule(model, c, f):
        return model.planted[c,f,'Spring'] + model.planted[c,f,'Summer'] <= model.field_area[f]
    model.rotation = Constraint(model.Crops, model.Fields, rule=rotation_rule)
    
    return model


# === Logistique - Cross-docking ===

def cross_docking_optimization():
    """Optimisation de plateforme de cross-docking"""
    model = ConcreteModel(name="Cross-Docking")
    
    # Fournisseurs, produits, clients
    model.Suppliers = Set(initialize=['S1', 'S2', 'S3'])
    model.Products = Set(initialize=['P1', 'P2', 'P3', 'P4'])
    model.Customers = Set(initialize=['C1', 'C2', 'C3', 'C4', 'C5'])
    model.TimeSlots = RangeSet(1, 12)  # Créneaux horaires
    
    # Quantité disponible par fournisseur
    supply_data = {
        ('S1', 'P1'): 100, ('S1', 'P2'): 80,
        ('S2', 'P2'): 120, ('S2', 'P3'): 100,
        ('S3', 'P3'): 90, ('S3', 'P4'): 110,
    }
    model.supply = Param(model.Suppliers, model.Products, 
                        initialize=supply_data, default=0)
    
    # Demande des clients
    demand_data = {
        ('C1', 'P1'): 30, ('C1', 'P2'): 20,
        ('C2', 'P2'): 40, ('C2', 'P3'): 25,
        ('C3', 'P1'): 25, ('C3', 'P4'): 35,
        ('C4', 'P3'): 30, ('C4', 'P4'): 40,
        ('C5', 'P1'): 20, ('C5', 'P4'): 30,
    }
    model.demand = Param(model.Customers, model.Products,
                        initialize=demand_data, default=0)
    
    # Coût de manutention
    model.handling_cost = Param(initialize=2)  # $ par unité
    
    # Capacité de la plateforme par créneau
    model.platform_capacity = Param(initialize=200)
    
    # Variables
    # Flux direct: fournisseur -> client (idéal pour cross-docking)
    model.direct_flow = Var(model.Suppliers, model.Customers, model.Products,
                           model.TimeSlots, within=NonNegativeReals)
    
    # Stockage temporaire (à minimiser)
    model.storage = Var(model.Products, model.TimeSlots, within=NonNegativeReals)
    
    # Objectif: minimiser coût de manutention
    def cost_rule(model):
        handling = sum(model.handling_cost * model.direct_flow[s,c,p,t]
                      for s in model.Suppliers for c in model.Customers
                      for p in model.Products for t in model.TimeSlots)
        storage_penalty = 5 * sum(model.storage[p,t] 
                                 for p in model.Products for t in model.TimeSlots)
        return handling + storage_penalty
    
    model.total_cost = Objective(rule=cost_rule, sense=minimize)
    
    # Contrainte: respecter l'offre
    def supply_rule(model, s, p):
        return sum(model.direct_flow[s,c,p,t] 
                  for c in model.Customers for t in model.TimeSlots) <= model.supply[s,p]
    model.supply_constraint = Constraint(model.Suppliers, model.Products,
                                        rule=supply_rule)
    
    # Contrainte: satisfaire la demande
    def demand_rule(model, c, p):
        return sum(model.direct_flow[s,c,p,t]
                  for s in model.Suppliers for t in model.TimeSlots) >= model.demand[c,p]
    model.demand_constraint = Constraint(model.Customers, model.Products,
                                        rule=demand_rule)
    
    # Contrainte: capacité de la plateforme
    def capacity_rule(model, t):
        return sum(model.direct_flow[s,c,p,t]
                  for s in model.Suppliers for c in model.Customers
                  for p in model.Products) <= model.platform_capacity
    model.capacity_constraint = Constraint(model.TimeSlots, rule=capacity_rule)
    
    return model


[OK] CALLBACKS ET ÉVÉNEMENTS

# === Callbacks pendant la résolution ===

from pyomo.opt import SolverFactory, TerminationCondition

class SolutionCallback:
    """Classe pour gérer les callbacks de solution"""
    
    def __init__(self):
        self.solutions = []
        self.iteration = 0
    
    def __call__(self, model, where=None):
        """Appelé pendant la résolution"""
        self.iteration += 1
        
        if where == 'mip_sol':
            # Nouvelle solution entière trouvée
            obj_value = model.obj()
            print(f"Iteration {self.iteration}: Nouvelle solution = {obj_value:.2f}")
            
            # Sauvegarder solution
            solution = {v.name: v.value for v in model.component_data_objects(Var)
                       if v.value is not None}
            self.solutions.append({
                'iteration': self.iteration,
                'objective': obj_value,
                'solution': solution
            })
        
        elif where == 'mip_node':
            # Nouveau nœud de branch-and-bound
            if self.iteration % 100 == 0:
                print(f"Nœud {self.iteration} exploré")

# Utilisation
callback = SolutionCallback()
# Note: Support de callbacks dépend du solveur


[OK] TECHNIQUES DE DÉCOMPOSITION

# === Décomposition de Dantzig-Wolfe ===

def dantzig_wolfe_decomposition():
    """Exemple de décomposition de Dantzig-Wolfe"""
    
    # Problème maître
    master = ConcreteModel(name="Master Problem")
    
    # Variables: poids des colonnes (solutions du sous-problème)
    master.Columns = Set(initialize=range(1, 10))  # Colonnes initiales
    master.lambda_var = Var(master.Columns, within=NonNegativeReals)
    
    # Contraintes de convexité
    master.convexity = Constraint(
        expr=sum(master.lambda_var[k] for k in master.Columns) == 1
    )
    
    # Sous-problème
    def solve_subproblem(dual_values):
        """Résoudre sous-problème avec valeurs duales du maître"""
        subproblem = ConcreteModel(name="Subproblem")
        
        # ... définir sous-problème ...
        # Objectif modifié avec valeurs duales
        
        solver = SolverFactory('glpk')
        results = solver.solve(subproblem)
        
        # Extraire nouvelle colonne
        new_column = extract_column(subproblem)
        return new_column
    
    # Algorithme de génération de colonnes
    iteration = 0
    max_iterations = 100
    
    while iteration < max_iterations:
        # Résoudre problème maître
        solver = SolverFactory('glpk')
        results = solver.solve(master)
        
        # Obtenir valeurs duales
        dual_values = get_dual_values(master)
        
        # Résoudre sous-problème
        new_column = solve_subproblem(dual_values)
        
        # Vérifier optimalité (reduced cost < 0)
        if new_column['reduced_cost'] >= -1e-6:
            print(f"Optimal trouvé après {iteration} itérations")
            break
        
        # Ajouter nouvelle colonne au maître
        add_column_to_master(master, new_column)
        
        iteration += 1

# === Décomposition de Benders ===

def benders_decomposition_algorithm(data):
    """Algorithme de décomposition de Benders"""
    
    # Problème maître
    master = ConcreteModel(name="Benders Master")
    master.I = Set(initialize=data['locations'])
    master.x = Var(master.I, within=Binary)  # Décisions de premier niveau
    master.eta = Var()  # Approximation coût second niveau
    
    # Objectif maître
    master.obj = Objective(
        expr=sum(data['fixed_cost'][i] * master.x[i] for i in master.I) + master.eta,
        sense=minimize
    )
    
    # Contraintes initiales
    master.budget = Constraint(
        expr=sum(master.x[i] for i in master.I) <= data['max_facilities']
    )
    
    # Liste des coupes
    master.optimality_cuts = ConstraintList()
    master.feasibility_cuts = ConstraintList()
    
    # Sous-problème
    def create_subproblem(x_values):
        """Créer sous-problème pour x fixé"""
        sub = ConcreteModel(name="Benders Subproblem")
        
        sub.I = Set(initialize=data['locations'])
        sub.J = Set(initialize=data['customers'])
        
        # Variables de second niveau
        sub.y = Var(sub.I, sub.J, within=NonNegativeReals)
        
        # Objectif
        def sub_obj_rule(model):
            return sum(data['transport_cost'][i,j] * model.y[i,j]
                      for i in model.I for j in model.J)
        sub.obj = Objective(rule=sub_obj_rule, sense=minimize)
        
        # Contraintes avec x fixé
        def capacity_rule(model, i):
            return sum(model.y[i,j] for j in model.J) <= data['capacity'][i] * x_values[i]
        sub.capacity = Constraint(sub.I, rule=capacity_rule)
        
        def demand_rule(model, j):
            return sum(model.y[i,j] for i in sub.I) >= data['demand'][j]
        sub.demand = Constraint(sub.J, rule=demand_rule)
        
        return sub
    
    # Algorithme principal
    iteration = 0
    max_iterations = 50
    tolerance = 1e-4
    
    lower_bound = -float('inf')
    upper_bound = float('inf')
    
    solver = SolverFactory('glpk')
    
    while iteration < max_iterations and upper_bound - lower_bound > tolerance:
        iteration += 1
        print(f"\n=== Iteration {iteration} ===")
        
        # Résoudre maître
        results = solver.solve(master)
        x_values = {i: master.x[i].value for i in master.I}
        lower_bound = master.obj()
        
        print(f"Lower bound: {lower_bound:.2f}")
        
        # Créer et résoudre sous-problème
        subproblem = create_subproblem(x_values)
        sub_results = solver.solve(subproblem)
        
        if sub_results.solver.termination_condition == TerminationCondition.optimal:
            # Sous-problème optimal
            sub_obj_value = subproblem.obj()
            current_obj = sum(data['fixed_cost'][i] * x_values[i] 
                            for i in master.I) + sub_obj_value
            
            upper_bound = min(upper_bound, current_obj)
            print(f"Upper bound: {upper_bound:.2f}")
            
            # Ajouter coupe d'optimalité
            # eta >= sub_obj + gradient * (x - x_current)
            # Simplifié ici
            
        else:
            # Sous-problème infaisable -> ajouter coupe de faisabilité
            print("Sous-problème infaisable")
            # Ajouter coupe de faisabilité
        
        print(f"Gap: {upper_bound - lower_bound:.2f}")
    
    print(f"\nConvergence après {iteration} itérations")
    print(f"Solution optimale: {upper_bound:.2f}")
    
    return master, x_values


[OK] OPTIMISATION SOUS INCERTITUDE

# === Programmation stochastique multi-scénarios ===

def multi_stage_stochastic_programming():
    """Programmation stochastique à plusieurs étages"""
    model = ConcreteModel(name="Multi-Stage Stochastic")
    
    # Étages temporels
    model.Stages = RangeSet(1, 3)
    
    # Scénarios avec probabilités
    model.Scenarios = Set(initialize=['Low', 'Medium', 'High'])
    model.prob = Param(model.Scenarios, initialize={
        'Low': 0.2, 'Medium': 0.5, 'High': 0.3
    })
    
    # Décisions de production par étage et scénario
    model.production = Var(model.Stages, model.Scenarios, within=NonNegativeReals)
    
    # Capacité installée (décision de premier étage - ici et maintenant)
    model.capacity = Var(within=NonNegativeReals)
    
    # Demande dépendante du scénario et de l'étage
    demand_data = {
        (1, 'Low'): 100, (1, 'Medium'): 120, (1, 'High'): 140,
        (2, 'Low'): 110, (2, 'Medium'): 140, (2, 'High'): 170,
        (3, 'Low'): 105, (3, 'Medium'): 135, (3, 'High'): 165,
    }
    model.demand = Param(model.Stages, model.Scenarios, initialize=demand_data)
    
    # Coûts
    model.capacity_cost = Param(initialize=1000)
    model.production_cost = Param(initialize=50)
    model.shortage_penalty = Param(initialize=200)
    
    # Variables de manque
    model.shortage = Var(model.Stages, model.Scenarios, within=NonNegativeReals)
    
    # Objectif: minimiser coût espéré total
    def cost_rule(model):
        capacity_cost = model.capacity_cost * model.capacity
        
        expected_operating_cost = sum(
            model.prob[s] * (
                sum(model.production_cost * model.production[t,s] for t in model.Stages) +
                sum(model.shortage_penalty * model.shortage[t,s] for t in model.Stages)
            )
            for s in model.Scenarios
        )
        
        return capacity_cost + expected_operating_cost
    
    model.total_cost = Objective(rule=cost_rule, sense=minimize)
    
    # Contrainte: non-anticipativité (décisions avant observation)
    # Les décisions du stage 1 doivent être identiques pour tous les scénarios
    def non_anticipativity_rule(model, s):
        if s == 'Low':
            return Constraint.Skip
        return model.production[1,s] == model.production[1,'Low']
    model.non_anticipativity = Constraint(model.Scenarios, rule=non_anticipativity_rule)
    
    # Contrainte: capacité
    def capacity_rule(model, t, s):
        return model.production[t,s] <= model.capacity
    model.capacity_constraint = Constraint(model.Stages, model.Scenarios, 
                                          rule=capacity_rule)
    
    # Contrainte: demande
    def demand_rule(model, t, s):
        return model.production[t,s] + model.shortage[t,s] >= model.demand[t,s]
    model.demand_constraint = Constraint(model.Stages, model.Scenarios,
                                        rule=demand_rule)
    
    return model

# === Optimisation robuste avec ensembles d'incertitude ===

def robust_optimization_uncertainty_sets():
    """Optimisation robuste avec ensembles d'incertitude"""
    model = ConcreteModel(name="Robust Optimization")
    
    # Décisions
    model.I = RangeSet(1, 5)
    model.x = Var(model.I, within=NonNegativeReals)
    
    # Paramètres nominaux
    model.demand_nominal = Param(model.I, initialize={
        1: 100, 2: 120, 3: 90, 4: 110, 5: 105
    })
    
    # Déviation maximale (intervalle d'incertitude)
    model.demand_deviation = Param(model.I, initialize={
        1: 20, 2: 25, 3: 15, 4: 20, 5: 18
    })
    
    # Budget d'incertitude Γ (contrôle conservatisme)
    # Γ = 0: solution nominale, Γ = |I|: très conservateur
    model.Gamma = Param(initialize=2.5)
    
    # Variables auxiliaires pour reformulation robuste
    model.z = Var(model.I, within=NonNegativeReals)
    model.p = Var(within=NonNegativeReals)
    
    # Coût
    model.cost = Param(model.I, initialize={i: 10 + i for i in model.I})
    
    # Objectif
    model.obj = Objective(
        expr=sum(model.cost[i] * model.x[i] for i in model.I),
        sense=minimize
    )
    
    # Contrainte robuste (reformulée avec dualité)
    # Original: x[i] >= demand[i] + u[i] * deviation[i], pour tout u tel que ||u||_1 <= Gamma
    # Reformulation robuste:
    def robust_constraint_rule(model, i):
        return model.x[i] >= model.demand_nominal[i] + model.z[i]
    model.robust_demand = Constraint(model.I, rule=robust_constraint_rule)
    
    # Contrainte du budget d'incertitude
    def uncertainty_budget_rule(model):
        return (sum(model.z[i] / model.demand_deviation[i] for i in model.I) + 
                model.p <= model.Gamma)
    model.uncertainty_budget = Constraint(rule=uncertainty_budget_rule)
    
    # Contrainte de liaison
    def linking_rule(model, i):
        return model.z[i] <= model.demand_deviation[i] + model.p
    model.linking = Constraint(model.I, rule=linking_rule)
    
    return model


[OK] ASTUCES ET BEST PRACTICES FINALES

# === Performance Tips ===

# 1. Utiliser sets filtrés pour éviter variables/contraintes inutiles
def filtered_sets_example():
    model = ConcreteModel()
    model.I = Set(initialize=range(1, 101))
    model.J = Set(initialize=range(1, 101))
    
    # Mauvais: toutes les paires
    # model.x = Var(model.I, model.J, within=Binary)  # 10,000 variables!
    
    # Bon: seulement paires valides
    model.ValidPairs = Set(within=model.I * model.J,
                          initialize=[(i,j) for i in model.I for j in model.J if i < j])
    model.x = Var(model.ValidPairs, within=Binary)  # ~5,000 variables

# 2. Précalculer expressions complexes
def precompute_expressions():
    model = ConcreteModel()
    model.I = RangeSet(1, 1000)
    model.x = Var(model.I, within=NonNegativeReals)
    
    # Calculer une seule fois
    model.total_x = Expression(expr=sum(model.x[i] for i in model.I))
    
    # Utiliser dans plusieurs contraintes
    model.c1 = Constraint(expr=model.total_x <= 1000)
    model.c2 = Constraint(expr=model.total_x >= 500)

# 3. Indexation efficace
def efficient_indexing():
    # Utiliser tuples pour indices multiples
    data = {(1, 'A'): 10, (1, 'B'): 15, (2, 'A'): 20}
    
    # Plus rapide que dictionnaires imbriqués
    model.param = Param(model.I, model.J, initialize=data, default=0)

# 4. Éviter boucles Python dans règles
def avoid_python_loops():
    model = ConcreteModel()
    model.I = RangeSet(1, 1000)
    model.x = Var(model.I)
    
    # Lent
    # def rule(model, i):
    #     result = 0
    #     for j in model.I:
    #         result += model.x[j]
    #     return result
    
    # Rapide
    def rule(model, i):
        return quicksum(model.x[j] for j in model.I)

# === Documentation du modèle ===

def document_model(model, filename='model_documentation.txt'):
    """Générer documentation complète du modèle"""
    with open(filename, 'w') as f:
        f.write("="*70 + "\n")
        f.write("DOCUMENTATION DU MODÈLE\n")
        f.write("="*70 + "\n\n")
        
        # Nom et description
        f.write(f"Nom: {model.name}\n\n")
        
        # Sets
        f.write("ENSEMBLES (SETS)\n")
        f.write("-"*70 + "\n")
        for s in model.component_objects(Set):
            f.write(f"{s.name}: {list(s)[:10]}")  # Premiers 10 éléments
            if len(s) > 10:
                f.write(f" ... ({len(s)} éléments)")
            f.write("\n")
        f.write("\n")
        
        # Parameters
        f.write("PARAMÈTRES\n")
        f.write("-"*70 + "\n")
        for p in model.component_objects(Param):
            f.write(f"{p.name}")
            if p.is_indexed():
                f.write(f" (indexé)")
            f.write("\n")
        f.write("\n")
        
        # Variables
        f.write("VARIABLES    print(f"\nStatut: {results.solver.status}")
    print(f"Condition de terminaison: {results.solver.termination_condition}")
    print(f"Valeur objectif: ${model.total_cost():.2f}")
    
    print("\n" + "="*60)
    print("ÉTAPE 4: Analyse des résultats")
    print("="*60)
    
    # Statistiques de production
    total_production = sum(model.production[p,t].value 
                          for p in model.Products for t in model.Periods)
    print(f"Production totale: {total_production:.0f} unités")
    
    # Statistiques d'inventaire
    max_inventory = max(model.inventory[p,t].value 
                       for p in model.Products for t in model.Periods)
    print(f"Inventaire maximum: {max_inventory:.0f} unités")
    
    print("\n" + "="*60)
    print("ÉTAPE 5: Export des résultats")
    print("="*60)
    
    # Export vers CSV
    production_data = []
    for p in model.Products:
        for t in model.Periods:
            production_data.append({
                'Product': p,
                'Period': t,
                'Production': model.production[p,t].value,
                'Inventory': model.inventory[p,t].value,
                'Demand': model.demand[p,t]
            })
    
    results_df = pd.DataFrame(production_data)
    results_df.to_csv('production_plan.csv', index=False)
    print("Résultats exportés vers 'production_plan.csv'")
    
    # Visualisation
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    for p in model.Products:
        periods = list(model.Periods)
        production = [model.production[p,t].value for t in periods]
        inventory = [model.inventory[p,t].value for t in periods]
        
        axes[0].plot(periods, production, marker='o', label=p)
        axes[1].plot(periods, inventory, marker='s', label=p)
    
    axes[0].set_title('Production par période')
    axes[0].set_xlabel('Période')
    axes[0].set_ylabel('Quantité')
    axes[0].legend()
    axes[0].grid(True)
    
    axes[1].set_title('Inventaire par période')
    axes[1].set_xlabel('Période')
    axes[1].set_ylabel('Quantité')
    axes[1].legend()
    axes[1].grid(True)
    
    plt.tight_layout()
    plt.savefig('production_analysis.png', dpi=300)
    print("Visualisation sauvegardée dans 'production_analysis.png'")
    
    print("\n" + "="*60)
    print("WORKFLOW TERMINÉ")
    print("="*60)
    
    return model, results_df


[OK] PATTERNS DE CONCEPTION AVANCÉS

# === Factory Pattern pour modèles ===

class ModelFactory:
    """Factory pour créer différents types de modèles"""
    
    @staticmethod
    def create_model(model_type, **kwargs):
        if model_type == 'production':
            return ProductionModel(**kwargs)
        elif model_type == 'transportation':
            return TransportationModel(**kwargs)
        elif model_type == 'scheduling':
            return SchedulingModel(**kwargs)
        else:
            raise ValueError(f"Unknown model type: {model_type}")

class BaseOptimizationModel:
    """Classe de base pour modèles d'optimisation"""
    
    def __init__(self, name):
        self.name = name
        self.model = None
        self.results = None
        self.solver = None
    
    def build(self):
        """Construire le modèle - à implémenter dans les sous-classes"""
        raise NotImplementedError
    
    def solve(self, solver_name='glpk', **options):
        """Résoudre le modèle"""
        if self.model is None:
            self.build()
        
        self.solver = SolverFactory(solver_name)
        for key, value in options.items():
            self.solver.options[key] = value
        
        self.results = self.solver.solve(self.model, tee=True)
        return self.results
    
    def get_solution(self):
        """Extraire la solution - à implémenter dans les sous-classes"""
        raise NotImplementedError
    
    def validate_solution(self):
        """Valider la solution"""
        if self.results is None:
            raise ValueError("No solution available")
        
        # Vérifier toutes les contraintes
        for c in self.model.component_data_objects(Constraint, active=True):
            if c.body is not None and c.upper is not None:
                if c.body() > c.upper + 1e-6:
                    print(f"Constraint {c.name} violated: {c.body()} > {c.upper}")
                    return False
            if c.body is not None and c.lower is not None:
                if c.body() < c.lower - 1e-6:
                    print(f"Constraint {c.name} violated: {c.body()} < {c.lower}")
                    return False
        
        return True

class ProductionModel(BaseOptimizationModel):
    """Modèle de production spécialisé"""
    
    def __init__(self, products, periods, **kwargs):
        super().__init__("Production Model")
        self.products = products
        self.periods = periods
        self.data = kwargs
    
    def build(self):
        self.model = ConcreteModel(name=self.name)
        
        # Sets
        self.model.Products = Set(initialize=self.products)
        self.model.Periods = Set(initialize=self.periods)
        
        # Parameters
        self.model.demand = Param(self.model.Products, self.model.Periods,
                                 initialize=self.data.get('demand', {}))
        
        # Variables
        self.model.production = Var(self.model.Products, self.model.Periods,
                                   within=NonNegativeReals)
        
        # Objective
        def cost_rule(model):
            return sum(model.production[p,t] for p in model.Products 
                      for t in model.Periods)
        self.model.objective = Objective(rule=cost_rule, sense=minimize)
        
        # Constraints
        def demand_rule(model, p, t):
            return model.production[p,t] >= model.demand[p,t]
        self.model.demand_constraint = Constraint(self.model.Products,
                                                  self.model.Periods,
                                                  rule=demand_rule)
    
    def get_solution(self):
        if self.results is None:
            raise ValueError("Model not solved")
        
        solution = {}
        for p in self.model.Products:
            for t in self.model.Periods:
                solution[(p,t)] = self.model.production[p,t].value
        
        return solution

# === Builder Pattern ===

class ModelBuilder:
    """Builder pour construction progressive de modèles"""
    
    def __init__(self, name):
        self.model = ConcreteModel(name=name)
        self._sets = {}
        self._params = {}
        self._vars = {}
    
    def add_set(self, name, elements):
        """Ajouter un ensemble"""
        setattr(self.model, name, Set(initialize=elements))
        self._sets[name] = elements
        return self
    
    def add_param(self, name, index_sets=None, values=None, **kwargs):
        """Ajouter un paramètre"""
        if index_sets is None:
            param = Param(initialize=values, **kwargs)
        elif isinstance(index_sets, list):
            sets = [getattr(self.model, s) for s in index_sets]
            param = Param(*sets, initialize=values, **kwargs)
        else:
            param = Param(getattr(self.model, index_sets), 
                         initialize=values, **kwargs)
        
        setattr(self.model, name, param)
        self._params[name] = param
        return self
    
    def add_var(self, name, index_sets=None, **kwargs):
        """Ajouter une variable"""
        if index_sets is None:
            var = Var(**kwargs)
        elif isinstance(index_sets, list):
            sets = [getattr(self.model, s) for s in index_sets]
            var = Var(*sets, **kwargs)
        else:
            var = Var(getattr(self.model, index_sets), **kwargs)
        
        setattr(self.model, name, var)
        self._vars[name] = var
        return self
    
    def add_objective(self, name, expr, sense=minimize):
        """Ajouter un objectif"""
        obj = Objective(expr=expr, sense=sense)
        setattr(self.model, name, obj)
        return self
    
    def add_constraint(self, name, rule, index_sets=None):
        """Ajouter une contrainte"""
        if index_sets is None:
            constraint = Constraint(rule=rule)
        elif isinstance(index_sets, list):
            sets = [getattr(self.model, s) for s in index_sets]
            constraint = Constraint(*sets, rule=rule)
        else:
            constraint = Constraint(getattr(self.model, index_sets), rule=rule)
        
        setattr(self.model, name, constraint)
        return self
    
    def build(self):
        """Retourner le modèle construit"""
        return self.model

# Utilisation du Builder
builder = ModelBuilder("MyModel")
model = (builder
         .add_set('Products', ['A', 'B', 'C'])
         .add_param('cost', 'Products', {'A': 10, 'B': 15, 'C': 12})
         .add_var('x', 'Products', within=NonNegativeReals)
         .add_objective('minimize_cost', 
                       expr=sum(builder.model.cost[p] * builder.model.x[p] 
                               for p in builder.model.Products))
         .build())

# === Strategy Pattern pour solveurs ===

class SolverStrategy:
    """Interface pour stratégies de résolution"""
    
    def solve(self, model):
        raise NotImplementedError

class GLPKStrategy(SolverStrategy):
    def __init__(self, **options):
        self.options = options
    
    def solve(self, model):
        solver = SolverFactory('glpk')
        for key, value in self.options.items():
            solver.options[key] = value
        return solver.solve(model, tee=True)

class CPLEXStrategy(SolverStrategy):
    def __init__(self, **options):
        self.options = options
    
    def solve(self, model):
        solver = SolverFactory('cplex')
        for key, value in self.options.items():
            solver.options[key] = value
        return solver.solve(model, tee=True)

class GurobiStrategy(SolverStrategy):
    def __init__(self, **options):
        self.options = options
    
    def solve(self, model):
        solver = SolverFactory('gurobi')
        for key, value in self.options.items():
            solver.options[key] = value
        return solver.solve(model, tee=True)

class OptimizationSolver:
    """Context pour résolution avec stratégie"""
    
    def __init__(self, strategy: SolverStrategy):
        self.strategy = strategy
    
    def set_strategy(self, strategy: SolverStrategy):
        self.strategy = strategy
    
    def solve(self, model):
        return self.strategy.solve(model)

# Utilisation
solver = OptimizationSolver(GLPKStrategy(tmlim=300))
results = solver.solve(model)

# Changer de stratégie
solver.set_strategy(CPLEXStrategy(timelimit=600))
results = solver.solve(model)


[OK] DEBUGGING ET TROUBLESHOOTING

# === Problèmes courants et solutions ===

# PROBLÈME 1: Modèle infaisable
def debug_infeasibility(model):
    """Diagnostiquer infaisabilité"""
    from pyomo.util.infeasible import log_infeasible_constraints
    
    print("\n=== ANALYSE D'INFAISABILITÉ ===\n")
    
    # Logger contraintes infaisables
    log_infeasible_constraints(model, log_expression=True, log_variables=True)
    
    # Relaxer modèle et trouver IIS (Irreducible Inconsistent Subsystem)
    relaxed = model.clone()
    
    # Ajouter variables de slack
    relaxed.slack_pos = Var(relaxed.component_objects(Constraint), 
                           within=NonNegativeReals)
    relaxed.slack_neg = Var(relaxed.component_objects(Constraint), 
                           within=NonNegativeReals)
    
    # Modifier contraintes avec slack
    for c in relaxed.component_objects(Constraint):
        for index in c:
            constraint = c[index]
            if constraint.body is not None:
                # Ajouter slack à la contrainte
                pass  # Implémentation simplifiée
    
    # Objectif: minimiser slack total
    def slack_obj_rule(model):
        return sum(model.slack_pos[c] + model.slack_neg[c] 
                  for c in model.component_objects(Constraint))
    relaxed.slack_obj = Objective(rule=slack_obj_rule, sense=minimize)
    
    # Résoudre modèle relaxé
    solver = SolverFactory('glpk')
    results = solver.solve(relaxed)
    
    # Identifier contraintes problématiques
    print("\nContraintes avec slack non nul:")
    for c in relaxed.component_objects(Constraint):
        if hasattr(relaxed, f'slack_pos_{c.name}'):
            slack_val = relaxed.component(f'slack_pos_{c.name}').value
            if slack_val is not None and slack_val > 1e-6:
                print(f"  {c.name}: slack = {slack_val}")

# PROBLÈME 2: Solution non bornée
def check_unbounded(model):
    """Vérifier si problème non borné"""
    print("\n=== VÉRIFICATION DES BORNES ===\n")
    
    # Vérifier bornes des variables
    unbounded_vars = []
    for v in model.component_data_objects(Var):
        if v.lb is None or v.ub is None:
            unbounded_vars.append(v.name)
    
    if unbounded_vars:
        print("Variables sans bornes:")
        for var_name in unbounded_vars[:10]:  # Afficher 10 premiers
            print(f"  {var_name}")
    
    # Vérifier si objectif peut devenir infini
    obj = list(model.component_objects(Objective, active=True))[0]
    print(f"\nSens de l'objectif: {obj.sense}")
    
    if obj.sense == minimize:
        print("Vérifier que l'objectif ne peut pas tendre vers -∞")
    else:
        print("Vérifier que l'objectif ne peut pas tendre vers +∞")

# PROBLÈME 3: Temps de résolution trop long
def optimize_model_performance(model):
    """Suggestions pour améliorer performance"""
    print("\n=== ANALYSE DE PERFORMANCE ===\n")
    
    # Compter composants
    n_vars = sum(1 for _ in model.component_data_objects(Var))
    n_constraints = sum(1 for _ in model.component_data_objects(Constraint, active=True))
    
    binary_vars = sum(1 for v in model.component_data_objects(Var) 
                     if v.domain == Binary)
    integer_vars = sum(1 for v in model.component_data_objects(Var) 
                      if v.domain in (Integers, NonNegativeIntegers))
    
    print(f"Variables totales: {n_vars}")
    print(f"  - Binaires: {binary_vars}")
    print(f"  - Entières: {integer_vars}")
    print(f"  - Continues: {n_vars - binary_vars - integer_vars}")
    print(f"Contraintes: {n_constraints}")
    
    # Suggestions
    print("\nSUGGESTIONS:")
    
    if binary_vars > 1000:
        print("  [ATTENTION] Nombre élevé de variables binaires")
        print("    -> Considérer agrégation ou décomposition")
    
    if n_constraints > 10000:
        print("  [ATTENTION] Nombre élevé de contraintes")
        print("    -> Vérifier contraintes redondantes")
    
    # Vérifier contraintes denses vs sparses
    print("\n  -> Utiliser solveurs spécialisés (Gurobi, CPLEX)")
    print("  -> Définir valeurs initiales (warm start)")
    print("  -> Réduire gap d'optimalité si solution approchée suffisante")

# PROBLÈME 4: Erreurs numériques
def check_numerical_issues(model):
    """Vérifier problèmes numériques"""
    print("\n=== VÉRIFICATION NUMÉRIQUE ===\n")
    
    # Vérifier échelles des coefficients
    param_values = []
    for p in model.component_data_objects(Param):
        if p.value is not None:
            param_values.append(abs(p.value))
    
    if param_values:
        min_val = min(param_values)
        max_val = max(param_values)
        ratio = max_val / min_val if min_val > 0 else float('inf')
        
        print(f"Plage des paramètres:")
        print(f"  Min: {min_val:.2e}")
        print(f"  Max: {max_val:.2e}")
        print(f"  Ratio: {ratio:.2e}")
        
        if ratio > 1e6:
            print("\n  [ATTENTION] Problème mal conditionné")
            print("    -> Considérer mise à l'échelle des variables/contraintes")
    
    # Vérifier contraintes presque redondantes
    print("\nVérifier contraintes similaires ou redondantes")

# PROBLÈME 5: Variables à valeurs fractionnaires non désirées
def enforce_integrality(model, tolerance=1e-6):
    """Forcer intégralité des variables"""
    print("\n=== VÉRIFICATION INTÉGRALITÉ ===\n")
    
    fractional_vars = []
    for v in model.component_data_objects(Var):
        if v.value is not None and v.domain in (Binary, Integers, NonNegativeIntegers):
            if abs(v.value - round(v.value)) > tolerance:
                fractional_vars.append((v.name, v.value))
    
    if fractional_vars:
        print(f"Trouvé {len(fractional_vars)} variables avec valeurs fractionnaires:")
        for var_name, val in fractional_vars[:10]:
            print(f"  {var_name} = {val}")
        
        print("\nSOLUTIONS:")
        print("  1. Augmenter limite de temps du solveur")
        print("  2. Réduire gap d'optimalité")
        print("  3. Arrondir manuellement si proche de l'entier")
    else:
        print("[OK] Toutes les variables entières ont des valeurs entières")


[OK] TECHNIQUES DE MODÉLISATION AVANCÉES

# === Symétrie et Coupes ===

def add_symmetry_breaking_constraints(model):
    """Ajouter contraintes de brisure de symétrie"""
    
    # Exemple: contraintes d'ordre pour variables binaires identiques
    # Si x[1], x[2], ..., x[n] sont symétriques, imposer x[1] >= x[2] >= ... >= x[n]
    
    n = len(model.x)
    model.symmetry_breaking = ConstraintList()
    
    for i in range(1, n):
        model.symmetry_breaking.add(model.x[i] >= model.x[i+1])

# === Coupes valides ===

def add_valid_cuts(model):
    """Ajouter coupes valides pour renforcer formulation"""
    
    model.valid_cuts = ConstraintList()
    
    # Exemple: inégalité de couverture
    # Si sum(x[i]) >= k requis, on peut ajouter des coupes plus fortes
    
    # Coupe 1: Au moins k variables parmi n doivent être à 1
    model.valid_cuts.add(
        sum(model.x[i] for i in model.I) >= model.k
    )
    
    # Coupe 2: Contrainte de clique
    # Si variables sont mutuellement exclusives par groupes
    for group in model.ConflictGroups:
        model.valid_cuts.add(
            sum(model.x[i] for i in group) <= 1
        )

# === Reformulations ===

def reformulate_product_terms(model):
    """Reformuler produits de variables binaires"""
    
    # Produit x * y où x, y binaires
    # Remplacer par variable z = x * y avec:
    # z <= x
    # z <= y
    # z >= x + y - 1
    
    model.z = Var(within=Binary)
    
    model.product_reform1 = Constraint(expr=model.z <= model.x)
    model.product_reform2 = Constraint(expr=model.z <= model.y)
    model.product_reform3 = Constraint(expr=model.z >= model.x + model.y - 1)

def reformulate_max_min(model):
    """Reformuler fonctions max/min"""
    
    # z = max(x[1], x[2], ..., x[n])
    # Remplacer par:
    # z >= x[i] pour tout i
    # z <= max_val
    
    model.z = Var(within=Reals)
    
    model.max_constraints = ConstraintList()
    for i in model.I:
        model.max_constraints.add(model.z >= model.x[i])
    
    # Minimiser z dans l'objectif force z = max(x[i])


[OK] RESSOURCES ET RÉFÉRENCES

# === Documentation officielle ===
# Pyomo: https://pyomo.readthedocs.io/
# Exemples: https://github.com/Pyomo/pyomo/tree/main/examples
# Forum: https://groups.google.com/g/pyomo-forum

# === Livres recommandés ===
# - "Pyomo - Optimization Modeling in Python" par Hart, Watson, Woodruff
# - "Model Building in Mathematical Programming" par H.P. Williams
# - "Integer Programming" par Wolsey

# === Solveurs ===
# Open-source:
# - GLPK: https://www.gnu.org/software/glpk/
# - CBC: https://github.com/coin-or/Cbc
# - IPOPT: https://coin-or.github.io/Ipopt/
# 
# Commerciaux (académique gratuit):
# - Gurobi: https://www.gurobi.com/
# - CPLEX: https://www.ibm.com/products/ilog-cplex-optimization-studio
# - FICO Xpress: https://www.fico.com/en/products/fico-xpress-optimization

# === Tutoriels et cours ===
# - Pyomo Workshop: https://github.com/Pyomo/PyomoWorkshops
# - OR-Tools: https://developers.google.com/optimization
# - NEOS Server (résolution en ligne): https://neos-server.org/


[OK] EXEMPLES PAR INDUSTRIE

# === Finance - Portfolio avec risque ===

def portfolio_with_risk_constraints():
    """Portfolio avec contraintes de risque VaR"""
    model = ConcreteModel(name="Portfolio Risk")
    
    # Actifs
    model.Assets = Set(initialize=['Stock1', 'Stock2', 'Bond1', 'Bond2'])
    
    # Rendements attendus
    model.returns = Param(model.Assets, initialize={
        'Stock1': 0.12, 'Stock2': 0.15, 'Bond1': 0.05, 'Bond2': 0.06
    })
    
    # Variance et covariance (matrice simplifiée)
    model.variance = Param(model.Assets, initialize={
        'Stock1': 0.04, 'Stock2': 0.06, 'Bond1': 0.01, 'Bond2': 0.01
    })
    
    # Poids du portfolio
    model.w = Var(model.Assets, bounds=(0, 1))
    
    # Budget
    model.budget = Constraint(expr=sum(model.w[a] for a in model.Assets) == 1)
    
    # Rendement minimum
    model.min_return = Constraint(
        expr=sum(model.returns[a] * model.w[a] for a in model.Assets) >= 0.08
    )
    
    # VaR constraint (simplifié)
    model.risk_limit = Constraint(
        expr=sum(model.variance[a] * model.w[a]**2 for a in model.Assets) <= 0.03
    )
    
    # Objectif: maximiser rendement
    model.obj = Objective(
        expr=sum(model.returns[a] * model.w[a] for a in model.Assets),
        sense=maximize
    )
    
    return model

# === Santé - Planification de personnel médical ===

def hospital_staff_scheduling():
    """Planification optimale du personnel hospitalier"""
    model = ConcreteModel(name="Hospital Scheduling")
    
    # Personnel et compétences
    model.Staff = Set(initialize=['Nurse1', 'Nurse2', 'Nurse3', 'Doctor1', 'Doctor2'])
    model.Shifts = Set(initialize=['Morning', 'Evening', 'Night'])
    model.Days = RangeSet(1, 7)
    model.Skills = Set(initialize=['ICU', 'ER', 'Surgery'])
    
    # Compétences du personnel
    staff_skills = {
        ('Nurse1', 'ICU'): 1, ('Nurse1', 'ER'): 1,
        ('Nurse2', 'ER'): 1, ('Nurse2', 'Surgery'): 1,
        ('Nurse3', 'ICU'): 1,
        ('Doctor1', 'ICU'): 1, ('Doctor1', 'ER'): 1, ('Doctor1', 'Surgery'): 1,
        ('Doctor2', 'Surgery'): 1, ('Doctor2', 'ER'): 1,
    }
    model.has_skill = Param(model.Staff, model.Skills, 
                           initialize=staff_skills, default=0)
    
    # Besoin en compétences par shift
    skill_demand = {}
    for s in model.Shifts:
        for d in model.Days:
            skill_demand[('ICU', s, d)] = 2
            skill_demand[('ER', s, d)] = 2
            skill_demand[('Surgery', s, d)] = 1
    
    model.demand = Param(model.Skills, model.Shifts, model.Days,
                        initialize=skill_demand)
    
    # Variables: assignation
    model.assigned = Var(model.Staff, model.Shifts, model.Days, within=Binary)
    
    # Objectif: minimiser coût (heures sup, nuit, etc.)
    def cost_rule(model):
        regular_cost = sum(model.assigned[st,sh,d] 
                          for st in model.Staff 
                          for sh in model.Shifts 
                          for d in model.Days)
        night_premium = 1.5 * sum(model.assigned[st,'Night',d]
                                 for st in model.Staff for d in model.Days)
        return regular_cost + night_premium
    
    model.total_cost = Objective(rule=cost_rule, sense=minimize)
    
    # Contraintes: satisfaction de la demande
    def demand_rule(model, sk, sh, d):
        return sum(model.has_skill[st,sk] * model.assigned[st,sh,d]
                  for st in model.Staff) >= model.demand[sk,sh,d]
    model.meet_demand = Constraint(model.Skills, model.Shifts, model.Days,
                                   rule=demand_rule)
    
    # Un shift max par jour
    def one_shift_rule(model, st, d):
        return sum(model.assigned[st,sh,d] for sh in model.Shifts) <= 1
    model.one_shift = Constraint(model.Staff, model.Days, rule=one_shift_rule)
    
    # Max 5 jours par semaine
    def max_days_rule(model, st):
        return sum(model.assigned[st,sh,d] 
                  for sh in model.Shifts for d in model.Days) <= 5
    model.max_days = Constraint(model.Staff, rule=max_days_rule)
    
    return model


# === Télécommunications - Placement d'antennes ===

def cell_tower_placement():
    """Placement optimal d'antennes de télécom"""
    model = ConcreteModel(name="Cell Tower Placement")
    
    # Emplacements potentiels et zones à couvrir
    model.Locations = RangeSet(1, 20)  # Emplacements possibles
    model.Zones = RangeSet(1, 50)  # Zones à couvrir
    
    # Coût d'installation
    import random
    random.seed(42)
    model.cost = Param(model.Locations, 
                      initialize={i: random.randint(50000, 150000) 
                                for i in model.Locations})
    
    # Couverture: quelle zone est couverte par quel emplacement
    coverage = {}
    for loc in model.Locations:
        for zone in model.Zones:
            # Emplacement couvre zones proches (simplifié)
            if abs(loc * 2.5 - zone) <= 5:
                coverage[(loc, zone)] = 1
            else:
                coverage[(loc, zone)] = 0
    
    model.covers = Param(model.Locations, model.Zones, initialize=coverage)
    
    # Variable binaire: installer antenne ou non
        # Multiple objectives (stored, not active)
    model.obj1 = Expression(expr=model.x[1]**2 + model.x[2]**2 + model.x[3]**2)
    model.obj2 = Expression(expr=2*model.x[1] + 3*model.x[2] + model.x[3])
    model.obj3 = Expression(expr=(model.x[1] - 5)**2 + (model.x[2] - 3)**2)
    
    # Weighted sum objective (scalarization)
    model.w1 = Param(initialize=0.5, mutable=True)
    model.w2 = Param(initialize=0.3, mutable=True)
    model.w3 = Param(initialize=0.2, mutable=True)
    
    model.weighted_obj = Objective(
        expr=model.w1 * model.obj1 + model.w2 * model.obj2 + model.w3 * model.obj3,
        sense=minimize
    )
    
    # Constraints
    model.c1 = Constraint(expr=model.x[1] + model.x[2] + model.x[3] >= 5)
    model.c2 = Constraint(expr=model.x[1] + 2*model.x[2] <= 15)
    
    return model

# Solve for Pareto frontier
def solve_pareto_frontier(model, solver, n_points=10):
    pareto_points = []
    
    for i in range(n_points):
        # Update weights
        w1 = i / n_points
        w2 = (n_points - i) / n_points / 2
        w3 = (n_points - i) / n_points / 2
        
        model.w1 = w1
        model.w2 = w2
        model.w3 = w3
        
        # Solve
        results = solver.solve(model)
        
        # Store Pareto point
        pareto_points.append({
            'weights': (w1, w2, w3),
            'obj1': value(model.obj1),
            'obj2': value(model.obj2),
            'obj3': value(model.obj3),
            'x': [model.x[j].value for j in [1,2,3]]
        })
    
    return pareto_points


[OK] OPTIMISATION DE RÉSEAUX

# === Max Flow Problem ===

def max_flow():
    model = ConcreteModel(name="Max Flow")
    
    # Nodes
    model.Nodes = Set(initialize=['S', 'A', 'B', 'C', 'D', 'T'])
    
    # Arcs with capacity
    arcs_capacity = {
        ('S','A'): 10, ('S','B'): 5,
        ('A','C'): 15, ('A','D'): 5,
        ('B','D'): 10,
        ('C','T'): 10, ('D','T'): 10
    }
    model.Arcs = Set(initialize=arcs_capacity.keys())
    model.capacity = Param(model.Arcs, initialize=arcs_capacity)
    
    # Source and sink
    model.source = Param(initialize='S')
    model.sink = Param(initialize='T')
    
    # Flow variables
    model.flow = Var(model.Arcs, within=NonNegativeReals)
    
    # Objective: maximize flow into sink
    def flow_rule(model):
        return sum(model.flow[i,j] for (i,j) in model.Arcs 
                  if j == model.sink)
    model.max_flow = Objective(rule=flow_rule, sense=maximize)
    
    # Flow conservation (except source and sink)
    def conservation_rule(model, node):
        if node == model.source or node == model.sink:
            return Constraint.Skip
        inflow = sum(model.flow[i,j] for (i,j) in model.Arcs if j == node)
        outflow = sum(model.flow[i,j] for (i,j) in model.Arcs if i == node)
        return inflow == outflow
    model.conservation = Constraint(model.Nodes, rule=conservation_rule)
    
    # Capacity constraints
    def capacity_rule(model, i, j):
        return model.flow[i,j] <= model.capacity[i,j]
    model.capacity_constraint = Constraint(model.Arcs, rule=capacity_rule)
    
    return model

# === Shortest Path Problem ===

def shortest_path():
    model = ConcreteModel(name="Shortest Path")
    
    # Nodes
    model.Nodes = Set(initialize=['A', 'B', 'C', 'D', 'E'])
    
    # Arcs with distance
    arcs_distance = {
        ('A','B'): 4, ('A','C'): 2,
        ('B','C'): 1, ('B','D'): 5,
        ('C','D'): 8, ('C','E'): 10,
        ('D','E'): 2
    }
    model.Arcs = Set(initialize=arcs_distance.keys())
    model.distance = Param(model.Arcs, initialize=arcs_distance)
    
    # Binary: arc used
    model.x = Var(model.Arcs, within=Binary)
    
    # Source and destination
    model.source = 'A'
    model.destination = 'E'
    
    # Objective: minimize total distance
    def distance_rule(model):
        return sum(model.distance[i,j] * model.x[i,j] for (i,j) in model.Arcs)
    model.total_distance = Objective(rule=distance_rule, sense=minimize)
    
    # Flow conservation
    def flow_rule(model, node):
        if node == model.source:
            supply = 1
        elif node == model.destination:
            supply = -1
        else:
            supply = 0
        
        inflow = sum(model.x[i,j] for (i,j) in model.Arcs if j == node)
        outflow = sum(model.x[i,j] for (i,j) in model.Arcs if i == node)
        return inflow - outflow == -supply
    model.flow_conservation = Constraint(model.Nodes, rule=flow_rule)
    
    return model

# === Minimum Spanning Tree ===

def minimum_spanning_tree():
    model = ConcreteModel(name="MST")
    
    # Nodes
    n_nodes = 5
    model.Nodes = RangeSet(1, n_nodes)
    
    # Edges with weights
    edges_weight = {
        (1,2): 2, (1,3): 3, (1,4): 5,
        (2,3): 4, (2,5): 7,
        (3,4): 1, (3,5): 6,
        (4,5): 3
    }
    model.Edges = Set(initialize=edges_weight.keys())
    model.weight = Param(model.Edges, initialize=edges_weight)
    
    # Binary: edge in MST
    model.x = Var(model.Edges, within=Binary)
    
    # Objective: minimize total weight
    model.total_weight = Objective(
        expr=sum(model.weight[i,j] * model.x[i,j] for (i,j) in model.Edges),
        sense=minimize
    )
    
    # Exactly n-1 edges in spanning tree
    model.n_edges = Constraint(
        expr=sum(model.x[i,j] for (i,j) in model.Edges) == n_nodes - 1
    )
    
    # Connectivity constraint (can be enforced via subtour elimination)
    # Note: Full implementation needs more complex constraints
    
    return model


[OK] SCHEDULING & PLANNING

# === Job Shop Scheduling ===

def job_shop_scheduling():
    model = ConcreteModel(name="Job Shop")
    
    # Jobs and machines
    model.Jobs = Set(initialize=['J1', 'J2', 'J3'])
    model.Machines = Set(initialize=['M1', 'M2', 'M3'])
    
    # Processing time
    process_time = {
        ('J1','M1'): 3, ('J1','M2'): 2, ('J1','M3'): 2,
        ('J2','M1'): 2, ('J2','M2'): 1, ('J2','M3'): 4,
        ('J3','M1'): 4, ('J3','M2'): 3, ('J3','M3'): 1,
    }
    model.proc_time = Param(model.Jobs, model.Machines, 
                            initialize=process_time)
    
    # Start time variables
    model.start = Var(model.Jobs, model.Machines, within=NonNegativeReals)
    
    # Makespan variable
    model.makespan = Var(within=NonNegativeReals)
    
    # Big M
    M = 1000
    
    # Binary: precedence between jobs on same machine
    model.y = Var(model.Jobs, model.Jobs, model.Machines, within=Binary)
    
    # Objective: minimize makespan
    model.minimize_makespan = Objective(expr=model.makespan, sense=minimize)
    
    # Makespan constraints
    def makespan_rule(model, j, m):
        return model.makespan >= model.start[j,m] + model.proc_time[j,m]
    model.makespan_constraint = Constraint(model.Jobs, model.Machines, 
                                           rule=makespan_rule)
    
    # Disjunctive constraints (one job at a time on each machine)
    def disjunctive_rule(model, j1, j2, m):
        if j1 >= j2:
            return Constraint.Skip
        return [
            model.start[j1,m] + model.proc_time[j1,m] <= 
            model.start[j2,m] + M * (1 - model.y[j1,j2,m]),
            
            model.start[j2,m] + model.proc_time[j2,m] <= 
            model.start[j1,m] + M * model.y[j1,j2,m]
        ]
    model.disjunctive = Constraint(model.Jobs, model.Jobs, model.Machines, 
                                   rule=disjunctive_rule)
    
    return model

# === Nurse Scheduling ===

def nurse_scheduling():
    model = ConcreteModel(name="Nurse Scheduling")
    
    # Nurses and days
    model.Nurses = Set(initialize=['N1', 'N2', 'N3', 'N4', 'N5'])
    model.Days = RangeSet(1, 7)
    model.Shifts = Set(initialize=['Morning', 'Evening', 'Night'])
    
    # Demand for each shift
    demand_data = {
        ('Morning',1): 2, ('Morning',2): 2, ('Morning',3): 2,
        ('Morning',4): 2, ('Morning',5): 3, ('Morning',6): 3, ('Morning',7): 2,
        ('Evening',1): 2, ('Evening',2): 2, ('Evening',3): 2,
        ('Evening',4): 2, ('Evening',5): 3, ('Evening',6): 3, ('Evening',7): 2,
        ('Night',1): 1, ('Night',2): 1, ('Night',3): 1,
        ('Night',4): 1, ('Night',5): 2, ('Night',6): 2, ('Night',7): 1,
    }
    model.demand = Param(model.Shifts, model.Days, initialize=demand_data)
    
    # Binary: nurse n works shift s on day d
    model.x = Var(model.Nurses, model.Shifts, model.Days, within=Binary)
    
    # Objective: minimize total assignments (minimize cost)
    model.total_assignments = Objective(
        expr=sum(model.x[n,s,d] 
                for n in model.Nurses 
                for s in model.Shifts 
                for d in model.Days),
        sense=minimize
    )
    
    # Meet demand
    def demand_rule(model, s, d):
        return sum(model.x[n,s,d] for n in model.Nurses) >= model.demand[s,d]
    model.meet_demand = Constraint(model.Shifts, model.Days, rule=demand_rule)
    
    # Max one shift per day per nurse
    def one_shift_rule(model, n, d):
        return sum(model.x[n,s,d] for s in model.Shifts) <= 1
    model.one_shift = Constraint(model.Nurses, model.Days, rule=one_shift_rule)
    
    # Max 5 days per week per nurse
    def max_days_rule(model, n):
        return sum(model.x[n,s,d] 
                  for s in model.Shifts for d in model.Days) <= 5
    model.max_days = Constraint(model.Nurses, rule=max_days_rule)
    
    return model

# === Production Planning with Inventory ===

def production_planning_inventory():
    model = ConcreteModel(name="Production Planning")
    
    # Time periods
    model.T = RangeSet(1, 12)  # 12 months
    
    # Products
    model.Products = Set(initialize=['P1', 'P2'])
    
    # Demand forecast
    demand_data = {
        ('P1',1): 100, ('P1',2): 120, ('P1',3): 110, ('P1',4): 130,
        ('P1',5): 150, ('P1',6): 140, ('P1',7): 160, ('P1',8): 150,
        ('P1',9): 130, ('P1',10): 120, ('P1',11): 110, ('P1',12): 100,
        ('P2',1): 80, ('P2',2): 90, ('P2',3): 85, ('P2',4): 95,
        ('P2',5): 110, ('P2',6): 100, ('P2',7): 120, ('P2',8): 110,
        ('P2',9): 95, ('P2',10): 90, ('P2',11): 85, ('P2',12): 80,
    }
    model.demand = Param(model.Products, model.T, initialize=demand_data)
    
    # Production cost
    model.prod_cost = Param(model.Products, initialize={'P1': 50, 'P2': 60})
    
    # Inventory cost
    model.inv_cost = Param(model.Products, initialize={'P1': 5, 'P2': 6})
    
    # Production capacity
    model.capacity = Param(initialize=300)
    
    # Decision variables
    model.production = Var(model.Products, model.T, within=NonNegativeReals)
    model.inventory = Var(model.Products, model.T, within=NonNegativeReals)
    
    # Initial inventory
    model.initial_inv = Param(model.Products, initialize={'P1': 50, 'P2': 40})
    
    # Objective: minimize total cost
    def cost_rule(model):
        prod_cost = sum(model.prod_cost[p] * model.production[p,t]
                       for p in model.Products for t in model.T)
        inv_cost = sum(model.inv_cost[p] * model.inventory[p,t]
                      for p in model.Products for t in model.T)
        return prod_cost + inv_cost
    model.total_cost = Objective(rule=cost_rule, sense=minimize)
    
    # Inventory balance
    def balance_rule(model, p, t):
        if t == 1:
            return (model.initial_inv[p] + model.production[p,t] - 
                   model.demand[p,t] == model.inventory[p,t])
        else:
            return (model.inventory[p,t-1] + model.production[p,t] - 
                   model.demand[p,t] == model.inventory[p,t])
    model.balance = Constraint(model.Products, model.T, rule=balance_rule)
    
    # Capacity constraint
    def capacity_rule(model, t):
        return sum(model.production[p,t] for p in model.Products) <= model.capacity
    model.capacity_constraint = Constraint(model.T, rule=capacity_rule)
    
    return model


[OK] TECHNIQUES AVANCÉES

# === Column Generation ===

def cutting_stock_master():
    """Master problem for cutting stock (column generation)"""
    model = ConcreteModel(name="Cutting Stock Master")
    
    # Items to cut
    model.Items = Set(initialize=[1, 2, 3])
    model.demand = Param(model.Items, initialize={1: 100, 2: 80, 3: 60})
    
    # Patterns (columns) - initially basic patterns
    model.Patterns = Set(initialize=[1, 2, 3])
    
    # Number of each item in each pattern
    pattern_data = {
        (1,1): 1, (1,2): 0, (1,3): 0,
        (2,1): 0, (2,2): 1, (2,3): 0,
        (3,1): 0, (3,2): 0, (3,3): 1,
    }
    model.pattern_count = Param(model.Items, model.Patterns, 
                                initialize=pattern_data, mutable=True)
    
    # Variables: number of times each pattern is used
    model.x = Var(model.Patterns, within=NonNegativeIntegers)
    
    # Objective: minimize number of rolls
    model.minimize_rolls = Objective(
        expr=sum(model.x[p] for p in model.Patterns),
        sense=minimize
    )
    
    # Demand satisfaction
    def demand_rule(model, i):
        return sum(model.pattern_count[i,p] * model.x[p] 
                  for p in model.Patterns) >= model.demand[i]
    model.demand_constraint = Constraint(model.Items, rule=demand_rule)
    
    return model

# === Benders Decomposition ===

def benders_master():
    """Master problem for Benders decomposition"""
    model = ConcreteModel(name="Benders Master")
    
    # First-stage variables
    model.I = RangeSet(1, 3)
    model.x = Var(model.I, within=Binary)
    model.eta = Var()  # Approximation of second-stage cost
    
    # First-stage cost
    model.first_cost = Param(model.I, initialize={1: 100, 2: 150, 3: 120})
    
    # Objective
    model.obj = Objective(
        expr=sum(model.first_cost[i] * model.x[i] for i in model.I) + model.eta,
        sense=minimize
    )
    
    # Constraints
    model.budget = Constraint(expr=sum(model.x[i] for i in model.I) <= 2)
    
    # Optimality cuts (added iteratively)
    model.optimality_cuts = ConstraintList()
    
    return model

# === Branch and Price ===

# Combine branch-and-bound with column generation
# Master problem solved at each node with column generation

# === Lagrangian Relaxation ===

def lagrangian_relaxation_subproblem(multipliers):
    """Subproblem after Lagrangian relaxation"""
    model = ConcreteModel(name="Lagrangian Subproblem")
    
    model.I = RangeSet(1, 5)
    model.x = Var(model.I, within=Binary)
    
    # Original objective coefficients
    model.c = Param(model.I, initialize={1: 10, 2: 15, 3: 12, 4: 18, 5: 14})
    
    # Lagrangian multipliers (for relaxed constraints)
    model.lambda_ = Param(model.I, initialize=multipliers, mutable=True)
    
    # Modified objective (original + penalty for relaxed constraints)
    def lagrangian_obj(model):
        return (sum(model.c[i] * model.x[i] for i in model.I) +
                sum(model.lambda_[i] * model.x[i] for i in model.I))
    model.obj = Objective(rule=lagrangian_obj, sense=minimize)
    
    # Only easy constraints remain
    model.simple_constraint = Constraint(expr=sum(model.x[i] for i in model.I) >= 2)
    
    return model

# === Rolling Horizon ===

def rolling_horizon_planning(horizon=12, window=4):
    """Solve problem with rolling horizon"""
    all_periods = range(1, horizon + 1)
    solutions = {}
    
    for start_period in range(1, horizon - window + 2):
        # Define window
        current_window = range(start_period, min(start_period + window, horizon + 1))
        
        # Create model for current window
        model = ConcreteModel(name=f"Window_{start_period}")
        model.T = Set(initialize=current_window)
        model.x = Var(model.T, within=NonNegativeReals)
        
        # ... define rest of model for window ...
        
        # Solve
        solver = SolverFactory('glpk')
        results = solver.solve(model)
        
        # Implement first period decision only
        solutions[start_period] = model.x[start_period].value
    
    return solutions


[OK] TRANSFORMATIONS ET PRÉTRAITEMENT

# === Transformations de modèle ===

from pyomo.core.base import TransformationFactory

# Relaxation LP d'un MILP
def relax_integrality(model):
    """Relax integer variables to continuous"""
    relaxed = model.clone()
    TransformationFactory('core.relax_integrality').apply_to(relaxed)
    return relaxed

# Big-M transformation pour disjonctions
def apply_bigm(model):
    """Apply Big-M transformation to disjunctions"""
    TransformationFactory('gdp.bigm').apply_to(model)
    return model

# Hull reformulation pour disjonctions
def apply_hull(model):
    """Apply convex hull reformulation"""
    TransformationFactory('gdp.hull').apply_to(model)
    return model

# === Prétraitement ===

# Fixer des variables
model.x[1].fix(5)
model.x[1].fixed = True

# Libérer des variables fixées
model.x[1].free()
model.x[1].fixed = False

# Désactiver/activer contraintes
model.my_constraint.deactivate()
model.my_constraint.activate()

# Désactiver/activer objectifs
model.obj1.deactivate()
model.obj2.activate()

# === Scaling ===

# Mise à l'échelle des variables et contraintes
from pyomo.util.infeasible import log_infeasible_constraints

# Identifier contraintes infaisables
log_infeasible_constraints(model)

# Mise à l'échelle manuelle
model.scaling_factor = Suffix(direction=Suffix.EXPORT)
model.scaling_factor[model.x] = 0.001  # Variable scaling
model.scaling_factor[model.my_constraint] = 100  # Constraint scaling


[OK] DÉBOGAGE ET DIAGNOSTICS

# === Identifier problèmes ===

from pyomo.util.infeasible import log_infeasible_constraints, log_infeasible_bounds
from pyomo.common.timing import TicTocTimer

# Logger les contraintes infaisables
log_infeasible_constraints(model)

# Logger les bornes infaisables
log_infeasible_bounds(model)

# === Temporisation ===

timer = TicTocTimer()
timer.tic("Building model")
# ... build model ...
timer.toc("Building model")

timer.tic("Solving")
results = solver.solve(model)
timer.toc("Solving")

# === Vérification du modèle ===

# Vérifier structure
def check_model_structure(model):
    print(f"Variables: {len(list(model.component_data_objects(Var)))}")
    print(f"Constraints: {len(list(model.component_data_objects(Constraint, active=True)))}")
    print(f"Objectives: {len(list(model.component_data_objects(Objective, active=True)))}")
    
    # Variables par type
    binary_vars = [v for v in model.component_data_objects(Var) 
                   if v.domain == Binary]
    integer_vars = [v for v in model.component_data_objects(Var) 
                    if v.domain in (Integers, NonNegativeIntegers)]
    continuous_vars = [v for v in model.component_data_objects(Var) 
                       if v.domain == Reals]
    
    print(f"Binary: {len(binary_vars)}, Integer: {len(integer_vars)}, " 
          f"Continuous: {len(continuous_vars)}")

check_model_structure(model)

# === Écrire modèle en format lisible ===

# Écrire modèle complet
with open('model_details.txt', 'w') as f:
    model.pprint(ostream=f)

# Écrire statistiques
def write_model_stats(model, filename):
    with open(filename, 'w') as f:
        f.write("MODEL STATISTICS\n")
        f.write("=" * 50 + "\n\n")
        
        # Count components
        n_vars = sum(1 for _ in model.component_data_objects(Var))
        n_constraints = sum(1 for _ in model.component_data_objects(Constraint, active=True))
        n_objectives = sum(1 for _ in model.component_data_objects(Objective, active=True))
        
        f.write(f"Variables: {n_vars}\n")
        f.write(f"Constraints: {n_constraints}\n")
        f.write(f"Objectives: {n_objectives}\n\n")
        
        # Variable bounds
        f.write("VARIABLE BOUNDS\n")
        f.write("-" * 50 + "\n")
        for v in model.component_data_objects(Var):
            f.write(f"{v.name}: [{v.lb}, {v.ub}]\n")

write_model_stats(model, 'model_stats.txt')


[OK] BONNES PRATIQUES

# === Organisation du code ===

# Structure modulaire recommandée
class OptimizationModel:
    def __init__(self, data):
        self.data = data
        self.model = self.create_model()
        self.results = None
    
    def create_model(self):
        model = ConcreteModel(name="My Model")
        self._create_sets(model)
        self._create_parameters(model)
        self._create_variables(model)
        self._create_constraints(model)
        self._create_objective(model)
        return model
    
    def _create_sets(self, model):
        model.I = Set(initialize=self.data['indices'])
    
    def _create_parameters(self, model):
        model.cost = Param(model.I, initialize=self.data['costs'])
    
    def _create_variables(self, model):
        model.x = Var(model.I, within=NonNegativeReals)
    
    def _create_constraints(self, model):
        def my_rule(model, i):
            return model.x[i] <= self.data['capacity'][i]
        model.my_constraint = Constraint(model.I, rule=my_rule)
    
    def _create_objective(self, model):
        def obj_rule(model):
            return sum(model.cost[i] * model.x[i] for i in model.I)
        model.obj = Objective(rule=obj_rule, sense=minimize)
    
    def solve(self, solver_name='glpk', **options):
        solver = SolverFactory(solver_name)
        for key, value in options.items():
            solver.options[key] = value
        self.results = solver.solve(self.model, tee=True)
        return self.results
    
    def get_solution(self):
        if self.results is None:
            raise ValueError("Model not solved yet")
        return {i: self.model.x[i].value for i in self.model.I}
    
    def display_solution(self):
        print(f"\nObjective: {self.model.obj():.2f}")
        for i in self.model.I:
            print(f"x[{i}] = {self.model.x[i].value:.2f}")

# Utilisation
data = {
    'indices': [1, 2, 3],
    'costs': {1: 10, 2: 15, 3: 12},
    'capacity': {1: 100, 2: 150, 3: 120}
}

opt_model = OptimizationModel(data)
opt_model.solve(solver_name='glpk')
solution = opt_model.get_solution()
opt_model.display_solution()

# === Validation des données ===

def validate_input_data(data):
    """Validate input data before model creation"""
    assert 'indices' in data, "Missing indices"
    assert 'costs' in data, "Missing costs"
    assert len(data['indices']) > 0, "Empty indices"
    
    for i in data['indices']:
        assert i in data['costs'], f"Cost missing for index {i}"
        assert data['costs'][i] >= 0, f"Negative cost for index {i}"

# === Gestion d'erreurs ===

def solve_with_error_handling(model, solver_name='glpk'):
    try:
        solver = SolverFactory(solver_name)
        if not solver.available():
            raise ValueError(f"Solver {solver_name} not available")
        
        results = solver.solve(model, tee=True)
        
        if results.solver.status != SolverStatus.ok:
            raise RuntimeError(f"Solver status: {results.solver.status}")
        
        if results.solver.termination_condition != TerminationCondition.optimal:
            print(f"Warning: {results.solver.termination_condition}")
        
        return results
        
    except Exception as e:
        print(f"Error during solving: {e}")
        raise

# === Tests unitaires ===

import unittest

class TestOptimizationModel(unittest.TestCase):
    def setUp(self):
        self.data = {'indices': [1, 2], 'costs': {1: 10, 2: 15}}
        self.model_builder = OptimizationModel(self.data)
    
    def test_model_creation(self):
        model = self.model_builder.model
        self.assertIsNotNone(model)
        self.assertEqual(len(model.I), 2)
    
    def test_solution_feasibility(self):
        self.model_builder.solve()
        solution = self.model_builder.get_solution()
        self.assertIsNotNone(solution)
        for val in solution.values():
            self.assertGreaterEqual(val, 0)


[OK] PERFORMANCE ET OPTIMISATION

# === Conseils de performance ===

# 1. Utiliser quicksum au lieu de sum pour grandes sommes
from pyomo.environ import quicksum

# Lent pour grandes dimensions
expr = sum(model.x[i] for i in model.I)

# Plus rapide
expr = quicksum(model.x[i] for i in model.I)

# 2. Précalculer les expressions réutilisées
# Mauvais: recalcule à chaque fois
def constraint_rule(model, i):
    total = sum(model.x[j] for j in model.J)
    return model.y[i] <= total

# Bon: utiliser Expression
model.total_x = Expression(expr=sum(model.x[j] for j in model.J))

def constraint_rule(model, i):
    return model.y[i] <= model.total_x

# 3. Éviter les boucles imbriquées profondes
# Mauvais
for i in model.I:
    for j in model.J:
        for k in model.K:
            model.x[i,j,k] = Var(within=NonNegativeReals)

# Bon: utiliser indexation multiple
model.x = Var(model.I, model.J, model.K, within=NonNegativeReals)

# 4. Filtrer les combinaisons inutiles
# Mauvais: crée variables pour toutes combinaisons
model.x = Var(model.I, model.J, within=Binary)

# Bon: seulement combinaisons valides
model.ValidPairs = Set(within=model.I * model.J, 
                       initialize=[(i,j) for i in model.I for j in model.J if valid(i,j)])
model.x = Var(model.ValidPairs, within=Binary)

# 5. Utiliser solveurs persistants pour résolutions multiples
solver = SolverFactory('cplex_persistent')
solver.set_instance(model)

for scenario in scenarios:
    model.demand.value = scenario['demand']
    solver.solve(save_results=False)  # Plus rapide

# === Profilage ===

import cProfile
import pstats

# Profiler la création du modèle
profiler = cProfile.Profile()
profiler.enable()

model = create_large_model()

profiler.disable()
stats = pstats.Stats(profiler)
stats.sort_stats('cumulative')
stats.print_stats(20)  # Top 20 fonctions

# === Parallélisation ===

from pyomo.opt.parallel import SolverManagerFactory
from multiprocessing import Pool

def solve_scenario(scenario_data):
    model = create_model(scenario_data)
    solver = SolverFactory('glpk')
    results = solver.solve(model)
    return extract_solution(model)

# Résoudre scenarios en parallèle
scenarios = [scenario1, scenario2, scenario3, scenario4]

with Pool(processes=4) as pool:
    solutions = pool.map(solve_scenario, scenarios)


[OK] EXEMPLES COMPLETS PAR DOMAINE

# === Supply Chain Network Design ===

def supply_chain_network():
    model = ConcreteModel(name="Supply Chain Network")
    
    # Sets
    model.Suppliers = Set(initialize=['S1', 'S2', 'S3'])
    model.Plants = Set(initialize=['P1', 'P2'])
    model.Warehouses = Set(initialize=['W1', 'W2', 'W3'])
    model.Customers = Set(initialize=['C1', 'C2', 'C3', 'C4'])
    model.Products = Set(initialize=['Prod1', 'Prod2'])
    
    # Capacities
    model.supplier_capacity = Param(model.Suppliers, model.Products, initialize={
        ('S1','Prod1'): 1000, ('S1','Prod2'): 800,
        ('S2','Prod1'): 1200, ('S2','Prod2'): 1000,
        ('S3','Prod1'): 900, ('S3','Prod2'): 1100,
    })
    
    model.plant_capacity = Param(model.Plants, initialize={
        'P1': 2000, 'P2': 2500
    })
    
    model.warehouse_capacity = Param(model.Warehouses, initialize={
        'W1': 1500, 'W2': 1800, 'W3': 1600
    })
    
    # Customer demand
    model.demand = Param(model.Customers, model.Products, initialize={
        ('C1','Prod1'): 200, ('C1','Prod2'): 150,
        ('C2','Prod1'): 250, ('C2','Prod2'): 200,
        ('C3','Prod1'): 180, ('C3','Prod2'): 220,
        ('C4','Prod1'): 300, ('C4','Prod2'): 180,
    })
    
    # Costs
    model.transport_cost_sp = Param(model.Suppliers, model.Plants, model.Products,
                                    initialize=lambda m,s,p,pr: 2)
    model.transport_cost_pw = Param(model.Plants, model.Warehouses,
                                    initialize=lambda m,p,w: 3)
    model.transport_cost_wc = Param(model.Warehouses, model.Customers, model.Products,
                                    initialize=lambda m,w,c,pr: 1.5)
    
    model.production_cost = Param(model.Plants, model.Products, initialize={
        ('P1','Prod1'): 10, ('P1','Prod2'): 12,
        ('P2','Prod1'): 11, ('P2','Prod2'): 11,
    })
    
    model.warehouse_fixed_cost = Param(model.Warehouses, initialize={
        'W1': 5000, 'W2': 6000, 'W3': 5500
    })
    
    # Variables
    model.flow_sp = Var(model.Suppliers, model.Plants, model.Products, 
                        within=NonNegativeReals)
    model.flow_pw = Var(model.Plants, model.Warehouses, 
                        within=NonNegativeReals)
    model.flow_wc = Var(model.Warehouses, model.Customers, model.Products,
                        within=NonNegativeReals)
    
    model.production = Var(model.Plants, model.Products, 
                          within=NonNegativeReals)
    
    model.warehouse_open = Var(model.Warehouses, within=Binary)
    
    # Objective
    def total_cost_rule(model):
        transport_sp = sum(model.transport_cost_sp[s,p,pr] * model.flow_sp[s,p,pr]
                          for s in model.Suppliers for p in model.Plants 
                          for pr in model.Products)
        transport_pw = sum(model.transport_cost_pw[p,w] * model.flow_pw[p,w]
                          for p in model.Plants for w in model.Warehouses)
        transport_wc = sum(model.transport_cost_wc[w,c,pr] * model.flow_wc[w,c,pr]
                          for w in model.Warehouses for c in model.Customers 
                          for pr in model.Products)
        prod_cost = sum(model.production_cost[p,pr] * model.production[p,pr]
                       for p in model.Plants for pr in model.Products)
        fixed_cost = sum(model.warehouse_fixed_cost[w] * model.warehouse_open[w]
                        for w in model.Warehouses)
        return transport_sp + transport_pw + transport_wc + prod_cost + fixed_cost
    
    model.total_cost = Objective(rule=total_cost_rule, sense=minimize)
    
    # Constraints
    # Supplier capacity
    def supplier_cap_rule(model, s, pr):
        return sum(model.flow_sp[s,p,pr] for p in model.Plants) <= model.supplier_capacity[s,pr]
    model.supplier_cap = Constraint(model.Suppliers, model.Products, rule=supplier_cap_rule)
    
    # Plant production balance
    def plant_balance_rule(model, p, pr):
        inflow = sum(model.flow_sp[s,p,pr] for s in model.Suppliers)
        return inflow == model.production[p,pr]
    model.plant_balance = Constraint(model.Plants, model.Products, rule=plant_balance_rule)
    
    # Plant capacity
    def plant_cap_rule(model, p):
        return sum(model.production[p,pr] for pr in model.Products) <= model.plant_capacity[p]
    model.plant_cap = Constraint(model.Plants, rule=plant_cap_rule)
    
    # Plant to warehouse flow
    def plant_warehouse_flow_rule(model, p):
        return sum(model.flow_pw[p,w] for w in model.Warehouses) == sum(model.production[p,pr] for pr in model.Products)
    model.pw_flow = Constraint(model.Plants, rule=plant_warehouse_flow_rule)
    
    # Warehouse balance
    def warehouse_balance_rule(model, w):
        inflow = sum(model.flow_pw[p,w] for p in model.Plants)
        outflow = sum(model.flow_wc[w,c,pr] for c in model.Customers for pr in model.Products)
        return inflow == outflow
    model.warehouse_balance = Constraint(model.Warehouses, rule=warehouse_balance_rule)
    
    # Warehouse capacity
    def warehouse_cap_rule(model, w):
        total_flow = sum(model.flow_wc[w,c,pr] for c in model.Customers for pr in model.Products)
        return total_flow <= model.warehouse_capacity[w] * model.warehouse_open[w]
    model.warehouse_cap = Constraint(model.Warehouses, rule=warehouse_cap_rule)
    
    # Customer demand
    def demand_rule(model, c, pr):
        return sum(model.flow_wc[w,c,pr] for w in model.Warehouses) >= model.demand[c,pr]
    model.demand_constraint = Constraint(model.Customers, model.Products, rule=demand_rule)
    
    return model


# === Energy System Optimization ===

def energy_system_optimization():
    model = ConcreteModel(name="Energy System")
    
    # Time periods (24 hours)
    model.T = RangeSet(1, 24)
    
    # Generation units
    model.Units = Set(initialize=['Coal', 'Gas', 'Wind', 'Solar'])
    
    # Demand per hour
    demand_profile = {
        1: 1000, 2: 950, 3: 900, 4: 880, 5: 900, 6: 1000,
        7: 1200, 8: 1400, 9: 1500, 10: 1550, 11: 1600, 12: 1580,
        13: 1550, 14: 1520, 15: 1500, 16: 1520, 17: 1600, 18: 1800,
        19: 1900, 20: 1850, 21: 1700, 22: 1500, 23: 1300, 24: 1100
    }
    model.demand = Param(model.T, initialize=demand_profile)
    
    # Generation capacity
    model.capacity = Param(model.Units, initialize={
        'Coal': 800, 'Gas': 600, 'Wind': 400, 'Solar': 300
    })
    
    # Marginal cost ($/MWh)
    model.marginal_cost = Param(model.Units, initialize={
        'Coal': 30, 'Gas': 50, 'Wind': 0, 'Solar': 0
    })
    
    # Availability factor (renewable availability)
    availability_data = {}
    for t in range(1, 25):
        availability_data[('Wind', t)] = 0.3 + 0.2 * (t % 6) / 6
        if 7 <= t <= 18:
            availability_data[('Solar', t)] = 0.8 * (1 - abs(t - 12.5) / 6)
        else:
            availability_data[('Solar', t)] = 0
        availability_data[('Coal', t)] = 1.0
        availability_data[('Gas', t)] = 1.0
    
    model.availability = Param(model.Units, model.T, initialize=availability_data)
    
    # Ramp rate limits (MW/hour)
    model.ramp_rate = Param(model.Units, initialize={
        'Coal': 200, 'Gas': 300, 'Wind': 1000, 'Solar': 1000
    })
    
    # CO2 emissions (ton/MWh)
    model.emissions = Param(model.Units, initialize={
        'Coal': 1.0, 'Gas': 0.5, 'Wind': 0, 'Solar': 0
    })
    
    # Variables
    model.generation = Var(model.Units, model.T, within=NonNegativeReals)
    model.unit_on = Var(model.Units, model.T, within=Binary)
    
    # Objective: minimize generation cost
    def cost_rule(model):
        return sum(model.marginal_cost[u] * model.generation[u,t]
                  for u in model.Units for t in model.T)
    model.total_cost = Objective(rule=cost_rule, sense=minimize)
    
    # Demand satisfaction
    def demand_rule(model, t):
        return sum(model.generation[u,t] for u in model.Units) >= model.demand[t]
    model.meet_demand = Constraint(model.T, rule=demand_rule)
    
    # Generation capacity
    def capacity_rule(model, u, t):
        return model.generation[u,t] <= model.capacity[u] * model.availability[u,t] * model.unit_on[u,t]
    model.gen_capacity = Constraint(model.Units, model.T, rule=capacity_rule)
    
    # Ramp rate constraints
    def ramp_up_rule(model, u, t):
        if t == 1:
            return Constraint.Skip
        return model.generation[u,t] - model.generation[u,t-1] <= model.ramp_rate[u]
    model.ramp_up = Constraint(model.Units, model.T, rule=ramp_up_rule)
    
    def ramp_down_rule(model, u, t):
        if t == 1:
            return Constraint.Skip
        return model.generation[u,t-1] - model.generation[u,t] <= model.ramp_rate[u]
    model.ramp_down = Constraint(model.Units, model.T, rule=ramp_down_rule)
    
    # Emissions limit (optional)
    model.emissions_limit = Param(initialize=30000)
    
    def emissions_rule(model):
        return sum(model.emissions[u] * model.generation[u,t]
                  for u in model.Units for t in model.T) <= model.emissions_limit
    model.emissions_constraint = Constraint(rule=emissions_rule)
    
    return model


# === Vehicle Routing Problem with Time Windows ===

def vehicle_routing_time_windows():
    model = ConcreteModel(name="VRPTW")
    
    # Locations (0 = depot)
    model.Locations = RangeSet(0, 10)
    model.Customers = RangeSet(1, 10)
    
    # Vehicles
    model.Vehicles = RangeSet(1, 3)
    
    # Distance matrix
    import random
    random.seed(42)
    distance_data = {}
    for i in model.Locations:
        for j in model.Locations:
            if i != j:
                distance_data[(i,j)] = random.randint(5, 30)
            else:
                distance_data[(i,j)] = 0
    model.distance = Param(model.Locations, model.Locations, initialize=distance_data)
    
    # Demand at each customer
    demand_data = {i: random.randint(5, 20) for i in model.Customers}
    demand_data[0] = 0  # Depot has no demand
    model.demand = Param(model.Locations, initialize=demand_data)
    
    # Vehicle capacity
    model.vehicle_capacity = Param(initialize=50)
    
    # Time windows [earliest, latest]
    time_window_data = {}
    time_window_data[0] = (0, 480)  # Depot: 0-8 hours
    for i in model.Customers:
        early = random.randint(0, 180)
        time_window_data[i] = (early, early + random.randint(60, 120))
    model.time_early = Param(model.Locations, initialize={i: tw[0] for i, tw in time_window_data.items()})
    model.time_late = Param(model.Locations, initialize={i: tw[1] for i, tw in time_window_data.items()})
    
    # Service time at each location
    model.service_time = Param(model.Customers, initialize=lambda m,i: 15)
    
    # Variables
    model.x = Var(model.Vehicles, model.Locations, model.Locations, within=Binary)
    model.arrival_time = Var(model.Vehicles, model.Locations, within=NonNegativeReals)
    model.load = Var(model.Vehicles, model.Locations, within=NonNegativeReals)
    
    # Objective: minimize total distance
    def distance_rule(model):
        return sum(model.distance[i,j] * model.x[v,i,j]
                  for v in model.Vehicles 
                  for i in model.Locations 
                  for j in model.Locations if i != j)
    model.total_distance = Objective(rule=distance_rule, sense=minimize)
    
    # Each customer visited exactly once
    def visit_rule(model, i):
        return sum(model.x[v,i,j] 
                  for v in model.Vehicles 
                  for j in model.Locations if i != j) == 1
    model.visit_once = Constraint(model.Customers, rule=visit_rule)
    
    # Flow conservation
    def flow_rule(model, v, i):
        if i == 0:
            return Constraint.Skip
        inflow = sum(model.x[v,j,i] for j in model.Locations if j != i)
        outflow = sum(model.x[v,i,j] for j in model.Locations if i != j)
        return inflow == outflow
    model.flow_conservation = Constraint(model.Vehicles, model.Locations, rule=flow_rule)
    
    # Each vehicle starts from depot
    def start_depot_rule(model, v):
        return sum(model.x[v,0,j] for j in model.Customers) <= 1
    model.start_depot = Constraint(model.Vehicles, rule=start_depot_rule)
    
    # Each vehicle returns to depot
    def return_depot_rule(model, v):
        return sum(model.x[v,i,0] for i in model.Customers) <= 1
    model.return_depot = Constraint(model.Vehicles, rule=return_depot_rule)
    
    # Time window constraints
    M = 1000  # Big M
    
    def time_constraint_rule(model, v, i, j):
        if i == j:
            return Constraint.Skip
        service = model.service_time[i] if i in model.Customers else 0
        return (model.arrival_time[v,i] + service + model.distance[i,j] - 
                model.arrival_time[v,j] <= M * (1 - model.x[v,i,j]))
    model.time_constraint = Constraint(model.Vehicles, model.Locations, model.Locations, 
                                       rule=time_constraint_rule)
    
    # Time window bounds
    def time_window_early_rule(model, v, i):
        return model.arrival_time[v,i] >= model.time_early[i]
    model.tw_early = Constraint(model.Vehicles, model.Locations, rule=time_window_early_rule)
    
    def time_window_late_rule(model, v, i):
        return model.arrival_time[v,i] <= model.time_late[i]
    model.tw_late = Constraint(model.Vehicles, model.Locations, rule=time_window_late_rule)
    
    # Capacity constraints
    def capacity_rule(model, v, i, j):
        if i == 0 or j == 0 or i == j:
            return Constraint.Skip
        return (model.load[v,i] + model.demand[j] - model.load[v,j] <= 
                M * (1 - model.x[v,i,j]))
    model.capacity_constraint = Constraint(model.Vehicles, model.Locations, model.Locations,
                                           rule=capacity_rule)
    
    def max_capacity_rule(model, v, i):
        return model.load[v,i] <= model.vehicle_capacity
    model.max_capacity = Constraint(model.Vehicles, model.Locations, rule=max_capacity_rule)
    
    return model


[OK] INTÉGRATION AVEC AUTRES OUTILS

# === Intégration avec Pandas ===

import pandas as pd

def load_data_from_csv():
    # Charger données
    costs_df = pd.read_csv('costs.csv')
    demand_df = pd.read_csv('demand.csv')
    
    # Convertir en dictionnaires pour Pyomo
    cost_dict = {(row['source'], row['dest']): row['cost'] 
                 for _, row in costs_df.iterrows()}
    demand_dict = {row['location']: row['quantity'] 
                   for _, row in demand_df.iterrows()}
    
    return cost_dict, demand_dict

def export_solution_to_csv(model):
    # Extraire solution
    solution_data = []
    for v in model.component_data_objects(Var):
        if v.value is not None and abs(v.value) > 1e-6:
            solution_data.append({
                'Variable': v.name,
                'Value': v.value
            })
    
    # Créer DataFrame
    df = pd.DataFrame(solution_data)
    df.to_csv('solution.csv', index=False)
    
    return df

# === Intégration avec NumPy ===

import numpy as np

def create_model_from_numpy(cost_matrix, supply_vector, demand_vector):
    m, n = cost_matrix.shape
    
    model = ConcreteModel()
    model.I = RangeSet(0, m-1)
    model.J = RangeSet(0, n-1)
    
    # Convertir arrays NumPy en dictionnaires
    cost_dict = {(i,j): cost_matrix[i,j] 
                 for i in range(m) for j in range(n)}
    supply_dict = {i: supply_vector[i] for i in range(m)}
    demand_dict = {j: demand_vector[j] for j in range(n)}
    
    model.cost = Param(model.I, model.J, initialize=cost_dict)
    model.supply = Param(model.I, initialize=supply_dict)
    model.demand = Param(model.J, initialize=demand_dict)
    
    # ... reste du modèle ...
    
    return model

# === Intégration avec NetworkX ===

import networkx as nx

def create_model_from_graph(G):
    """Créer modèle Pyomo depuis graphe NetworkX"""
    model = ConcreteModel(name="Network Model")
    
    # Nœuds et arcs depuis le graphe
    model.Nodes = Set(initialize=G.nodes())
    model.Arcs = Set(initialize=G.edges())
    
    # Attributs du graphe comme paramètres
    cost_dict = nx.get_edge_attributes(G, 'weight')
    model.cost = Param(model.Arcs, initialize=cost_dict)
    
    return model

def export_solution_to_graph(model, G):
    """Exporter solution vers graphe NetworkX"""
    # Colorier arcs utilisés
    edge_colors = []
    edge_widths = []
    
    for (i,j) in G.edges():
        if hasattr(model, 'x') and (i,j) in model.x:
            if model.x[i,j].value is not None and model.x[i,j].value > 0.5:
                edge_colors.append('red')
                edge_widths.append(3)
            else:
                edge_colors.append('gray')
                edge_widths.append(1)
    
    return edge_colors, edge_widths

# === Intégration avec Matplotlib ===

import matplotlib.pyplot as plt

def plot_solution(model):
    """Visualiser solution"""
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # Graphique 1: Valeurs des variables
    var_names = []
    var_values = []
    for v in model.component_data_objects(Var):
        if v.value is not None and abs(v.value) > 1e-6:
            var_names.append(v.name)
            var_values.append(v.value)
    
    axes[0,0].barh(var_names[:10], var_values[:10])
    axes[0,0].set_xlabel('Value')
    axes[0,0].set_title('Top 10 Variables')
    
    # Graphique 2: Valeurs d'objectif par itération (si disponible)
    # ...
    
    plt.tight_layout()
    plt.savefig('solution_visualization.png', dpi=300)
    plt.show()

# === Intégration avec Plotly (interactif) ===

import plotly.graph_objects as go
import plotly.express as px

def create_interactive_visualization(model):
    """Créer visualisation interactive de la solution"""
    # Extraire données
    solution_data = []
    for v in model.component_data_objects(Var):
        if v.value is not None:
            solution_data.append({
                'Variable': v.name,
                'Value': v.value,
                'Type': str(v.domain)
            })
    
    df = pd.DataFrame(solution_data)
    
    # Créer figure interactive
    fig = px.bar(df.head(20), x='Variable', y='Value', color='Type',
                 title='Solution Variables')
    fig.update_layout(height=600)
    fig.write_html('solution_interactive.html')
    
    return fig


[OK] EXEMPLES DE WORKFLOWS COMPLETS

# === Workflow complet d'optimisation ===

def complete_optimization_workflow():
    """Workflow complet du chargement des données à l'export"""
    
    print("="*60)
    print("ÉTAPE 1: Chargement des données")
    print("="*60)
    
    # Charger depuis CSV
    costs_df = pd.read_csv('costs.csv')
    demand_df = pd.read_csv('demand.csv')
    
    # Validation des données
    assert not costs_df.isnull().any().any(), "Missing values in costs"
    assert (costs_df['cost'] >= 0).all(), "Negative costs found"
    
    print(f"Loaded {len(costs_df)} cost entries")
    print(f"Loaded {len(demand_df)} demand entries")
    
    print("\n" + "="*60)
    print("ÉTAPE 2: Construction du modèle")
    print("="*60)
    
    model = ConcreteModel(name="Production Planning")
    
    # Créer sets
    model.Products = Set(initialize=demand_df['product'].unique())
    model.Periods = Set(initialize=demand_df['period'].unique())
    
    # Créer paramètres
    demand_dict = {(row['product'], row['period']): row['demand']
                   for _, row in demand_df.iterrows()}
    model.demand = Param(model.Products, model.Periods, initialize=demand_dict)
    
    cost_dict = {row['product']: row['cost']
                 for _, row in costs_df.iterrows()}
    model.cost = Param(model.Products, initialize=cost_dict)
    
    # Variables
    model.production = Var(model.Products, model.Periods, within=NonNegativeReals)
    model.inventory = Var(model.Products, model.Periods, within=NonNegativeReals)
    
    # Objectif
    def cost_rule(model):
        return sum(model.cost[p] * model.production[p,t]
                  for p in model.Products for t in model.Periods)
    model.total_cost = Objective(rule=cost_rule, sense=minimize)
    
    # Contraintes
    def balance_rule(model, p, t):
        if t == model.Periods.first():
            return model.production[p,t] == model.demand[p,t] + model.inventory[p,t]
        else:
            prev_t = model.Periods.prev(t)
            return (model.inventory[p,prev_t] + model.production[p,t] ==
                   model.demand[p,t] + model.inventory[p,t])
    model.balance = Constraint(model.Products, model.Periods, rule=balance_rule)
    
    print(f"Modèle créé avec {len(list(model.component_data_objects(Var)))} variables")
    print(f"et {len(list(model.component_data_objects(Constraint, active=True)))} contraintes")
    
    print("\n" + "="*60)
    print("ÉTAPE 3: Résolution")
    print("="*60)
    
    solver = SolverFactory('glpk')
    solver.options['tmlim'] = 300
    
    results = solver.solve(model, tee=True)
    
    print(f"\nStatut: {results.solver.status# Fichier: python_cheats/cheatsheets/pyomo.txt
# Cheatsheet Pyomo - Guide Complet d'Optimisation Mathématique


[OK] INTRODUCTION & INSTALLATION

# Pyomo (Python Optimization Modeling Objects)
# Framework pour modélisation et résolution de problèmes d'optimisation

# Installation de base
pip install pyomo

# Installation complète avec extensions
pip install pyomo[optional]

# Installer solveurs open-source
conda install -c conda-forge glpk          # Solveur linéaire
conda install -c conda-forge ipopt         # Solveur non-linéaire
conda install -c conda-forge coincbc       # Solveur MILP

# Installation via Conda (recommandé pour data science)
conda install -c conda-forge pyomo
conda install -c conda-forge pyomo.extras

# Vérifier installation
pyomo help --solvers

# Importer Pyomo
from pyomo.environ import *
import pyomo.environ as pyo


[OK] TYPES DE MODÈLES

# Modèle Concret (ConcreteModel)
# Données définies lors de la création du modèle
model = ConcreteModel()
model = ConcreteModel(name="Mon Modèle")

# Modèle Abstrait (AbstractModel)
# Séparation modèle/données - données chargées après
model = AbstractModel()
model = AbstractModel(name="Modèle Générique")


[OK] COMPOSANTS DE BASE - SETS (ENSEMBLES)

# === Sets simples ===

# Set avec liste explicite
model.I = Set(initialize=[1, 2, 3, 4, 5])
model.Cities = Set(initialize=['Paris', 'Lyon', 'Marseille'])

# Set avec range
model.T = Set(initialize=range(1, 11))  # 1 à 10
model.Hours = RangeSet(1, 24)           # RangeSet optimisé

# Set vide (rempli plus tard)
model.Nodes = Set()

# Set ordonné
model.Days = Set(initialize=['Lun', 'Mar', 'Mer'], ordered=True)

# Set avec domaine
model.Positive = Set(within=PositiveIntegers)
model.Real = Set(within=Reals)

# === Sets multidimensionnels ===

# Set de tuples (produit cartésien)
model.Arcs = Set(initialize=[
    ('A', 'B'), ('B', 'C'), ('A', 'C')
])

# Set 2D avec within (sous-ensemble de produit cartésien)
model.I = Set(initialize=[1, 2, 3])
model.J = Set(initialize=[1, 2, 3])
model.Pairs = Set(within=model.I * model.J)

# Set avec dimension explicite
model.Coords = Set(dimen=2, initialize=[
    (1, 2), (2, 3), (3, 4)
])

# Set avec filtrage
model.ValidPairs = Set(initialize=model.I * model.J, 
                       filter=lambda model, i, j: i != j)

# === Sets avec règles ===

def init_arcs(model):
    return [(i, j) for i in model.I for j in model.J if i < j]

model.Arcs = Set(initialize=init_arcs, dimen=2)

# Set avec validation
def validate_set(model, value):
    return value > 0

model.PositiveSet = Set(initialize=[1, 2, 3], validate=validate_set)


[OK] PARAMÈTRES (PARAMETERS)

# === Paramètres scalaires ===

# Paramètre simple
model.Cost = Param(initialize=100)
model.Rate = Param(default=0.05)

# Paramètre mutable (modifiable après création)
model.Demand = Param(initialize=50, mutable=True)

# Paramètre avec domaine
model.Capacity = Param(within=PositiveReals, initialize=1000)

# === Paramètres indexés ===

# Paramètre indexé par un set
model.I = Set(initialize=[1, 2, 3])
model.Cost = Param(model.I, initialize={1: 10, 2: 20, 3: 15})

# Paramètre avec valeur par défaut
model.Supply = Param(model.I, default=0)

# Paramètre multi-indexé
model.Distance = Param(model.I, model.J, initialize={
    (1, 1): 0, (1, 2): 10, (1, 3): 20,
    (2, 1): 10, (2, 2): 0, (2, 3): 15,
    (3, 1): 20, (3, 2): 15, (3, 3): 0
})

# === Paramètres avec règles ===

def distance_rule(model, i, j):
    if i == j:
        return 0
    return abs(i - j) * 10

model.Dist = Param(model.I, model.J, initialize=distance_rule)

# === Paramètres calculés ===

def total_demand_rule(model):
    return sum(model.Demand[i] for i in model.I)

model.TotalDemand = Param(initialize=total_demand_rule)


[OK] VARIABLES (VARIABLES DE DÉCISION)

# === Variables continues ===

# Variable continue sans bornes
model.x = Var()

# Variable avec domaine
model.y = Var(within=Reals)
model.z = Var(domain=NonNegativeReals)

# Variable avec bornes
model.production = Var(bounds=(0, 100))
model.temp = Var(bounds=(-50, 150))

# Variable avec borne inférieure seulement
model.positive = Var(bounds=(0, None))

# === Variables entières ===

model.n_workers = Var(within=NonNegativeIntegers)
model.selection = Var(domain=Integers)
model.binary_choice = Var(within=Binary)  # 0 ou 1

# === Variables booléennes ===

model.is_active = Var(domain=Boolean)
model.flag = Var(within=Binary)

# === Variables indexées ===

# Variable indexée par un set
model.I = Set(initialize=[1, 2, 3, 4])
model.x = Var(model.I, within=NonNegativeReals)

# Variable multi-indexée
model.flow = Var(model.Arcs, within=NonNegativeReals, bounds=(0, 100))

# Variable avec bornes fonction d'index
def bounds_rule(model, i):
    return (0, model.Capacity[i])

model.y = Var(model.I, bounds=bounds_rule)

# === Variables avec initialisation ===

model.x = Var(initialize=10)  # Valeur initiale
model.y = Var(model.I, initialize={1: 5, 2: 10, 3: 15})

# === Variables avec règles complexes ===

def var_bounds_rule(model, i, j):
    if i == j:
        return (0, 0)  # Fixé à 0
    return (0, model.Capacity)

model.transport = Var(model.I, model.J, bounds=var_bounds_rule)


[OK] EXPRESSIONS

# === Expressions simples ===

# Expression linéaire
model.total_cost = Expression(expr=sum(model.Cost[i] * model.x[i] 
                                       for i in model.I))

# Expression avec règle
def profit_rule(model):
    return model.Revenue - model.Cost

model.profit = Expression(rule=profit_rule)

# === Expressions indexées ===

def margin_rule(model, i):
    return model.Price[i] - model.Cost[i]

model.margin = Expression(model.I, rule=margin_rule)

# === Opérations mathématiques ===

# Somme
total = sum(model.x[i] for i in model.I)
weighted_sum = sum(model.Weight[i] * model.x[i] for i in model.I)

# Produit
product = prod(model.x[i] for i in model.I)

# Maximum/Minimum (linéarisable)
from pyomo.core.expr import Expr_if
max_expr = max(model.x[i] for i in model.I)

# Valeur absolue (pour contraintes)
abs_value = abs(model.x[1] - model.x[2])

# Puissance
squared = model.x[1]**2
power = model.x[1]**model.n

# Fonctions mathématiques
from pyomo.environ import sin, cos, exp, log, sqrt

model.nonlinear = Expression(expr=exp(model.x) + log(model.y))
model.trig = Expression(expr=sin(model.theta) + cos(model.phi))


[OK] FONCTION OBJECTIF (OBJECTIVE)

# === Objectif simple ===

# Minimisation
model.obj = Objective(expr=sum(model.Cost[i] * model.x[i] 
                               for i in model.I), 
                      sense=minimize)

# Maximisation
model.obj = Objective(expr=sum(model.Revenue[i] * model.x[i] 
                               for i in model.I), 
                      sense=maximize)

# === Objectif avec règle ===

def objective_rule(model):
    production_cost = sum(model.ProdCost[i] * model.x[i] for i in model.I)
    transport_cost = sum(model.TransCost[i,j] * model.y[i,j] 
                        for i in model.I for j in model.J)
    return production_cost + transport_cost

model.total_cost = Objective(rule=objective_rule, sense=minimize)

# === Objectifs multiples (hiérarchie) ===

# Objectif prioritaire
model.primary_obj = Objective(expr=model.profit, sense=maximize)

# Objectif secondaire (désactivé initialement)
model.secondary_obj = Objective(expr=model.emissions, 
                                sense=minimize)
model.secondary_obj.deactivate()

# === Objectifs complexes ===

# Objectif quadratique
model.quad_obj = Objective(
    expr=sum(model.x[i]**2 for i in model.I) + 
         sum(model.x[i] * model.x[j] for i in model.I for j in model.J),
    sense=minimize
)

# Objectif non-linéaire
model.nonlin_obj = Objective(
    expr=sum(model.a[i] * exp(model.x[i]) for i in model.I),
    sense=minimize
)

# Objectif avec pénalité
model.penalized_obj = Objective(
    expr=model.profit - model.penalty_weight * model.violations,
    sense=maximize
)


[OK] CONTRAINTES (CONSTRAINTS)

# === Contraintes simples ===

# Contrainte d'égalité
model.balance = Constraint(expr=model.supply == model.demand)

# Contrainte d'inégalité
model.capacity_limit = Constraint(expr=model.production <= model.capacity)
model.min_production = Constraint(expr=model.production >= model.min_level)

# Contrainte composée
model.range_constraint = Constraint(
    expr=inequality(model.lower, model.x, model.upper)
)

# === Contraintes indexées ===

# Contrainte pour chaque élément d'un set
def demand_rule(model, i):
    return model.supply[i] >= model.demand[i]

model.demand_constraint = Constraint(model.I, rule=demand_rule)

# Contrainte multi-indexée
def flow_rule(model, i, j):
    return model.flow[i, j] <= model.capacity[i, j]

model.flow_limit = Constraint(model.Arcs, rule=flow_rule)

# === Contraintes de conservation ===

# Conservation de flux (nœud)
def flow_balance_rule(model, node):
    inflow = sum(model.flow[i, node] for i in model.I if (i, node) in model.Arcs)
    outflow = sum(model.flow[node, j] for j in model.J if (node, j) in model.Arcs)
    return inflow - outflow == model.demand[node]

model.flow_balance = Constraint(model.Nodes, rule=flow_balance_rule)

# === Contraintes conditionnelles ===

# Contrainte avec condition
def conditional_rule(model, i):
    if i in model.SpecialSet:
        return model.x[i] >= model.MinValue
    return Constraint.Skip  # Pas de contrainte

model.conditional_const = Constraint(model.I, rule=conditional_rule)

# Contrainte avec implication logique (Big-M)
M = 1000  # Big-M suffisamment grand
model.logical = Constraint(
    expr=model.x <= M * model.binary_var
)

# === Contraintes linéaires par morceaux ===

from pyomo.environ import Piecewise

model.pwl = Piecewise(
    model.y,           # Variable dépendante
    model.x,           # Variable indépendante
    pw_pts=[0, 1, 2, 3, 4],
    pw_constr_type='EQ',
    f_rule={0: 0, 1: 1, 2: 3, 3: 4, 4: 5}
)

# === Contraintes non-linéaires ===

# Contrainte quadratique
model.quad_const = Constraint(
    expr=sum(model.x[i]**2 for i in model.I) <= model.limit
)

# Contrainte avec fonctions transcendantes
model.nonlin_const = Constraint(
    expr=exp(model.x) + log(model.y) <= model.bound
)

# === Contraintes disjonctives ===

from pyomo.gdp import Disjunct, Disjunction

# Définir disjonctions
model.d1 = Disjunct()
model.d1.constraint = Constraint(expr=model.x >= 10)

model.d2 = Disjunct()
model.d2.constraint = Constraint(expr=model.x <= 5)

# Créer disjonction (l'une OU l'autre)
model.disjunction = Disjunction(expr=[model.d1, model.d2])

# === Contraintes avec indicateurs ===

# Contrainte activée par variable binaire
model.indicator_const = Constraint(
    expr=model.x >= 10,
    indicator=model.binary_var
)

# === Contraintes SOC (Second Order Cone) ===

from pyomo.environ import quicksum

# Contrainte de type ||x|| <= t
model.soc = Constraint(
    expr=quicksum(model.x[i]**2 for i in model.I) <= model.t**2
)


[OK] CONTRAINTES SPÉCIALISÉES

# === Contraintes SOS (Special Ordered Sets) ===

from pyomo.environ import SOSConstraint

# SOS Type 1: au plus une variable non-nulle
model.sos1 = SOSConstraint(var=model.x, sos=1)

# SOS Type 2: au plus deux variables consécutives non-nulles
model.sos2 = SOSConstraint(var=model.y, sos=2)

# === Contraintes de complémentarité ===

from pyomo.mpec import Complementarity

model.complementarity = Complementarity(
    expr=complements(model.x >= 0, model.y >= 0)
)

# === Contraintes de cardinalité ===

# Limiter nombre de variables non-nulles
def cardinality_rule(model):
    return sum(model.is_active[i] for i in model.I) <= model.max_active

model.cardinality = Constraint(rule=cardinality_rule)


[OK] CONSTRUCTION DE MODÈLES

# === Modèle simple complet ===

def create_simple_model():
    model = ConcreteModel(name="Simple LP")
    
    # Sets
    model.I = Set(initialize=[1, 2, 3])
    
    # Parameters
    model.c = Param(model.I, initialize={1: 3, 2: 4, 3: 5})
    model.b = Param(initialize=100)
    
    # Variables
    model.x = Var(model.I, within=NonNegativeReals)
    
    # Objective
    model.obj = Objective(
        expr=sum(model.c[i] * model.x[i] for i in model.I),
        sense=maximize
    )
    
    # Constraint
    model.resource = Constraint(
        expr=sum(model.x[i] for i in model.I) <= model.b
    )
    
    return model

# === Modèle avec fonctions de construction ===

def create_transportation_model(sources, destinations, supply, demand, cost):
    model = ConcreteModel(name="Transportation")
    
    # Sets
    model.I = Set(initialize=sources)
    model.J = Set(initialize=destinations)
    
    # Parameters
    model.supply = Param(model.I, initialize=supply)
    model.demand = Param(model.J, initialize=demand)
    model.cost = Param(model.I, model.J, initialize=cost)
    
    # Variables
    model.x = Var(model.I, model.J, within=NonNegativeReals)
    
    # Objective
    def obj_rule(model):
        return sum(model.cost[i,j] * model.x[i,j] 
                  for i in model.I for j in model.J)
    model.obj = Objective(rule=obj_rule, sense=minimize)
    
    # Supply constraints
    def supply_rule(model, i):
        return sum(model.x[i,j] for j in model.J) <= model.supply[i]
    model.supply_constraint = Constraint(model.I, rule=supply_rule)
    
    # Demand constraints
    def demand_rule(model, j):
        return sum(model.x[i,j] for i in model.I) >= model.demand[j]
    model.demand_constraint = Constraint(model.J, rule=demand_rule)
    
    return model


[OK] CHARGEMENT DE DONNÉES

# === Depuis dictionnaires Python ===

data = {
    'I': [1, 2, 3],
    'cost': {1: 10, 2: 20, 3: 15},
    'capacity': {1: 100, 2: 150, 3: 200}
}

model.I = Set(initialize=data['I'])
model.cost = Param(model.I, initialize=data['cost'])

# === Depuis fichiers DAT (format Pyomo) ===

# data.dat:
# set I := 1 2 3 ;
# param cost := 1 10  2 20  3 15 ;

model = AbstractModel()
# ... définir modèle ...
instance = model.create_instance('data.dat')

# === Depuis JSON ===

import json

with open('data.json', 'r') as f:
    data = json.load(f)

model.I = Set(initialize=data['sources'])
model.cost = Param(model.I, initialize=data['costs'])

# === Depuis CSV ===

import pandas as pd

# Charger données
df = pd.read_csv('costs.csv')

# Créer dictionnaire pour Pyomo
cost_dict = {(row['source'], row['dest']): row['cost'] 
             for _, row in df.iterrows()}

model.cost = Param(model.I, model.J, initialize=cost_dict)

# === Depuis Excel ===

df = pd.read_excel('data.xlsx', sheet_name='Costs')
cost_dict = df.set_index(['source', 'dest'])['cost'].to_dict()
model.cost = Param(model.I, model.J, initialize=cost_dict)

# === Depuis base de données ===

import sqlite3

conn = sqlite3.connect('database.db')
df = pd.read_sql_query("SELECT * FROM costs", conn)
cost_dict = {(row['i'], row['j']): row['cost'] 
             for _, row in df.iterrows()}

# === DataPortal (chargement avancé) ===

from pyomo.dataportal import DataPortal

data = DataPortal()
data.load(filename='data.dat')
data.load(filename='sets.dat', model=model)

instance = model.create_instance(data)


[OK] RÉSOLUTION DE MODÈLES

# === Solveurs disponibles ===

# Vérifier solveurs installés
from pyomo.opt import SolverFactory

available_solvers = []
for solver_name in ['glpk', 'cplex', 'gurobi', 'ipopt', 'bonmin']:
    solver = SolverFactory(solver_name)
    if solver.available():
        available_solvers.append(solver_name)

# === Résolution de base ===

# Créer solveur
solver = SolverFactory('glpk')

# Résoudre
results = solver.solve(model)

# Résoudre avec affichage
results = solver.solve(model, tee=True)

# === Options de solveur ===

# Options générales
solver = SolverFactory('glpk')
solver.options['mipgap'] = 0.01      # Gap d'optimalité 1%
solver.options['timelimit'] = 3600   # Limite de temps en secondes

# Options GLPK
solver.options['tmlim'] = 300
solver.options['mipgap'] = 0.05

# Options CPLEX
solver = SolverFactory('cplex')
solver.options['mip_tolerances_mipgap'] = 0.01
solver.options['timelimit'] = 1800
solver.options['threads'] = 4

# Options Gurobi
solver = SolverFactory('gurobi')
solver.options['MIPGap'] = 0.01
solver.options['TimeLimit'] = 3600
solver.options['Threads'] = 4
solver.options['Method'] = 2  # Barrier method

# Options IPOPT (non-linéaire)
solver = SolverFactory('ipopt')
solver.options['max_iter'] = 3000
solver.options['tol'] = 1e-6
solver.options['acceptable_tol'] = 1e-4

# === Résolution avec warm start ===

# Définir valeurs initiales
for i in model.I:
    model.x[i] = initial_values[i]

# Résoudre avec warm start
results = solver.solve(model, warmstart=True)

# === Résolution persistante (pour résolutions multiples) ===

from pyomo.solvers.plugins.solvers.persistent_solver import PersistentSolver

solver = SolverFactory('cplex_persistent')
solver.set_instance(model)

# Première résolution
solver.solve()

# Modifier paramètre
model.demand[1] = 150

# Résoudre à nouveau (plus rapide)
solver.solve()

# === Résolution avec callbacks ===

def callback(model, where):
    if where == 'mip_sol':
        print(f"Solution trouvée avec objectif: {model.obj()}")

results = solver.solve(model, callback=callback)

# === Résolution parallèle ===

from pyomo.opt import SolverManagerFactory

# Résolveur parallèle
solver_manager = SolverManagerFactory('serial')

# Soumettre plusieurs modèles
action_handles = []
for scenario in scenarios:
    model_instance = create_model(scenario)
    action_handle = solver_manager.queue(model_instance, 
                                         opt=solver, 
                                         tee=False)
    action_handles.append(action_handle)

# Récupérer résultats
for action_handle in action_handles:
    results = solver_manager.wait_for(action_handle)


[OK] ANALYSE DES RÉSULTATS

# === Status de résolution ===

from pyomo.opt import TerminationCondition, SolverStatus

# Vérifier si optimal
if results.solver.status == SolverStatus.ok:
    if results.solver.termination_condition == TerminationCondition.optimal:
        print("Solution optimale trouvée!")
    elif results.solver.termination_condition == TerminationCondition.infeasible:
        print("Problème infaisable")
    elif results.solver.termination_condition == TerminationCondition.unbounded:
        print("Problème non borné")

# === Récupérer valeurs des variables ===

# Valeur d'une variable
x_value = model.x[1].value
print(f"x[1] = {x_value}")

# Valeurs de toutes les variables indexées
for i in model.I:
    print(f"x[{i}] = {model.x[i].value}")

# Valeur de l'objectif
obj_value = model.obj()
print(f"Objectif = {obj_value}")

# === Récupérer informations supplémentaires ===

# Variables avec leurs valeurs
solution = {}
for v in model.component_data_objects(Var):
    if v.value is not None:
        solution[v.name] = v.value

# Variables non-nulles seulement
nonzero_vars = {v.name: v.value 
                for v in model.component_data_objects(Var) 
                if v.value is not None and abs(v.value) > 1e-6}

# === Analyser contraintes ===

# Slack des contraintes
for c in model.component_data_objects(Constraint, active=True):
    if c.body is not None:
        slack = c.upper - c.body()
        print(f"{c.name}: slack = {slack}")

# Contraintes saturées (actives)
active_constraints = []
for c in model.component_data_objects(Constraint, active=True):
    if c.body is not None and c.upper is not None:
        if abs(c.upper - c.body()) < 1e-6:
            active_constraints.append(c.name)

# === Informations de dualité (LP) ===

# Charger informations duales
model.dual = Suffix(direction=Suffix.IMPORT)
model.rc = Suffix(direction=Suffix.IMPORT)  # Reduced costs

results = solver.solve(model, load_solutions=False)
model.solutions.load_from(results)

# Variables duales (prix fictifs)
for c in model.component_data_objects(Constraint):
    if c in model.dual:
        print(f"{c.name}: dual = {model.dual[c]}")

# Coûts réduits
for v in model.component_data_objects(Var):
    if v in model.rc:
        print(f"{v.name}: reduced cost = {model.rc[v]}")

# === Analyse de sensibilité ===

# Sensibilité des paramètres
base_cost = model.cost[1].value
perturbations = [0.9, 0.95, 1.0, 1.05, 1.1]

sensitivity_results = {}
for factor in perturbations:
    model.cost[1] = base_cost * factor
    results = solver.solve(model)
    sensitivity_results[factor] = model.obj()

# === Statistiques du solveur ===

print(f"Temps de résolution: {results.solver.time:.2f}s")
print(f"Nombre d'itérations: {results.solver.statistics.branch_and_bound.number_of_created_subproblems}")
print(f"Gap d'optimalité: {results.solver.gap:.4f}")


[OK] AFFICHAGE ET EXPORT

# === Afficher modèle ===

# Affichage complet
model.pprint()

# Afficher composants spécifiques
model.x.pprint()
model.obj.pprint()
model.demand_constraint.pprint()

# Affichage compact
print(model.x.display())

# === Afficher solution ===

# Solution formatée
def display_solution(model):
    print("\n=== SOLUTION ===")
    print(f"Objectif: {model.obj():.2f}\n")
    
    print("Variables:")
    for v in model.component_data_objects(Var):
        if v.value is not None and abs(v.value) > 1e-6:
            print(f"  {v.name:20s} = {v.value:10.2f}")

display_solution(model)

# === Exporter vers fichiers ===

# Exporter vers LP format
from pyomo.opt import WriterFactory

writer = WriterFactory('lp')
writer(model, 'model.lp')

# Exporter vers MPS format
writer = WriterFactory('mps')
writer(model, 'model.mps')

# Exporter vers NL format (AMPL)
writer = WriterFactory('nl')
writer(model, 'model.nl')

# === Exporter solution vers CSV ===

import csv

with open('solution.csv', 'w', newline='') as f:
    writer = csv.writer(f)
    writer.writerow(['Variable', 'Value'])
    for v in model.component_data_objects(Var):
        if v.value is not None:
            writer.writerow([v.name, v.value])

# === Exporter vers DataFrame ===

import pandas as pd

# Variables
var_data = []
for v in model.component_data_objects(Var):
    if v.value is not None:
        var_data.append({'Variable': v.name, 'Value': v.value})

df_vars = pd.DataFrame(var_data)
df_vars.to_csv('variables.csv', index=False)
df_vars.to_excel('variables.xlsx', index=False)

# === Visualisation ===

import matplotlib.pyplot as plt

# Graphique des variables
values = [model.x[i].value for i in model.I]
plt.bar(model.I, values)
plt.xlabel('Index')
plt.ylabel('Value')
plt.title('Variable Values')
plt.savefig('solution.png')
plt.close()


[OK] PROGRAMMATION LINÉAIRE (LP)

# === Problème de production ===

def production_planning():
    model = ConcreteModel(name="Production Planning")
    
    # Products
    model.Products = Set(initialize=['A', 'B', 'C'])
    
    # Profit per unit
    model.profit = Param(model.Products, initialize={
        'A': 40, 'B': 30, 'C': 50
    })
    
    # Resource consumption
    model.labor_hours = Param(model.Products, initialize={
        'A': 2, 'B': 1, 'C': 3
    })
    model.material = Param(model.Products, initialize={
        'A': 3, 'B': 2, 'C': 4
    })
    
    # Available resources
    model.max_labor = Param(initialize=100)
    model.max_material = Param(initialize=150)
    
    # Decision variables: quantity to produce
    model.x = Var(model.Products, within=NonNegativeReals)
    
    # Objective: maximize profit
    def profit_rule(model):
        return sum(model.profit[p] * model.x[p] for p in model.Products)
    model.objective = Objective(rule=profit_rule, sense=maximize)
    
    # Constraints
    def labor_rule(model):
        return sum(model.labor_hours[p] * model.x[p] 
                  for p in model.Products) <= model.max_labor
    model.labor_constraint = Constraint(rule=labor_rule)
    
    def material_rule(model):
        return sum(model.material[p] * model.x[p] 
                  for p in model.Products) <= model.max_material
    model.material_constraint = Constraint(rule=material_rule)
    
    return model

# === Problème de transport ===

def transportation_problem():
    model = ConcreteModel(name="Transportation")
    
    # Sets
    model.Warehouses = Set(initialize=['W1', 'W2', 'W3'])
    model.Customers = Set(initialize=['C1', 'C2', 'C3', 'C4'])
    
    # Supply and demand
    model.supply = Param(model.Warehouses, initialize={
        'W1': 100, 'W2': 150, 'W3': 200
    })
    model.demand = Param(model.Customers, initialize={
        'C1': 80, 'C2': 90, 'C3': 120, 'C4': 150
    })
    
    # Transportation costs
    cost_data = {
        ('W1','C1'): 4, ('W1','C2'): 6, ('W1','C3'): 8, ('W1','C4'): 5,
        ('W2','C1'): 5, ('W2','C2'): 4, ('W2','C3'): 7, ('W2','C4'): 6,
        ('W3','C1'): 6, ('W3','C2'): 5, ('W3','C3'): 4, ('W3','C4'): 3,
    }
    model.cost = Param(model.Warehouses, model.Customers, 
                       initialize=cost_data)
    
    # Variables: quantity shipped
    model.x = Var(model.Warehouses, model.Customers, 
                  within=NonNegativeReals)
    
    # Objective: minimize cost
    def cost_rule(model):
        return sum(model.cost[w,c] * model.x[w,c] 
                  for w in model.Warehouses for c in model.Customers)
    model.total_cost = Objective(rule=cost_rule, sense=minimize)
    
    # Supply constraints
    def supply_rule(model, w):
        return sum(model.x[w,c] for c in model.Customers) <= model.supply[w]
    model.supply_constraint = Constraint(model.Warehouses, rule=supply_rule)
    
    # Demand constraints
    def demand_rule(model, c):
        return sum(model.x[w,c] for w in model.Warehouses) >= model.demand[c]
    model.demand_constraint = Constraint(model.Customers, rule=demand_rule)
    
    return model

# === Problème de mélange (Blending) ===

def blending_problem():
    model = ConcreteModel(name="Blending")
    
    # Ingredients
    model.Ingredients = Set(initialize=['I1', 'I2', 'I3'])
    
    # Cost per unit
    model.cost = Param(model.Ingredients, initialize={
        'I1': 10, 'I2': 15, 'I3': 12
    })
    
    # Nutritional content (protein %)
    model.protein = Param(model.Ingredients, initialize={
        'I1': 0.20, 'I2': 0.15, 'I3': 0.25
    })
    
    # Decision variables
    model.x = Var(model.Ingredients, within=NonNegativeReals)
    
    # Total quantity
    model.total_qty = Param(initialize=1000)
    
    # Objective: minimize cost
    def cost_rule(model):
        return sum(model.cost[i] * model.x[i] for i in model.Ingredients)
    model.objective = Objective(rule=cost_rule, sense=minimize)
    
    # Total quantity constraint
    def quantity_rule(model):
        return sum(model.x[i] for i in model.Ingredients) == model.total_qty
    model.quantity_constraint = Constraint(rule=quantity_rule)
    
    # Minimum protein requirement (18%)
    def protein_rule(model):
        return (sum(model.protein[i] * model.x[i] for i in model.Ingredients) 
                >= 0.18 * model.total_qty)
    model.protein_constraint = Constraint(rule=protein_rule)
    
    return model


[OK] PROGRAMMATION LINÉAIRE EN NOMBRES ENTIERS (MILP)

# === Problème du sac à dos (Knapsack) ===

def knapsack_problem():
    model = ConcreteModel(name="Knapsack")
    
    # Items
    model.Items = Set(initialize=range(1, 11))
    
    # Value and weight
    model.value = Param(model.Items, initialize={
        1: 10, 2: 15, 3: 8, 4: 20, 5: 12,
        6: 18, 7: 9, 8: 14, 9: 11, 10: 16
    })
    model.weight = Param(model.Items, initialize={
        1: 5, 2: 7, 3: 4, 4: 9, 5: 6,
        6: 8, 7: 3, 8: 7, 9: 5, 10: 8
    })
    
    # Capacity
    model.capacity = Param(initialize=40)
    
    # Binary decision variables
    model.x = Var(model.Items, within=Binary)
    
    # Objective: maximize value
    def value_rule(model):
        return sum(model.value[i] * model.x[i] for i in model.Items)
    model.objective = Objective(rule=value_rule, sense=maximize)
    
    # Capacity constraint
    def capacity_rule(model):
        return sum(model.weight[i] * model.x[i] 
                  for i in model.Items) <= model.capacity
    model.capacity_constraint = Constraint(rule=capacity_rule)
    
    return model

# === Problème d'affectation ===

def assignment_problem():
    model = ConcreteModel(name="Assignment")
    
    # Workers and tasks
    model.Workers = Set(initialize=['W1', 'W2', 'W3', 'W4'])
    model.Tasks = Set(initialize=['T1', 'T2', 'T3', 'T4'])
    
    # Cost matrix
    cost_data = {
        ('W1','T1'): 9, ('W1','T2'): 2, ('W1','T3'): 7, ('W1','T4'): 8,
        ('W2','T1'): 6, ('W2','T2'): 4, ('W2','T3'): 3, ('W2','T4'): 7,
        ('W3','T1'): 5, ('W3','T2'): 8, ('W3','T3'): 1, ('W3','T4'): 8,
        ('W4','T1'): 7, ('W4','T2'): 6, ('W4','T3'): 9, ('W4','T4'): 4,
    }
    model.cost = Param(model.Workers, model.Tasks, initialize=cost_data)
    
    # Binary assignment variables
    model.x = Var(model.Workers, model.Tasks, within=Binary)
    
    # Objective
    def cost_rule(model):
        return sum(model.cost[w,t] * model.x[w,t] 
                  for w in model.Workers for t in model.Tasks)
    model.total_cost = Objective(rule=cost_rule, sense=minimize)
    
    # Each worker assigned to exactly one task
    def worker_rule(model, w):
        return sum(model.x[w,t] for t in model.Tasks) == 1
    model.worker_constraint = Constraint(model.Workers, rule=worker_rule)
    
    # Each task assigned to exactly one worker
    def task_rule(model, t):
        return sum(model.x[w,t] for w in model.Workers) == 1
    model.task_constraint = Constraint(model.Tasks, rule=task_rule)
    
    return model

# === Facility Location Problem ===

def facility_location():
    model = ConcreteModel(name="Facility Location")
    
    # Sets
    model.Facilities = Set(initialize=['F1', 'F2', 'F3'])
    model.Customers = Set(initialize=['C1', 'C2', 'C3', 'C4', 'C5'])
    
    # Fixed cost to open facility
    model.fixed_cost = Param(model.Facilities, initialize={
        'F1': 1000, 'F2': 1200, 'F3': 900
    })
    
    # Transportation cost
    transport_cost = {
        ('F1','C1'): 4, ('F1','C2'): 6, ('F1','C3'): 9, ('F1','C4'): 8, ('F1','C5'): 7,
        ('F2','C1'): 5, ('F2','C2'): 4, ('F2','C3'): 7, ('F2','C4'): 5, ('F2','C5'): 6,
        ('F3','C1'): 6, ('F3','C2'): 7, ('F3','C3'): 4, ('F3','C4'): 5, ('F3','C5'): 3,
    }
    model.transport_cost = Param(model.Facilities, model.Customers, 
                                 initialize=transport_cost)
    
    # Capacity
    model.capacity = Param(model.Facilities, initialize={
        'F1': 100, 'F2': 150, 'F3': 120
    })
    
    # Demand
    model.demand = Param(model.Customers, initialize={
        'C1': 30, 'C2': 40, 'C3': 35, 'C4': 45, 'C5': 25
    })
    
    # Binary: is facility open?
    model.y = Var(model.Facilities, within=Binary)
    
    # Continuous: quantity shipped
    model.x = Var(model.Facilities, model.Customers, 
                  within=NonNegativeReals)
    
    # Objective
    def cost_rule(model):
        fixed = sum(model.fixed_cost[f] * model.y[f] 
                   for f in model.Facilities)
        transport = sum(model.transport_cost[f,c] * model.x[f,c] 
                       for f in model.Facilities for c in model.Customers)
        return fixed + transport
    model.total_cost = Objective(rule=cost_rule, sense=minimize)
    
    # Demand satisfaction
    def demand_rule(model, c):
        return sum(model.x[f,c] for f in model.Facilities) >= model.demand[c]
    model.demand_constraint = Constraint(model.Customers, rule=demand_rule)
    
    # Capacity constraint
    def capacity_rule(model, f):
        return sum(model.x[f,c] for c in model.Customers) <= model.capacity[f] * model.y[f]
    model.capacity_constraint = Constraint(model.Facilities, rule=capacity_rule)
    
    return model

# === Traveling Salesman Problem (TSP) ===

def tsp_problem():
    model = ConcreteModel(name="TSP")
    
    # Cities
    n_cities = 5
    model.Cities = RangeSet(1, n_cities)
    
    # Distance matrix
    distances = {
        (1,2): 10, (1,3): 15, (1,4): 20, (1,5): 25,
        (2,1): 10, (2,3): 35, (2,4): 25, (2,5): 30,
        (3,1): 15, (3,2): 35, (3,4): 30, (3,5): 20,
        (4,1): 20, (4,2): 25, (4,3): 30, (4,5): 16,
        (5,1): 25, (5,2): 30, (5,3): 20, (5,4): 16,
    }
    model.distance = Param(model.Cities, model.Cities, 
                           initialize=distances, default=0)
    
    # Binary: edge used in tour
    model.x = Var(model.Cities, model.Cities, within=Binary)
    
    # Subtour elimination variables (MTZ formulation)
    model.u = Var(model.Cities, within=NonNegativeIntegers, bounds=(1, n_cities))
    
    # Objective: minimize total distance
    def distance_rule(model):
        return sum(model.distance[i,j] * model.x[i,j] 
                  for i in model.Cities for j in model.Cities if i != j)
    model.total_distance = Objective(rule=distance_rule, sense=minimize)
    
    # Each city has exactly one incoming edge
    def in_degree_rule(model, j):
        return sum(model.x[i,j] for i in model.Cities if i != j) == 1
    model.in_degree = Constraint(model.Cities, rule=in_degree_rule)
    
    # Each city has exactly one outgoing edge
    def out_degree_rule(model, i):
        return sum(model.x[i,j] for j in model.Cities if i != j) == 1
    model.out_degree = Constraint(model.Cities, rule=out_degree_rule)
    
    # Subtour elimination (MTZ constraints)
    def subtour_rule(model, i, j):
        if i != j and (i != 1 and j != 1):
            return model.u[i] - model.u[j] + n_cities * model.x[i,j] <= n_cities - 1
        return Constraint.Skip
    model.subtour = Constraint(model.Cities, model.Cities, rule=subtour_rule)
    
    return model

# === Bin Packing Problem ===

def bin_packing():
    model = ConcreteModel(name="Bin Packing")
    
    # Items
    model.Items = RangeSet(1, 10)
    
    # Bins (maximum needed = number of items)
    model.Bins = RangeSet(1, 10)
    
    # Item sizes
    model.size = Param(model.Items, initialize={
        1: 4, 2: 7, 3: 3, 4: 6, 5: 5,
        6: 8, 7: 2, 8: 5, 9: 4, 10: 6
    })
    
    # Bin capacity
    model.capacity = Param(initialize=15)
    
    # Binary: item i in bin j
    model.x = Var(model.Items, model.Bins, within=Binary)
    
    # Binary: bin j is used
    model.y = Var(model.Bins, within=Binary)
    
    # Objective: minimize number of bins
    model.objective = Objective(expr=sum(model.y[j] for j in model.Bins), 
                                sense=minimize)
    
    # Each item in exactly one bin
    def assignment_rule(model, i):
        return sum(model.x[i,j] for j in model.Bins) == 1
    model.assignment = Constraint(model.Items, rule=assignment_rule)
    
    # Bin capacity constraint
    def capacity_rule(model, j):
        return sum(model.size[i] * model.x[i,j] 
                  for i in model.Items) <= model.capacity * model.y[j]
    model.bin_capacity = Constraint(model.Bins, rule=capacity_rule)
    
    return model


[OK] PROGRAMMATION NON-LINÉAIRE (NLP)

# === Portfolio Optimization (Markowitz) ===

def portfolio_optimization():
    model = ConcreteModel(name="Portfolio")
    
    # Assets
    model.Assets = Set(initialize=['A1', 'A2', 'A3', 'A4'])
    
    # Expected returns
    model.returns = Param(model.Assets, initialize={
        'A1': 0.10, 'A2': 0.15, 'A3': 0.12, 'A4': 0.08
    })
    
    # Covariance matrix (risk)
    cov_data = {
        ('A1','A1'): 0.01, ('A1','A2'): 0.005, ('A1','A3'): 0.003, ('A1','A4'): 0.002,
        ('A2','A1'): 0.005, ('A2','A2'): 0.02, ('A2','A3'): 0.008, ('A2','A4'): 0.004,
        ('A3','A1'): 0.003, ('A3','A2'): 0.008, ('A3','A3'): 0.015, ('A3','A4'): 0.003,
        ('A4','A1'): 0.002, ('A4','A2'): 0.004, ('A4','A3'): 0.003, ('A4','A4'): 0.005,
    }
    model.covariance = Param(model.Assets, model.Assets, initialize=cov_data)
    
    # Decision variables: portfolio weights
    model.w = Var(model.Assets, bounds=(0, 1))
    
    # Risk aversion parameter
    model.lambda_risk = Param(initialize=0.5)
    
    # Objective: maximize return - risk penalty
    def portfolio_rule(model):
        expected_return = sum(model.returns[i] * model.w[i] 
                             for i in model.Assets)
        portfolio_risk = sum(model.w[i] * model.covariance[i,j] * model.w[j]
                            for i in model.Assets for j in model.Assets)
        return expected_return - model.lambda_risk * portfolio_risk
    model.objective = Objective(rule=portfolio_rule, sense=maximize)
    
    # Budget constraint: weights sum to 1
    def budget_rule(model):
        return sum(model.w[i] for i in model.Assets) == 1
    model.budget = Constraint(rule=budget_rule)
    
    return model

# === Nonlinear Regression ===

def nonlinear_regression(x_data, y_data):
    model = ConcreteModel(name="Nonlinear Regression")
    
    # Data points
    model.N = RangeSet(0, len(x_data)-1)
    model.x_data = Param(model.N, initialize=dict(enumerate(x_data)))
    model.y_data = Param(model.N, initialize=dict(enumerate(y_data)))
    
    # Parameters to estimate: y = a * exp(b * x) + c
    model.a = Var(initialize=1.0)
    model.b = Var(initialize=0.1)
    model.c = Var(initialize=0.0)
    
    # Objective: minimize sum of squared errors
    def sse_rule(model):
        return sum((model.y_data[i] - 
                   (model.a * exp(model.b * model.x_data[i]) + model.c))**2
                  for i in model.N)
    model.sse = Objective(rule=sse_rule, sense=minimize)
    
    return model

# === Optimal Control Problem ===

def optimal_control():
    model = ConcreteModel(name="Optimal Control")
    
    # Time discretization
    model.T = RangeSet(0, 100)
    dt = 0.1
    
    # State variables
    model.x = Var(model.T, bounds=(-10, 10), initialize=0)
    model.v = Var(model.T, bounds=(-5, 5), initialize=0)
    
    # Control variable
    model.u = Var(model.T, bounds=(-1, 1))
    
    # Parameters
    model.mass = Param(initialize=1.0)
    model.drag = Param(initialize=0.1)
    
    # Objective: minimize control effort and final distance to target
    def cost_rule(model):
        control_cost = sum(model.u[t]**2 for t in model.T)
        final_cost = 10 * (model.x[100] - 10)**2  # Target position = 10
        return control_cost + final_cost
    model.cost = Objective(rule=cost_rule, sense=minimize)
    
    # Dynamics: velocity
    def velocity_dynamics(model, t):
        if t == 0:
            return Constraint.Skip
        return model.v[t] == model.v[t-1] + dt * (model.u[t-1] - model.drag * model.v[t-1]) / model.mass
    model.velocity_constraint = Constraint(model.T, rule=velocity_dynamics)
    
    # Dynamics: position
    def position_dynamics(model, t):
        if t == 0:
            return Constraint.Skip
        return model.x[t] == model.x[t-1] + dt * model.v[t-1]
    model.position_constraint = Constraint(model.T, rule=position_dynamics)
    
    # Initial conditions
    model.x[0].fix(0)
    model.v[0].fix(0)
    
    return model


[OK] PROGRAMMATION STOCHASTIQUE

# === Two-Stage Stochastic Programming ===

def two_stage_stochastic():
    model = ConcreteModel(name="Two-Stage Stochastic")
    
    # Scenarios
    model.Scenarios = Set(initialize=['S1', 'S2', 'S3'])
    model.prob = Param(model.Scenarios, initialize={
        'S1': 0.3, 'S2': 0.5, 'S3': 0.2
    })
    
    # Products
    model.Products = Set(initialize=['P1', 'P2'])
    
    # First-stage decision: production quantity
    model.x = Var(model.Products, within=NonNegativeReals)
    
    # First-stage cost
    model.prod_cost = Param(model.Products, initialize={'P1': 100, 'P2': 120})
    
    # Second-stage decision: recourse (shortage/excess)
    model.y_short = Var(model.Products, model.Scenarios, within=NonNegativeReals)
    model.y_excess = Var(model.Products, model.Scenarios, within=NonNegativeReals)
    
    # Scenario-dependent demand
    demand_data = {
        ('P1','S1'): 80, ('P1','S2'): 100, ('P1','S3'): 120,
        ('P2','S1'): 60, ('P2','S2'): 80, ('P2','S3'): 100,
    }
    model.demand = Param(model.Products, model.Scenarios, initialize=demand_data)
    
    # Recourse costs
    model.shortage_cost = Param(model.Products, initialize={'P1': 50, 'P2': 60})
    model.excess_cost = Param(model.Products, initialize={'P1': 10, 'P2': 15})
    
    # Objective: minimize expected total cost
    def cost_rule(model):
        first_stage = sum(model.prod_cost[p] * model.x[p] for p in model.Products)
        second_stage = sum(model.prob[s] * (
            sum(model.shortage_cost[p] * model.y_short[p,s] + 
                model.excess_cost[p] * model.y_excess[p,s]
                for p in model.Products)
        ) for s in model.Scenarios)
        return first_stage + second_stage
    model.total_cost = Objective(rule=cost_rule, sense=minimize)
    
    # Balance constraint for each scenario
    def balance_rule(model, p, s):
        return (model.x[p] + model.y_short[p,s] - model.y_excess[p,s] 
                == model.demand[p,s])
    model.balance = Constraint(model.Products, model.Scenarios, 
                               rule=balance_rule)
    
    return model


[OK] PROGRAMMATION PAR OBJECTIFS (GOAL PROGRAMMING)

def goal_programming():
    model = ConcreteModel(name="Goal Programming")
    
    # Decision variables
    model.x1 = Var(within=NonNegativeReals)
    model.x2 = Var(within=NonNegativeReals)
    
    # Deviation variables for goals
    model.d1_pos = Var(within=NonNegativeReals)  # Overachievement
    model.d1_neg = Var(within=NonNegativeReals)  # Underachievement
    model.d2_pos = Var(within=NonNegativeReals)
    model.d2_neg = Var(within=NonNegativeReals)
    
    # Goal targets
    model.target1 = Param(initialize=100)  # Profit goal
    model.target2 = Param(initialize=50)   # Production goal
    
    # Priority weights
    model.w1 = Param(initialize=2)  # Profit priority
    model.w2 = Param(initialize=1)  # Production priority
    
    # Objective: minimize weighted deviations
    model.obj = Objective(
        expr=model.w1 * (model.d1_pos + model.d1_neg) + 
             model.w2 * (model.d2_pos + model.d2_neg),
        sense=minimize
    )
    
    # Goal constraints
    model.goal1 = Constraint(
        expr=5*model.x1 + 3*model.x2 + model.d1_neg - model.d1_pos == model.target1
    )
    model.goal2 = Constraint(
        expr=model.x1 + model.x2 + model.d2_neg - model.d2_pos == model.target2
    )
    
    # Hard constraints
    model.resource1 = Constraint(expr=model.x1 + 2*model.x2 <= 80)
    model.resource2 = Constraint(expr=3*model.x1 + model.x2 <= 100)
    
    return model


[OK] OPTIMISATION ROBUSTE

def robust_optimization():
    model = ConcreteModel(name="Robust Optimization")
    
    # Products
    model.Products = Set(initialize=['P1', 'P2', 'P3'])
    
    # Nominal demand
    model.demand_nominal = Param(model.Products, initialize={
        'P1': 100, 'P2': 80, 'P3': 120
    })
    
    # Uncertainty deviation
    model.demand_dev = Param(model.Products, initialize={
        'P1': 20, 'P2': 15, 'P3': 25
    })
    
    # Budget of uncertainty (controls conservatism)
    model.Gamma = Param(initialize=1.5)
    
    # Decision variables
    model.x = Var(model.Products, within=NonNegativeReals)
    
    # Auxiliary variables for robust counterpart
    model.z = Var(model.Products, within=NonNegativeReals)
    model.xi = Var(within=NonNegativeReals)
    
    # Cost
    model.cost = Param(model.Products, initialize={
        'P1': 10, 'P2': 12, 'P3': 11
    })
    
    # Objective
    model.obj = Objective(
        expr=sum(model.cost[p] * model.x[p] for p in model.Products),
        sense=minimize
    )
    
    # Robust constraint (worst-case demand satisfaction)
    def robust_demand_rule(model, p):
        return model.x[p] >= model.demand_nominal[p] + model.z[p]
    model.robust_demand = Constraint(model.Products, rule=robust_demand_rule)
    
    # Uncertainty budget constraint
    def uncertainty_budget_rule(model):
        return sum(model.z[p] / model.demand_dev[p] 
                  for p in model.Products) <= model.Gamma
    model.uncertainty_budget = Constraint(rule=uncertainty_budget_rule)
    
    # Linking constraints
    def linking_rule(model, p):
        return model.z[p] <= model.demand_dev[p]
    model.linking = Constraint(model.Products, rule=linking_rule)
    
    return model


[OK] OPTIMISATION MULTI-OBJECTIFS (PARETO)

def multi_objective_optimization():
    model = ConcreteModel(name="Multi-Objective")
    
    # Decision variables
    model.x = Var(RangeSet(1, 3), within=NonNegativeReals, bounds=(0, 10))
    
    # Multiple objectives (stored, not active)