Apprendimento di operatori tramite Physics-Informed DeepONet Implementiamolo da zero

'Apprendimento operatori tramite Physics-Informed DeepONet'

Un’immersione profonda nelle DeepONets, reti neurali informate dalla fisica e DeepONets informate dalla fisica

Le equazioni differenziali ordinarie e parziali (ODE/PDE) sono il fondamento di molte discipline scientifiche ed ingegneristiche, dalla fisica alla biologia, dall’economia alla climatologia. Sono strumenti fondamentali utilizzati per descrivere sistemi e processi fisici, catturando il cambiamento continuo delle quantità nel tempo e nello spazio.

Tuttavia, una caratteristica unica di molte di queste equazioni è che non prendono solo valori singoli come input, ma prendono funzioni. Ad esempio, considera il caso della previsione delle vibrazioni in un edificio a causa di un terremoto. La scossa del terreno, che varia nel tempo, può essere rappresentata come una funzione che agisce come input per l’equazione differenziale che descrive il movimento dell’edificio. Allo stesso modo, nel caso delle onde sonore che si propagano in una sala da concerto, le onde sonore prodotte da uno strumento musicale possono essere una funzione di input con volume e tono variabili nel tempo. Queste funzioni di input variabili influenzano fondamentalmente le funzioni di output risultanti – le vibrazioni dell’edificio e il campo acustico nella sala, rispettivamente.

Tradizionalmente, queste ODE/PDE sono affrontate utilizzando risolutori numerici come i metodi delle differenze finite o degli elementi finiti. Tuttavia, questi metodi presentano un limite: per ogni nuova funzione di input, il risolutore deve essere eseguito nuovamente. Questo processo può essere computazionalmente intensivo e lento, specialmente per sistemi complessi o input ad alta dimensione.

Per affrontare questa sfida, è stato introdotto nel 2019 un nuovo framework da Lu et al.: la Deep Operator Network, o DeepONet. Le DeepONet mirano a imparare l’operatore che mappa le funzioni di input alle funzioni di output, imparando essenzialmente a prevedere l’output di queste ODE/PDE per qualsiasi funzione di input senza dover eseguire nuovamente un risolutore numerico ogni volta.

• Migliora la moderazione dei contenuti con l’API di Moderazione di OpenAI
• Cosa ho imparato da oltre 50 interviste di Machine Learning (come intervistatore)
• Iniziare con JAX

Ma le DeepONet, sebbene potenti, ereditano i problemi comuni dei metodi basati sui dati: come possiamo garantire che le previsioni della rete siano in linea con le leggi fisiche conosciute racchiuse nelle equazioni governanti?

Entra in gioco l’apprendimento informato dalla fisica.

L’apprendimento informato dalla fisica è un ramo in rapida evoluzione del machine learning che combina principi fisici con la scienza dei dati per migliorare la modellazione e la comprensione di sistemi fisici complessi. Coinvolge l’utilizzo della conoscenza specifica del dominio e delle leggi fisiche per guidare il processo di apprendimento e migliorare l’accuratezza, la generalizzazione e l’interpretabilità dei modelli di machine learning.

In questo contesto, nel 2021, Wang et al. hanno introdotto una nuova variante delle DeepONet: le Physics-Informed DeepONet. Questo approccio innovativo si basa sulla fondazione delle DeepONet incorporando la nostra comprensione delle leggi fisiche nel processo di apprendimento. Non stiamo più solo chiedendo al nostro modello di imparare dai dati; lo stiamo guidando con principi derivati da secoli di indagine scientifica.

Sembra essere un approccio molto promettente! Ma come dovremmo implementarlo nella pratica? Proprio questo sarà l’argomento che esploreremo oggi 🤗

In questo blog, discuteremo la teoria dietro le Physics-Informed DeepONet e mostreremo come implementarle da zero. Metteremo anche in azione il nostro modello sviluppato e ne dimostreremo la potenza attraverso uno studio di caso pratico.

Iniziamo!

Indice · 1. Studio di caso · 2. Physics-informed DeepONet ∘ 2.1 DeepONet: Una panoramica ∘ 2.2 Physics-Informed Neural Networks (PINNs) ∘ 2.3 Physics-Informed DeepONet · 3. Implementazione di Physics-Informed DeepONet ∘ 3.1 Definizione dell’architettura ∘ 3.2 Definizione della perdita ODE ∘ 3.3 Definizione del passo del gradiente discendente · 4. Generazione e organizzazione dei dati ∘ 4.1 Generazione dei profili u(·) ∘ 4.2 Generazione del dataset ∘ 4.3 Organizzazione del dataset · 5. Addestramento di Physics-Informed DeepONet · 6. Discussione dei risultati · 7. Conclusioni · Riferimenti

1. Studio di caso

Basiamo la nostra discussione su un esempio concreto. In questo blog, riprodurremo il primo studio di caso considerato nel paper originale di Wang et al., ovvero un problema di valore iniziale descritto dalla seguente equazione differenziale ordinaria (ODE):

con una condizione iniziale s(0) = 0.

In questa equazione, u( t ) è la funzione di input che varia nel tempo, e s( t ) è lo stato del sistema al tempo t di cui siamo interessati a prevedere. In uno scenario fisico, u( t ) potrebbe rappresentare una forza applicata a un sistema, e s( t ) potrebbe rappresentare la risposta del sistema, come il suo spostamento o velocità, a seconda del contesto. Il nostro obiettivo qui è apprendere la mappatura tra il termine di forzatura u( t ) e la soluzione dell’ODE s( t ).

Metodi numerici tradizionali come il metodo di Eulero o i metodi di Runge-Kutta possono risolvere efficacemente questa equazione. Tuttavia, notiamo che il termine di forzatura u( t ) può assumere vari profili, come mostrato dalla seguente figura:

Di conseguenza, ogni volta che u( t ) cambia, dovremmo rieseguire l’intera simulazione per ottenere il corrispondente s( t ) (come mostrato nella Figura 3), il che può essere computazionalmente intensivo e inefficiente.

Allora, come possiamo affrontare questo tipo di problema in modo più efficiente?

2. DeepONet informato dalla fisica

Come accennato nell’introduzione, il DeepONet informato dalla fisica costituisce una soluzione promettente per il nostro obiettivo. In questa sezione, analizzeremo i suoi concetti fondamentali per renderli più comprensibili.

Discuteremo prima i principi che sono alla base del DeepONet originale. Successivamente, esploreremo il concetto di reti neurali informate dalla fisica e come questo porti un’ulteriore dimensione al tavolo della risoluzione dei problemi. Infine, dimostreremo come possiamo integrare in modo fluido queste due idee per costruire il DeepONet informato dalla fisica.

2.1 DeepONet: Una panoramica

DeepONet, abbreviazione di Deep Operator Network, rappresenta una nuova frontiera nell’apprendimento profondo. A differenza dei metodi tradizionali di apprendimento automatico che mappano un insieme di valori di input a valori di output, DeepONet è progettato per mappare intere funzioni ad altre funzioni. Ciò rende DeepONet particolarmente potente quando si tratta di problemi che coinvolgono naturalmente input e output funzionali. Quindi come riesce esattamente ad ottenere tale obiettivo?

Per formulare ciò che vogliamo raggiungere in modo simbolico:

A sinistra, abbiamo l’operatore G che mappa da una funzione di input u(·) a una funzione di output s(·). A destra, vorremmo utilizzare una rete neurale per approssimare l’operatore G. Una volta che ciò può essere raggiunto, potremmo utilizzare la rete neurale addestrata per eseguire un calcolo rapido di s(·) dato qualsiasi u(·).

Per lo studio di caso attuale, sia la funzione di input u(·) che la funzione di output s(·) prendono come unico argomento il tempo t. Pertanto, l'”input” e l'”output” della rete neurale che intendiamo costruire dovrebbero essere così:

In sostanza, la nostra rete neurale dovrebbe accettare l’intero profilo di u( t ) come primo input, così come un’istanza di tempo specifica t come secondo input. Successivamente, dovrebbe restituire l’output target della funzione s(·) valutata nell’istanza di tempo t, cioè s( t ).

Per comprendere meglio questa configurazione, riconosciamo che il valore di s( t ) dipende innanzitutto dal profilo di s(·), che a sua volta dipende da u(·), e in secondo luogo dipende dal momento in cui viene valutato s(·). Questo è anche il motivo per cui è necessario che il tempo t sia tra gli input.

Ci sono due domande che dobbiamo chiarire al momento: innanzitutto, come dovremmo inserire un profilo continuo di u(·) nella rete? E in secondo luogo, come dovremmo concatenare i due input, cioè t e u(·)?

1️⃣ Come dovremmo inserire un profilo continuo di u(·)?

Bene, in realtà non lo facciamo. Una soluzione semplice consiste nel rappresentare la funzione u(·) in modo discreto. Più specificamente, valutiamo semplicemente i valori di u(·) in un numero sufficiente ma finito di posizioni e successivamente alimentiamo quei valori discreti di u(·) nella rete neurale:

Queste posizioni sono indicate come sensors nel documento originale di DeepONet.

2️⃣ Come dovremmo concatenare l’input t e u(·)?

A prima vista, potremmo volerli concatenare direttamente nel layer di input. Tuttavia, si scopre che questo approccio ingenuo non solo impone un vincolo su quali tipi di reti neurali possiamo utilizzare, ma porta anche a una precisione di previsione subottimale nella pratica. Esiste però un modo migliore. È il momento di introdurre DeepONet.

In poche parole, DeepONet ha proposto un’architettura di rete neurale per eseguire l’apprendimento di operatori: è composta da due componenti principali: una rete branch e una rete trunk. La rete branch prende in input i valori discreti della funzione e li trasforma in un vettore di feature. Nel frattempo, la rete trunk prende le coordinate (nel nostro caso di studio attuale, la coordinata è solo t. Per le equazioni differenziali alle derivate parziali, includerà le coordinate temporali e spaziali) e le converte anche in un vettore di feature con le stesse dimensioni. Questi due vettori di feature vengono quindi uniti tramite un prodotto scalare e il risultato finale viene utilizzato come previsione di s(·) valutata alla coordinata di input.

Nel paper originale di DeepONet, gli autori hanno affermato che questa strategia “divide-et-impera”, esemplificata in reti separate “branch” e “trunk”, è ispirata al teorema di approssimazione universale per operatori e serve per introdurre un forte bias induttivo specificamente per l’apprendimento degli operatori. Questo è anche il punto chiave che rende DeepONet una soluzione efficace, come affermato dagli autori.

Se sei curioso di approfondire le basi teoriche di DeepONet, consulta l’Appendice A nel paper originale.

