4.7 Redes neuronales informadas por la física - Machine Learning en las Geociencias

1. Introducción¶

Las redes neuronales informadas por la física (Physics-Informed Neural Networks, PINN) incorporan leyes físicas, expresadas como ecuaciones diferenciales ordinarias o en derivadas parciales (EDO/EDP), directamente en el entrenamiento de una red neuronal. La física entra a través de la función de pérdida: además de ajustar los datos, la red se penaliza cuando sus predicciones violan la ecuación de gobierno. Las PINN pueden:

resolver ecuaciones diferenciales con o sin datos,
trabajar con pocas muestras etiquetadas, porque la física restringe la solución,
incorporar conocimiento previo a un flujo de trabajo de aprendizaje automático.

Las geociencias están llenas de procesos gobernados por ecuaciones bien conocidas: la difusión del calor en la corteza, la propagación de ondas sísmicas, el flujo de agua subterránea (ley de Darcy) y la dinámica de fluidos (Navier-Stokes). Las PINN fueron introducidas por Raissi et al. (2019) y revisadas para la comunidad física en general por Karniadakis et al. (2021); vea las referencias al final de este cuaderno. Ejemplos de uso en geociencias:

Difusión de temperatura en la corteza: inferir la estructura térmica y el potencial geotérmico a partir de perfiles de temperatura escasos en pozos.
Propagación de ondas sísmicas: resolver ecuaciones de onda en medios heterogéneos sin mallar el dominio.
Flujo de agua subterránea: resolver las ecuaciones de Darcy y de advección-difusión para modelar la respuesta de un acuífero y el transporte de contaminantes.

En los tres casos el atractivo es el mismo: los datos de campo son escasos y ruidosos, pero la EDP de gobierno se conoce. La física rellena donde los datos no llegan.

🖥️ Diapositivas de clase — Sesión 26 (mié 2 dic, contexto de enriquecimiento)

2. Marco matemático¶

Una PINN tiene tres componentes:

Una red neuronal $u_\theta(x, t)$ con parámetros θ que aproxima la solución $u(x, t)$ .
Una ecuación de gobierno escrita en forma de residuo $\mathcal{L}(u, x, t) = 0$ , donde $\mathcal{L}$ es un operador diferencial. Las derivadas de $u_\theta$ que aparecen en $\mathcal{L}$ se calculan de manera exacta con diferenciación automática, la misma maquinaria que usa la retropropagación.
Una función de pérdida compuesta

\text{Loss} = \text{Data loss} + \text{Physics loss} + \text{Initial/boundary condition loss},

(1)

donde

la pérdida de datos mide el desajuste respecto de las observaciones,
la pérdida física es el residuo cuadrático medio de la EDP $\mathcal{L}(u_\theta)$ evaluado en puntos de colocación muestreados dentro del dominio,
la pérdida de condiciones iniciales y de contorno fija la solución allí donde se conoce.

El entrenamiento minimiza la pérdida combinada con un optimizador basado en gradientes como Adam. Los pesos relativos de los tres términos son hiperparámetros, y equilibrarlos es una de las dificultades prácticas de las PINN (más sobre esto en la sección 5).

import functools
import time

import matplotlib.pyplot as plt
import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim

device = torch.device("cuda" if torch.cuda.is_available()
                      else "mps" if torch.backends.mps.is_available()
                      else "cpu")
# PINN training differentiates through gradients (double-backward autograd),
# which is most reliable on CPU, and these models are tiny: we train on CPU.
DEVICE = torch.device("cpu")
torch.set_num_threads(1)  # single thread is fastest for such small tensors

torch.manual_seed(42)
np.random.seed(10)

3. Calentamiento: la ley de enfriamiento de Newton¶

Antes de abordar una EDP, empezamos con una EDO, donde la diferenciación automática solo necesita primeras derivadas. Considere una muestra de roca que se enfría hacia la temperatura ambiente $T_{\text{env}}$ . La ley de enfriamiento de Newton establece que

\frac{dT}{dt} = R \, (T_{\text{env}} - T),

(2)

donde $R$ es una constante de tasa de enfriamiento. Con temperatura inicial $T(0) = T_0$ , la solución analítica es

T(t) = T_{\text{env}} + (T_0 - T_{\text{env}}) \, e^{-Rt}.

(3)

Haremos una comparación de tres vías sobre las mismas 10 mediciones ruidosas, todas tomadas temprano en la historia de enfriamiento:

una red neuronal simple (modelo de referencia),
la misma red con regularización L2 de los pesos,
la misma red con una pérdida física (una PINN).

La pregunta es cuál de los modelos extrapola correctamente más allá del rango temporal que cubren los datos.

def cooling_law(time, Tenv, T0, R):
    T = Tenv + (T0 - Tenv) * np.exp(-R * time)
    return T

Generamos datos sintéticos ruidosos: la curva verdadera abarca 1000 s, pero las 10 muestras de entrenamiento solo cubren los primeros 300 s.

Tenv = 25
T0 = 100
R = 0.005
times = np.linspace(0, 1000, 1000)
eq = functools.partial(cooling_law, Tenv=Tenv, T0=T0, R=R)
temps = eq(times)

# Make training data
t = np.linspace(0, 300, 10)
T = eq(t) + 2 * np.random.randn(10)

plt.figure(figsize=(6, 4))
plt.plot(times, temps)
plt.plot(t, T, 'o')
plt.legend(['Equation', 'Training data'])
plt.ylabel('Temperature (C)')
plt.xlabel('Time (s)')
plt.title('Newton cooling: truth and noisy samples')
plt.show()

3.1 Una red con una segunda pérdida intercambiable¶

La clase Net de abajo es una red pequeña totalmente conectada. Su método fit minimiza el desajuste con los datos más una segunda pérdida opcional loss2, ponderada por loss2_weight. Reutilizaremos la misma clase en los tres experimentos y solo intercambiaremos loss2: None para el modelo de referencia, una penalización L2 para la red regularizada y un residuo físico para la PINN.

El ayudante grad calcula la derivada de las salidas de la red con respecto a las entradas usando torch.autograd.grad con create_graph=True, de modo que el resultado pueda a su vez derivarse durante la retropropagación. Este es el truco central de las PINN.

