Corrigés Optimisation sous Contraintes

import numpy as np
import matplotlib.pyplot as plt

def rosenbrock(x):
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

def grad_rosenbrock(x):
    # Dérivées partielles
    df_dx = -2*(1 - x[0]) - 400*x[0]*(x[1] - x[0]**2)
    df_dy = 200*(x[1] - x[0]**2)
    return np.array([df_dx, df_dy])

def gradient_descent_rosenbrock(x0, alpha=0.001, tol=1e-6, max_iter=10000):
    x = x0.copy()
    trajectory = [x.copy()]
    
    for i in range(max_iter):
        grad = grad_rosenbrock(x)
        if np.linalg.norm(grad) < tol:
            break
        x = x - alpha * grad
        trajectory.append(x.copy())
        
        # Adaptation du pas (line search simplifié)
        if i > 0 and rosenbrock(x) > rosenbrock(trajectory[-2]):
            alpha *= 0.5  # Réduire le pas si divergence
    
    return x, np.array(trajectory), i

# Test
x0 = np.array([-1.0, 1.0])
x_opt, traj, n_iter = gradient_descent_rosenbrock(x0)

print(f"Solution: {x_opt}")  # Proche de [1, 1]
print(f"Itérations: {n_iter}")  # ~5000-8000 itérations

def newton_barrier(x0, mu=0.1, tol=1e-8):
    """
    Minimise -sum(log(1-x_i^2)) par méthode de Newton
    Domaine: ]-1,1[^n
    """
    x = x0.copy()
    n = len(x)
    
    for k in range(100):
        # Gradient: 2x/(1-x^2)
        grad = 2*x / (1 - x**2)
        
        # Hessienne: diag(2(1+x^2)/(1-x^2)^2)
        hess_diag = 2*(1 + x**2) / (1 - x**2)**2
        
        # Direction de Newton: d = -H^{-1}g
        d = -grad / hess_diag
        
        if np.linalg.norm(grad) < tol:
            break
        
        # Pas unitaire (quadratique)
        x = x + d
        
        # Projection si sortie du domaine (sécurité)
        x = np.clip(x, -0.999, 0.999)
    
    return x, k

# Test avec n=10
x0 = np.random.uniform(-0.5, 0.5, 10)
x_opt, iters = newton_barrier(x0)
print(f"Convergé en {iters} itérations vs ~500 pour gradient")

from scipy.optimize import linprog
import numpy as np

# Données
capacities = [100, 200, 150]
demands = [80, 120, 100, 150]
costs = np.array([
    [2, 3, 1, 4],
    [3, 2, 4, 2],
    [4, 1, 2, 3]
])

# Formulation: variable x_ij pour i=3 usines, j=4 magasins
# Ordonnées: x_00, x_01, x_02, x_03, x_10, x_11, ...
c = costs.flatten()

# Contraintes de capacité (lignes)
A_eq = []
b_eq = []

# Capacités usines: sum_j x_ij = capacity_i
for i in range(3):
    row = np.zeros(12)
    row[i*4:(i+1)*4] = 1
    A_eq.append(row)
    b_eq.append(capacities[i])

# Demandes magasins: sum_i x_ij = demand_j
for j in range(4):
    row = np.zeros(12)
    for i 3):
        row[i*4 + j] = 1
    A_eq.append(row)
    b_eq.append(demands[j])

result = linprog(c, A_eq=np.array(A_eq), b_eq=np.array(b_eq), 
                 bounds=(0, None), method='highs')

print(f"Coût optimal: {result.fun}")  # 940
print(f"Flux:\n{result.x.reshape(3,4)}")

def water_filling(alpha, P_total, tol=1e-9):
    """
    Algorithme water-filling pour allocation de puissance
    """
    n = len(alpha)
    
    # Tri pour efficacité (optionnel mais pratique)
    idx = np.argsort(alpha)[::-1]
    alpha_sorted = alpha[idx]
    
    # Recherche du niveau d'eau 1/lambda
    def total_power(level):
        return np.sum(np.maximum(0, level - 1/alpha_sorted))
    
    # Binary search sur level
    low, high = 0, P_total + np.max(1/alpha)
    while high - low > tol:
        mid = (low + high) / 2
        if total_power(mid) > P_total:
            high = mid
        else:
            low = mid
    
    level = (low + high) / 2
    p = np.maximum(0, level - 1/alpha)
    return p

# Exemple
alpha = np.array([1.0, 0.5, 2.0, 0.8])
P = 10
p_opt = water_filling(alpha, P)
print(f"Puissances optimales: {p_opt}")
print(f"Total utilisé: {np.sum(p_opt)}")

import numpy as np

class ACO:
    def __init__(self, n_ants, n_iterations, alpha=1, beta=2, rho=0.5):
        self.n_ants = n_ants
        self.n_iterations = n_iterations
        self.alpha = alpha  # poids phéromone
        self.beta = beta    # poids heuristique (1/distance)
        self.rho = rho      # évaporation
    
    def solve(self, dist_matrix):
        n = len(dist_matrix)
        pheromone = np.ones((n, n))
        best_path = None
        best_length = np.inf
        
        for iteration in range(self.n_iterations):
            all_paths = []
            all_lengths = []
            
            for ant in range(self.n_ants):
                path = self._construct_path(pheromone, dist_matrix)
                length = self._path_length(path, dist_matrix)
                all_paths.append(path)
                all_lengths.append(length)
                
                if length < best_length:
                    best_length = length
                    best_path = path
            
            # Mise à jour phéromones
            pheromone *= (1 - self.rho)
            for path, length in zip(all_paths, all_lengths):
                for i in range(len(path)-1):
                    pheromone[path[i], path[i+1]] += 1.0 / length
        
        return best_path, best_length
    
    def _construct_path(self, pheromone, dist):
        n = len(dist)
        start = np.random.randint(n)
        path = [start]
        visited = {start}
        
        while len(path) < n:
            current = path[-1]
            unvisited = [i for i in range(n) if i not in visited]
            
            # Probabilités
            probs = np.array([
                (pheromone[current][j]**self.alpha) * 
                ((1.0/dist[current][j])**self.beta)
                for j in unvisited
            ])
            probs = probs / np.sum(probs)
            
            next_city = np.random.choice(unvisited, p=probs)
            path.append(next_city)
            visited.add(next_city)
        
        path.append(path[0])  # Retour départ
        return path

