Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

No9 - PINNs Pt2

Physics-Informed Neural Networks (PINNs)

Fecha: 11/05/2026

Tres tipos de algoritmos para incorporar ecuaciones diferenciales

Cuando queremos estimar o aproximar una ecuación diferencial con modelos de aprendizaje automático, hay tres enfoques principales:

  1. Emuladores: se estima la ecuación diferencial con otro modelo de ML. Por ejemplo, cuando hay ecuaciones costosas (Navier-Stokes, etc.)

  2. Restricciones suaves (vía Lagrangiano): el modelo se entrena minimizando una función de costo empírica Lemp\mathcal{L}_{\text{emp}} más un término de penalización D[x(θ)]D[x(\theta)] que incentiva satisfacer la ecuación diferencial. Por ejemplo, las PINNs:

min⁡θ, xLemp(y,x,θ)+λ∥D[x(θ)]∥.\min_{\theta,\, x} \mathcal{L}_{\text{emp}}(y, x, \theta) + \lambda \| D[x(\theta)] \|.
(1)

  1. Restricciones fuertes (se satisfacen estrictamente): el modelo debe satisfacer la ecuación diferencial D[x(θ)]=0D[x(\theta)] = 0 en todo momento. Son las NODE y UDE:

min⁡θ, xLemp(y,x(θ),θ)s.a. D[x(θ)]=0.\min_{\theta,\, x} \mathcal{L}_{\text{emp}}(y, x(\theta), \theta) \quad \text{s.a. } D[x(\theta)] = 0.
(2)

Dualidad Lagrangiana

Un resultado fundamental de la optimización con restricciones establece que, si relajamos la restricción fuerte D[x(θ)]=0D[x(\theta)] = 0 por una restricción aproximada ∥D[x(θ)]∥≤ε\| D[x(\theta)] \| \leq \varepsilon, existe un λ(ε)\lambda(\varepsilon) tal que la solución del problema con restricción suave (penalizado con λ(ε)\lambda(\varepsilon)) coincide con la solución del problema con restricción fuerte relajada Boyd & Vandenberghe (2004):

λ(ε)=λ⟺restriccioˊn suave=restriccioˊn fuerte relajada\lambda(\varepsilon) = \lambda \quad \Longleftrightarrow \quad \text{restricción suave} = \text{restricción fuerte relajada}
(3)

Dificultades al entrenar PINNs

En la práctica, entrenar PINNs tiene problemas importantes de optimización. La función de pérdida puede tener múltiples términos: Lemp\mathcal{L}_{\text{emp}} (ajuste a los datos) y uno o mas términos de penalización λi∥Di[x(θ)]∥\lambda_i \|D_i[x(\theta)]\|, uno por cada restricción (condición inicial, condición de borde, residuo de la ecuación diferencial, etc.).

LPINN=Lemp+∑i=1nλi∥Di[x(θ)]∥\mathcal{L}_{\text{PINN}} = \mathcal{L}_{\text{emp}} + \sum_{i=1}^{n} \lambda_i \|D_i[x(\theta)]\|
(4)

Cuando λ\lambda es muy grande, el problema se vuelve mal condicionado: los gradientes del término de penalización dominan, entonces el optimizador tiene problemas para elegir la dirección óptima.

Curvas de nivel de \mathcal{L}_{\text{PINN}} (rosa) alrededor de la variedad donde D[x(\cdot,\theta)]=0 (negro). En un entorno de D[x(\cdot,\theta)]=0 es probable que la PINN esté mal condicionada.

Curvas de nivel de LPINN\mathcal{L}_{\text{PINN}} (rosa) alrededor de la variedad donde D[x(⋅,θ)]=0D[x(\cdot,\theta)]=0 (negro). En un entorno de D[x(⋅,θ)]=0D[x(\cdot,\theta)]=0 es probable que la PINN esté mal condicionada.

Selección de λ\lambda

En general, el vector de hiperparámetros λ=(λ1,…,λn)\lambda = (\lambda_1, \ldots, \lambda_n) controla cuánto pesa cada restricción. Hay tres estrategias para elegirlos:

Lempk≈λ1kLCI≈λ2kL∂≈⋯\mathcal{L}_{\text{emp}}^k \approx \lambda_1^k \mathcal{L}_{CI} \approx \lambda_2^k \mathcal{L}_{\partial} \approx \cdots
(6)
∥∇θLempk∥≈∥λ1k∇θLCI∥≈∥λ2k∇θL∂∥≈⋯\| \nabla_\theta \mathcal{L}_{\text{emp}}^k \| \approx \| \lambda_1^k \nabla_\theta \mathcal{L}_{CI} \| \approx \| \lambda_2^k \nabla_\theta \mathcal{L}_{\partial} \| \approx \cdots
(8)

Implementación

Problema directo

Vemos una implementación de PINN directo.

t0, T = tspan
scale_layer_inp = WrappedFunction(t -> (t .- t0) ./ (T - t0))
scale_layer_out = WrappedFunction(x -> 10.0 .* x)

Definimos la red neuronal:

nn = Chain(
    scale_layer_inp,
    Dense(1 => 5, tanh),
    Dense(5 => 10, tanh),
    Dense(10 => 10, tanh),
    Dense(10 => 2),
    scale_layer_out
)

Puntos de colocación: En una PINN hay que elegir puntos en el dominio temporal donde evaluar el residuo de la ecuación diferencial D[x(θ)]D[x(\theta)].

N_col = 100
# t_col = collect(range(tspan[1], tspan[2], length=N_col))
# t_col = 10.0.^collect(range(log10(tspan[1] + 1e-10), log10(tspan[2]), length = N_col))
t_col = [
    10.0.^collect(range(log10(tspan[1] + 1e-10), 
	log10(tspan[2]), length = div(N_col, 2)));
    collect(range(tspan[1], tspan[2], length=div(N_col, 2)))
]
Las cruces son los puntos de colocación: donde la ecuación diferencial se está evaluando.

Las cruces son los puntos de colocación: donde la ecuación diferencial se está evaluando.

Este es un ejemplo que claramente no convergió bien, probablemente el resultado mejoraría ejecutando el código por mas tiempo, pero.. ¿Qué mas podemos hacer?

Técnicas para mejorar la convergencia

Puntos de colocación: La ubicación de los puntos de colocación no es arbitraria. En este ejemplo, exploramos solo dos opciones:

Forzar la condición inicial exactamente: Vemos cómo implementamos la función de costo:

function pinn_loss_components(ps)
    u_ic    = nn([tspan[1]], ps, st)[1]
    ic_loss = sum((u_ic .- u0) .^ 2)

    phys_residuals = map(t_col) do t_i
        u    = nn([t_i], ps, st)[1]
        dudt = nn_dudt(t_i, ps, st)
        rhs  = lv_rhs(u)
        sum((dudt .- rhs) .^ 2)
    end
    phys_loss_weighted = λ_physics * mean(phys_residuals)

    return ic_loss, phys_loss_weighted
end

pinn_loss(ps) = sum(pinn_loss_components(ps))

Estamos definiendo a la condición inicial como una restricción suave. No es ideal. ¿Cómo podemos imponerla como restricción fuerte? Una reparametrización directa es:

u(t)=NNθ(t)(t−t0)+u0u(t) = \text{NN}_\theta(t)(t - t_0) + u_0
(9)

En t=t0t = t_0 se tiene u(t0)=NNθ(t0)⋅0+u0=u0u(t_0) = \text{NN}_\theta(t_0) \cdot 0 + u_0 = u_0 independientemente de los parámetros. Sin embargo, el factor (t−t0)(t - t_0) crece linealmente, obligando a la red a desaprender un trend lineal.

Una mejor alternativa es usar una función ϕ\phi acotada:

uθ(t)=u0+ϕ(t) NNθ(t)u_\theta(t) = u_0 + \phi(t)\, \text{NN}_\theta(t)
(10)

donde ϕ(t0)=0\phi(t_0) = 0. Por ejemplo, ϕ(t)=t−t0t−T0\phi(t) = \frac{t - t_0}{t - T_0} con T0T_0 cualquier número finito satisface ϕ(t0)=0\phi(t_0) = 0 y ϕ(t)→1\phi(t) \to 1 cuando t→∞t \to \infty.

Podemos aplicar esto al ejemplo:

function u_ansatz(t_val, ps, st)
    raw = nn([t_val], ps, st)[1]   # NN(t) ∈ ℝ²
    return u0 .+ t_val .* raw      # u₀ + t·NN(t)
