Volver al blog
General

La Ecuación de Hamilton-Jacobi-Bellman en Control Óptimo

La Ecuación de Hamilton-Jacobi-Bellman (HJB) es una piedra angular en la teoría del control óptimo, proporcionando un marco fundamental para resolver problemas ...

30 de marzo de 202614 min read

La Ecuación de Hamilton-Jacobi-Bellman (HJB) es una piedra angular en la teoría del control óptimo, proporcionando un marco fundamental para resolver problemas de decisión secuencial tanto en entornos deterministas como estocásticos. Su relevancia se extiende más allá del control tradicional, encontrando aplicaciones profundas en campos emergentes como el aprendizaje por refuerzo continuo y, más recientemente, en la modelización generativa a través de los modelos de difusión. Este artículo explorará la HJB, su formulación, y su impacto unificador en estas disciplinas.

La Ecuación de Hamilton-Jacobi-Bellman en Control Óptimo

La HJB es una ecuación diferencial parcial (PDE) que describe la función de valor óptimo de un sistema dinámico. Derivada del Principio de Optimalidad de Bellman, establece que una política óptima tiene la propiedad de que, independientemente de los estados y decisiones anteriores, las decisiones restantes deben constituir una política óptima con respecto al estado resultante de esas decisiones anteriores.

Consideremos un problema de control óptimo estocástico en tiempo continuo. El objetivo es minimizar un costo esperado a lo largo de una trayectoria:

$$ J(x_0, t_0) = \mathbb{E}\left[ \int_{t_0}^T L(x_t, u_t, t) dt + \Phi(x_T) \right] $$

donde $x_t$ es el estado del sistema en el tiempo $t$, $u_t$ es la acción de control, $L(\cdot)$ es la función de costo instantáneo, $\Phi(\cdot)$ es el costo terminal, y el sistema evoluciona según una ecuación diferencial estocástica (SDE):

$$ dx_t = f(x_t, u_t, t) dt + g(x_t, t) dW_t $$

Aquí, $f(\cdot)$ es la función de deriva, $g(\cdot)$ es la función de difusión que escala el ruido Wiener $dW_t$. La función de valor óptimo, $V(x, t)$, se define como el costo mínimo esperado desde el estado $x$ en el tiempo $t$ hasta el tiempo final $T$:

$$ V(x, t) = \min_{u} \mathbb{E}\left[ \int_t^T L(x_s, u_s, s) ds + \Phi(x_T) \mid x_t = x \right] $$

Aplicando el Principio de Optimalidad de Bellman en un intervalo de tiempo infinitesimal $dt$, y utilizando el Lema de Itô para el diferencial de la función de valor $dV(x,t)$, se obtiene la Ecuación de Hamilton-Jacobi-Bellman:

$$ -\frac{\partial V}{\partial t} = \min_u \left{ L(x, u, t) + \nabla_x V \cdot f(x, u, t) + \frac{1}{2} \text{Tr}\left( g(x, t) g(x, t)^T \nabla_x^2 V \right) \right} $$

con la condición de contorno $V(x, T) = \Phi(x)$.

En esta ecuación:

  • $\frac{\partial V}{\partial t}$ representa la tasa de cambio del valor óptimo con el tiempo.
  • $L(x, u, t)$ es el costo instantáneo de la acción $u$ en el estado $x$ en el tiempo $t$.
  • $\nabla_x V \cdot f(x, u, t)$ es la contribución del cambio determinista en el estado a la función de valor, ponderada por el gradiente de $V$.
  • $\frac{1}{2} \text{Tr}\left( g(x, t) g(x, t)^T \nabla_x^2 V \right)$ es el término de difusión, que captura el impacto de la estocasticidad del sistema en la función de valor, involucrando el Hessiano $\nabla_x^2 V$ (matriz de segundas derivadas parciales).

La solución de la HJB proporciona la función de valor óptimo $V(x, t)$, a partir de la cual se puede derivar la política de control óptima $\pi^*(x, t)$ obteniendo el argumento que minimiza el término entre llaves:

$$ \pi^*(x, t) = \arg\min_u \left{ L(x, u, t) + \nabla_x V \cdot f(x, u, t) + \frac{1}{2} \text{Tr}\left( g(x, t) g(x, t)^T \nabla_x^2 V \right) \right} $$

La HJB es una PDE no lineal que rara vez tiene soluciones analíticas cerradas, excepto para casos muy específicos (por ejemplo, sistemas lineales con costos cuadráticos, donde se reduce a la ecuación de Riccati). Esto impulsa la necesidad de métodos numéricos o de aproximación.

HJB en Aprendizaje por Refuerzo Continuo

El Aprendizaje por Refuerzo (RL) se ocupa de cómo un agente debe tomar decisiones en un entorno para maximizar una noción acumulativa de recompensa. En el contexto de espacios de estados y acciones continuos, RL comparte una conexión profunda con el control óptimo y, por extensión, con la HJB.

En RL, el objetivo es maximizar la recompensa total esperada, lo que es equivalente a minimizar un costo total esperado (si las recompensas son negativas de costos). La función de valor $V(s)$ y la función de acción-valor $Q(s, a)$ son centrales, representando la recompensa esperada desde un estado $s$ o después de tomar una acción $a$ en el estado $s$, respectivamente.

Para un problema de RL con estados y acciones continuos, podemos reformular la ecuación de Bellman en tiempo continuo, que es la contraparte discreta de la HJB. Si consideramos un modelo estocástico con dinámica $dx_t = f(x_t, u_t) dt + g(x_t) dW_t$ y una función de recompensa $R(x_t, u_t)$, la función de valor óptimo $V(x)$ satisface:

$$ \rho V(x) = \max_u \left{ R(x, u) + \nabla_x V \cdot f(x, u) + \frac{1}{2} \text{Tr}\left( g(x) g(x)^T \nabla_x^2 V \right) \right} $$

donde $\rho$ es una tasa de descuento o un término que define el objetivo del valor promedio a largo plazo.

La principal dificultad en la resolución directa de la HJB en RL radica en la "maldición de la dimensionalidad" y en la necesidad de conocer explícitamente las dinámicas del sistema ($f, g$) y la función de recompensa $R$. Los algoritmos de RL a menudo operan en un régimen "model-free", donde estas dinámicas son desconocidas y deben aprenderse a través de la interacción.

Los enfoques de Deep Reinforcement Learning (DRL) para el control continuo, como Deep Deterministic Policy Gradient (DDPG), Twin Delayed DDPG (TD3) y Soft Actor-Critic (SAC), se inspiran indirectamente en la estructura de la HJB. Estos algoritmos utilizan redes neuronales para aproximar las funciones de valor (crítico) y las políticas (actor).

Consideremos el marco Actor-Critic:

  • El crítico aprende la función de valor $Q(s, a)$ (o $V(s)$), que se puede interpretar como una aproximación de la solución de la HJB. Su objetivo es minimizar el error de Bellman: $$ L_{critic} = \mathbb{E} \left[ (Q(s, a) - (R(s, a) + \gamma Q(s', \pi(s'))))^2 \right] $$ En tiempo continuo, esto se traduce en minimizar la discrepancia con respecto a la ecuación de Bellman en tiempo continuo.
  • El actor aprende la política $\pi(a|s)$ (o $a = \mu(s)$ para políticas deterministas) que maximiza la función de valor. Esto es análogo a resolver la operación $\min_u$ (o $\max_u$ para recompensas) en la HJB. Para una política determinista $\mu(s)$: $$ L_{actor} = - \mathbb{E} \left[ Q(s, \mu(s)) \right] $$

En SAC, la función objetivo incorpora un término de entropía para fomentar la exploración, lo que modifica la forma de la HJB subyacente, llevando a una función de valor suave y más robusta. La ecuación de Bellman para SAC es:

$$ Q(s,a) = \mathbb{E}{s' \sim P} \left[ R(s,a) + \gamma \left( \mathbb{E}{a' \sim \pi} [Q(s', a')] - \alpha \log \pi(a'|s') \right) \right] $$

Este enfoque resalta cómo la HJB, aunque no se resuelve explícitamente, sirve como el principio teórico que subyace a la formulación de las funciones objetivo de los algoritmos DRL para aproximar las soluciones de control óptimo. Los términos de gradiente y hessiano en la HJB se capturan implícitamente a través de las actualizaciones de la red neuronal y la propagación de errores.

Fragmento de Código Conceptual para un Crítico en DRL

import torch
import torch.nn as nn
import torch.optim as optim