def np_to_th(x):
    """Convert a numpy array to a float32 torch tensor of shape (n, -1)."""
    n_samples = len(x)
    return torch.from_numpy(x).to(torch.float).to(DEVICE).reshape(n_samples, -1)


def grad(outputs, inputs):
    """Partial derivative of outputs with respect to inputs.

    Args:
        outputs: (N, 1) tensor
        inputs: (N, D) tensor
    """
    return torch.autograd.grad(
        outputs, inputs, grad_outputs=torch.ones_like(outputs), create_graph=True
    )


class Net(nn.Module):
    def __init__(
        self,
        input_dim,
        output_dim,
        n_units=100,
        epochs=1000,
        loss=nn.MSELoss(),
        lr=1e-3,
        loss2=None,
        loss2_weight=0.1,
    ) -> None:
        super().__init__()

        self.epochs = epochs
        self.loss = loss
        self.loss2 = loss2
        self.loss2_weight = loss2_weight
        self.lr = lr
        self.n_units = n_units

        self.layers = nn.Sequential(
            nn.Linear(input_dim, self.n_units),
            nn.ReLU(),
            nn.Linear(self.n_units, self.n_units),
            nn.ReLU(),
            nn.Linear(self.n_units, self.n_units),
            nn.ReLU(),
            nn.Linear(self.n_units, self.n_units),
            nn.ReLU(),
        )
        self.out = nn.Linear(self.n_units, output_dim)

    def forward(self, x):
        h = self.layers(x)
        out = self.out(h)
        return out

    def fit(self, X, y):
        Xt = np_to_th(X)
        yt = np_to_th(y)

        optimiser = optim.Adam(self.parameters(), lr=self.lr)
        self.train()
        losses = []
        for ep in range(self.epochs):
            optimiser.zero_grad()
            outputs = self.forward(Xt)
            loss = self.loss(yt, outputs)
            if self.loss2:
                loss += self.loss2_weight * self.loss2(self)
            loss.backward()
            optimiser.step()
            losses.append(loss.item())
            if ep % int(self.epochs / 10) == 0:
                print(f"Epoch {ep}/{self.epochs}, loss: {losses[-1]:.2f}")
        return losses

    def predict(self, X):
        self.eval()
        out = self.forward(np_to_th(X))
        return out.detach().cpu().numpy()

3.2 Red de referencia¶

Entrene sobre las 10 muestras ruidosas sin ningún término de pérdida adicional.

net = Net(1, 1, loss2=None, epochs=2000, lr=1e-4).to(DEVICE)

losses = net.fit(t, T)

plt.figure(figsize=(6, 3))
plt.plot(losses)
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('Baseline network: training loss')
plt.show()

Epoch 0/2000, loss: 4713.88
Epoch 200/2000, loss: 2527.58
Epoch 400/2000, loss: 2299.78

Epoch 600/2000, loss: 1233.46
Epoch 800/2000, loss: 74.98
Epoch 1000/2000, loss: 3.04

Epoch 1200/2000, loss: 2.11
Epoch 1400/2000, loss: 1.84
Epoch 1600/2000, loss: 1.69

Epoch 1800/2000, loss: 1.57

prediction_temp_baseline = net.predict(times)

plt.figure(figsize=(6, 4))
plt.plot(times, prediction_temp_baseline, alpha=0.8)
plt.plot(t, T, 'o')
plt.legend(['Prediction', 'Training data'])
plt.ylabel('Temperature (C)')
plt.xlabel('Time (s)')
plt.title('Baseline network prediction')
plt.show()

El modelo de referencia ajusta las muestras que vio y hace lo que se le antoja después de los 300 s: nada restringe la extrapolación.

3.3 Red regularizada con L2¶

Un remedio clásico contra el sobreajuste es penalizar los pesos grandes. Pasamos la norma L2 de los parámetros como loss2.

# Second loss: L2 norm of the network weights
def l2_reg(model: torch.nn.Module):
    return torch.sum(sum([p.pow(2.) for p in model.parameters()]))

netreg = Net(1, 1, loss2=l2_reg, epochs=20000, lr=1e-4, loss2_weight=1).to(DEVICE)

losses = netreg.fit(t, T)

plt.figure(figsize=(6, 3))
plt.plot(losses)
plt.yscale('log')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('L2-regularized network: training loss')
plt.show()

Epoch 0/20000, loss: 11219.99

Epoch 2000/20000, loss: 3914.45

Epoch 4000/20000, loss: 2424.52

Epoch 6000/20000, loss: 1626.11

Epoch 8000/20000, loss: 1206.96

Epoch 10000/20000, loss: 1010.50

Epoch 12000/20000, loss: 912.78

Epoch 14000/20000, loss: 840.27

Epoch 16000/20000, loss: 775.95

Epoch 18000/20000, loss: 715.29

prediction_temp_regularized = netreg.predict(times)

plt.figure(figsize=(6, 4))
plt.plot(times, temps, alpha=0.8)
plt.plot(t, T, 'o')
plt.plot(times, prediction_temp_regularized, alpha=0.8)
plt.plot(times, prediction_temp_baseline, alpha=0.8)
plt.legend(labels=['Equation', 'Training data', 'Regularized', 'Baseline'])
plt.ylabel('Temperature (C)')
plt.xlabel('Time (s)')
plt.title('Baseline vs. L2 regularization')
plt.show()

La regularización suaviza la predicción, pero no tiene razón alguna para seguir el decaimiento exponencial hacia $T_{\text{env}}$ : simplemente mantiene pequeños los pesos.

3.4 PINN¶

Ahora reemplace la penalización L2 por una pérdida física. Evaluamos la red en 1000 tiempos de colocación que abarcan la ventana completa de 1000 s (mucho más allá de los datos), derivamos la salida con respecto al tiempo usando grad y penalizamos el residuo cuadrático medio de la ley de enfriamiento

r(t) = R\,(T_{\text{env}} - T_\theta(t)) - \frac{dT_\theta}{dt}.

(4)

Note que la pérdida física no usa ninguna observación de temperatura, solo la ecuación.