end

# Derivada temporal del ansatz por diferencias finitas centrales
function u_ansatz_dudt(t_val, ps, st; h=1e-4)
    u_fwd = u_ansatz(t_val + h, ps, st)
    u_bwd = u_ansatz(t_val - h, ps, st)
    return (u_fwd .- u_bwd) ./ (2h)
end
function pinn_loss(ps)
    phys_residuals = map(t_col) do t_i
        u    = u_ansatz(t_i, ps, st)
        dudt = u_ansatz_dudt(t_i, ps, st)
        rhs  = lv_rhs(u)
        sum((dudt .- rhs) .^ 2)
    end
    return mean(phys_residuals)
end
Mejora... pero no es ideal.

Mejora... pero no es ideal.

Aún podría ser mejorable este resultado. Pero no es lo que nos interesa, no buscamos que una PINN sea excelente resolviendo una ecuación diferencial. Donde verdaderamente brillan las PINN es en el problema inverso.

Problema inverso

En el problema inverso no solo se conoce la condición inicial sino también observaciones de la trayectoria. El objetivo deja de ser únicamente ajustar u(t0)=u0u(t_0) = u_0 y pasa a ser ajustar la trayectoria completa. En este caso se aplica a Lotka-Volterra sin conocer los parámetros del sistema.

Generación de datos observados

Se resuelve la ODE con los parámetros verdaderos (desconocidos para la PINN) y se submuestran M=40M = 40 puntos uniformes:

const α_true = 1.0;   const β_true = 0.1
const δ_true = 0.075; const γ_true = 1.5
const u0    = [10.0, 5.0]
const tspan = (0.0, 15.0)

prob_true = ODEProblem(lotka_volterra!, u0, tspan, [α_true, β_true, δ_true, γ_true])
sol_true  = solve(prob_true, Tsit5(), saveat=0.01)

M_obs = 40
t_obs = collect(range(tspan[1], tspan[2], length=M_obs))
u_obs = hcat([sol_true(t) for t in t_obs]...)   # (2, M_obs)
Parámetros entrenables

Los parámetros de la ODE se inicializan lejos de los valores verdaderos y se agrupan junto con los pesos de la red en un único NamedTuple que Optimisers.jl puede diferenciar recursivamente:

p_ode_init = [0.5, 0.05, 0.04, 1.0]   # [α̂, β̂, δ̂, γ̂]

theta = (nn = nn_ps, ode = copy(p_ode_init))
Función de pérdida
const λ_phys = 1.0

function inverse_loss_components(theta)
    nn_ps = theta.nn
    α̂, β̂, δ̂, γ̂ = theta.ode

    data_residuals = map(1:M_obs) do j
        u_pred = nn([t_obs[j]], nn_ps, st)[1]
        sum((u_pred .- u_obs[:, j]) .^ 2)
    end
    data_loss = mean(data_residuals)

    phys_residuals = map(t_col) do t_i
        u    = nn([t_i], nn_ps, st)[1]
        dudt = nn_dudt(t_i, nn_ps, st)
        x, y = u[1], u[2]
        rhs  = [α̂*x - β̂*x*y,
                δ̂*x*y - γ̂*y]
        sum((dudt .- rhs) .^ 2)
    end
    phys_loss = λ_phys * mean(phys_residuals)

    return data_loss, phys_loss
end

inverse_loss(theta) = sum(inverse_loss_components(theta))
Resultados
Convergencia de los parámetros estimados \hat{p} hacia los valores verdaderos a lo largo del entrenamiento.

Convergencia de los parámetros estimados p^\hat{p} hacia los valores verdaderos a lo largo del entrenamiento.

El error puntual muestra que la PINN ajusta exactamente en los puntos de colocación y en los puntos de observación, pero el error crece en las regiones intermedias.

References
  1. Boyd, S., & Vandenberghe, L. (2004). Convex optimization. Cambridge university press.
  2. Lu, L., Pestourie, R., Yao, W., Wang, Z., Verdugo, F., & Johnson, S. G. (2021). Physics-Informed Neural Networks with Hard Constraints for Inverse Design. SIAM Journal on Scientific Computing, 43(6), B1105–B1132. 10.1137/21m1397908
Clases
No8 - PINNs
Clases
No10 - PINNs Pt2