Lagrange Multipliers for Constraints¶

da PowerShell -- jupyter nbconvert --to html notebook.ipynb

Il metodo dei moltiplicatori di Lagrange è una tecnica di ottimizzazione utilizzata per trovare i massimi e i minimi di una funzione soggetta a uno o più vincoli di uguaglianza.

Problema di ottimizzazione vincolata¶

Sia $$ f : \mathbb{R}^n \to \mathbb{R} $$ una funzione da ottimizzare, soggetta al vincolo $$ g(x_1, x_2, \dots, x_n) = 0. $$

L'obiettivo è determinare i punti in cui $f$ assume valori estremi rispettando il vincolo imposto da $g$.

Idea geometrica¶

Geometricamente, i punti di estremo vincolato si trovano quando il gradiente di $f$ è parallelo al gradiente del vincolo $g$, cioè: $$ \nabla f = \lambda \nabla g, $$ dove $\lambda \in \mathbb{R}$ è detto $moltiplicatore di Lagrange$.

Funzione Lagrangiana¶

Si definisce la funzione Lagrangiana: $$ \mathcal{L}(x_1, \dots, x_n, \lambda) = f(x_1, \dots, x_n) + \lambda g(x_1, \dots, x_n). $$

I punti critici vincolati si ottengono risolvendo il sistema: $$ \begin{cases} \frac{\partial \mathcal{L}}{\partial x_1} = 0 \\ \frac{\partial \mathcal{L}}{\partial x_2} = 0 \\ \vdots \\ \frac{\partial \mathcal{L}}{\partial x_n} = 0 \\ \frac{\partial \mathcal{L}}{\partial \lambda} = g(x_1, \dots, x_n) = 0 \end{cases} $$

Esempio¶

Massimizzare la funzione $$ f(x,y) = xy $$ soggetta al vincolo $$ g(x,y) = x^2 + y^2 - 1 = 0. $$

La Lagrangiana è: $$ \mathcal{L}(x,y,\lambda) = xy + \lambda(x^2 + y^2 - 1). $$

Il sistema delle derivate parziali è: $$ \begin{cases} y + 2\lambda x = 0 \\ x + 2\lambda y = 0 \\ x^2 + y^2 = 1 \end{cases} $$

Risolvendo il sistema si ottengono i punti: $$ \left(\pm \frac{1}{\sqrt{2}}, \pm \frac{1}{\sqrt{2}}\right), $$ che corrispondono a massimi e minimi vincolati.

Estensione a più vincoli¶

Se sono presenti più vincoli $$ g_1(x) = 0, \quad g_2(x) = 0, \dots, g_m(x) = 0, $$ la condizione diventa: $$ \nabla f = \lambda_1 \nabla g_1 + \cdots + \lambda_m \nabla g_m, $$ con una Lagrangiana: $$ \mathcal{L}(x, \lambda_1, \dots, \lambda_m) = f(x) + \sum_{i=1}^{m} \lambda_i g_i(x). $$

