Identificazione parametrica e meccanica respiratoria

Nel seguente report viene analizzato un modello lineare di meccanica polmonare. Il modello corrisponde ad un circuito RLC con parametri incogniti. L’identificazione parametrica permette di risalire ai coefficienti di governo del sistema una volta stimato un modello matematico. Diversi approcci permettono di minimizzare la differenza tra i parametri stimati e le misure reali. Viene affrontato un approccio privo di differenziazione numerica, tramite algoritmo di Nelder-Mead e un approccio mediante il metodo di Gauss-Newton rilassato. Entrambi i metodi portano alla stima corretta dei parametri ma alcune considerazioni permettono di ottimizzarne le prestazioni.

Il seguente articolo è un estratto di un progetto realizzato per il corso di Modellazione e Simulazione di Sistemi Fisiologici. Tenuto per Ingegneria Medica presso l’Università degli Studi di Roma Tor Vergata. 

Introduzione

La meccanica respiratoria viene descritta in letteratura da numerosi modelli. Tra questi, un modello valido in determinate condizioni è rappresentato da un circuito RLC. Il modello permette di descrivere la meccanica respiratoria di un individuo sano.

A partire dalle misure dell’indivuduo è possibile stimare i parametri per rappresentarlo tramite il modello stimato. Questo processo, sotto il nome di ottimizzazione parametrica, permette di adattare il sistema generale al caso specifico. Nella seguente analisi viene affrontata la tecnica di identificazione parametrica tramite minimizzazione della differenza tra la risposta misurata e quella stimata.

Differenti tecniche di ottimizzazione permettono di stimare i parametri che permettono di descrivere il sistema considerato. Viene descritto sia il metodo di Nelder-Mead sia il metodo di Gauss-Newton, analizzandone punti in comune e differenze.

Vengono fatte anche alcune considerazioni sulla convergenza del metodo e sul corretto condizionamento del problema.

Background

L’identificazione di sistemi è un insieme di tecniche che si prefigge lo scopo di identificare i parametri e il sistema tale da descrivere un certo fenomeno. La prima grande distinzione vede la separazione di due tecniche: identificazione parametrica e non parametrica. La prima parte stimando un modello del fenomeno e punta ad identificare i parametri che lo rappresentano. L’identificazione non parametrica, invece, procede senza conoscere il modello del fenomeno, al più ne stima alcune caratteristiche, e punta ad identificare il sistema che meglio lo descrive.

In questo caso si parte da un modello ben noto, un modello di meccanica polmonare lineare descritto da un circuito RLC, e se ne identificano i parametri descrittivi.

Sistema RLC

Il circuito RLC (fig. 2) permette, seguendo l’analogia elettrica, di descrivere un modello approssimato della meccanica respiratoria. Il modello prevedere una serie di resistenza, induttore e capacità. Sulla base dell’analogia tra la corrente elettrica e il flusso polmonare è possibile rappresentare, in un modello a parametri concentrati, le componenti di compliance e resistenza elastica dei polmoni [1].

Questo modello permette di affrontare i casi fisiologici in cui il flusso entrante nei polmoni è predominante rispetto al flusso bloccato (flusso che invece aumenta notevolmente in caso di patologie ostruttive). L’analogia vede l’associazione della differenza di pressione alla differenza di potenziale, la corrente al flusso d’aria e l’induttanza come un’inertanza, ovvero la differenza di pressione richiesta per causare una variazione unitaria nel tasso di variazione del flusso nel tempo.

In un caso paziente specifico sarebbero necessarie misure sperimentali così da ottenere il modello reale sul quale effettuare l’identificazione. In questo caso, non avendo a disposizione dati sperimentali, si utilizzano dei dati noti in letteratura [2].

\mathrm{R}=0.1\left[\frac{\mathrm{cmH} 2 \mathrm{Os}}{\mathrm{L}}\right]\\ \:\\
\mathrm{L}=0.01\left[\frac{\mathrm{cmH} 2 \mathrm{O} \mathrm{s}^{2}}{\mathrm{~L}}\right]\\\: \\
 \mathrm{C}=0.1\left[\frac{\mathrm{L}}{\mathrm{cmH} 2 \mathrm{O}}\right]\\

Si considerano un’unica resistenza e un’unica compliance che vanno a rappresentare la complessiva e la capacità di accumulare aria del sistema respiratorio nel suo complesso. Quindi la resisitenza R porta con se il contributo resistivo delle vie aeree, dei tessuti polmonari e della parete toracica. La capacità C rappresenta il compliance dei tessuti e della parete toracica.

Legge di Kirchoff

L’obiettivo è la predizione della pressione alveolare 𝑝𝐴 e della sua risposta dinamica a differenti forme d’onda di pressione applicare all’apertura delle vie aeree 𝑝𝑎𝑜.

Applicando la prima legge di Kirchoff al circuito in fig. 2 si ottiene:

\begin{equation}
p_{a O}=L \frac{d i}{d t}+R i+p_{a}
\end{equation}

Dove vale il legame:

\begin{equation}
i=C \frac{d p_{a}}{d t}
\end{equation}

Da cui si ricava la funzione di trasferimento:

\begin{equation}
H(s)=\frac{P_{a}(s)}{P_{a o}(s)}=\frac{1}{1+s R C+s^{2} L C}
\end{equation}

Questo mostra come 𝑝𝑎 è interamente dipendente da 𝑝𝑎𝑜 e il sistema è un sistema con configurazione open-loop.

Identificazione dei parametri

Non avendo a disposizione dati paziente specifici sulla meccanica respiratoria i dati verranno generati tramite la procedura rlc_fun.m che restituisce la risposta del circuito RLC fornendogli in ingressi i parametri del sistema. A questo può essere aggiunto anche del rumore.

Algoritmo di Nelder-Mead

Una volta assunto come modello del sistema la funzione di trasferimento di un circuito RLC, per stimare i parametri è necessario definire una funzione obiettivo. Tale funzione rappresenta la differenza tra la risposta misurata (vera) del sistema e quella stimata:

\begin{equation}
\underline{e}=\underline{y}^{\mathrm{pred}}-\underline{y}^{\mathrm{mis}}
\end{equation}

E la ricerca dei parametri corretti punta alla minimizzazione della funzione obiettivo definita come:

\begin{equation}
\mathrm{E}=\frac{1}{2}\sum \|\underline{e}\|^{2}
\end{equation}

Uno dei possibili modi per affrontare questo problema consiste nell’affidarsi alla routine di Matlab nota come fminsearch(). Tale procedura si occupa di restituire i parametri ottimali. In particolare, questa funzione sfrutta l’algoritmo di Nelder Mead [3]. Tale metodo non fa uso di derivate e si basa sul concetto di un simplesso approssimando il punto di ottimo locale in 𝑛 variabili.

Metodo di Gauss-Newton

Un’alternativa è il metodo di Gauss-Newton. Il metodo, diversamente dal precedente, richiede di portare in conto anche le derivate. La minimizzazione dell’eq. (5) è un problema non lineare ai minimi quadrati. La soluzione numerica sfrutta il metodo di Gauss Newton. Come il metodo di Newton, si impone la stazionarietà, condizione che si riflette sulla jma componente:

\begin{equation}
(\nabla \mathrm{E})_{j}=\sum_{k} e_{k} \cdot \frac{\partial e_{k}}{\partial p_{j}}
\end{equation}

Ovvero tale differenziazione corrisponde a premoltiplicare per la matrice di sensibilità J:

\begin{equation}
\nabla \mathrm{E}=\underline{\underline{J^{}}}^{T} \:\underline{e}
\end{equation}

Utilizzando quindi l’algoritmo di Newton l’iterazione al passo 𝑘 + 1 cercherà i parametri (𝑝) come:

\begin{equation}
\underline{p}^{(k+1)}=\underline{p}^{(k)}-\left.[\nabla(\nabla \mathrm{E})]^{-1}(\nabla \mathrm{E})\right|_{\underline{p}^{k}}
\end{equation}

Per l’approssimazione introdotta da Gauss [5] è possibile trascurare le componenti quadratiche dovute alla derivata seconda:

\begin{equation}
\begin{aligned}
\left(\nabla^{2} \mathrm{E}\right)_{i j}=\frac{\partial \mathrm{E}}{\partial p_{j} \partial p_{j}} &=\frac{\partial}{\partial p_{i}}\left(e_{k} \frac{\partial e_{k}}{\partial p_{j}}\right) \\
&=\frac{\partial e_{k}}{\partial p_{i}} \frac{\partial e_{k}}{\partial p_{j}}+\cancel{e_{k} \frac{\partial^{2} e_{k}}{\partial p_{j} \partial p_{i}}}
\end{aligned}
\end{equation}

Ovvero considerare l’iterazione semplicemente come:

\begin{equation}
\underline{p}^{(k+1)}=\underline{p}^{(k)}+\underline{h}^{(k)}
\end{equation}

Dove l’incremento sarà:

\begin{equation}
\underline{h}^{(k)}=-\left.\left[\underline{\underline{J^{T}}} \underline{J}\right]^{-1} \cdot \underline{\underline{J}} \underline{\underline{e}}\right|_{\underline{p}}(k)
\end{equation}

Metodi