def physics_loss(model: torch.nn.Module):
    ts = torch.linspace(0, 1000, steps=1000).view(-1, 1).requires_grad_(True).to(DEVICE)
    temps = model(ts)
    dT = grad(temps, ts)[0]
    pde = R * (Tenv - temps) - dT
    return torch.mean(pde**2)

net_PINN = Net(1, 1, loss2=physics_loss, epochs=10000, loss2_weight=1, lr=1e-4).to(DEVICE)

losses = net_PINN.fit(t, T)

plt.figure(figsize=(6, 3))
plt.plot(losses)
plt.yscale('log')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('PINN: training loss')
plt.show()

Epoch 0/10000, loss: 4775.84

Epoch 1000/10000, loss: 1.04

Epoch 2000/10000, loss: 0.63

Epoch 3000/10000, loss: 1.32

Epoch 4000/10000, loss: 0.63

Epoch 5000/10000, loss: 2.21

Epoch 6000/10000, loss: 0.35

Epoch 7000/10000, loss: 0.58

Epoch 8000/10000, loss: 0.76

Epoch 9000/10000, loss: 0.66

prediction_temp_pinn = net_PINN.predict(times)

plt.figure(figsize=(6, 4))
plt.plot(times, temps, alpha=0.8)
plt.plot(t, T, 'o')
plt.plot(times, prediction_temp_baseline, alpha=0.8)
plt.plot(times, prediction_temp_regularized, alpha=0.8)
plt.plot(times, prediction_temp_pinn, alpha=0.8)
plt.legend(labels=['Equation', 'Training data', 'Baseline', 'Regularized', 'PINN'])
plt.ylabel('Temperature (C)')
plt.xlabel('Time (s)')
plt.title('Three-way comparison')
plt.show()

La PINN sigue el decaimiento exponencial verdadero mucho más allá del último dato, porque la pérdida física restringe la solución en todas partes hasta donde llegan los puntos de colocación. Los datos anclan la amplitud; la ecuación aporta la forma. Esa es toda la idea.

El entrenamiento aquí es deliberadamente corto. En su propia máquina puede subir el número de épocas para obtener ajustes más suaves.

4. Una EDP de verdad: difusión de calor 1D en el subsuelo somero¶

La ley de enfriamiento era una EDO. Ahora resolvemos una EDP real con una PINN: la difusión de calor,

\frac{\partial T}{\partial t} = k \, \frac{\partial^2 T}{\partial x^2},

(5)

donde $T(x,t)$ es la temperatura y $k$ es la difusividad térmica.

Contexto físico. Piense en una anomalía de temperatura en los primeros $L = 10$ m de suelo o roca, la capa activa de un permafrost o un perfil geotérmico somero. Una anomalía cálida alcanza su máximo a media columna y se relaja por difusión, mientras que la superficie ( $x=0$ ) y la base ( $x=L$ ) permanecen fijas en la temperatura media anual, que tomamos como la referencia de 0 °C. Una difusividad térmica típica de roca o suelo es $k \approx 10^{-6}\ \mathrm{m^2/s}$ .

Adimensionalización. Los valores crudos en unidades SI ( $L=10$ , $k=10^{-6}$ , $t$ en unidades de 10⁷ s) producen números feos para una red neuronal, así que reescalamos: $\hat{x} = x/L$ y $\hat{t} = k t / L^2$ . Una unidad de $\hat{t}$ corresponde a $L^2/k = 10^8\ \mathrm{s} \approx 3.2$ años; resolvemos en $\hat{t} \in [0, 0.3]$ , aproximadamente un año. La temperatura se mantiene en °C. En estas unidades la EDP queda como $\partial T/\partial \hat{t} = \partial^2 T/\partial \hat{x}^2$ (difusividad 1).

Referencia analítica. Un único modo senoidal difunde sin cambiar de forma:

T(\hat{x}, \hat{t}) = T_0 \, \sin(\pi \hat{x}) \, e^{-\pi^2 \hat{t}},

(6)

con $T_0 = 5$ °C. Satisface la EDP y las condiciones de contorno $T=0$ en $\hat{x}=0$ y $\hat{x}=1$ , así que podemos medir el error de la PINN de manera exacta.

Red. $T_\theta(\hat{x}, \hat{t})$ es un MLP pequeño: 2 entradas, tres capas ocultas de 64 unidades, 1 salida. Usamos activaciones tanh en lugar de ReLU porque el residuo de la EDP necesita segundas derivadas, y la segunda derivada de ReLU es cero en casi todas partes.

T0_heat = 5.0   # amplitude of the initial anomaly, deg C
t_max = 0.3     # nondimensional time horizon (~1 year)


def analytic_T(x, t):
    """Analytic solution: single diffusing sine mode (nondimensional x, t)."""
    return T0_heat * np.sin(np.pi * x) * np.exp(-np.pi**2 * t)


class HeatPINN(nn.Module):
    """T_theta(x, t): 2 inputs -> 3 x 64 tanh -> 1 output."""

    def __init__(self, n_units=64):
        super().__init__()
        self.layers = nn.Sequential(
            nn.Linear(2, n_units), nn.Tanh(),
            nn.Linear(n_units, n_units), nn.Tanh(),
            nn.Linear(n_units, n_units), nn.Tanh(),
            nn.Linear(n_units, 1),
        )

    def forward(self, xt):
        return self.layers(xt)

4.1 Puntos de entrenamiento¶

Tres conjuntos de puntos, uno por término de pérdida:

Puntos de datos: 40 muestras ruidosas de temperatura en ubicaciones $(\hat{x}, \hat{t})$ aleatorias, como lecturas dispersas de sensores, con ruido gaussiano de 0.2 °C.
Puntos de colocación: 1000 puntos interiores aleatorios donde imponemos el residuo de la EDP. No llevan asociado ningún valor de temperatura.
Puntos iniciales y de contorno: 200 puntos sobre la recta $\hat{t}=0$ (perfil inicial) y 200 sobre cada contorno $\hat{x}=0$ y $\hat{x}=1$ (temperatura mantenida en 0 °C).

# Data points: sparse noisy samples of the true solution
n_data = 40
x_d = np.random.rand(n_data)
t_d = np.random.rand(n_data) * t_max
T_d = analytic_T(x_d, t_d) + 0.2 * np.random.randn(n_data)
X_data = torch.tensor(np.stack([x_d, t_d], axis=1), dtype=torch.float32)
y_data = torch.tensor(T_d[:, None], dtype=torch.float32)

