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

1. Il Problema

L'articolo affronta il problema inverso della stima dei parametri per equazioni differenziali ordinarie (ODE) e ritardate (DDE), un compito fondamentale nell'ingegneria chimica e nei sistemi dinamici. L'obiettivo è trovare i valori dei parametri che minimizzano la discrepanza tra le previsioni del modello e i dati sperimentali.

Le sfide principali identificate sono:

• Costo Computazionale: I metodi tradizionali come il single-shooting (ottimizzazione sequenziale) richiedono l'integrazione numerica delle ODE ad ogni iterazione e il calcolo delle sensibilità (derivate degli stati rispetto ai parametri), che è computazionalmente oneroso.
• Convergenza: I metodi sequenziali faticano a convergere per traiettorie lunghe, stime iniziali povere o dinamiche altamente non lineari e complesse.
• Instabilità: I sistemi rigidi (stiff) o caotici possono portare a instabilità numerica durante l'integrazione, specialmente se i parametri iniziali sono lontani dal valore ottimo.
• Limitazioni dei Metodi Esistenti: I metodi basati sull'interpolazione precedenti erano limitati a parametri che apparivano linearmente nel modello.

2. Metodologia Proposta

Gli autori propongono un framework di ottimizzazione a due livelli (bi-level) che combina l'interpolazione dei dati con la struttura matematica delle equazioni differenziali.

A. Formulazione del Problema

Il problema è formulato come un'ottimizzazione gerarchica:

1. Livello Interno (Inner Problem): Risolve un problema di ottimizzazione convessa rispetto ai parametri lineari ($p$), dato un valore fisso dei parametri non lineari ($\phi$).
   • Gli stati del sistema $X(t)$ sono approssimati utilizzando tecniche di interpolazione sui dati misurati ($\Psi(X)(t)$), eliminando la necessità di integrare le ODE in ogni passo.
   • L'integrale delle dinamiche viene approssimato separando i parametri lineari.
2. Livello Esterno (Outer Problem): Ottimizza i parametri non lineari ($\phi$) minimizzando l'errore quadratico medio tra i dati e il modello.

B. Calcolo dei Gradienti e Teorema della Funzione Implicita

Per collegare i due livelli e permettere l'aggiornamento dei parametri non lineari, è necessario calcolare il gradiente della soluzione ottima del livello interno rispetto ai parametri esterni.

• Gli autori utilizzano il Teorema della Funzione Implicita applicato alle condizioni di Karush-Kuhn-Tucker (KKT) del problema interno convesso.
• Questo permette di derivare analiticamente le sensibilità ($dp^*/d\phi$) senza dover ricorrere a differenze finite o a un'integrazione numerica diretta delle sensibilità attraverso l'ODE.
• Viene utilizzata la differenziazione automatica (AD) (tramite la libreria JAX) per calcolare efficientemente i gradienti e gli Hessiani, sfruttando la fattorizzazione della matrice KKT del risolutore interno.

C. Vantaggi dell'Interpolazione

Utilizzando l'interpolazione sui dati misurati invece di integrare le ODE forward:

• Si evitano le traiettorie instabili che possono derivare da stime iniziali errate dei parametri.
• Si riduce drasticamente il costo computazionale del calcolo delle sensibilità, specialmente per sistemi rigidi.
• Il metodo diventa robusto anche per equazioni differenziali con ritardi (DDE) e sistemi con stati parzialmente osservati.

3. Contributi Chiave

1. Estensione ai Parametri Non Lineari: A differenza dei lavori precedenti basati sull'interpolazione (limitati a parametri lineari), questo metodo gestisce genericamente parametri che appaiono in modo non lineare nelle equazioni.
2. Framework Bi-Level Convesso: La decomposizione del problema in un sottoproblema interno convesso e un problema esterno non convesso garantisce una migliore convergenza verso l'ottimo globale rispetto ai metodi sequenziali.
3. Generalità del Metodo: La metodologia è stata dimostrata efficace su una vasta gamma di problemi, inclusi:
   • ODE standard (dinamiche oscillanti, modelli biochimici).
   • Equazioni differenziali con ritardi (modelli epidemici).
   • Sistemi rigidi (reazione di Belousov).
   • Sistemi con stati parzialmente osservati (reattori CSTR).
   • Scoperta di Modelli (Model Discovery): Estensione del metodo per apprendere sia la struttura del modello (selezione delle basi polinomiali) che i parametri, utilizzando algoritmi di thresholding sequenziale.

4. Risultati Sperimentali

Gli autori hanno confrontato il metodo proposto con il single-shooting e il collocamento ortogonale su diversi benchmark (es. ioni Calcio, fermentazione etanolo, reazione di Belousov).

• Convergenza: Il metodo bi-level ha converguto all'ottimo globale in tutti i casi di test. Al contrario, il single-shooting e il collocamento ortogonale hanno fallito la convergenza o sono rimasti bloccati in minimi locali in diversi scenari (specialmente per sistemi rigidi o con dinamiche complesse).
• Accuratezza: Gli errori quadratici medi (MSE) tra le traiettorie simulate e quelle reali sono stati significativamente inferiori con il metodo proposto.
• Iterazioni vs. Tempo: Sebbene il numero di iterazioni richieste dal metodo bi-level sia inferiore (grazie alla convessità interna), il tempo di esecuzione (wall time) è generalmente più alto rispetto al single-shooting. Questo è dovuto al costo di risolvere il problema interno convesso e di invertire la matrice KKT ad ogni passo esterno. Tuttavia, la maggiore affidabilità e la capacità di risolvere problemi altrimenti intrattabili giustificano questo costo.
• Robustezza: Il metodo ha dimostrato di funzionare bene anche con sistemi rigidi (dove i solver espliciti falliscono nel single-shooting) e con stati non osservati, utilizzando l'interpolazione per stimare le traiettorie latenti.

5. Significato e Conclusioni

Questo lavoro rappresenta un avanzamento significativo nella stima dei parametri per sistemi dinamici complessi.

• Superamento delle limitazioni numeriche: Offrendo un'alternativa stabile all'integrazione numerica diretta, il metodo mitiga i problemi di instabilità tipici dei sistemi rigidi e caotici.
• Scalabilità verso la scoperta di modelli: La capacità di integrare la selezione della struttura del modello (model discovery) all'interno dello stesso framework di ottimizzazione apre nuove possibilità per l'analisi di dati complessi in ingegneria chimica e biologia.
• Implicazioni Pratiche: Sebbene il costo computazionale per iterazione sia superiore, la capacità di trovare soluzioni globali affidabili senza richiedere stime iniziali estremamente precise rende questo approccio superiore per applicazioni reali dove i dati sono abbondanti ma i modelli sono complessi e non lineari.

In sintesi, l'approccio proposto trasforma un problema di ottimizzazione non convesso e difficile in una sequenza di problemi più gestibili (uno convesso interno e uno di ottimizzazione dei parametri non lineari), sfruttando l'interpolazione per bypassare le difficoltà numeriche dell'integrazione diretta delle ODE.