13  Métricas del modelo

Conforme se ha referido previamente, el desarrollo del modelo predictivo se realizó utilizando el framework DeepXDE (versión 1.10.1) con backend de TensorFlow 1.x (configurado mediante tensorflow.compat.v1). Para asegurar reproducibilidad, se fijó la semilla aleatoria en 123 a nivel de DeepXDE, TensorFlow y NumPy. La red neuronal se implementó como un DeepONetCartesianProd con la siguiente estructura especializada:

  • Rama (branch)
    • Capa de entrada: (num_sensors + 1)² = 49 neuronas.
    • 3 capas ocultas de 20 neuronas cada una.
  • Tronco (trunk)
    • Capa de entrada: 3 neuronas (coordenadas espaciotemporales x, y, t).
    • Misma configuración de capas ocultas que la rama.
  • Hiperparámetros clave
    • Función de activación: ELU (Exponential Linear Unit).
    • Inicialización de pesos: Glorot normal.
    • Optimizador: ADAM con tasa de aprendizaje inicial de 2×10⁻³ y decaimiento exponencial (decay_rate=0.05 cada 500 pasos).
Código 
import deepxde as dde
import numpy as np
import tensorflow as tf

# -------------------------------------------------------------------------
# Constantes y Parametros
# -------------------------------------------------------------------------

# Backend y semilla
dde.backend.set_default_backend("tensorflow.compat.v1")
dde.config.set_random_seed(123)

# Parametros fisicos
p = 1050
c = 3639
keff = 5
final_time = 1800
L0 = 0.05
cb = 3825
Q = 0
TM = 45
Ta = 37
alpha = p * c / keff

# Coeficientes adimencionales
a1 = final_time / (alpha * L0**2)
a2 = final_time * cb / (p * c)
a3 = (final_time * Q) / (p * c * (TM - Ta))

# Dominio de las fronteras
x_initial, x_boundary = 0.0, 1.0
y_initial, y_boundary = 0.0, 1.0
t_initial, t_final = 0.0, 1.0

# Configuracion del numero de datos
pts_dom = 45
pts_bc = 30
pts_ic = 20
num_test = 25

# Malla de sensores y espacio de funciones
num_sensors = 6
num_function = 25
size_cov_matrix = 50

# Arquitectura de la red
width_net = 20
len_net = 3
AF = "elu"
k_initializer = "Glorot normal"


# Parámetros de entrenamiento
num_iterations = 20000
learning_rate = 2e-3
decay_rate = 0.05
decay_steps = 500

# -------------------------------------------------------------------------
# Dominio espacial y temporal
# -------------------------------------------------------------------------

spatial_domain = dde.geometry.Rectangle([x_initial, y_initial],
                                        [x_boundary, y_boundary])
time_domain = dde.geometry.TimeDomain(t_initial, t_final)
geomtime = dde.geometry.GeometryXTime(spatial_domain, time_domain)

# -------------------------------------------------------------------------
# EDP, CI y CF 
# -------------------------------------------------------------------------

def initial_condition(X):
    X = np.asarray(X)
    if X.ndim == 1:
        return 0.0
    return np.zeros((X.shape[0], 1))

def heat_equation(x, u, coords):
    u_t = dde.grad.jacobian(u, x, i=0, j=2)
    u_xx = dde.grad.hessian(u, x, i=0, j=0)
    u_yy = dde.grad.hessian(u, x, i=1, j=1)
    return u_t - a1*(u_xx + u_yy) + a2*u

def zero_value(X):
    return 0

def time_value(X):
    return X[:, 2]

def is_on_vertex(x):
    vertices = np.array([[x_initial, y_initial],
                         [x_boundary, y_initial],
                         [x_initial, y_boundary],
                         [x_boundary, y_boundary]])
    return any(np.allclose(x, v) for v in vertices)

def is_initial(X, on_initial):
    return on_initial and np.isclose(X[2], t_initial)

def left_boundary(X, on_boundary):
    spatial = X[0:2]
    t = X[2]
    return (
        on_boundary 
        and np.isclose(spatial[0], x_initial) 
        and not np.isclose(t, t_initial) 
        and not is_on_vertex(spatial)
    )

def right_boundary(X, on_boundary):
    spatial = X[0:2]
    t = X[2]
    return (
        on_boundary 
        and np.isclose(spatial[0], x_boundary) 
        and not np.isclose(t, t_initial) 
        and not is_on_vertex(spatial)
    )

def up_low_boundary(X, on_boundary):
    spatial = X[0:2]
    t = X[2]
    return (on_boundary 
    and (np.isclose(spatial[1], y_initial) 
    or np.isclose(spatial[1], y_boundary)) 
    and not np.isclose(t, t_initial) 
    and not is_on_vertex(spatial)
    )

# Condiciones iniciales y de frontera
ic = dde.icbc.IC(geomtime, initial_condition, is_initial)
left_bc = dde.icbc.DirichletBC(geomtime, 
                                zero_value, left_boundary)
right_bc = dde.icbc.NeumannBC(geomtime,
                                time_value, right_boundary)
up_low_bc = dde.icbc.NeumannBC(geomtime, 
                                zero_value, up_low_boundary)

# -------------------------------------------------------------------------
# Construccion de los datos
# -------------------------------------------------------------------------

pde_data = dde.data.TimePDE(
    geomtime,
    heat_equation,
    [ic, left_bc, right_bc, up_low_bc],
    num_domain=pts_dom,
    num_boundary=pts_bc,
    num_initial=pts_ic 
)

# -------------------------------------------------------------------------
# Sensores y espacio de funciones
# -------------------------------------------------------------------------

side = np.linspace(x_initial, x_boundary, num_sensors + 1)
x, y = np.meshgrid(side, side, indexing='xy')
sensor_pts = np.stack([x.ravel(), y.ravel()], axis=1)

fs = dde.data.function_spaces.GRF2D(N=size_cov_matrix, 
                                    interp="linear")

data = dde.data.PDEOperatorCartesianProd(
    pde_data,
    fs,
    sensor_pts,
    num_function=num_function,
    function_variables=[0, 1],
    num_test=num_test
)

# -------------------------------------------------------------------------
# Definicion de la red
# -------------------------------------------------------------------------

branch_layers = [(num_sensors + 1)**2] + len_net * [width_net]
trunk_layers = [3] + len_net * [width_net]

net = dde.nn.DeepONetCartesianProd(
    branch_layers,
    trunk_layers,
    activation=AF,
    kernel_initializer=k_initializer
)

# -------------------------------------------------------------------------
# Compilacion y entrenamiento del modelo
# -------------------------------------------------------------------------

model = dde.Model(data, net)
model.compile("adam", lr=learning_rate,
            decay=("inverse time", decay_steps, decay_rate))
losshistory, train_state = model.train(iterations=num_iterations,
display_every=decay_steps)

# Refinamiento con el optimizador L-BFGS
model.compile("L-BFGS")
losshistory, train_state = model.train()

13.1 Gráficas de pérdida del modelo

El proceso de entrenamiento del modelo se monitoreó mediante el seguimiento detallado de cinco componentes de pérdida, cada una asociada a restricciones físicas y matemáticas específicas del problema:

  1. Pérdida residual de la EDP
    • Función: Mide el cumplimiento de la ecuación de Bio-Calor en el dominio interior.
    • Importancia: Garantiza que la solución aprendida satisfaga la física subyacente.
    • Comportamiento esperado: Debe converger a valores cercanos a cero (típicamente < 1e-3).
  2. Pérdida de condición inicial
    • Función: Controla la precisión en t=0.
    • Importancia: Asegura coherencia con el estado inicial del sistema.
    • Patrón típico: Suele ser la primera en converger por su carácter puntual.
  3. Pérdida de frontera izquierda (Dirichlet)
    • Función: Evalúa el cumplimiento de condiciones de valor prescrito.
    • Relevancia: Mantiene valores fijos en bordes específicos.
    • Convergencia: Normalmente rápida por ser restrictiva.
  4. Pérdida de frontera derecha (Neumann)
    • Función: Verifica gradientes normales en esta frontera
    • Dificultad característica: Puede mostrar oscilaciones iniciales
  5. Pérdida de fronteras superior/inferior (Neumann)
    • Función: Controla condiciones de flujo en estos bordes
    • Complejidad: En problemas 2D/3D suele ser la última en estabilizarse

13.1.1 Perdida para el conjunto de entrenamiento

La tendencia decreciente de la pérdida durante el entrenamiento evidencia que el modelo neuronal está aprendiendo progresivamente a satisfacer las restricciones físicas impuestas por la ecuación de Bio-Calor y sus condiciones asociadas. Este comportamiento indica una correcta adaptación de los parámetros de la red, permitiendo reducir consistentemente el error en las distintas componentes de la función objetivo y aproximando con mayor precisión la dinámica térmica en el dominio.

Código 
import plotly.graph_objects as go

# Nombres de las componentes del loss
loss_labels = [
    "Pérdida residual PDE",
    "Pérdida de condición inicial",
    "Pérdida de frontera izquierda (Dirichlet)",
    "Pérdida de frontera derecha (Neumann)",
    "Pérdida de fronteras superior/inferior (Neumann)"
]

# Extraer pasos y pérdida de entrenamiento
steps = losshistory.steps
train_loss = np.array(losshistory.loss_train)

# Crear figura
fig_train = go.Figure()

for i in range(train_loss.shape[1]):
    fig_train.add_trace(go.Scatter(
        x=steps,
        y=train_loss[:, i],
        mode='lines',
        name=loss_labels[i]
    ))