# Collocation points: interior points where the PDE residual is enforced
n_col = 1000
x_c = torch.rand(n_col, 1)
t_c = torch.rand(n_col, 1) * t_max
X_col = torch.cat([x_c, t_c], dim=1).requires_grad_(True)

# Initial condition points (t = 0) and boundary points (x = 0 and x = 1)
n_b = 200
x_ic = torch.rand(n_b, 1)
X_ic = torch.cat([x_ic, torch.zeros(n_b, 1)], dim=1)
y_ic = T0_heat * torch.sin(np.pi * x_ic)
t_b = torch.rand(n_b, 1) * t_max
X_b0 = torch.cat([torch.zeros(n_b, 1), t_b], dim=1)
X_b1 = torch.cat([torch.ones(n_b, 1), t_b], dim=1)

print(f"data: {X_data.shape}, collocation: {X_col.shape}, "
      f"IC: {X_ic.shape}, BC: {X_b0.shape} + {X_b1.shape}")

data: torch.Size([40, 2]), collocation: torch.Size([1000, 2]), IC: torch.Size([200, 2]), BC: torch.Size([200, 2]) + torch.Size([200, 2])

4.2 Bucle de entrenamiento¶

Cada paso calcula tres términos de pérdida y los suma con pesos iguales:

Pérdida de datos: MSE entre $T_\theta$ y las muestras ruidosas.
Pérdida de la EDP: en los puntos de colocación, grad entrega $\partial T_\theta / \partial \hat{x}$ y $\partial T_\theta / \partial \hat{t}$ en una sola llamada (las columnas del gradiente); una segunda llamada a grad sobre la primera columna entrega $\partial^2 T_\theta / \partial \hat{x}^2$ . El residuo es $\partial T_\theta/\partial \hat{t} - \partial^2 T_\theta/\partial \hat{x}^2$ .
Pérdida de condiciones iniciales y de contorno: MSE contra el perfil senoidal inicial más la temperatura al cuadrado en ambos contornos.

Unos pocos miles de pasos de Adam bastan para este problema suave.

model_heat = HeatPINN().to(DEVICE)
optimizer = optim.Adam(model_heat.parameters(), lr=2e-3)
mse = nn.MSELoss()

n_steps = 2000
history = {"data": [], "pde": [], "icbc": []}

tic = time.perf_counter()
for step in range(n_steps):
    optimizer.zero_grad()

    # Data loss
    loss_data = mse(model_heat(X_data), y_data)

    # PDE residual loss at collocation points
    T_c = model_heat(X_col)
    g = grad(T_c, X_col)[0]           # (n_col, 2): columns are dT/dx, dT/dt
    T_x, T_t = g[:, 0:1], g[:, 1:2]
    T_xx = grad(T_x, X_col)[0][:, 0:1]
    loss_pde = torch.mean((T_t - T_xx) ** 2)

    # Initial and boundary condition loss
    loss_icbc = (mse(model_heat(X_ic), y_ic)
                 + torch.mean(model_heat(X_b0) ** 2)
                 + torch.mean(model_heat(X_b1) ** 2))

    loss = loss_data + loss_pde + loss_icbc
    loss.backward()
    optimizer.step()

    history["data"].append(loss_data.item())
    history["pde"].append(loss_pde.item())
    history["icbc"].append(loss_icbc.item())
    if step % 400 == 0 or step == n_steps - 1:
        print(f"step {step:4d}  data {loss_data.item():.4f}  "
              f"pde {loss_pde.item():.5f}  ic/bc {loss_icbc.item():.5f}")
pinn_train_seconds = time.perf_counter() - tic
print(f"wall-clock training time: {pinn_train_seconds:.1f} s")

step    0  data 2.7585  pde 0.00039  ic/bc 12.18578

step  400  data 0.0475  pde 0.02406  ic/bc 0.01504

step  800  data 0.0425  pde 0.01021  ic/bc 0.00182

step 1200  data 0.0416  pde 0.00178  ic/bc 0.00086

step 1600  data 0.0413  pde 0.00164  ic/bc 0.00100

step 1999  data 0.0414  pde 0.01566  ic/bc 0.00158
wall-clock training time: 13.6 s

4.3 Resultados¶

Primero, los perfiles de temperatura predichos frente a los analíticos en cuatro instantes. El eje de profundidad vuelve a estar en metros ( $x = L\hat{x}$ ).

L = 10.0  # m
x_grid = np.linspace(0, 1, 200)
plot_times = [0.0, 0.05, 0.15, 0.3]
years_per_that = 1e8 / (365.25 * 86400)  # one unit of t_hat in years

plt.figure(figsize=(7, 4.5))
colors = plt.cm.viridis(np.linspace(0, 0.85, len(plot_times)))
for that, c in zip(plot_times, colors):
    X_plot = torch.tensor(np.stack([x_grid, np.full_like(x_grid, that)], axis=1),
                          dtype=torch.float32)
    T_pred = model_heat(X_plot).detach().numpy().ravel()
    label_t = f"t = {that * years_per_that:.2f} yr"
    plt.plot(x_grid * L, analytic_T(x_grid, that), '-', color=c,
             label=f'analytic, {label_t}')
    plt.plot(x_grid * L, T_pred, '--', color=c, label=f'PINN, {label_t}')
plt.xlabel('Depth x (m)')
plt.ylabel('Temperature anomaly (C)')
plt.title('Heat diffusion: PINN vs. analytic solution')
plt.legend(fontsize=8)
plt.show()

# Quantify the error (kept in a dict for the comparison in Section 4.4)
pinn_max_err = {}
for that in plot_times:
    X_plot = torch.tensor(np.stack([x_grid, np.full_like(x_grid, that)], axis=1),
                          dtype=torch.float32)
    T_pred = model_heat(X_plot).detach().numpy().ravel()
    pinn_max_err[that] = np.abs(T_pred - analytic_T(x_grid, that)).max()
    print(f"t_hat = {that:.2f}: max abs error = {pinn_max_err[that]:.3f} C")