class CriticNetwork(nn.Module):
    def __init__(self, state_dim, action_dim, hidden_size=256):
        super(CriticNetwork, self).__init__()
        self.fc1 = nn.Linear(state_dim + action_dim, hidden_size)
        self.fc2 = nn.Linear(hidden_size, hidden_size)
        self.fc3 = nn.Linear(hidden_size, 1)

    def forward(self, state, action):
        x = torch.cat([state, action], dim=1)
        x = torch.relu(self.fc1(x))
        x = torch.relu(self.fc2(x))
        return self.fc3(x)

class ActorNetwork(nn.Module):
    def __init__(self, state_dim, action_dim, hidden_size=256):
        super(ActorNetwork, self).__init__()
        self.fc1 = nn.Linear(state_dim, hidden_size)
        self.fc2 = nn.Linear(hidden_size, hidden_size)
        self.mu = nn.Linear(hidden_size, action_dim)
        self.log_std = nn.Linear(hidden_size, action_dim) # Para políticas estocásticas

    def forward(self, state):
        x = torch.relu(self.fc1(state))
        x = torch.relu(self.fc2(x))
        mu = torch.tanh(self.mu(x)) # Asegura que las acciones estén en un rango
        log_std = self.log_std(x).clamp(-20, 2) # Limitar para estabilidad
        std = torch.exp(log_std)
        return mu, std

# Simulación de un paso de entrenamiento de Q-network (Critic)
def train_critic_step(critic_net, target_critic_net, actor_net, states, actions, rewards, next_states, dones, gamma, alpha, critic_optimizer):
    # Calcular next_actions de la política actual (o de la target_actor para TD3/DDPG)
    with torch.no_grad():
        next_actions_mu, next_actions_std = actor_net(next_states)
        # Para SAC, muestrear desde la política y calcular log_prob
        policy_dist = torch.distributions.Normal(next_actions_mu, next_actions_std)
        next_actions_sampled = policy_dist.sample()
        next_log_probs = policy_dist.log_prob(next_actions_sampled).sum(axis=-1, keepdim=True)

        # Calcular target Q-value
        target_q_values = target_critic_net(next_states, next_actions_sampled)
        # Incorporar el término de entropía para SAC
        target_q_values_with_entropy = target_q_values - alpha * next_log_probs
        y_q = rewards + (1 - dones) * gamma * target_q_values_with_entropy

    # Calcular Q-value actual
    q_values = critic_net(states, actions)

    # Calcular pérdida y actualizar la red
    critic_loss = nn.functional.mse_loss(q_values, y_q)
    critic_optimizer.zero_grad()
    critic_loss.backward()
    critic_optimizer.step()

    return critic_loss.item()

Este pseudocódigo ilustra cómo el crítico en un algoritmo como SAC aprende a aproximar la solución a una forma modificada de la Ecuación de Bellman, que en el límite de tiempo continuo se relaciona con la HJB.

HJB y Modelos de Difusión

Los modelos de difusión estocástica han surgido como una clase potente de modelos generativos, capaces de sintetizar datos de alta calidad (imágenes, audio). Sorprendentemente, su formulación tiene conexiones profundas con el control óptimo y la HJB.

Un modelo de difusión funciona en dos fases:

  1. Proceso Forward (Difusión): Se añade ruido gaussiano de forma iterativa a los datos de forma controlada, transformando progresivamente una distribución de datos complejos $x_0 \sim q(x_0)$ en una simple distribución de ruido $x_T \sim \mathcal{N}(0, I)$. Esto se modela con una SDE: $$ dx_t = f(x_t, t) dt + g(t) dW_t $$ donde $f(x_t, t)$ es un término de deriva (drift) y $g(t)$ es un término de difusión que escala el ruido. Un ejemplo común es el proceso de Ornstein-Uhlenbeck: $dx_t = -\frac{1}{2}\beta(t)x_t dt + \sqrt{\beta(t)} dW_t$.

  2. Proceso Reverse (Generación): El objetivo es revertir el proceso de difusión, comenzando con ruido y transformándolo de nuevo en datos reales. La SDE del proceso inverso, que va de $T$ a $0$, tiene la forma: $$ dx_t = \left[ f(x_t, t) - g(t)^2 \nabla_{x_t} \log p_t(x_t) \right] dt + g(t) d\bar{W}_t $$ donde $d\bar{W}t$ es un proceso Wiener en tiempo invertido, y $\nabla{x_t} \log p_t(x_t)$ es la función de score de la distribución de los datos ruidosos en el tiempo $t$. Esta función de score guía el proceso de reversión, indicando la dirección en la que se debe mover el estado $x_t$ para alcanzar densidades de probabilidad más altas.

La conexión con el control óptimo se vuelve evidente cuando interpretamos el término $-\frac{1}{2} g(t)^2 \nabla_{x_t} \log p_t(x_t)$ como una política de control óptima. El problema de generar datos a partir de ruido puede verse como un problema de control estocástico donde el objetivo es "conducir" el proceso de $x_T$ (ruido) a $x_0$ (datos) de manera eficiente, típicamente minimizando alguna distancia entre la distribución generada y la distribución de datos reales.

En esta perspectiva, la HJB no se resuelve directamente para encontrar una política de control. Más bien, la función de valor óptimo (o una función relacionada) y la política óptima están intrínsecamente ligadas a la distribución de probabilidad $p_t(x_t)$ del proceso forward. La función de score, $\nabla_{x_t} \log p_t(x_t)$, es el gradiente logarítmico de esta densidad, y se aprende típicamente con redes neuronales mediante el denoising score matching.

El puente formal entre los modelos de difusión y la HJB se puede establecer a través de la teoría de control estocástico óptimo de transporte o problemas de puente de Schrödinger. En estos problemas, se busca encontrar el proceso estocástico más probable que conecta dos distribuciones de probabilidad marginales dadas en dos instantes de tiempo. La solución a estos problemas a menudo implica resolver una HJB acoplada con una Ecuación de Fokker-Planck (que describe la evolución temporal de la densidad de probabilidad del estado).

La Ecuación de Fokker-Planck para la densidad $p(x,t)$ asociada a la SDE forward es: $$ \frac{\partial p}{\partial t} = -\nabla_x \cdot (f(x, t) p(x, t)) + \frac{1}{2} \nabla_x^2 : (g(x, t) g(x, t)^T p(x, t)) $$

Si se considera un problema de control óptimo donde la función de costo está relacionada con la Divergencia KL entre la trayectoria del sistema y una trayectoria de referencia, o si se busca el control que transforma una densidad inicial en una final, la HJB y la Fokker-Planck se acoplan. El gradiente de la función de valor de la HJB para ciertos costos (por ejemplo, costos cuadráticos de control) puede estar directamente relacionado con el gradiente logarítmico de la densidad de probabilidad.

Ejemplo de Simulación de SDE y Concepto de Score

import numpy as np
import matplotlib.pyplot as plt

# Parámetros del proceso de Ornstein-Uhlenbeck (ejemplo de SDE forward)
beta_0 = 0.1
beta_1 = 20.0
T_max = 1.0
num_steps = 1000
dt = T_max / num_steps

# Función de beta (p.ej., lineal)
def beta(t):
    return beta_0 + t * (beta_1 - beta_0) / T_max

# Simulación del proceso forward SDE
def simulate_forward_sde(x0, T, num_steps):
    xs = [x0]
    x_t = x0
    for i in range(num_steps):
        t = i * dt
        # dx_t = -0.5 * beta(t) * x_t dt + sqrt(beta(t)) dW_t
        drift = -0.5 * beta(t) * x_t
        diffusion = np.sqrt(beta(t))
        dW = np.random.normal(0, np.sqrt(dt))
        x_t += drift * dt + diffusion * dW
        xs.append(x_t)
    return np.array(xs)

# Simular algunos datos iniciales (p.ej., de una distribución gaussiana)
x_initial = np.random.normal(0, 1, 100) # 100 puntos de datos iniciales

# Simular el proceso forward para un punto
example_x0 = 0.5
forward_path = simulate_forward_sde(example_x0, T_max, num_steps)

# Plot de la trayectoria forward (ruido incremental)
# plt.figure(figsize=(10, 6))
# plt.plot(np.linspace(0, T_max, num_steps + 1), forward_path)
# plt.title('Simulación de Proceso de Difusión Forward (Ornstein-Uhlenbeck)')
# plt.xlabel('Tiempo')
# plt.ylabel('Estado $x_t$')
# plt.grid(True)
# plt.show()

# Concepto de Score Function:
# En la práctica, la función de score $\nabla_{x_t} \log p_t(x_t)$ no se conoce analíticamente
# y se aprende con una red neuronal (Score Network) entrenada para predecir el ruido añadido.

