Optimisation sous Contraintes

# Environnement D.I.A.M Optimization
pip install numpy scipy matplotlib cvxpy
pip install ipywidgets  # Pour visualisations interactives

# Pour l'optimisation symbolique (vérification des gradients)
pip install sympy

import numpy as np

def gradient_descent(f, grad_f, x0, alpha=0.1, tol=1e-6, max_iter=1000):
    """
    Minimisation de f par descente de gradient avec backtracking line search
    
    Args:
        f: fonction objectif (pour le suivi)
        grad_f: gradient de f (doit retourner un np.array)
        x0: point initial
        alpha: pas d'apprentissage initial
    """
    x = x0.copy()
    history = [x.copy()]
    f_history = [f(x)]
    
    for k in range(max_iter):
        gradient = grad_f(x)
        
        # Critère d'arrêt sur la norme du gradient
        if np.linalg.norm(gradient) < tol:
            print(f"Convergence atteinte en {k} itérations")
            break
        
        # Backtracking line search (condition d'Armijo)
        t = alpha
        while f(x - t * gradient) > f(x) - 0.5 * t * np.linalg.norm(gradient)**2:
            t *= 0.5
        
        # Mise à jour
        x = x - t * gradient
        history.append(x.copy())
        f_history.append(f(x))
    
    return x, np.array(history), f_history

# Exemple : f(x,y) = x^2 + 2y^2 (conditionnement 2)
f = lambda x: x[0]**2 + 2*x[1]**2
grad_f = lambda x: np.array([2*x[0], 4*x[1]])

x_opt, hist, f_hist = gradient_descent(f, grad_f, x0=np.array([5.0, 5.0]), alpha=0.1)
print(f"Minimum trouvé: {x_opt}")  # [0, 0]

import sympy as sp

# Résolution analytique: min x^2 + y^2 sous x + y = 1
x, y, lam = sp.symbols('x y lambda', real=True)

# Lagrangien
L = x**2 + y**2 + lam*(x + y - 1)

# Système: gradient L = 0
eq1 = sp.diff(L, x)  # 2x + lambda = 0
eq2 = sp.diff(L, y)  # 2y + lambda = 0
eq3 = sp.diff(L, lam) # x + y - 1 = 0

sol = sp.solve([eq1, eq2, eq3], [x, y, lam])
print(sol)  # x=0.5, y=0.5, lambda=-1

from scipy.optimize import linprog

# Problème: Diet Problem (minimiser coût nutritionnel)
# Variables: x1 = pain, x2 = viande, x3 = légumes

# Coûts à minimiser (€ par unité)
c = [2.0, 10.0, 3.0]

# Contraintes nutritionnelles (Ax >= b devient -Ax <= -b)
A_ub = [
    [-100, -200, -50],   # Calories >= 2000
    [-10, -50, -5]       # Protéines >= 100g
]
b_ub = [-2000, -100]

# Bornes (x >= 0)
bounds = [(0, None), (0, None), (0, None)]

result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

if result.success:
    print(f"Coût optimal: {result.fun:.2f}€")
    print(f"Quantités: {result.x}")
else:
    print("Pas de solution réalisable")

import cvxpy as cp
import numpy as np

# Problème: QP (Quadratic Programming)
# min 0.5*x^T*P*x + q^T*x sous Gx <= h

n = 2
x = cp.Variable(n)

# Données
P = np.array([[4, 1], [1, 2]])  # Définie positive
q = np.array([1, 1])
G = np.array([[-1, 0], [0, -1], [1, 1], [-1, 2]])
h = np.array([0, 0, 1, 1])

# Formulation
objective = cp.Minimize(0.5 * cp.quad_form(x, P) + q.T @ x)
constraints = [G @ x <= h]
prob = cp.Problem(objective, constraints)

result = prob.solve()
print(f"Valeur optimale: {result}")
print(f"Solution: {x.value}")
print(f"Multiplicateurs duaux: {[c.dual_value for c in constraints]}")

from scipy.optimize import minimize

def objective(x):
    return x[0]**2 + x[1]**2 + x[2]**2

def constraint1(x):
    return x[0] + x[1] + x[2] - 1  # = 0

def constraint2(x):
    return 0.5 - x[0]  # >= 0 donc x[0] <= 0.5

constraints = [
    {'type': 'eq', 'fun': constraint1},
    {'type': 'ineq', 'fun': constraint2}
]

x0 = np.array([0.5, 0.5, 0.0])

sol = minimize(objective, x0, method='SLSQP', 
               constraints=constraints,
               options={'ftol': 1e-9, 'disp': True})

print(f"Solution: {sol.x}")
print(f"Valeur f: {sol.fun}")

def simulated_annealing(f, x0, T_init=1000, cooling=0.95, n_iter=1000):
    """
    Minimisation par recuit simulé
    Accepte les mauvaises solutions avec proba exp(-deltaE/T)
    """
    x = x0.copy()
    x_best, f_best = x.copy(), f(x)
    T = T_init
    
    for i in range(n_iter):
        # Voisin aléatoire
        x_new = x + np.random.normal(0, 0.1, size=x.shape)
        f_new = f(x_new)
        delta = f_new - f(x)
        
        # Critère de Metropolis
        if delta < 0 or np.random.random() < np.exp(-delta / T):
            x = x_new
            if f_new < f_best:
                x_best, f_best = x_new, f_new
        
        T *= cooling
    
    return x_best, f_best