t_hat = 0.00: max abs error = 0.029 C
t_hat = 0.05: max abs error = 0.023 C
t_hat = 0.15: max abs error = 0.023 C
t_hat = 0.30: max abs error = 0.024 C

plt.figure(figsize=(7, 4))
plt.plot(history["data"], label='data loss')
plt.plot(history["pde"], label='PDE residual loss')
plt.plot(history["icbc"], label='IC/BC loss')
plt.yscale('log')
plt.xlabel('Adam iteration')
plt.ylabel('Loss')
plt.title('Loss components during training')
plt.legend()
plt.show()

La PINN reproduce los perfiles analíticos con un error de unas pocas centésimas de grado a partir de 40 mediciones puntuales ruidosas, porque la EDP y las condiciones de contorno cargan con la mayor parte de la información. La pérdida de datos se estanca cerca del piso de ruido (la red no debería ajustar el ruido de 0.2 °C), mientras que las pérdidas de la EDP y de las condiciones iniciales y de contorno caen tres o cuatro órdenes de magnitud. Los picos periódicos en esas curvas son Adam pasándose brevemente del mínimo y recuperándose, algo común en el entrenamiento de PINN.

El entrenamiento se mantiene corto a propósito; en su propia máquina, suba n_steps y n_col para obtener residuos más ajustados.

Ejercicio. Ponga en cero el peso de la pérdida de la EDP (reemplace loss_pde por 0 * loss_pde en la suma) y reentrene. Compare los perfiles en $\hat{t} = 0.3$ , donde hay pocos datos y la anomalía es pequeña.

Solución

Sin el término de la EDP, el modelo se convierte en una regresión simple sobre 40 puntos ruidosos más las restricciones de condiciones iniciales y de contorno. Cerca del perfil inicial, donde los datos son densos en relación con la señal, todavía se ve bien. En tiempos posteriores el perfil predicho se aleja del decaimiento exponencial: nada obliga al interior del dominio a comportarse de manera difusiva, así que la red interpola el ruido en su lugar. El error en t_hat = 0.3 típicamente crece varias veces. Esto refleja la ablación de la ley de enfriamiento de la sección 3: el término físico es lo que vuelve confiable la extrapolación.

4.4 El modelo de referencia: 15 líneas de diferencias finitas¶

La regla de este libro es modelos de referencia antes que modelos, y el capítulo de física no está exento. El problema directo que acabamos de resolver, con difusividad conocida, condiciones iniciales y de contorno conocidas y una solución suave, es territorio de manual para el análisis numérico clásico, así que la PINN tiene que enfrentarse al método clásico antes de que la elogiemos. El esquema más simple es FTCS (forward-time, centered-space): ponga $\hat{x}$ sobre una malla, reemplace $\partial^2 T/\partial \hat{x}^2$ por la segunda diferencia centrada y avance en el tiempo con Euler explícito:

T_i^{n+1} = T_i^n + r\,(T_{i+1}^n - 2T_i^n + T_{i-1}^n), \qquad r = \frac{\Delta \hat{t}}{\Delta \hat{x}^2}.

(7)

El esquema solo es estable para $r \le 1/2$ ; usamos $r = 0.4$ . Esa restricción es el precio de un método explícito (Crank–Nicolson la elimina a cambio de resolver un sistema tridiagonal), y en este problema no cuesta casi nada. El solucionador completo, cronometrado:

# FTCS finite differences: the entire solver is the 8 lines between tic and toc
nx = 101
dx = 1.0 / (nx - 1)
dt = 0.4 * dx**2                      # explicit stability requires dt <= dx^2 / 2
nt = int(round(t_max / dt))
x_fd = np.linspace(0.0, 1.0, nx)
r = dt / dx**2
save_steps = {int(round(that / dt)): that for that in plot_times if that > 0}

tic = time.perf_counter()
T_fd = T0_heat * np.sin(np.pi * x_fd)             # initial condition
fd_snapshots = {0.0: T_fd.copy()}
for n in range(1, nt + 1):
    T_fd[1:-1] += r * (T_fd[2:] - 2 * T_fd[1:-1] + T_fd[:-2])
    if n in save_steps:                           # endpoints never updated: T = 0 C
        fd_snapshots[save_steps[n]] = T_fd.copy()
fd_seconds = time.perf_counter() - tic

print(f"FTCS: {nx} grid points, {nt} time steps, "
      f"wall-clock {fd_seconds * 1e3:.0f} ms")
print(f"\n{'t_hat':>6} {'FD max err (C)':>16} {'PINN max err (C)':>18}")
for that in plot_times:
    fd_err = np.abs(fd_snapshots[that] - analytic_T(x_fd, that)).max()
    print(f"{that:6.2f} {fd_err:16.2e} {pinn_max_err[that]:18.3f}")
print(f"\nwall-clock: FD {fd_seconds * 1e3:.0f} ms vs. PINN training "
      f"{pinn_train_seconds:.0f} s -> the PINN is "
      f"{pinn_train_seconds / fd_seconds:.0f}x slower")

FTCS: 101 grid points, 7500 time steps, wall-clock 20 ms

 t_hat   FD max err (C)   PINN max err (C)
  0.00         0.00e+00              0.029
  0.05         1.73e-04              0.023
  0.15         1.94e-04              0.023
  0.30         8.83e-05              0.024

wall-clock: FD 20 ms vs. PINN training 14 s -> the PINN is 698x slower

El veredicto no es reñido. El solucionador de diferencias finitas termina en unos 20 ms y coincide con la solución analítica hasta unos $2 \times 10^{-4}$ °C; la PINN necesitó 14 s de entrenamiento para quedar dentro de 0.03 °C. En este problema directo el solucionador clásico gana por un factor de aproximadamente 750 en velocidad y por dos órdenes de magnitud en exactitud, sin tasa de aprendizaje, sin arquitectura y sin semilla. Esa es la lectura honesta de la sección 4.3, y se generaliza: cuando los coeficientes, la condición inicial y las condiciones de contorno se conocen y existe una discretización estándar, use la discretización estándar. El argumento a favor de la PINN hay que construirlo en otra parte.

