Bi-Level optimization for interpolation-based parameter estimation of differential equations

Resumen Técnico: Optimización de Dos Niveles para la Estimación de Parámetros Basada en Interpolación de Ecuaciones Diferenciales

1. Planteamiento del Problema

La estimación de parámetros (o problema inverso) en ecuaciones diferenciales ordinarias (EDO) es fundamental en ingeniería química para ajustar modelos a datos experimentales. El enfoque tradicional, conocido como optimización secuencial o "single-shooting", integra numéricamente las EDO en cada iteración del optimizador. Sin embargo, este método presenta desafíos significativos:

• Costo Computacional: El cálculo de gradientes requiere la integración de ecuaciones de sensibilidad (derivadas de los estados respecto a los parámetros), lo cual es costoso.
• Convergencia: Es difícil de converger para trayectorias largas, estimaciones iniciales pobres o dinámicas altamente no lineales y complejas.
• Inestabilidad: En sistemas rígidos (stiff) o con caos, pequeñas variaciones en los parámetros pueden desviar la trayectoria simulada de la solución real, llevando a mínimos locales o fallos de convergencia.

Métodos alternativos como el "multiple-shooting" o la colocación ortogonal aumentan la dimensionalidad del problema o requieren suposiciones estrictas sobre las condiciones iniciales. Además, métodos previos basados en interpolación estaban limitados a parámetros que aparecen linealmente en las ecuaciones.

2. Metodología Propuesta

Los autores proponen un marco de optimización de dos niveles (bi-level) que combina interpolación de datos con técnicas de diferenciación automática y teorema de la función implícita.

• Formulación del Problema:
  Se divide el problema en dos niveles:

  1. Nivel Interno (Convexo): Se asume que los parámetros lineales ($p$) pueden estimarse mediante un problema de optimización convexa dado un valor fijo de los parámetros no lineales ($\phi$). Los estados del sistema ($X$) se aproximan mediante interpolación de los datos medidos ($\hat{X}$), eliminando la necesidad de integración numérica explícita en este paso.
  2. Nivel Externo (No Convexo): Se optimiza sobre los parámetros no lineales ($\phi$).

• Cálculo de Gradientes (Teorema de la Función Implícita):
  Para conectar los niveles, se derivan las sensibilidades de la solución óptima del nivel interno ($p^*$) respecto a los parámetros no lineales ($\phi$). Utilizando las condiciones de Karush-Kuhn-Tucker (KKT) del problema interno y el teorema de la función implícita, se calculan los gradientes del nivel externo sin necesidad de diferenciar a través de un integrador numérico inestable.

  • Se utiliza diferenciación automática (AD) en modo forward para calcular estas sensibilidades.
  • Se aprovecha la estructura de la matriz Hessiana del Lagrangiano (o matriz KKT) para acelerar los cálculos, reutilizando factorizaciones cuando es posible.

• Ventaja de la Interpolación:
  Al aproximar los estados mediante interpolación en lugar de integrar las EDO, se evitan las trayectorias inestables que suelen ocurrir en el método single-shooting cuando los parámetros iniciales son pobres. Esto estabiliza el cálculo de sensibilidades, especialmente en sistemas rígidos.

3. Contribuciones Clave

1. Generalización de Métodos Basados en Interpolación: Extiende el uso de interpolación más allá de los parámetros lineales, permitiendo su aplicación a una clase amplia de EDO donde los parámetros pueden aparecer de forma no lineal.
2. Marco Bi-Level Eficiente: Propone una estructura que resuelve un problema convexo interno (rápido y garantizado) y un problema no convexo externo, utilizando derivadas exactas obtenidas mediante el teorema de la función implícita.
3. Versatilidad: El método se demuestra aplicable a:
   • Estimación de parámetros en EDO estándar.
   • Descubrimiento de modelos (identificación de la estructura de la ecuación y sus parámetros simultáneamente).
   • Ecuaciones diferenciales con retardo (Delay Differential Equations).
   • Ecuaciones diferenciales rígidas (stiff).
   • Sistemas con estados parcialmente observados.
4. Implementación Abierta: El código está disponible públicamente, utilizando librerías modernas como JAX (para diferenciación automática), IPOPT (optimización no lineal) y SciPy.

4. Resultados Experimentales

El método se evaluó en siete casos de estudio de referencia, incluyendo dinámicas de iones de calcio, modelos bioquímicos no lineales, fermentación, modelos de enfermedades (Kermack-McKendrick), la reacción de Belousov (rígida) y un reactor de tanque agitado (CSTR) con estados no observados. También se aplicó al descubrimiento de modelos en la esterificación de ácidos carboxílicos.

• Comparativa: Se comparó contra Single-Shooting y Colocación Ortogonal.
• Convergencia: El método propuesto convergió a la solución óptima global en todos los casos de prueba. En contraste, Single-Shooting y la colocación ortogonal fallaron en converger o se estancaron en mínimos locales en varios problemas (especialmente en sistemas rígidos como la reacción de Belousov).
• Iteraciones vs. Tiempo:
  • El método propuesto requirió significativamente menos iteraciones que Single-Shooting debido a la convexidad del subproblema interno.
  • Sin embargo, el tiempo de ejecución total (wall time) fue mayor en el método propuesto. Esto se debe al costo computacional de resolver el problema interno convexo y calcular las derivadas a través de la solución óptima (inversión de matrices KKT) en cada iteración externa.
• Robustez: El método demostró ser robusto frente a estimaciones iniciales pobres y capaz de manejar sistemas rígidos sin necesidad de solvers de EDO rígidos durante la optimización, gracias a la interpolación de los datos medidos.

5. Significado y Conclusión

Este trabajo representa un avance significativo en la estimación de parámetros para sistemas dinámicos complejos. Al descomponer el problema en niveles y utilizar interpolación, logra superar las limitaciones de inestabilidad numérica y convergencia local que plagan a los métodos tradicionales de shooting.

Aunque el costo computacional por iteración es mayor, la capacidad de garantizar la convergencia a la solución global y manejar problemas difíciles (rígidos, con retardo, parcialmente observados) lo convierte en una herramienta valiosa, especialmente cuando los datos experimentales son abundantes y de alta frecuencia. El enfoque también abre nuevas vías para el descubrimiento de modelos dinámicos basados en datos, permitiendo identificar tanto la estructura como los parámetros de sistemas complejos de manera robusta.

Limitaciones: El método depende de la calidad de la interpolación; datos ruidosos o con muestreo muy escaso pueden afectar la precisión de las aproximaciones integrales, aunque técnicas de suavizado pueden mitigar este problema.