Una delle principali forze di DeepONet è la sua efficienza. Una volta addestrato, un DeepONet può inferire la funzione di output per una nuova funzione di input in tempo reale, senza la necessità di ulteriori addestramenti, purché la nuova funzione di input sia all’interno del range delle funzioni di input su cui è stato addestrato. Questo rende DeepONet uno strumento potente in applicazioni che richiedono inferenze in tempo reale.

Un’altra nota forza di DeepONet risiede nella sua flessibilità e versatilità. Mentre la scelta più comune per le reti trunk e branch potrebbe essere costituita da layer completamente connessi, il framework di DeepONet consente un alto livello di personalizzazione dell’architettura. A seconda delle caratteristiche della funzione di input u(·) e delle coordinate, possono essere impiegate diverse architetture di reti neurali come CNN, RNN, ecc. Questa adattabilità rende DeepONet uno strumento altamente versatile.

Tuttavia, nonostante queste forze, le limitazioni di DeepONet sono anche evidenti: essendo un metodo puramente basato sui dati, DeepONet non può garantire che le sue previsioni seguiranno conoscenze precedenti o equazioni che descrivono il sistema fisico preso in considerazione. Di conseguenza, DeepONet potrebbe non generalizzare bene, specialmente quando si trova di fronte a funzioni di input che si trovano al di fuori della distribuzione dei dati di addestramento, definite come input out-of-distribution (OOD). Un rimedio comune a ciò è semplicemente preparare una grande quantità di dati per l’addestramento, il che potrebbe non sempre essere fattibile in pratica, specialmente nei campi scientifici e ingegneristici in cui la raccolta dei dati può essere costosa o richiedere molto tempo.

Quindi come dovremmo affrontare queste limitazioni? È ora di parlare di apprendimento informato dalla fisica, e più specificamente, di reti neurali informate dalla fisica (PINN).

2.2 Reti neurali informate dalla fisica (PINN)

Nei modelli di apprendimento automatico tradizionali, ci affidiamo principalmente ai dati per apprendere i pattern sottostanti. Tuttavia, in molti campi scientifici e ingegneristici, sono disponibili equazioni governanti (ODE/PDE) che racchiudono le nostre conoscenze precedenti sul sistema dinamico, e rappresentano un’altra fonte di informazioni su cui possiamo fare affidamento oltre ai dati osservati. Questa conoscenza aggiuntiva, se incorporata correttamente, potrebbe migliorare potenzialmente le prestazioni del modello e la sua capacità di generalizzazione, specialmente quando si tratta di dati limitati o rumorosi. Qui entra in gioco l’apprendimento informato dalla fisica.

Quando uniamo il concetto di apprendimento informato dalla fisica e le reti neurali, otteniamo le reti neurali informate dalla fisica (PINN).

Le PINN sono un tipo di rete neurale in cui la rete viene addestrata non solo per adattarsi ai dati, ma anche per rispettare le leggi fisiche conosciute descritte dalle equazioni differenziali. Questo viene realizzato introducendo una perdita di ODE/PDE, che misura il grado di violazione delle equazioni differenziali governanti. In questo modo, iniettiamo le leggi fisiche nel processo di addestramento della rete e lo rendiamo informato dal punto di vista fisico.

Nonostante le PINN si siano dimostrate efficaci in molte applicazioni, non sono senza limitazioni. Le PINN vengono tipicamente addestrate per parametri di input specifici (ad esempio, condizioni al contorno e iniziali, forzanti esterne, ecc.). Di conseguenza, ogni volta che i parametri di input cambiano, è necessario riaddestrare la PINN. Pertanto, non sono particolarmente efficaci per l’inferenza in tempo reale in diverse condizioni operative.

Ricordi ancora quale metodo è specificamente progettato per gestire parametri di input variabili? Esatto, è il DeepONet! È ora di combinare l’idea di apprendimento informato dalla fisica con DeepONet.

2.3 Physics-informed DeepONet

L’idea principale dietro il Physics-informed DeepONet è quella di combinare i punti di forza sia dei DeepONet che dei PINN. Proprio come un DeepONet, un Physics-informed DeepONet è in grado di prendere una funzione in input e produrre una funzione in output. Ciò lo rende altamente efficiente per l’inferenza in tempo reale di nuove funzioni di input, senza la necessità di un nuovo addestramento.

D’altra parte, come un PINN, un Physics-informed DeepONet incorpora le leggi fisiche conosciute nel suo processo di apprendimento. Queste leggi vengono introdotte come vincoli aggiuntivi nella funzione di perdita durante l’addestramento. Questo approccio consente al modello di fare previsioni coerenti dal punto di vista fisico, anche quando si lavora con dati limitati o rumorosi.

Come raggiungiamo questa integrazione? Similmente ai PINN, aggiungiamo un contributo di perdita extra per misurare quanto bene le previsioni del modello aderiscono all’equazione differenziale conosciuta. Ottimizzando questa funzione di perdita, il modello impara a fare previsioni che sono sia coerenti con i dati (se i dati di misurazione sono forniti durante l’addestramento) che coerenti dal punto di vista fisico.

In sintesi, il Physics-informed DeepONet è uno strumento potente che combina il meglio di entrambi i mondi: l’efficienza del DeepONet e l’accuratezza dell’apprendimento informato dalla fisica. Rappresenta un approccio promettente per risolvere problemi complessi in campi in cui sia l’inferenza in tempo reale che la coerenza fisica sono cruciali.

Nella prossima sezione, inizieremo a lavorare sul nostro caso di studio e trasformeremo la teoria in codice effettivo.

3. Implementazione del Physics-informed DeepONet