Hay dos cosas, sin embargo, que el solucionador de diferencias finitas no hizo. Nunca tocó las 40 mediciones interiores ruidosas, porque un solucionador directo no tiene ranura para observaciones dispersas. Y exigió $k$ como entrada. Ambas carencias apuntan al mismo lugar, y hacia allá vamos.

4.5 El problema inverso: recuperar la difusividad¶

La sección 5 argumentará que los problemas inversos son el primer régimen donde las PINN siguen siendo la herramienta correcta. Aquí está la demostración, sobre los datos que ya tenemos. Suponga que la difusividad es desconocida. En unidades escaladas, el valor del generador es exactamente $\hat{k} = 1$ , porque la adimensionalización usó el $k = 10^{-6}\ \mathrm{m^2/s}$ verdadero; al modelo no se le dice eso. Inicializamos un $\hat{k}$ entrenable en 0.2, cinco veces demasiado bajo, y dejamos que el optimizador lo recupere a partir de las mismas 40 muestras ruidosas.

El cambio respecto del bucle de entrenamiento de la sección 4.2 es de cinco líneas, marcadas (1)–(4) abajo: $\hat{k}$ se convierte en un nn.Parameter (almacenado como $\log \hat{k}$ para que se mantenga positivo), se suma a la lista de parámetros del optimizador y el residuo pasa a ser $\partial T_\theta/\partial \hat{t} - \hat{k}\, \partial^2 T_\theta/\partial \hat{x}^2$ . Todo lo demás queda intacto. La información fluye así: la condición inicial fija la forma de la anomalía, las 40 muestras fijan qué tan rápido decae, y la EDP ata esa tasa de decaimiento a $\hat{k}$ .

model_inv = HeatPINN().to(DEVICE)
log_k = nn.Parameter(torch.log(torch.tensor(0.2)))    # (1) trainable k_hat, start 5x too low
optimizer_inv = optim.Adam(list(model_inv.parameters()) + [log_k],  # (2) k joins the optimizer
                           lr=2e-3)

k_history = []
tic = time.perf_counter()
for step in range(n_steps):
    optimizer_inv.zero_grad()
    loss_data = mse(model_inv(X_data), y_data)
    T_c = model_inv(X_col)
    g = grad(T_c, X_col)[0]
    T_x, T_t = g[:, 0:1], g[:, 1:2]
    T_xx = grad(T_x, X_col)[0][:, 0:1]
    loss_pde = torch.mean((T_t - torch.exp(log_k) * T_xx) ** 2)     # (3) k in the residual
    loss_icbc = (mse(model_inv(X_ic), y_ic)
                 + torch.mean(model_inv(X_b0) ** 2)
                 + torch.mean(model_inv(X_b1) ** 2))
    (loss_data + loss_pde + loss_icbc).backward()
    optimizer_inv.step()
    k_history.append(torch.exp(log_k).item())                       # (4) track the estimate
    if step % 400 == 0 or step == n_steps - 1:
        print(f"step {step:4d}  k_hat = {k_history[-1]:.3f}")
inv_seconds = time.perf_counter() - tic

k_rec = k_history[-1]
print(f"\nrecovered k_hat = {k_rec:.3f}  (generator truth 1.0, "
      f"error {abs(k_rec - 1.0) * 100:.1f}%)")
print(f"physical units: k = {k_rec * 1e-6:.2e} m^2/s  (truth 1.00e-06)")
print(f"wall-clock: {inv_seconds:.0f} s")

step    0  k_hat = 0.200

step  400  k_hat = 0.460

step  800  k_hat = 0.728

step 1200  k_hat = 0.875

step 1600  k_hat = 0.947

step 1999  k_hat = 0.981

recovered k_hat = 0.981  (generator truth 1.0, error 1.9%)
physical units: k = 9.81e-07 m^2/s  (truth 1.00e-06)
wall-clock: 14 s

plt.figure(figsize=(6, 3.5))
plt.plot(k_history)
plt.axhline(1.0, color='k', ls=':', label='generator truth')
plt.xlabel('Adam iteration')
plt.ylabel(r'$\hat{k}$ estimate')
plt.title('Diffusivity recovered jointly with the temperature field')
plt.legend()
plt.show()

La difusividad recuperada queda dentro del 2 % del valor verdadero del generador, a partir de 40 mediciones puntuales ruidosas y de nada que esté en malla. Seamos claros sobre qué es esta corrida: una calibración, también llamada ajuste histórico (history matching). El solucionador de diferencias finitas de la sección 4.4 no puede hacerla por sí solo, porque necesita $k$ como entrada. La vía clásica envuelve el solucionador en un bucle externo de optimización que lo vuelve a correr para cada $k$ candidato; en la práctica hidrogeológica eso es lo que hace PEST alrededor de un modelo MODFLOW al calibrar la conductividad hidráulica contra registros escasos de pozos. La PINN fusiona los dos bucles: una sola corrida de entrenamiento ajusta el campo de temperatura e invierte el parámetro al mismo tiempo, y el mismo truco se extiende a muchas incógnitas (una $k(\hat{x})$ variable en el espacio puede ser una segunda red pequeña). Datos puntuales escasos y ruidosos, más una ecuación confiable, más coeficientes desconocidos: este es el régimen donde la maquinaria adicional se paga sola.

4.6 Romper la PINN: desbalance de los pesos de la pérdida¶

Hasta ahora la pérdida compuesta se ha portado bien porque sus tres términos resultaron estar en escalas comparables; los pesos iguales simplemente funcionaron. Eso es suerte, no una ley, y la advertencia de la sección 5 sobre el equilibrio de la pérdida merece el mismo trato que recibieron los entrenamientos rotos del cuaderno 4.5: mostrar la falla, leer los síntomas, nombrar el arreglo.

El ejercicio de la sección 4.3 quitó el término físico y observó cómo el modelo degeneraba en una regresión ruidosa. Aquí rompemos el equilibrio en la dirección contraria: el término físico ahoga a los datos y a la condición inicial. La corrida de abajo repite la sección 4.2 con un solo cambio, w_pde = 1e4 sobre el residuo de la EDP. Observe el total impreso: sigue cayendo, como si el entrenamiento fuera bien.

model_bad = HeatPINN().to(DEVICE)
optimizer_bad = optim.Adam(model_bad.parameters(), lr=2e-3)
w_pde = 1e4                       # the only change from Section 4.2