import cvxpy as cp
import numpy as np

class SmartGridOptimizer:
    def __init__(self, T=24, B_max=100, R_max=20):
        self.T = T  # Horizon de planification (heures)
        self.B_max = B_max  # Capacité batterie MWh
        self.R_max = R_max  # Taux de charge max MW
        
    def optimize_deterministic(self, demand, renewable, costs, initial_soc=0.5):
        """
        Optimisation QP déterministe (moyenne des renouvelables)
        """
        # Variables
        p_thermal = cp.Variable(self.T, nonneg=True)
        p_hydro = cp.Variable(self.T)
        b_soc = cp.Variable(self.T, nonneg=True)  # State of charge
        s_market = cp.Variable(self.T)  # Positif=achat, négatif=vente
        
        # Paramètres
        d = np.array(demand)
        r = np.array(renewable)
        c_g, c_co2, c_m = costs['thermal'], costs['co2'], costs['market']
        
        # Objectif
        cost_thermal = c_g * cp.sum(p_thermal) + c_co2 * cp.sum(p_thermal)
        cost_market = cp.sum(c_m * s_market)
        cost_degrad = 0.01 * cp.sum(cp.abs(b_soc[1:] - b_soc[:-1]))  # Usure batterie
        
        objective = cp.Minimize(cost_thermal + cost_market + cost_degrad)
        
        # Contraintes
        constraints = []
        
        # Équilibre offre-demande
        for t in range(self.T):
            constraints.append(
                p_thermal[t] + p_hydro[t] + r[t] + s_market[t] == 
                d[t] + (b_soc[t] - (b_soc[t-1] if t > 0 else initial_soc*self.B_max))
            )
        
        # Limites batterie
        constraints += [b_soc <= 0.9*self.B_max, b_soc >= 0.2*self.B_max]
        constraints += [cp.abs(b_soc[1:] - b_soc[:-1]) <= self.R_max]
        
        # Limites thermique et hydro
        constraints += [p_thermal <= 80]  # MW max
        constraints += [cp.abs(p_hydro) <= 50]
        
        # Résolution
        prob = cp.Problem(objective, constraints)
        result = prob.solve(solver=cp.OSQP)
        
        return {
            'cost': result,
            'thermal': p_thermal.value,
            'hydro': p_hydro.value,
            'battery': b_soc.value,
            'market': s_market.value
        }
    
    def optimize_robust(self, demand, renewable_mean, renewable_std, costs, gamma=2.0):
        """
        Optimisation robuste (protection contre l'incertitude)
        Approche: budget d'incertitude (Bertsimas & Sim)
        """
        p_thermal = cp.Variable(self.T, nonneg=True)
        p_hydro = cp.Variable(self.T)
        b_soc = cp.Variable(self.T, nonneg=True)
        s_market = cp.Variable(self.T)
        
        # Variables d'ajustement pour le robuste
        z = cp.Variable(self.T, nonneg=True)  # Dépassement dû à l'incertitude
        
        d = np.array(demand)
        r_bar = np.array(renewable_mean)
        sigma = np.array(renewable_std)
        
        # Objectif avec coût de réserve
        cost = (costs['thermal'] + costs['co2']) * cp.sum(p_thermal) + \
               cp.sum(costs['market'] * s_market) + \
               0.1 * cp.sum(z)  # Pénalité robustesse
        
        constraints = []
        
        for t in range(self.T):
            # Contrainte robuste: l'incertitude ne doit pas causer de blackout
            # Worst case: renewable = r_bar - gamma*sigma (scénario défavorable)
            r_worst = r_bar[t] - gamma * sigma[t]
            
            constraints.append(
                p_thermal[t] + p_hydro[t] + r_worst + s_market[t] >= 
                d[t] + (b_soc[t] - (b_soc[t-1] if t > 0 else 50)) - z[t]
            )
        
        # ... autres contraintes identiques ...
        
        prob = cp.Problem(cp.Minimize(cost), constraints)
        prob.solve()
        
        return {'cost': prob.value, 'thermal': p_thermal.value, ...}

# Scénario de test
optimizer = SmartGridOptimizer(T=24)
demand = 50 + 30 * np.sin(np.linspace(0, 2*np.pi, 24))  # Charge journalière
solar = np.maximum(0, 40 * np.sin(np.linspace(0, np.pi, 24)))  # Production solaire

result = optimizer.optimize_deterministic(
    demand, solar, 
    {'thermal': 50, 'co2': 20, 'market': np.random.uniform(30, 100, 24)}
)

print(f"Coût total: {result['cost']:.2f} €")
print(f"Émissions CO2 évitées: {np.sum(result['thermal'])*0.5:.2f} tonnes")