In questa sezione, vedremo come definire un modello Physics-informed DeepONet per affrontare il nostro caso di studio di riferimento. Lo implementeremo in TensorFlow. Iniziamo importando le librerie necessarie:

import numpy as npimport matplotlib.pyplot as pltimport tensorflow as tffrom tensorflow import kerastf.random.set_seed(42)

3.1 Definire l’Architettura

Come discusso in precedenza, il Physics-informed DeepONet condivide la stessa architettura del DeepONet originale. La seguente funzione definisce l’architettura per il DeepONet:

def create_model(mean, var, verbose=False):    """Definizione di un DeepONet con rami e tronchi completamente connessi.        Args:    ----    mean: dizionario, valori medi degli input    var: dizionario, valori di varianza degli input    verbose: booleano, indica se mostrare il riepilogo del modello        Outputs:    --------    model: il modello DeepONet    """        # Rete del ramo    branch_input = tf.keras.Input(shape=(len(mean['forcing'])), name="forcing")    branch = tf.keras.layers.Normalization(mean=mean['forcing'], variance=var['forcing'])(branch_input)    for i in range(3):        branch = tf.keras.layers.Dense(50, activation="tanh")(branch)        # Rete del tronco    trunk_input = tf.keras.Input(shape=(len(mean['time'])), name="time")    trunk = tf.keras.layers.Normalization(mean=mean['time'], variance=var['time'])(trunk_input)       for i in range(3):        trunk = tf.keras.layers.Dense(50, activation="tanh")(trunk)        # Calcola il prodotto scalare tra la rete del ramo e la rete del tronco    dot_product = tf.reduce_sum(tf.multiply(branch, trunk), axis=1, keepdims=True)        # Aggiungi il bias    output = BiasLayer()(dot_product)        # Crea il modello    model = tf.keras.models.Model(inputs=[branch_input, trunk_input], outputs=output)        if verbose:        model.summary()            return model   

Nel codice sopra:

1. Assumiamo che sia la rete del tronco che quella del ramo siano reti completamente connesse, con 3 strati nascosti, ognuno contenente 50 neuroni e con funzione di attivazione tangente iperbolica. Questa architettura è stata scelta in base a test preliminari e dovrebbe servire come un buon punto di partenza per questo problema. Tuttavia, è facile sostituirla con altre architetture (ad esempio, CNN, RNN, ecc.) e altri iperparametri di layer.
2. Le uscite delle reti del tronco e del ramo vengono unite tramite un prodotto scalare. Come suggerito nell’articolo originale del DeepONet, aggiungiamo un termine di bias per migliorare l’accuratezza delle previsioni. La classe BiasLayer() è una classe definita dall’utente per raggiungere questo obiettivo:

class BiasLayer(tf.keras.layers.Layer):    def build(self, input_shape):        self.bias = self.add_weight(shape=(1,),                                    initializer=tf.keras.initializers.Zeros(),                                    trainable=True)    def call(self, inputs):        return inputs + self.bias

3.2 Definire la perdita ODE

Successivamente, definiamo una funzione per calcolare la perdita ODE. Ricordiamo che la nostra ODE di destinazione è:

Pertanto, possiamo definire la funzione come segue:

@tf.functiondef Calcolatore_residual_ODE(t, u, u_t, model):    """Calcolo del residuo ODE.        Argomenti:    ----    t: coordinata temporale    u: funzione di input valutata in coordinate temporali discrete    u_t: funzione di input valutata in t    model: modello DeepONet        Output:    --------    ODE_residual: residuo dell'ODE governante    """        with tf.GradientTape() as tape:        tape.watch(t)        s = model({"forcing": u, "time": t})        # Calcola i gradienti    ds_dt = tape.gradient(s, t)        # Residuo ODE    ODE_residual = ds_dt - u_t        return ODE_residual

Nel codice sopra:

1. Abbiamo utilizzato tf.GradientTape() per calcolare il gradiente di s(·) rispetto a t . Notare che in TensorFlow, tf.GradientTape() viene utilizzato come un gestore di contesto, e tutte le operazioni eseguite all’interno del contesto del nastro verranno registrate dal nastro. Qui, abbiamo esplicitamente monitorato la variabile t . Di conseguenza, TensorFlow traccerà automaticamente tutte le operazioni che coinvolgono t , che in questo caso è un’esecuzione in avanti del modello DeepONet. Successivamente, utilizziamo il metodo gradient() del nastro per calcolare il gradiente di s(·) rispetto a t .
2. Abbiamo incluso un ulteriore argomento di input u_t , che indica il valore della funzione di input u(·) valutata in t . Questo costituisce il termine a destra della nostra ODE di destinazione ed è necessario per il calcolo della perdita residua ODE.
3. Abbiamo utilizzato il decoratore @tf.function per convertire la normale funzione Python appena definita in un grafo TensorFlow. È utile farlo poiché il calcolo del gradiente può essere piuttosto costoso ed eseguirlo in modalità grafica può accelerare significativamente i calcoli.

3.3 Definire il passo di discesa del gradiente

Successivamente, definiamo la funzione per compilare la funzione di perdita totale e calcolare i gradienti della perdita totale rispetto ai parametri del modello di rete:

@tf.functiondef passo_discesa_gradiente(X, X_init, peso_IC, peso_ODE, model):    """Calcola i gradienti della perdita totale rispetto ai parametri del modello di rete.        Argomenti:    ----    X: set di dati di addestramento per valutare i residui ODE    X_init: set di dati di addestramento per valutare le condizioni iniziali    peso_IC: peso per la perdita delle condizioni iniziali    peso_ODE: peso per la perdita ODE    model: modello DeepONet        Output:    --------    ODE_loss: perdita ODE calcolata    IC_loss: perdita delle condizioni iniziali calcolata    total_loss: somma ponderata della perdita ODE e della perdita delle condizioni iniziali    gradients: gradienti della perdita totale rispetto ai parametri del modello di rete.    """    with tf.GradientTape() as tape:        tape.watch(model.trainable_weights)        # Predizione delle condizioni iniziali        y_pred_IC = model({"forcing": X_init[:, 1:-1], "time": X_init[:, :1]})        # Residuo dell'equazione        ODE_residual = Calcolatore_residual_ODE(t=X[:, :1], u=X[:, 1:-1], u_t=X[:, -1:], model=model)        # Calcola la perdita        IC_loss = tf.reduce_mean(keras.losses.mean_squared_error(0, y_pred_IC))        ODE_loss = tf.reduce_mean(tf.square(ODE_residual))                # Perdita totale        total_loss = IC_loss*peso_IC + ODE_loss*peso_ODE    gradients = tape.gradient(total_loss, model.trainable_variables)    return ODE_loss, IC_loss, total_loss, gradients

Nel codice sopra:

1. Consideriamo solo due termini di perdita: la perdita associata alle condizioni iniziali IC_loss e la perdita residua ODE ODE_loss . La IC_loss viene calcolata confrontando il valore predetto dal modello s( t =0) con il valore iniziale noto di 0, e la ODE_loss viene calcolata chiamando la nostra funzione Calcolatore_residual_ODE precedentemente definita. La perdita dei dati può anche essere calcolata e aggiunta alla perdita totale se sono disponibili i valori misurati di s( t ) (non implementati nel codice sopra).
2. In generale, la perdita totale è una somma ponderata di IC_loss e ODE_loss , dove i pesi controllano quanto enfasi o priorità viene data a quei singoli termini di perdita durante il processo di addestramento. Nel nostro caso di studio, è sufficiente impostare entrambi peso_IC e peso_ODE a 1.
3. Similmente al calcolo della ODE_loss , abbiamo adottato anche tf.GradientTape() come gestore di contesto per calcolare i gradienti. Tuttavia, qui calcoliamo i gradienti della perdita totale rispetto ai parametri del modello di rete, che sono necessari per eseguire la discesa del gradiente.

Prima di procedere, riassumiamo brevemente ciò che abbiamo sviluppato finora:

1️⃣ Possiamo inizializzare un modello DeepONet con la funzione create_model().

2️⃣ Possiamo calcolare i residui ODE per valutare quanto bene le previsioni del modello si adattano all’ODE governante. Questo viene realizzato con la funzione ODE_residual_calculator.

3️⃣ Possiamo calcolare la perdita totale e i suoi gradienti rispetto ai parametri del modello di rete con train_step.

Ora il lavoro di preparazione è a metà 🚀 Nella prossima sezione, discuteremo la generazione dei dati e le questioni di organizzazione dei dati (il codice strano X[:, :1] nel codice sopra diventerà sperabilmente chiaro). Dopo di che, possiamo finalmente allenare il modello e vedere come si comporta.

4. Generazione e organizzazione dei dati

In questa sezione, discutiamo la generazione di dati sintetici e come organizzarli per addestrare il modello DeepONet informato dalla fisica.

4.1 Generazione di profili u(·)

I dati utilizzati per l’addestramento, la convalida e il test saranno generati in modo sintetico. La motivazione dietro questo approccio è duplice: è comodo ma consente anche un controllo completo sulle caratteristiche dei dati.

Nel contesto del nostro caso di studio, genereremo la funzione di input u(·) utilizzando un Gaussian Process a media zero, con un kernel di funzione di base radiale (RBF).

Un Gaussian Process è un potente framework matematico comunemente usato nell’apprendimento automatico per modellare funzioni. Il kernel RBF è una scelta popolare per catturare la similarità tra i punti di input. Utilizzando il kernel RBF all’interno del Gaussian Process, garantiamo che i dati sintetici generati presentino un modello liscio e continuo, che è spesso desiderabile in diverse applicazioni. Per saperne di più sui Gaussian Process, sentiti libero di consultare il mio precedente blog.

In scikit-learn, ciò può essere ottenuto in poche righe di codice:

from sklearn.gaussian_process import GaussianProcessRegressorfrom sklearn.gaussian_process.kernels import RBFdef create_samples(length_scale, sample_num):    """Crea dati sintetici per u(·)        Args:    ----    length_scale: float, scala del kernel RBF    sample_num: numero di profili u(·) da generare        Outputs:    --------    u_sample: profili u(·) generati    """    # Definisci il kernel con la scala di lunghezza data    kernel = RBF(length_scale)    # Crea regressore del processo gaussiano    gp = GaussianProcessRegressor(kernel=kernel)    # Posizioni dei punti di collocamento    X_sample = np.linspace(0, 1, 100).reshape(-1, 1)         # Crea campioni    u_sample = np.zeros((sample_num, 100))    for i in range(sample_num):        # campionamento dal priore direttamente        n = np.random.randint(0, 10000)        u_sample[i, :] = gp.sample_y(X_sample, random_state=n).flatten()              return u_sample

Nel codice sopra:

1. Utilizziamo length_scale per controllare la forma della funzione generata. Per un kernel RBF, la Figura 11 mostra il profilo u(·) dato diversi valori di lunghezza del kernel.
2. Ricordiamo che è necessario discretizzare u(·) prima di alimentarlo al DeepONet. Questo viene fatto specificando una variabile X_sample, che alloca 100 punti distribuiti uniformemente nel nostro dominio temporale di interesse.
3. In scikit-learn, l’oggetto GaussianProcessRegressor espone un metodo sample_y per consentire il disegno di campioni casuali dal processo gaussiano con il kernel specificato dalla scala di lunghezza. Nota che non abbiamo chiamato .fit() prima di utilizzare l’oggetto GaussianProcessRegressor, a differenza di quanto facciamo normalmente con altri regressori di scikit-learn. Questo è intenzionale poiché vogliamo che GaussianProcessRegressor utilizzi l’esatta length_scale che abbiamo fornito. Se chiami .fit(), la length_scale sarà ottimizzata per un altro valore per adattarsi meglio ai dati forniti.
4. L’output u_sample è una matrice con una dimensione di sample_num * 100. Ogni riga di u_sample rappresenta un profilo di u(·), che consiste in 100 valori discreti.

4.2 Generazione del Dataset

Ora che abbiamo generato i profili di u(·), concentriamoci su come organizzare il dataset in modo che possa essere alimentato nel modello DeepONet.

Ricordiamo che il modello DeepONet che abbiamo sviluppato nella sezione precedente richiede 3 input:

1. la coordinata temporale t, che è uno scalare compreso tra 0 e 1 (per il momento non consideriamo la dimensione del batch);
2. il profilo di u(·), che è un vettore composto dai valori di u(·) valutati in coordinate temporali predefinite e fisse tra 0 e 1;
3. il valore di u(t), che è nuovamente uno scalare. Questo valore di u(t) viene utilizzato per calcolare la perdita ODE alla coordinata temporale t.

Pertanto, possiamo formulare un singolo campione in questo modo:

Ovviamente, per ogni profilo di u(·) (marcato in verde nell’illustrazione sopra), dovremmo considerare più t (e i corrispondenti u(t)) per valutare la perdita ODE al fine di imporre meglio i vincoli fisici. In teoria, t può assumere qualsiasi valore all’interno del dominio temporale considerato (cioè tra 0 e 1 per il nostro caso di studio). Tuttavia, per semplificare le cose, considereremo t negli stessi punti temporali in cui il profilo di u(·) è discretizzato. Di conseguenza, il nostro dataset aggiornato sarà così:

Si osservi che la discussione precedente considera solo un singolo profilo di u(·). Se tenessimo conto di tutti i profili di u(·), il nostro dataset finale sarebbe così:

dove N indica il numero di profili di u(·). Ora tenendo presente ciò, vediamo un po’ di codice:

from tqdm import tqdmfrom scipy.integrate import solve_ivpdef generate_dataset(N, length_scale, ODE_solve=False):    """Genera il dataset per l'addestramento di DeepONet informato dalla fisica.        Argomenti:    ----    N: int, numero di profili di u(·)    length_scale: float, scala di lunghezza per il kernel RNF    ODE_solve: booleano, indica se calcolare il corrispondente s(·)        Output:    --------    X: il dataset per t, i profili di u(·) e u(t)    y: il dataset per la soluzione ODE corrispondente s(·)    """        # Creazione dei campi casuali    random_field = create_samples(length_scale, N)        # Compilazione del dataset    X = np.zeros((N*100, 100+2))    y = np.zeros((N*100, 1))    for i in tqdm(range(N)):        u = np.tile(random_field[i, :], (100, 1))        t = np.linspace(0, 1, 100).reshape(-1, 1)        # u(·) valutato in t        u_t = np.diag(u).reshape(-1, 1)        # Aggiornamento della matrice complessiva        X[i*100:(i+1)*100, :] = np.concatenate((t, u, u_t), axis=1)        # Risoluzione dell'ODE        if ODE_solve:            sol = solve_ivp(lambda var_t, var_s: np.interp(var_t, t.flatten(), random_field[i, :]),                             t_span=[0, 1], y0=[0], t_eval=t.flatten(), method='RK45')            y[i*100:(i+1)*100, :] = sol.y[0].reshape(-1, 1)            return X, y

Nel codice sopra, aggiungiamo una opzione per calcolare il corrispondente s(·) per un dato profilo u(·). Anche se non useremo i valori di s(·) nell’addestramento, avremmo ancora bisogno di loro per testare le prestazioni del modello. Il calcolo di s(·) viene ottenuto utilizzando scipy.integrate.solve_ivp, che è un risolutore di equazioni differenziali ordinarie di SciPy appositamente progettato per risolvere problemi di valore iniziale.

Ora possiamo generare il dataset di addestramento, di validazione e di test. Nota che per questo caso di studio, utilizzeremo una scala di lunghezza di 0.4 per generare i profili u(·) e addestrare il DeepONet con informazioni fisiche.

# Creare dataset di addestramentoN_train = 2000length_scale_train = 0.4X_train, y_train = generate_dataset(N_train, length_scale_train)# Creare dataset di validazioneN_val = 100length_scale_test = 0.4X_val, y_val = generate_dataset(N_val, length_scale_test)# Creare dataset di testN_test = 100length_scale_test = 0.4X_test, y_test = generate_dataset(N_test, length_scale_test, ODE_solve=True)

4.3 Organizzazione del dataset

Infine, convertiamo l’array NumPy negli oggetti dataset TensorFlow per facilitare l’ingestione dei dati.