history_bad = {"data": [], "pde": [], "icbc": []}
for step in range(n_steps):
    optimizer_bad.zero_grad()
    loss_data = mse(model_bad(X_data), y_data)
    T_c = model_bad(X_col)
    g = grad(T_c, X_col)[0]
    T_x, T_t = g[:, 0:1], g[:, 1:2]
    T_xx = grad(T_x, X_col)[0][:, 0:1]
    loss_pde = torch.mean((T_t - T_xx) ** 2)
    loss_icbc = (mse(model_bad(X_ic), y_ic)
                 + torch.mean(model_bad(X_b0) ** 2)
                 + torch.mean(model_bad(X_b1) ** 2))
    total = loss_data + w_pde * loss_pde + loss_icbc
    total.backward()
    optimizer_bad.step()
    for key, val in zip(("data", "pde", "icbc"), (loss_data, loss_pde, loss_icbc)):
        history_bad[key].append(val.item())
    if step % 400 == 0 or step == n_steps - 1:
        print(f"step {step:4d}  total {total.item():8.3f}  data {loss_data.item():.3f}  "
              f"pde {loss_pde.item():.5f}  ic/bc {loss_icbc.item():.3f}")

step    0  total   96.174  data 2.569  pde 0.00819  ic/bc 11.750

step  400  total    9.546  data 1.123  pde 0.00000  ic/bc 8.419

step  800  total    9.476  data 1.111  pde 0.00000  ic/bc 8.361

step 1200  total    9.482  data 1.114  pde 0.00000  ic/bc 8.365

step 1600  total    9.499  data 1.118  pde 0.00000  ic/bc 8.378

step 1999  total    9.215  data 1.066  pde 0.00000  ic/bc 8.134

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

axes[0].plot(history_bad["data"], label='data loss')
axes[0].plot(history_bad["pde"], label='PDE residual loss')
axes[0].plot(history_bad["icbc"], label='IC/BC loss')
axes[0].set_yscale('log')
axes[0].set_xlabel('Adam iteration')
axes[0].set_ylabel('Loss (unweighted)')
axes[0].set_title(f'Loss components, w_pde = {w_pde:.0e}')
axes[0].legend()

X_plot0 = torch.tensor(np.stack([x_grid, np.zeros_like(x_grid)], axis=1),
                       dtype=torch.float32)
axes[1].plot(x_grid * L, analytic_T(x_grid, 0.0), 'k-', label='analytic, t = 0')
axes[1].plot(x_grid * L, model_heat(X_plot0).detach().numpy().ravel(), '--',
             label='balanced PINN (Sec. 4.2)')
axes[1].plot(x_grid * L, model_bad(X_plot0).detach().numpy().ravel(), '--',
             label=f'w_pde = {w_pde:.0e}')
axes[1].set_xlabel('Depth x (m)')
axes[1].set_ylabel('Temperature anomaly (C)')
axes[1].set_title('Initial profile: balanced vs. imbalanced')
axes[1].legend()
plt.tight_layout()
plt.show()

err_bad = np.abs(model_bad(X_plot0).detach().numpy().ravel()
                 - analytic_T(x_grid, 0.0)).max()
print(f"max error at t_hat = 0: broken run {err_bad:.2f} C, "
      f"balanced run {pinn_max_err[0.0]:.3f} C  (signal amplitude {T0_heat} C)")

max error at t_hat = 0: broken run 3.85 C, balanced run 0.029 C  (signal amplitude 5.0 C)

La pérdida total cae y la respuesta es basura. Esa combinación es la firma del desbalance de la pérdida, y resulta invisible a menos que usted registre las componentes por separado. Léalas a ellas, no al total: después de unos cientos de pasos el residuo de la EDP queda fijado por debajo de 10^-5, mientras que la pérdida de datos se sitúa cerca de 1 (muy por encima del piso de ruido de 0.04 de la sección 4.2) y la pérdida de condiciones iniciales y de contorno se estanca cerca de 8, del mismo orden que el 12.5 que obtendría una red que produjera cero en todas partes (la media de $(5\sin\pi\hat{x})^2$ ). La red encontró una solución casi trivial: $T \approx 0$ satisface exactamente la ecuación de calor y las condiciones de contorno, así que el objetivo ponderado premia quedarse cerca de ella, y cualquier movimiento hacia el perfil inicial verdadero se castiga a través del término dominante antes de que pueda rendir frutos. El perfil en $\hat{t}=0$ se queda unos 4 °C corto respecto de la anomalía de 5 °C.

El arreglo es el equilibrio que ya teníamos: con pesos iguales, la misma arquitectura y los mismos 2000 pasos llegaron a 0.03 °C en la sección 4.2. Cuando ninguna ponderación fija funciona, porque los términos tienen escalas genuinamente distintas, las jugadas estándar son normalizar cada término por su valor inicial, ajustar los pesos contra una métrica de validación reservada, o usar un esquema de ponderación adaptativa (Karniadakis et al., 2021); entrenar más tiempo y ensanchar las capas solo ayuda cuando el desbalance es leve. El hábito que hay que llevarse: grafique siempre las componentes de la pérdida por separado. Una sola curva de pérdida total esconde exactamente esta falla.

Ejercicio (sesgo espectral). El otro modo de falla que nombra la sección 5 es el sesgo espectral. Agregue un modo de frecuencia más alta a la condición inicial:

T(\hat{x}, 0) = 5 \sin(\pi \hat{x}) + 2 \sin(7\pi \hat{x}),

(8)

cuya solución analítica suma $2\sin(7\pi\hat{x})\,e^{-49\pi^2 \hat{t}}$ a la referencia de un solo modo. Regenere y_ic y las 40 muestras de datos a partir de la solución de dos modos, reentrene la PINN de la sección 4.2 (con pesos iguales) y vuelva a correr el solucionador FTCS de la sección 4.4 con el nuevo perfil inicial. Compare ambos con la solución analítica en $\hat{t} = 0$ .

Solución