Osservazioni¶

  • Il metodo fornisce condizioni necessarie, non sempre sufficienti.
  • È ampiamente utilizzato in fisica, economia, machine learning e ottimizzazione numerica.
  • Costituisce la base teorica dei problemi di ottimizzazione convessa e del metodo KKT.
  • Il Metodo KKT (Karush-Kuhn-Tucker) fornisce condizioni necessarie (e sufficienti nell'ottimizzazione convessa) per la soluzione di problemi vincolati, generalizzando i moltiplicatori di Lagrange e introducendo le condizioni di primalità, dualità e regolarità (qualificazione dei vincoli) per trovare i punti ottimali.
In [22]:
import numpy as np
from scipy.optimize import minimize 

def augmented_lagrangian(L, constraints, q, q_dot, t, lambda_vals):
    """
    Construct the augmented Lagrangian incorporating constraints
    using Lagrange multipliers.
        :param L: Original Lagrangian function, L(q, q_dot, t).
        :param constraints: List of constraint functions.
        :param q: Generalized coordinates.
        :param q_dot: Derivatives of generalized coordinates.
        : par am t: Time.
        :param lambda_vals: List of Lagrange multipliers.
    :return: 
        Augmented Lagrangian.
    """
    constraint_terms = sum(lmbda * f(*q, t) for lmbda, f in zip(lambda_vals, constraints))
    return L(q, q_dot, t) + constraint_terms 
    
# Example Lagrangian for a free particle
def lagrangian(q, q_dot, t):
    x, y, z = q
    x_dot, y_dot, z_dot = q_dot
    return 0.5 * (x_dot**2 + y_dot**2 + z_dot**2) 
    
# Constraint for a particle on a sphere of radius R
def sphere_constraint(x, y, z, t, R=1):
    return x**2 + y**2 + z**2 - R**2
    
# Solving the Euler-Lagrange equation with constraints
def euler_lagrange_with_constraints(L, constraints, q_initial, q_dot_initial, t_initial, lambda_initial):
    """
    Solve the Eulei—Lagrange equations incorporating constraints.
        :param L: Lagrangian function.
        :param constraints: Constraint functions.
        :param q_initial: Initial generalized coordinates.
        :param q_dot_initial: Initial velocities.
        :param t_initial: Initial time.
        :param lambda_initial: Initial estimates for Lagrange multip liers.
    :return: 
        Optimal coordinates and multipliers.
    """
    def objective(params):
        q, q_dot, lambdas = np.split(params, indices_or_sections=[3, 6])
        return augmented_lagrangian(L, constraints, q, q_dot, t_initial, lambdas) 
    
    # Initial guess combining coordinates, velocities, and multipliers
    initial_guess = np.concatenate((q_initial, q_dot_initial, lambda_initial))

    # Optimization to satisfy Euler-Lagrange with constraints
    result = minimize(objective, initial_guess)
    return result.x

# Initial conditions
q0 = np.array([1, 0, 0]) # Initial position on the sphere
q_dot0 = np.array([0, 1, 0]) # Initial velocity
lambda0 = np.array([1]) # Initial guess for lambda

Lagrangian = lagrangian(q0, q_dot0, 0)

# Solving for the particle constrained on a sphere
opt_params = euler_lagrange_with_constraints(
    Lagrangian,
    [sphere_constraint],
    q_initial=q0,
    q_dot_initial=q_dot0,
    t_initial=0,
    lambda_initial=lambda0
)

q_opt, q_dot_opt, lambda_opt = np.split(opt_params, indices_or_sections=[3, 6]) 
print("Optimal position:", q_opt)
print("Optimal velocity:", q_dot_opt)
print("Optimal Lagrange multipliers:", lambda_opt)

Final position: [0.40808206 0.91294525 0.        ]
Norm = 0.9999999999994466
In [28]:
import numpy as np
from scipy.optimize import minimize

def augmented_lagrangian(L, constraints, q, q_dot, t, lambda_vals):
    constraint_terms = sum(lmbda * f(*q, t) for lmbda, f in zip(lambda_vals, constraints))
    return L(q, q_dot, t) + constraint_terms

def lagrangian(q, q_dot, t):
    x, y, z = q
    x_dot, y_dot, z_dot = q_dot
    return 0.5 * (x_dot**2 + y_dot**2 + z_dot**2)

def sphere_constraint(x, y, z, t, R=1):
    return x**2 + y**2 + z**2 - R**2

def euler_lagrange_with_constraints(L, constraints, q_initial, q_dot_initial, t_initial, lambda_initial):

    def objective(params):
        q, q_dot, lambdas = np.split(params, indices_or_sections=[3, 6])
        return augmented_lagrangian(L, constraints, q, q_dot, t_initial, lambdas)

    initial_guess = np.concatenate((q_initial, q_dot_initial, lambda_initial))
    result = minimize(objective, initial_guess)
    return result.x

# -------------------------------------------------------------------
# ✔️ INPUT RAGIONEVOLE PER FARLO GIRARE
# -------------------------------------------------------------------

q0 = np.array([1.0, 0.0, 0.0])      # punto sulla sfera
q_dot0 = np.array([0.1, 0.2, 0.0])  # una piccola velocità tangenziale
lambda0 = np.array([0.0])          # valore iniziale del moltiplicatore

# Nota: qui si passa *lagrangian*, non il suo valore
opt_params = euler_lagrange_with_constraints(
    lagrangian,
    [sphere_constraint],
    q_initial=q0,
    q_dot_initial=q_dot0,
    t_initial=0.0,
    lambda_initial=lambda0
)

q_opt, q_dot_opt, lambda_opt = np.split(opt_params, indices_or_sections=[3, 6]) 
print("Optimal position:", q_opt)
print("Optimal velocity:", q_dot_opt)
print("Optimal Lagrange multipliers:", lambda_opt)

Optimal position: [1. 0. 0.]
Optimal velocity: [-7.31088221e-09 -7.40401446e-09 -7.45058060e-09]
Optimal Lagrange multipliers: [0.]
In [ ]: