Optimización bayesiana: surrogate models para CFD costoso

El problema del CFD costoso

Una simulación CFD típica tarda entre 20 minutos y varias horas. Los algoritmos genéticos necesitan miles de evaluaciones para converger. Si cada evaluación cuesta 1 hora y tu algoritmo necesita 1000 evaluaciones, tu optimización tarda 42 días. No es viable para plazos de proyecto reales.

La optimización bayesiana resuelve esto: construye un modelo sustituto (surrogate) de la función objetivo a partir de las evaluaciones ya realizadas, y usa ese modelo para decidir inteligentemente qué punto evaluar a continuación. Encontrarás diseños competitivos con 50-100 evaluaciones en lugar de miles.

Qué es un Gaussian Process

Un Gaussian Process (GP) es una distribución de probabilidad sobre funciones, no sobre vectores. Se define completamente mediante:

• Función media m(x): tendencia central de la función (típicamente m(x) = 0 después de estandarizar las salidas).
• Kernel de covarianza k(x, x'): cómo de correladas están las salidas de dos puntos x y x'. Codifica el supuesto de suavidad.

La propiedad fundamental: cualquier conjunto finito de puntos evaluados tiene distribución conjunta Gaussiana con media y covarianza definidas por el kernel. Esto permite hacer predicciones con incertidumbre.

Kernels comunes en ingeniería:

• RBF / Squared Exponential: k(x, x') = sigma^2 * exp(-||x - x'||^2 / (2 * l^2)). Supone funciones infinitamente diferenciables — a menudo demasiado suave para funciones de ingeniería.
• Matérn 5/2: asume solo 2 derivadas continuas, más realista para superficies de respuesta en CFD y FEA. Es la recomendación por defecto.

import numpy as np
from scipy.linalg import cholesky, solve_triangular

def matern52_kernel(X1, X2, length_scale=1.0, sigma=1.0):
    """Kernel Matern 5/2 con un solo length scale."""
    dist = np.sqrt(np.sum((X1[:, None] - X2[None, :]) ** 2, axis=-1))
    sqrt5 = np.sqrt(5.0)
    r = sqrt5 * dist / length_scale
    return sigma**2 * (1.0 + r + r**2 / 3.0) * np.exp(-r)

def gp_predict(X_train, y_train, X_test, kernel, noise=1e-6):
    """Prediccion GP con media y varianza posterior."""
    K = kernel(X_train, X_train)
    K += noise * np.eye(len(X_train))
    K_s = kernel(X_train, X_test)
    K_ss = kernel(X_test, X_test)
    
    L = cholesky(K, lower=True)
    alpha = solve_triangular(L.T, solve_triangular(L, y_train, lower=True))
    
    mu = K_s.T @ alpha
    v = solve_triangular(L, K_s, lower=True)
    sigma = np.diag(K_ss) - np.sum(v**2, axis=0)
    return mu, np.maximum(sigma, 0)

Función de adquisición: decidiendo qué evaluar

El GP entrenado nos da para cada punto x una predicción con media mu(x) y desviación estándar sigma(x). Pero no nos dice qué punto evaluar a continuación. La función de adquisición (acquisition function) cuantifica numéricamente lo prometedor que es un punto, equilibrando exploración (zonas con alta incertidumbre) y explotación (zonas con baja media).

Expected Improvement (EI): La más usada. Calcula la mejora esperada sobre el mejor valor observado f_best:

EI(x) = E[max(f_best - f(x), 0)]

Fórmula analítica: EI(x) = (f_best - mu(x)) * Phi(Z) + sigma(x) * phi(Z), donde Z = (f_best - mu(x)) / sigma(x), Phi es la CDF normal y phi es la PDF normal.

Upper Confidence Bound (UCB): Más simple: UCB(x) = mu(x) - kappa * sigma(x) en minimización. El parámetro kappa controla el balance exploración-explotación.

Pipeline completo de optimización bayesiana

1. DoE inicial: Generar 10-20 puntos con Latin Hypercube Sampling para cubrir el espacio de diseño de forma eficiente.
2. Evaluar CFD en todos los puntos iniciales. Añadir al dataset D.
3. Entrenar GP con todos los datos en D. Ajustar hiperparámetros del kernel (length scales, sigma) por máxima verosimilitud marginal.
4. Optimizar función de adquisición: Encontrar el punto x_next que maximiza EI(x) (o minimiza UCB). Se usa L-BFGS o búsqueda aleatoria multistart.
5. Evaluar CFD en x_next y añadir al dataset D.
6. Repetir pasos 3-5 hasta agotar presupuesto de evaluaciones o alcanzar convergencia.
7. Devolver el mejor punto evaluado (no el mínimo del GP — siempre validar con simulación real).

Caso práctico: optimización de un difusor 2D con CFD

Un difusor es un conducto que se expande, reduciendo velocidad y recuperando presión. Optimizamos uno con 5 parámetros geométricos: ancho de entrada, ancho de salida, longitud, ángulo de expansión y curvatura de pared.

Dos objetivos (convertidos a single-objective vía restricciones): maximizar el coeficiente de recuperación de presión Cp, sujeto a que las pérdidas de presión total no excedan un umbral. Cada evaluación: OpenFOAM simpleFoam, malla de ~200k celdas, ~20 minutos en 8 núcleos.

• Latin Hypercube: 15 puntos iniciales
• GP con kernel Matérn 5/2
• Adquisición: Expected Improvement
• 60 evaluaciones totales (15 iniciales + 45 iteraciones bayesianas)

Resultado: calidad del óptimo equivalente a NSGA-II con 500 evaluaciones. Una reducción del 88% en tiempo de cómputo.

Ventajas sobre algoritmos genéticos

• Eficiencia muestral: 10-100x menos evaluaciones. Ideal cuando cada evaluación es costosa.
• Cuantificación de incertidumbre: Los GPs proporcionan intervalos de confianza de forma natural.
• Sin parámetros de población: No hay que ajustar tamaño de población ni operadores genéticos.
• Ruido tolerante: Funciona con observaciones ruidosas, como datos experimentales (no solo simulación).

Desventajas: limitado a ~20 variables de diseño (los GPs escalan mal con dimensionalidad) y naturaleza secuencial (más difícil de paralelizar).

Optimización multi-objetivo con modelos Bayesianos

• ParEGO: Escalariza múltiples objetivos con pesos aleatorios en cada iteración, luego optimiza con BO mono-objetivo estándar.
• Expected Hypervolume Improvement (EHVI): Función de adquisición directa multi-objetivo. Evalúa cuánto hipervolumen añadiría un nuevo punto al frente de Pareto. Caro computacionalmente pero el gold standard.
• Multi-task GPs: Comparten información entre objetivos relacionados mediante kernels multitarea, aprovechando correlaciones entre objetivos.

Recomendaciones prácticas

• Kernel Matérn 5/2 como opción por defecto para ingeniería (suavidad realista para superficies de respuesta físicas).
• Normalizar entradas a [0, 1] y estandarizar salidas a media cero y varianza unitaria.
• Multistart al optimizar la función de adquisición: varios reinicios aleatorios evitan óptimos locales del propio surrogate.
• Validación periódica del GP: separar un conjunto de validación para comprobar que el surrogate no se aleja de la realidad.
• Nugget de 1e-6 en la diagonal del kernel si el GP está mal calibrado (incertidumbre subestimada).
• Paralelización: En lugar de una sola recomendación, generar k candidatos por iteración (q-EI, fantasizing) y lanzar varias simulaciones en paralelo.