Dunque il problema di identificazione parametrica viene affrontato come un problema di minimizzazione che porta a trovare i parametri ottimali. Sono quindi riportate le due differenti analisi.

Minimizzazione

Nel primo metodo si utilizza la routine di Matlab fminsearch(). La risposta vera del sistema viene calcolata tramite la routine rlc_fun.m che si occupa, una volta fornitegli i coefficienti 𝑅, 𝐿, 𝐶 e l’ingresso, la risposta del sistema dinamica. Una volta generati i dati corrispondenti alla risposta vera del sistema è possibile procedere all’ottimizzazione.

Si richiama allora la procedura fminsearch() fornendogli la funzione obiettivo, il punto iniziale e le variabile necessarie per il calcolo della funzione obiettivo (ingresso e discretizzazione temporale). La funzione obiettivo è presente all’interno della routine obj_fun.m e si occupa di fare la differenza tra l’uscita vera e quella predetta. L’uscita predetta viene calcolata ad ogni iterazione sfruttando la procedura di rlc_fun.m andando a calcolare la risposta del sistema RLC con i parametri aggiornati di volta in volta mediante l’algoritmo di Nelder-Mead.

Effettuando un primo test, con rumore gaussiano, è possibile stimare i parametri.

Inoltre, è possibile estrarre alcune informazioni dal processo di ottimizzazione. In particolare, si vede come nonostante la risposta del sistema viene stimata correttamente (fig. 4), l’errore sui parametri è piuttosto alto (Tab. 1).

Non identificabilità strutturale

Il problema, per come è stato posto, presenta una formulazione sbagliata. Infatti, andando a calcolare la matrice di sensibilità si vedono alcuni problemi. La matrice di sensibilità non è altro che l’insieme delle derivate della misura nma rispetto al jmo parametro. Calcolando quindi numericamente lo Jacobiano è possibile fare alcune considerazioni. Per calcolarlo è sufficiente perturbare rispetto al jmo parametro e calcolarne la derivata come rapporto incrementale tra la risposta del sistema con ingresso perturbato e lo stato non perturbato rispetto la perturbazione:

\begin{equation}
J_{i}=\frac{y_{\mathrm{pert}}-y_{\mathrm{ref}}}{\mid \text { pert }_{\mathrm{i}} \mid}
\end{equation}

I risultati numerici mostrano come la matrice ha un numero di condizionamento molto alto 1.856.107 ma sembrerebbe avere rango massimo. Tuttavia, andando ad indagare i valori singolari si vede come uno dei tre è in realtà molto piccolo, abbastanza da poter essere considerato numericamente nullo. Nonostante il software percepisce un numero diverso da zero questo è in realtà nullo ed indica come nel problema c’è una dipendenza tra i parametri.

Questo è legato al fatto che nella risposta di trasferimento in realtà i parametri sono soltanto due, dati dal prodotto dei tre valori numerici legati ai componenti circuitali:

\begin{equation}
\begin{aligned}
&\theta_{1}=L C \\
&\theta_{2}=R C
\end{aligned}
\end{equation}

Parametri ottimi

Considerando quindi il problema nella sua versione ben condizionata è possibile procedere a stimare i parametri 𝜃1 e 𝜃2. In particolare è possibile analizzare il risultato ottenuto considerando diversi ingressi, mostrati in fig. 3.

Per farlo viene riorganizzato il codice in modo da racchiudere la generazione dei dati e l’ottimizzazione in due procedure che possono essere richiamate di volta in volta con i differenti ingressi. Con tutti gli ingressi testati l’errore rimane sotto il 2% pur mantenendo una certa variabilità tra le differenti iterazioni.

Gauss Newton

Un approccio differente prevede l’utilizzo dell’algoritmo di Gauss Newton. Si cerca la minimizzazione della funzione obiettivo per il tramite dell’iterazione in eq. (5).

È necessario calcolare, per procedere lungo le iterazioni, lo Jacobiano. Questo viene calcolato tramite la routine jacobian_fun.m. Tale routine prende in ingresso il vettore dei parametri, l’ingresso e il vettore dei tempi e si occupa di calcolare la derivata numerica delle misure (risposta stimata) rispetto al parametro. Per farlo va a perturbare la jma componente della risposta facendo riferimento al jmo parametro e ne calcola il rapporto incrementale. La procedura è analogo a quanto espresso per la matrice di sensibilità nell’eq. (12).

È possibile sfruttare gli stessi ingressi del caso precedente (fig. 3). I risultanti seguenti fanno riferimento ad un ingresso tipo gradino.

Per aiutare il metodo a convergere si aggiunge un parametro di rilassamento (𝛼) che permette di modulare l’incremento tra un’iterazione e la successiva:

\begin{equation}
\underline{p}^{(k+1)}=\underline{p}^{(k)}+\alpha^{(k)} \underline{h}^{(k)}
\end{equation}

Regolarizzazione

Introducendo invece parametro di regolarizzazione 𝜆 > 0 si ottiene sempre l’inversione di una matrice non singolare [5] e si evitano problemi di cattivo condizionamento.

\begin{equation}
\left[{\underline{\underline{J}}^T \underline{\underline{J}}}\right]^{-1} \longrightarrow\left[\underline{\underline{J^{T}}} \underline{\underline{J}}+\lambda^{(k)} \underline{I}\right]^{-1}
\end{equation}

Il problema ricade sulla scelta del parametro ottimo, se troppo grande può rovinare la velocità di convergenza del metodo mentre se troppo vicino allo zero può essere inefficiente e portare verso le singolarità dello Jacobiano.

Fissando il numero di iterazioni massimo è possibile confrontare come varia la convergenza (tempo di convergenza) al variare dei parametri di rilassamento e di regolarizzazione (fig. 6). Si considera il tempo di convergenza come il numero di passi necessari ad avere una soluzione valida, soluzioni tali da non dare il risultato entro il numero di passi fissati vengono considerate non convergenti.

Per il parametro di regolarizzazione non si ottengono miglioramenti in quanto il problema è già ben condizionato. Per il parametro di rilassamento invece si vede fortemente la sua influenza sulla velocità di convergenza del metodo.

Guess iniziale

Una notevole importanza è ricoperta anche dal punto iniziale. Le tecniche utilizzate finora considerano il punto iniziale come una quantità ottenuta scaldando i parametri reali con i seguenti fattori:

\begin{equation}
\begin{aligned}
&\theta_{1}^{0}=1.5 \cdot \theta_{1} \\
&\theta_{2}^{0}=0.8 \cdot \theta_{2}
\end{aligned}
\end{equation}

È possibile scalare ulteriormente questi valori e osservarne l’effetto. In particolare, scalando entrambi i parametri è evidente come allontanandosi troppo dalla soluzione vera il metodo non converge. I parametri vengono scalati di 0.01 ÷ 3× rispetto il valore iniziale precedentemente considerato.

Conclusioni

Entrambi i metodi sono validi e permettono di raggiungere il risultato stimando i parametri. È importante considerare il problema nella forma corretta, considerando i parametri effettivamente linearmente indipendenti e non i differenti elementi circuitali.

Il metodo di Nelder-Mead permette di evitare la derivazione ma porta ad errori medi più elevati. Il metodo di Gauss-Newton si mostra più preciso anche se può risultare computazionalmente più costoso.

Essendo entrambi metodi non lineari, risulta molto importante scegliere il punto iniziale altrimenti il metodo potrebbe non convergere.

Inoltre, per il metodo di Gauss Newton risulta fondamentale selezionare correttamente i parametri di rilassamento e regolarizzazione per evitare una divergenza del metodo o un costo computazionale più elevato di quanto necessario.

Disponibilità dei dati

Il codice è disponibile alla repository GitHub del progetto.

• Il codice per l’ottimizzazione con NelderMead a tre parametri è presente nel file par_sys_id.m
• Il codice per l’ottimizzazione con NelderMead a due parametri è presente nel file par_sys_id_two_param.m
• Il codice per l’ottimizzazione con Gauss-Newton è presente nel file par_sys_id_gauss_newton.m

Riferimenti

• [1] Pardis Ghafarian, Hamidreza Jamaati, and Seyed Mohammadreza Hashemian. “A Review on Human Respiratory Modeling”. en. In: (), p. 9.
• [2] Michael C. K. Khoo. Physiological control systems: analysis, simulation, and estimation. en. Second editon. IEEE Press series in biomedical engineering. Piscataway, NJ : Hoboken, New Jersey: IEEE Press ; Wiley, 2018. ISBN: 978-1-119-05533-4.
• [3] Jeffrey C. Lagarias et al. “Convergence Properties of the Nelder–Mead Simplex Method in Low Dimensions”. en. In: SIAM Journal on Optimization 9.1 (Jan. 1998), pp. 112147. ISSN: 1052-6234, 1095-7189. DOI: 10 . 1137 / S1052623496303470. URL: http : / / epubs . siam . org / doi / 10 . 1137 / S1052623496303470 (visited on 04/01/2022).
• [4] MathWorks Matlab. Generate input signals to support system identification. URL: https:// it . mathworks . com / help / ident / ref / idinput.html. [5] Jorge Nocedal and Stephen J. Wright. Numerical optimization. 2nd ed. Springer series in operations research. OCLC: ocm68629100. New York: Springer, 2006. ISBN: 978-0-38730303-1.