# Determinare la dimensione del batchini_batch_size = int(2000/100)col_batch_size = 2000# Creare oggetto dataset (condizioni iniziali)X_train_ini = tf.convert_to_tensor(X_train[X_train[:, 0]==0], dtype=tf.float32)ini_ds = tf.data.Dataset.from_tensor_slices((X_train_ini))ini_ds = ini_ds.shuffle(5000).batch(ini_batch_size)# Creare oggetto dataset (punti di collocation)X_train = tf.convert_to_tensor(X_train, dtype=tf.float32)train_ds = tf.data.Dataset.from_tensor_slices((X_train))train_ds = train_ds.shuffle(100000).batch(col_batch_size)# Media della scala media = {    'forcing': np.mean(X_train[:, 1:-1], axis=0),    'time': np.mean(X_train[:, :1], axis=0)}var = {    'forcing': np.var(X_train[:, 1:-1], axis=0),    'time': np.var(X_train[:, :1], axis=0)}

Nel codice sopra, creiamo due dataset distinti: uno per valutare la perdita dell’ODE (train_ds), e l’altro per valutare la perdita della condizione iniziale (ini_ds). Abbiamo anche calcolato in anticipo i valori di media e varianza per t e u(·). Questi valori saranno utilizzati per standardizzare gli input.

Questo è tutto per la generazione e l’organizzazione dei dati. Successivamente, inizieremo l’addestramento del modello e vedremo come si comporta.

5. Addestramento di DeepONet con informazioni fisiche

Come primo passo, creiamo una classe personalizzata per tracciare l’evoluzione della perdita:

from collections import defaultdictclass LossTracking:    def __init__(self):        self.mean_total_loss = keras.metrics.Mean()        self.mean_IC_loss = keras.metrics.Mean()        self.mean_ODE_loss = keras.metrics.Mean()        self.loss_history = defaultdict(list)    def update(self, total_loss, IC_loss, ODE_loss):        self.mean_total_loss(total_loss)        self.mean_IC_loss(IC_loss)        self.mean_ODE_loss(ODE_loss)    def reset(self):        self.mean_total_loss.reset_states()        self.mean_IC_loss.reset_states()        self.mean_ODE_loss.reset_states()    def print(self):        print(f"IC={self.mean_IC_loss.result().numpy():.4e}, \              ODE={self.mean_ODE_loss.result().numpy():.4e}, \              total_loss={self.mean_total_loss.result().numpy():.4e}")            def history(self):        self.loss_history['total_loss'].append(self.mean_total_loss.result().numpy())        self.loss_history['IC_loss'].append(self.mean_IC_loss.result().numpy())        self.loss_history['ODE_loss'].append(self.mean_ODE_loss.result().numpy())

Poi, definiamo la logica principale di addestramento/validazione:

# Impostare le configurazioni di addestramenton_epochs = 300IC_weight= tf.constant(1.0, dtype=tf.float32)   ODE_weight= tf.constant(1.0, dtype=tf.float32)loss_tracker = LossTracking()val_loss_hist = []# Impostare l'ottimizzatoreoptimizer = keras.optimizers.Adam(learning_rate=1e-3)# Istanziare il modello PINNPI_DeepONet = create_model(mean, var)PI_DeepONet.compile(optimizer=optimizer)  # Configurare i callbackscallbacks = [keras.callbacks.ReduceLROnPlateau(factor=0.5, patience=30),             tf.keras.callbacks.ModelCheckpoint('NN_model.h5', monitor='val_loss', save_best_only=True)]callbacks = tf.keras.callbacks.CallbackList(_callbacks, add_history=False, model=PI_DeepONet)# Iniziare il processo di addestramentofor epoch in range(1, n_epochs + 1):      print(f"Epoch {epoch}:")    for X_init, X in zip(ini_ds, train_ds):                # Calcolare i gradienti        ODE_loss, IC_loss, total_loss, gradients = train_step(X, X_init,                                                             IC_weight, ODE_weight,                                                            PI_DeepONet)        # Discesa del gradiente        PI_DeepONet.optimizer.apply_gradients(zip(gradients, PI_DeepONet.trainable_variables))                # Tracciamento della perdita        loss_tracker.update(total_loss, IC_loss, ODE_loss)    # Riepilogo della perdita    loss_tracker.history()    loss_tracker.print()    loss_tracker.reset()    ####### Validazione    val_res = ODE_residual_calculator(X_val[:, :1], X_val[:, 1:-1], X_val[:, -1:], PI_DeepONet)    val_ODE = tf.cast(tf.reduce_mean(tf.square(val_res)), tf.float32)        X_val_ini = X_val[X_val[:, 0]==0]    pred_ini_valid = PI_DeepONet.predict({"forcing": X_val_ini[:, 1:-1], "time": X_val_ini[:, :1]}, batch_size=12800)    val_IC = tf.reduce_mean(keras.losses.mean_squared_error(0, pred_ini_valid))    print(f"val_IC: {val_IC.numpy():.4e}, val_ODE: {val_ODE.numpy():.4e}, lr: {PI_DeepONet.optimizer.lr.numpy():.2e}")        # Callback alla fine dell'epoca    callbacks.on_epoch_end(epoch, logs={'val_loss': val_IC+val_ODE})    val_loss_hist.append(val_IC+val_ODE)        # Rimescolamento del dataset    ini_ds = tf.data.Dataset.from_tensor_slices((X_train_ini))    ini_ds = ini_ds.shuffle(5000).batch(ini_batch_size)    train_ds = tf.data.Dataset.from_tensor_slices((X_train))    train_ds = train_ds.shuffle(100000).batch(col_batch_size)