El solucionador FTCS maneja el nuevo modo sin cambiar una línea de código: con nx = 101 el séptimo modo tiene unos 28 puntos de malla por longitud de onda. La PINN no. Después de los mismos 2000 pasos, su perfil inicial es el primer modo suave con las ondulaciones aplanadas: la pérdida de la condición inicial se estanca cerca de 2 (el valor cuadrático medio del modo faltante, 2^2/2), y el error máximo en t_hat = 0 es de unos 2 grados, la amplitud completa del modo que la red se negó a aprender. Eso es el sesgo espectral: los MLP con tanh ajustan primero las frecuencias bajas y las frecuencias altas lentamente, si es que las ajustan. Remedios, en orden creciente de esfuerzo: entrenar mucho más tiempo, ensanchar las capas o darle a la red entradas de características de Fourier (incrustaciones sin/cos de x_hat), la cura estándar. Para la difusión lo que está en juego es poco: el séptimo modo decae como e^(-49 pi^2 t_hat) y físicamente desaparece en días. Para las ecuaciones de onda no hay tal clemencia, porque las frecuencias altas son la señal.

5. Dónde están las PINN en 2026¶

Las PINN ya no son una novedad, y sus modos de falla están bien documentados; dos de ellos están ahora mismo en este cuaderno. El equilibrio de la pérdida es un problema de ajuste persistente: un mal equilibrio hace que el entrenamiento converja a una solución que satisface un término e ignora los demás, y la sección 4.6 lo mostró en nuestro propio problema, donde un factor de 10⁴ sobre el término físico produjo una red atorada cerca de la solución trivial cero, que se quedó unos 4 °C corta respecto de la anomalía de 5 °C mientras su pérdida total seguía cayendo. Los MLP estándar también arrastran un sesgo espectral hacia funciones suaves y de baja frecuencia: aprenden rápido las soluciones suaves y lentamente la estructura oscilatoria o de escala fina, si es que la aprenden (Karniadakis et al., 2021); el ejercicio de 4.6 hace que un MLP con tanh aplane una ondulación de séptimo modo que el solucionador de diferencias finitas resuelve sin cambio alguno. Las EDP rígidas, los frentes abruptos y las escalas temporales muy separadas agravan ambos problemas, porque entonces el residuo varía en órdenes de magnitud a lo largo del dominio y ninguna ponderación fija es correcta en todas partes.

La primera regla de uso es negativa: no resuelva un problema directo limpio con una PINN. La sección 4.4 le puso números a esto. Quince líneas de FTCS le ganaron a la PINN entrenada por un factor de aproximadamente 750 en tiempo de reloj y por dos órdenes de magnitud en exactitud, sin ningún hiperparámetro que ajustar. Cuando los coeficientes y las condiciones de contorno se conocen y existe una discretización estándar, gana el solucionador clásico. Y para el problema vecino de construir sustitutos rápidos que se evalúan muchas veces, los operadores neuronales se han quedado en gran medida con el terreno: el Fourier Neural Operator (Li et al., 2021) y DeepONet (Lu et al., 2021) aprenden el operador de solución, el mapa que va de una condición inicial, una condición de contorno o un campo de coeficientes a la solución, de modo que, tras entrenar sobre muchos pares de simulaciones, un caso nuevo cuesta milisegundos, allí donde una PINN tendría que reentrenar desde cero. Para estudios paramétricos, cuantificación de incertidumbre y pronóstico operativo, esa es la herramienta que se usa.

Las PINN todavía se ganan su lugar en dos regímenes, y ambos quedan ahora demostrados en vez de afirmados. El primero son los problemas inversos: la sección 4.5 convirtió la difusividad en un parámetro entrenable y la recuperó con un error menor al 2 % a partir de las 40 muestras ruidosas, con la EDP actuando como modelo directo dentro de una regresión (Raissi et al., 2019). Esto es inversión conjunta, el patrón de PEST-alrededor-de-MODFLOW de la calibración hidrogeológica comprimido en un solo bucle de entrenamiento. El segundo es el régimen de datos escasos que representan esos mismos 40 puntos: observaciones dispersas y ruidosas más una ecuación de gobierno confiable y ninguna biblioteca de simulaciones sobre la cual entrenar un operador. Un solucionador directo no tiene ranura para tales datos; la PINN los trata como un término más de la pérdida. Ambas situaciones son comunes en geociencias, y por eso el método sigue siendo parte de la caja de herramientas (Karniadakis et al., 2021).

Referencias¶

Raissi, M., Perdikaris, P., and Karniadakis, G. E. (2019). Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378, 686-707.
Karniadakis, G. E., Kevrekidis, I. G., Lu, L., Perdikaris, P., Wang, S., and Yang, L. (2021). Physics-informed machine learning. Nature Reviews Physics, 3, 422-440.
Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A. (2021). Fourier neural operator for parametric partial differential equations. International Conference on Learning Representations (ICLR).
Lu, L., Jin, P., Pang, G., Zhang, Z., and Karniadakis, G. E. (2021). Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3, 218-229.

Resumen¶

Una PINN agrega a la pérdida de entrenamiento el residuo de una ecuación de gobierno, evaluado por diferenciación automática en puntos de colocación.
En la ablación de la ley de enfriamiento, solo el término físico produjo una extrapolación correcta más allá de los datos; la regularización L2 no.
En la ecuación de calor 1D, un MLP de 3 capas con tanh recuperó el modo difusivo analítico a partir de 40 muestras ruidosas más la EDP y las condiciones de contorno.
El modelo de referencia clásico ganó el problema directo de manera rotunda: 15 líneas de FTCS fueron unas 750 veces más rápidas y dos órdenes de magnitud más exactas que la PINN (sección 4.4).
Convertir la difusividad en un parámetro entrenable transformó el mismo bucle en una inversión conjunta que recuperó $\hat{k}$ con un error menor al 2 % a partir de las 40 muestras ruidosas, una calibración que el solucionador directo no puede hacer solo (sección 4.5).
Un peso de 10⁴ sobre la pérdida física dio una pérdida total decreciente y una respuesta equivocada; la falla solo es visible en las curvas de pérdida por término, así que grafique siempre las componentes (sección 4.6).
En 2026: solucionadores clásicos para problemas directos limpios, operadores neuronales (FNO, DeepONet) para sustitutos paramétricos, PINN para problemas inversos y de datos escasos.