# Red neuronal para aproximar el score function (conceptual)
class ScoreNetwork(nn.Module):
    def __init__(self, input_dim, hidden_size=256):
        super(ScoreNetwork, self).__init__()
        self.fc1 = nn.Linear(input_dim + 1, hidden_size) # input_dim para x, +1 para t
        self.fc2 = nn.Linear(hidden_size, hidden_size)
        self.fc3 = nn.Linear(hidden_size, input_dim) # Output es el score vector

    def forward(self, x, t):
        # t es un escalar, convertir a tensor para concatenar
        t_tensor = torch.full((x.shape[0], 1), t, dtype=x.dtype, device=x.device)
        input_vec = torch.cat([x, t_tensor], dim=1)
        h = torch.relu(self.fc1(input_vec))
        h = torch.relu(self.fc2(h))
        return self.fc3(h) # score_t(x_t)

# Función para simular el proceso inverso usando la ScoreNetwork (conceptual)
def simulate_reverse_sde_conceptual(score_net, x_T, T, num_steps):
    dt_rev = -T / num_steps # Paso de tiempo negativo para ir hacia atrás
    xs_rev = [x_T]
    x_t = x_T
    for i in range(num_steps):
        t = T + i * dt_rev # Tiempo que disminuye de T a 0
        t_tensor = torch.tensor(t, dtype=torch.float32).unsqueeze(0)
        x_tensor = torch.tensor(x_t, dtype=torch.float32).unsqueeze(0)

        with torch.no_grad():
            predicted_score = score_net(x_tensor, t_tensor).squeeze(0).numpy()

        drift_rev = -0.5 * beta(t) * x_t - beta(t) * predicted_score # Termino de difusión g(t)^2 = beta(t)
        diffusion_rev = np.sqrt(beta(t))
        dW_rev = np.random.normal(0, np.sqrt(abs(dt_rev))) # dW_bar en el tiempo invertido

        x_t += drift_rev * dt_rev + diffusion_rev * dW_rev
        xs_rev.append(x_t)
    return np.array(xs_rev)

# Esto solo es conceptual, una red real requiere entrenamiento intensivo
# Para una dimensión 1D, el score analítico para N(0, sigma^2) es -x/sigma^2
# p_t(x_t) se vuelve más complicado en el proceso forward

Este fragmento ilustra cómo la función de score, crucial en los modelos de difusión, está directamente relacionada con el término de control que guía el proceso inverso. La HJB ofrece una perspectiva de cómo estos "controles" emergen de la minimización de un costo o la maximización de una recompensa en un entorno estocástico.

Conclusiones

La Ecuación de Hamilton-Jacobi-Bellman es un concepto matemáticamente profundo y fundamental que unifica diversas áreas de la inteligencia artificial y la ingeniería. Desde el control óptimo clásico hasta el aprendizaje por refuerzo moderno y, más recientemente, los modelos generativos de difusión, la HJB proporciona el marco teórico para comprender cómo se toman decisiones óptimas en entornos dinámicos y estocásticos.

En el aprendizaje por refuerzo continuo, la HJB sirve como el principio subyacente que inspira la formulación de los objetivos de las redes neuronales Actor-Critic, permitiendo la aproximación de funciones de valor y políticas óptimas sin la necesidad de resolver explícitamente la PDE. En los modelos de difusión, la HJB se manifiesta a través de la interpretación del proceso inverso como un problema de control estocástico, donde la función de score actúa como la política que guía la transformación del ruido en datos.

La dificultad intrínseca de resolver analíticamente la HJB, especialmente en alta dimensionalidad, ha impulsado la innovación en métodos computacionales, incluyendo el uso de redes neuronales y algoritmos de muestreo estocástico. A medida que continuamos explorando la complejidad de los sistemas autónomos y la generación de datos realistas, la HJB y sus extensiones seguirán siendo herramientas esenciales para formular y resolver estos desafíos.

Para aquellos interesados en profundizar en las aplicaciones prácticas de la ingeniería de datos y la inteligencia artificial, incluyendo el control óptimo, el aprendizaje por refuerzo y los modelos generativos avanzados, les invitamos a visitar https://www.mgatc.com para conocer nuestros servicios de consultoría especializada.

Escrito por Mariano Gobea Alcoba

← Otros artículosConsultoría →