È un pezzo di codice piuttosto lungo, ma dovrebbe essere autoesplicativo poiché abbiamo già coperto tutte le parti importanti.

Per visualizzare le prestazioni di addestramento, possiamo tracciare le curve di convergenza della perdita:

# Historyfig, ax = plt.subplots(1, 3, figsize=(12, 4))ax[0].plot(range(n_epochs), loss_tracker.loss_history['IC_loss'])ax[1].plot(range(n_epochs), loss_tracker.loss_history['ODE_loss'])ax[2].plot(range(n_epochs), val_loss_hist)ax[0].set_title('IC Loss')ax[1].set_title('ODE Loss')ax[2].set_title('Val Loss')for axs in ax:    axs.set_yscale('log')

I risultati dell’addestramento appaiono così:

Inoltre, possiamo anche vedere come l’accuratezza delle previsioni per un determinato target s(·) si evolve durante l’addestramento:

All’inizio dell’addestramento, possiamo vedere una discrepanza visibile tra la previsione del modello e il vero valore. Tuttavia, verso la fine dell’addestramento, la previsione di s(·) converge al vero valore. Ciò indica che il nostro DeepONet informato dalla fisica apprende correttamente.

6. Discussione dei risultati

Una volta completato l’addestramento, possiamo ricaricare i pesi salvati e valutare le prestazioni.

In questo caso abbiamo selezionato casualmente tre profili u(·) dal dataset di testing e confrontato i rispettivi s(·) predetti dal nostro DeepONet informato dalla fisica con quelli calcolati dal risolutore numerico delle equazioni differenziali ordinarie. Possiamo vedere che le previsioni e il vero valore sono quasi indistinguibili.

Questi risultati sono incredibili, considerando il fatto che non abbiamo nemmeno utilizzato alcun dato osservativo di s(·) (tranne la condizione iniziale) per addestrare il DeepONet. Questo dimostra che l’ODE governante stesso ha fornito un segnale di “supervisione” sufficiente per il modello per effettuare previsioni accurate.

Un’altra cosa interessante da valutare è la cosiddetta capacità di predizione “fuori distribuzione”. Poiché abbiamo imposto l’equazione governante durante l’addestramento del DeepONet, possiamo aspettarci che il DeepONet informato dalla fisica addestrato sia comunque in grado di effettuare previsioni decenti quando i profili u(·) si trovano al di fuori della distribuzione degli u(·) di addestramento.

Per testarlo, possiamo generare profili u(·) utilizzando una scala di lunghezza diversa. I seguenti risultati mostrano tre profili u(·) generati con una scala di lunghezza pari a 0,6, così come le s(·) previste. Questi risultati sono molto buoni, considerando che il DeepONet informato dalla fisica è stato addestrato con una scala di lunghezza pari a 0,4.

Tuttavia, se riduciamo ulteriormente la scala di lunghezza a 0,2, noteremo che iniziano a comparire discrepanze visibili. Ciò indica che c’è un limite alla capacità di generalizzazione del DeepONet informato dalla fisica addestrato.

Scale di lunghezza più piccole in generale portano a profili u(·) più complessi, che sarebbero molto diversi dai profili u(·) utilizzati per l’addestramento. Ciò potrebbe spiegare perché il modello addestrato ha incontrato difficoltà nel fare previsioni accurate nelle regioni a scala di lunghezza più piccola.

Complessivamente, potremmo dire che il DeepONet informato dalla fisica sviluppato può imparare correttamente la dinamica del sistema e mappare dalla funzione di input alla funzione di output dati solo i vincoli delle equazioni differenziali ordinarie. Inoltre, il DeepONet informato dalla fisica mostra un certo livello di capacità di gestire previsioni “fuori distribuzione”, indicando che l’addestramento del modello per allinearsi con l’equazione differenziale governante migliora la capacità di generalizzazione del modello.

7. Take-away

Siamo arrivati ​​lontano nella nostra esplorazione di Physics-Informed DeepONet. Dalla comprensione dei concetti fondamentali di DeepONet e dell’apprendimento informato dalla fisica, alla loro applicazione attraverso l’implementazione del codice, abbiamo coperto molto su questo potente metodo per risolvere equazioni differenziali.

Ecco alcuni punti chiave da ricordare:

1️⃣ DeepONet è un framework potente per eseguire l’apprendimento degli operatori, grazie alla sua architettura innovativa di reti a rami e tronco.

2️⃣ Apprendimento informato dalla fisica incorpora esplicitamente le equazioni differenziali governanti del sistema dinamico nel processo di apprendimento, possedendo così il potenziale di migliorare l’interpretabilità e la capacità di generalizzazione del modello.

3️⃣ DeepONet informato dalla fisica combina i punti di forza di DeepONet e dell’apprendimento informato dalla fisica, presentandosi come uno strumento promettente per apprendere mappe funzionali nel rispetto delle equazioni governanti associate.

Spero che tu abbia apprezzato questa immersione profonda in Physics-Informed DeepONet. Prossimamente, passeremo a risolvere problemi inversi con Physics-Informed DeepONet. Resta sintonizzato🤗

Puoi trovare il notebook di accompagnamento con il codice completo qui 💻

Referenze

[1] Lu et al., DeepONet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv, 2019.

[2] Wang et al., Learning the solution operator of parametric partial differential equations with physics-informed DeepOnets. arXiv, 2021.