fig_train.update_layout(
    title="Historial de pérdida en el entrenamiento",
    xaxis=dict(title="Iteración", tickformat=".1e"),
    yaxis=dict(title="Pérdida (log)", type="log", tickformat=".1e"),
    template="plotly_white",
    legend=dict(x=0.99, y=0.99),
    font=dict(size=14)
)
Figura 13.1: Historial de pérdida en el entrenamiento de la red neuronal.

13.1.2 Pérdida para el conjunto de prueba

La disminución del error en el conjunto de prueba confirma que el modelo no solo memoriza los datos de entrenamiento, sino que logra generalizar a situaciones no vistas. Esto constituye un resultado favorable, pues asegura que la red neuronal mantiene su capacidad predictiva fuera de los escenarios empleados para el ajuste, garantizando robustez y confiabilidad en aplicaciones biomédicas donde la precisión en la estimación térmica es fundamental.

Código 
import plotly.graph_objects as go

# Nombres de las componentes del loss
loss_labels = [
    "Pérdida residual PDE",
    "Pérdida de condición inicial",
    "Pérdida de frontera izquierda (Dirichlet)",
    "Pérdida de frontera derecha (Neumann)",
    "Pérdida de fronteras superior/inferior (Neumann)"
]

# Extraer pasos y pérdida de entrenamiento
steps = losshistory.steps
test_loss = np.array(losshistory.loss_test)

# Crear figura
fig_test = go.Figure()

for i in range(test_loss.shape[1]):
    fig_test.add_trace(go.Scatter(
        x=steps,
        y=test_loss[:, i],
        mode='lines',
        name=loss_labels[i]
    ))

fig_test.update_layout(
    title="Historial de pérdida en el conjunto de prueba",
    xaxis=dict(title="Iteración", tickformat=".1e"),
    yaxis=dict(title="Pérdida (log)", type="log", tickformat=".1e"),
    template="plotly_white",
    legend=dict(x=0.99, y=0.99),
    font=dict(size=14)
)
Figura 13.2: Historial de pérdida en el conjunto de prueba de la red neuronal.

13.2 Guardado de datos

Para permitir la comparación cuantitativa con el método de Crank-Nicolson y facilitar la generación de visualizaciones consistentes, se exportaron las predicciones del modelo neuronal en formato CSV. El proceso consistió en:

  1. Generación de la malla de evaluación:
    • Dominio espacial: Cuadrado unitario [0,1] × [0,1].
    • Discretización: 26 segmentos equiespaciados en cada eje (x, y).
    • Puntos totales: 676 (26 × 26).
    • Tiempos evaluados: t = [0.0, 0.25, 0.50, 0.75, 1.0].
  2. Estructura del archivo:
    • Coordenadas espacio-temporales (t, x, y) para cada punto de la grilla 26×26.
    • Valores de la solución en los tiempos de interés.
Listado 13.1: Guardado de los datos de la red neuronal.
Código 
import pandas as pd
# Lista de tiempos
times = [0.0, 0.25, 0.5, 0.75, 1.0]

# Crear la malla (x, y)
num_points = 26
x = np.linspace(0, 1, num_points)
y = np.linspace(0, 1, num_points)
X, Y = np.meshgrid(x, y)

# Lista para almacenar resultados
results = []

for t_val in times:
    # Crear entrada trunk: (num_points^2, 3)
    points = np.vstack((X.flatten(), Y.flatten(),
                         t_val * np.ones_like(X.flatten()))).T

    # Crear entrada branch: condición inicial constante cero
    branch_input = np.zeros((1, sensor_pts.shape[0]))

    # Predecir
    predicted = model.predict((branch_input, points)).flatten()

    # Agregar los datos al resultado
    for xi, yi, thetai in zip(points[:, 0], points[:, 1], predicted):
        results.append([t_val, xi, yi, thetai])

# Crear el DataFrame
df = pd.DataFrame(results, columns=["time", "X", "Y", "Theta"])

# Obtener la ruta del script actual y guardar el archivo CSV
ruta = r"data/model_DoN.csv"
df.to_csv(ruta, index=False)

13.2.1 Visualización interactiva

Para facilitar el análisis exploratorio de las predicciones generadas por el modelo DeepONet, se implementó una tabla dinámica interactiva mediante la librería itables (versión 1.5.2), que extiende las funcionalidades de Pandas para su visualización (Solo en versión en línea).

Código 
from itables import show, options

options.maxBytes = 0

# Mostrar la tabla interactiva
show(df, paging=True,
        ordering=True,
        searching = False,
        scrollY="350px",
        buttons=["pageLength", "csv", "excel"],
        lengthMenu=[10, 25, 50, 100],
        classes="display nowrap cell-border"
        )
Tabla 13.1: Predicciones de la red neuronal.
time X Y Theta
Loading ITables v2.3.0 from the internet... (need help?)