A linear model of lung mechanics is analyzed in the following report. The model corresponds to an RLC circuit with unknown parameters. Parametric identification allows to trace the system governance coefficients once a mathematical model has been estimated. Different approaches allow to minimize the difference between the estimated parameters and the real measures. An approach without numerical differentiation is addressed, using the Nelder-Mead algorithm and an approach using the relaxed Gauss-Newton method. Both methods lead to the correct estimation of the parameters but some considerations allow to optimize their performance.
The following report is an extract from a project carried out for the course of Modeling and Simulation of Physiological Systems, held for Medical Engineering at the University of Rome Tor Vergata.
Introduction
Respiratory mechanics is described in the literature by numerous models. Among these, a valid model under certain conditions is represented by an RLC circuit. The model allows to describe the respiratory mechanics of a healthy individual.
Starting from the measurements of the individual, it is possible to estimate the parameters to represent it through the estimated model. This process, under the name of parametric optimization, allows to adapt the general system to the specific case. The following analysis deals with the parametric identification technique by minimizing the difference between the measured response and the estimated one.
Different optimization techniques allow to estimate the parameters that allow to describe the considered system. Both the Nelder-Mead method and the Gauss-Newton method are described, analyzing their commonalities and differences. Some considerations are also made on the convergence of the method and on the correct conditioning of the problem.
Background
The identification of systems is a set of techniques that aims to identify the parameters and the system to describe a certain phenomenon. Therefore the first major distinction sees the separation of two techniques: parametric and non-parametric identification. The first part by estimating a model of the phenomenon and aims to identify the parameters that represent it. Non-parametric identification, on the other hand, proceeds without knowing the model of the phenomenon, at most it estimates some of its characteristics, and aims to identify the system that best describes it.
In this case we start from a well known model, a linear lung mechanics model described by an RLC circuit, and the descriptive parameters are identified.
RLC system
The RLC circuit (fig. 2) allows, by following the electrical analogy, to describe an approximate model of respiratory mechanics. The model predict a series of resistance, inductor and capacitance. Based on the analogy between electric current and pulmonary airflow, it is possible to represent the compliance and elastic resistance components of the lungs into a lumped parameter model [1].
This model allows to face the physiological cases in which the flow entering the lungs is predominant compared to the blocked flow (flow which instead increases considerably in the case of obstructive diseases). The analogy sees the association of the pressure difference to the potential difference, the current to the air flow, and the inductance as an inertance, that is, the pressure difference required to cause a unit change in the rate of change of flow in the time.
In a specific patient case, experimental measures would be necessary in order to obtain the real model on which to carry out the identification. In this case, since no experimental data are available, data known in the literature are used [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]\\A single resistance and a single compliance are considered, which represent the overall properties and the ability to accumulate air in the respiratory system as a whole. Therefore the resistance R brings with it the resistive contribution of the airways, lung tissues and chest wall. Capacity C represents tissue and chest wall compliance.
Kirchoff law
The goal is to predict the alveolar pressure 𝑝𝐴 and its dynamic response to different pressure waveforms applied to the opening of the airways 𝑝𝑎𝑜.
Applying Kirchoff’s first law to the circuit in fig. 2 we obtain:
\begin{equation}
p_{a O}=L \frac{d i}{d t}+R i+p_{a}
\end{equation}Where is present the following relation:
\begin{equation}
i=C \frac{d p_{a}}{d t}
\end{equation}From which we derive the transfer function:
\begin{equation}
H(s)=\frac{P_{a}(s)}{P_{a o}(s)}=\frac{1}{1+s R C+s^{2} L C}
\end{equation}This shows how 𝑝𝑎 is entirely dependent on 𝑝𝑎𝑜 and the system is an open-loop configuration system.
Parameters identification
Since specific patient data on respiratory mechanics are not available, the data will be generated through the rlc_fun.m procedure which returns the response of the RLC circuit providing the system parameters in inputs. Noise can also be added to this.
Nelder-Mead algorithm
Once the transfer function of an RLC circuit has been assumed as a model of the system, to estimate the parameters it is necessary to define an objective function. This function represents the difference between the measured (true) response of the system and the estimated one:
\begin{equation}
\underline{e}=\underline{y}^{\mathrm{pred}}-\underline{y}^{\mathrm{mis}}
\end{equation}And the search for the correct parameters aims at minimizing the objective function defined as:
\begin{equation}
\mathrm{E}=\frac{1}{2}\sum \|\underline{e}\|^{2}
\end{equation}One of the possible ways to deal with this problem is to rely on the Matlab routine known as fminsearch(). This procedure takes care of returning the optimal parameters. In particular, this function uses the Nelder Mead algorithm [3]. This method does not use derivatives and is based on the concept of a simplex by approximating the local optimal point in 𝑛 variables.
Gauss-Newton method
An alternative is the Gauss-Newton method. The method, unlike the previous one, requires the derivatives to be taken into account as well. The minimization of eq. (5) is a nonlinear least squares problem. The numerical solution uses the Gauss Newton method. Like Newton’s method, stationarity is imposed, a condition that is reflected in the jth component:
\begin{equation}
(\nabla \mathrm{E})_{j}=\sum_{k} e_{k} \cdot \frac{\partial e_{k}}{\partial p_{j}}
\end{equation}That is, this differentiation corresponds to premultiplying by the sensitivity matrix J:
\begin{equation}
\nabla \mathrm{E}=\underline{\underline{J^{}}}^{T} \:\underline{e}
\end{equation}Using Newton’s algorithm, the iteration at step 𝑘 + 1 will search for parameters (𝑝) such as:
\begin{equation}
\underline{p}^{(k+1)}=\underline{p}^{(k)}-\left.[\nabla(\nabla \mathrm{E})]^{-1}(\nabla \mathrm{E})\right|_{\underline{p}^{k}}
\end{equation}For the approximation introduced by Gauss [5] it is possible to neglect the quadratic components due to the second derivative:
\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}That is, consider the iteration simply as:
\begin{equation}
\underline{p}^{(k+1)}=\underline{p}^{(k)}+\underline{h}^{(k)}
\end{equation}Where the step will be:
\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}Methods
Therefore the parametric identification problem is approached as a minimization problem that leads to finding the optimal parameters. The two different analyzes are then reported.
Minimization
The first method uses the Matlab fminsearch() routine. The true response of the system is calculated using the rlc_fun.m routine which deals with the dynamic system response once the coefficients 𝑅, 𝐿, 𝐶 and the input have been provided. Once the data corresponding to the true response of the system has been generated, it is possible to proceed with optimization.
The procedure fminsearch() is then called, providing it with the objective function, the starting point and the variables necessary for the calculation of the objective function (input and time discretization). An so the objective function is present within the obj_fun.m routine and takes care of making the difference between the true output and the predicted one. The predicted output is calculated at each iteration using the rlc_fun.m procedure by calculating the response of the RLC system with the parameters updated from time to time using the Nelder-Mead algorithm.
By carrying out a first test, with Gaussian noise, it is possible to estimate the parameters.
In addition, some information can be extracted from the optimization process. In particular, it can be seen that despite the system response being correctly estimated (fig. 4), the error on the parameters is quite high (Tab. 1).
| Parameter | True | Initial | Estimate | Relative error |
|---|---|---|---|---|
| R | 0.100 | 0.150 | 0.119 | 0.197 |
| L | 0.010 | 0.008 | 0.012 | 0.206 |
| C | 0.100 | 0.280 | 0.083 | 0.169 |
Structural non-identifiability
The problem, as it has been posed, is incorrectly formulated. In fact, when calculating the sensitivity matrix, we see some problems. The sensitivity matrix is nothing more than the set of derivatives of the measure nth with respect to the jth parameter. Therefore, by calculating the Jacobian numerically, it is possible to make some considerations. To calculate it, it is sufficient to perturb with respect to the jmo parameter and calculate its derivative as an incremental ratio between the response of the system with perturbed input and the state not perturbed with respect to the perturbation:
\begin{equation}
J_{i}=\frac{y_{\mathrm{pert}}-y_{\mathrm{ref}}}{\mid \text { pert }_{\mathrm{i}} \mid}
\end{equation}The numerical results show how the matrix has a very high conditioning number 1.856.107 but it would appear to have maximum rank. However, investigating the singular values we see that one of the three is actually very small, enough to be considered numerically null. Although the software perceives a number other than zero, this is actually null and indicates that there is a dependence between the parameters in the problem.
This is linked to the fact that in the transfer response there are actually only two parameters, given by the product of the three numerical values linked to the circuit components:
\begin{equation}
\begin{aligned}
&\theta_{1}=L C \\
&\theta_{2}=R C
\end{aligned}
\end{equation}
Optimal parameters
Considering therefore the problem in its well conditioned version it is possible to proceed to estimate the parameters 𝜃1 and 𝜃2. In particular, it is possible to analyze the result obtained by considering different inputs, shown in fig. 3.
To do this, the code is reorganized so as to enclose the data generation and optimization in two procedures that can be recalled from time to time with the different inputs. With all the inputs tested, the error remains below 2% while maintaining a certain variability between the different iterations.
Gauss Newton
A different approach involves the use of the Gauss Newton algorithm. We seek the minimization of the objective function through the iteration in eq. (5).
It is necessary to calculate the Jacobian to proceed through the iterations. This is calculated through the jacobian_fun.m routine. This routine takes in input the vector of the parameters, the input and the vector of the times and takes care of calculating the numerical derivative of the measures (estimated response) with respect to the parameter. To do this, it perturbs the jma component of the response by referring to the jmo parameter and calculates its incremental ratio. The procedure is similar to what is expressed for the sensitivity matrix in Eq. (12).
It is possible to use the same inputs as in the previous case (fig. 3). The following resultants refer to a step type input.
To help the method converge, a relaxation parameter (𝛼) is added that allows you to modulate the increase between one iteration and the next:
\begin{equation}
\underline{p}^{(k+1)}=\underline{p}^{(k)}+\alpha^{(k)} \underline{h}^{(k)}
\end{equation}Regularization
By introducing instead the regularization parameter 𝜆> 0, the inversion of a non-singular matrix [5] is always obtained and problems of bad conditioning are avoided.
\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}The problem falls on the choice of the optimal parameter, if too large it can ruin the convergence speed of the method while if it is too close to zero it can be inefficient and lead to the Jacobian singularities.
By setting the maximum number of iterations it is possible to compare how the convergence (convergence time) varies as the relaxation and regularization parameters vary (fig. 6). The convergence time is considered as the number of steps necessary to have a valid solution, solutions that do not give the result within the number of steps set are considered non-convergent.
No improvements are obtained for the regularization parameter as the problem is already well conditioned. For the relaxation parameter, on the other hand, its influence on the convergence speed of the method is strongly seen.
Initial guess
Considerable importance is also covered by the starting point. The techniques used so far consider the starting point as a quantity obtained by heating the real parameters with the following factors:
\begin{equation}
\begin{aligned}
&\theta_{1}^{0}=1.5 \cdot \theta_{1} \\
&\theta_{2}^{0}=0.8 \cdot \theta_{2}
\end{aligned}
\end{equation}It is possible to further scale these values and observe their effect. In particular, by scaling both parameters it is evident that by moving too far from the true solution the method does not converge. The parameters are scaled by 0.01 ÷ 3 × with respect to the initial value previously considered.
Conclusions
Both methods are valid and allow you to achieve the result by estimating the parameters. It is important to consider the problem in the correct form, considering the parameters that are actually linearly independent and not the different circuit elements.
The Nelder-Mead method avoids the derivation but leads to higher average errors. The Gauss-Newton method is more accurate even if it can be computationally more expensive.
As both methods are non-linear, it is very important to choose the starting point otherwise the method may not converge.
Furthermore, for the Gauss Newton method it is essential to correctly select the relaxation and regularization parameters to avoid a divergence of the method or a higher computational cost than necessary.
Data availability
The code is available in the project’s GitHub repository.
References
- [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 identifica
