Modeling thermal systems with fractional models: human bronchus application

System thermal modeling allows heat and temperature simulations for many applications, such as refrigeration design, heat dissipation in power electronics, melting processes and bio-heat transfers. Sufficiently accurate models are especially needed in open-heart surgery where lung thermal modeling will prevent pulmonary cell dying. For simplicity purposes, simple RC circuits are often used, but such models are too simple and lack of precision in dynamical terms. A more complete description of conductive heat transfer can be obtained from the heat equation by means of a two-port network. The analytical expressions obtained from such circuit models are complex and nonlinear in the frequency ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document}. This complexity in Laplace domain is difficult to handle when it comes to control applications and more specifically during surgery, as heat transfer and temperature control of a tissue may help in reducing necrosis and preserving a greater amount of a given organ. Therefore, a frequency-domain analysis of the series and shunt impedances will be presented and different techniques of approximations will be explored in order to obtain simple but sufficiently precise linear fractional transfer function models. Several approximations are proposed to model heat transfers of a human middle bronchus and will be quantified by the absolute errors.


Introduction
Thermal modeling of systems is of particular interest in applications where temperature might be critical. This includes air conditioning, industrial refrigeration, electronic device cooling, biological tissue heat losses, etc. The complexity of heat transfer usually implies the use of finite element methods to solve the heat equation in a chosen region of space. However, finite element models tend to be complex to calculate and in many cases a simpler model that only takes into account temperatures around particular points of interest are precise enough for specific applications. Diffusive phenomena as well as the complex geometries found in biology can be approached through fractional calculus, which has become an important tool in bioengineering [18]. A typical simple model to take into account heat dissipation and its dynamics is the RC circuit. It is usually used in the domain of power electronics [16], building simulation [8,28] and even to model human heat losses [13,14]. RC circuit models are also present in other similar applications, such as the measurement of bioimpedances [5] or lithium-ion battery models [36]. In most of these applications, the system input is assumed a low-frequency signal. In order to widen the frequency range and especially in high frequency, the thermal twoport network was introduced in [20] and such a model is directly derived from the heat equation.
The thermal two-port network models heat conduction in a single direction as a T circuit (see Fig. 2) with two series impedances Z 1 (s) and Z 2 (s) and a shunt impedance Z 3 (s). These impedance expressions are complex and nonlinear in ω which do not allow obtaining rational transfer functions, the latter being more suited models for control design. It can be shown in high frequency that the thermal impedance of a plane wall is given by a half-order integrator [2,22]: Such an expression shows that the heat transfer in high frequency may be modeled by using constantphase elements (CPE) [15,26,34]. Constant-phase elements have particularly been used in biological and medical applications such as modeling intestine tissue [12], porous films [9], cardiac tissue [19] or even lung mechanics [6], thus showing the growing importance of fractional calculus in modeling applications.
During surgery, heat transfer and temperature control of a tissue may help in reducing necrosis and preserving a greater amount of a given organ. Patients can suffer from lung injuries during cardiopulmonary bypass, such as ischemia [3]. Perfusion techniques involving some type of indirect temperature control exist [21]. Some lung injuries may be treated by applying mild hypothermia [33]. The precision required to keep a fairly precise temperature in the tissue justifies our interest in enhancing the knowledge regarding the dynamics of thermal models for physiological scenarios.
The contributions of this paper lie on at giving different approximation models for both the series and shunt impedances. The series impedances will be approximated by means of an asymptotic model (taking into account only low and high-frequency behav- ior), zero-pole cells and the use of a fractional Butterworth impedance model. On the other hand, the shunt impedance will be approximated by a capacitance (the most basic model) and by filtered capacitance models: fractional slope filter and multiple pole filters. The poles for the last filter will be analyzed with and without a recursive factor imposed between the parameters. All these model approximations for the series and shunt impedances will be tested on a biological system: the heat transfer in a human bronchus.
This paper is organized as follows. The presentation of the thermal two-port network and its frequencydomain analysis will be presented in Sect. 2. The approximation propositions will be explored in Sect. 3 by providing simulation results. A biological application by modeling the human bronchus will be presented in Sect. 4. Finally, Sect. 5 presents conclusions, final remarks and research perspectives.

The analytical thermal two-port network
Let us consider the heat conduction on a simple plane wall in its longitudinal direction x, as shown in Fig. 1, where k, ρ and c, respectively, stand for the medium thermal conductivity, density and heat capacity. The total length is L, and the cross section is S w . The heat equation is expressed as: and as conduction only goes in x direction, it can be written as: The Laplace transform of equation (3), with null initial conditions, leads to: where s denotes the Laplace variable.  Fig. 2 Thermal two-port network By considering a heat flux input at x = 0 and a heat flux output at x = L, it then comes: or even, put under a matrix form: where with δ = s a and a = k ρc . Note that δ involves the presence of a fractional (noninteger) operator, as its time-domain representation is a half-order derivative. This matrix equation (6) can be represented by a T circuit model (see Fig. 2).
The series impedances can then be expressed by: and the shunt impedance by: where coth and csch, respectively, are the hyperbolic cotangent and cosecant functions.

Low-frequency behavior
Processes in thermal applications are usually assumed to be quasi-static, which means that the impedance behavior is only considered in low frequency, thus leading to: and lim ω→0 Z 3 ( jω) = 1 ρcL S w jω = 1 jωC t (11) where C t = ρcL S w . As can be seen, low-frequency heat transfers lead to a network consisting only of thermal resistances for the series impedances (which represent energy dissipation) and a capacitance for the shunt impedance (which models energy storage in the medium). This is coherent with the classic RC circuit model used to characterize heat transfer dynamics. However, it should be noted that, depending on the thermal properties of the medium, it is possible to have situations in which the transfer occurs at middle range frequencies, which limits the validity of the RC model.

High-frequency behavior
Now, by considering the limits in high frequency, the impedance analysis of the thermal quadrupole, by decoupling them into gain and argument, gives: The impedance Z 1 clearly exhibits a low-pass filter in high frequency, as shown by its zero gain in high frequency. This suggests a capacitive behavior, but the argument has to be analyzed with more care. Indeed, the high-frequency argument of this impedance is not −90 • , but its half. This allows us to deduce that the high-frequency behavior of the series impedance Z 1 is a constant-phase element of phase −45 • or even: At first, the shunt impedance Z 3 seems coherent with the capacitance model in low frequency, as it is normal for a capacitance to exhibit zero gain in high frequency. The argument, on the other hand, suggests that there is an additional filtering that appears in high frequency, which explains why this argument is not the expected The high-frequency model of Z 3 can be a combination of a capacitance and a filter, such as: Remark 1 Adding a correction filter will provide a wider validity domain for the approximation of Z 3 , but it is actually impossible to design a filter that will allow an infinite frequency range: Theoretically, such a filter would have an infinite number of poles in order to obtain the argument shown in Eq. (13).

Impedance approximations
It is impossible to well approximate the quadrupole by a model valid for any frequency just by using transfer function models, unless having an infinite number of parameters. However, a truncated transfer function model may be a fairly useful approximation with a wider validity compared to RC circuits. Consequently, a model going from low frequency (or even a static case) to a high but finite upper frequency will be developed. An academic example will be first used to illustrate the main ideas of this section, and a more complex example will be described in Sect. 4. The academic model parameters are presented in Table 1.
The proposed approximations will require the determination of different parameters in order to be as close as possible to the true impedance behavior. A common quadratic error criterion will be used by taking the error of the gain in d B defined as: where ω = [ω 1 , . . . , ω N ] is a frequency vector where ω i ∈ [0.1rad/s, 100 rad/s] for i = 1, . . . , N and λ is a weighting coefficient (0 < λ < 1) that provides a weight for low or high frequencies. If λ = 1, all frequencies have the same weight. Unless indicated, it will be taken as unity. This criterion (16) measures the error in the gain curve given in d B.
The series impedance approximations will be optimized by using the flower pollination algorithm (FPA) (see [35] for the complete description), which is a novel meta-heuristic optimization method. The parameters to be tuned in order to apply FPA are pretty low, which makes it faster to tune and test when compared to other algorithms. The main parameters to be tuned are the population size n (normally between 10 and 25) and a switching probability p which states the possibility of doing one of the two principal operations of the algorithm in an iteration: global pollination or local pollination. At each iteration, a random number between 0 and 1 is calculated and compared to switching probability p. One of the two following corrections is performed: -for a global pollination, the next population of possible solutions is corrected through the difference between the solution x t i and the current global best solution g best where stands for Levy flight distribution; -for a local pollination, the next generation of solutions n is obtained by using a correction term that takes the difference of two random solutions x t j and x t k multiplied by a random coefficient The most critical parameter is the switching probability, but applications found in the literature tend to use a value of 0.8 [35]. A mathematical study [17] stated that even though the value of p does not degrade the results when it comes to unimodal functions, a range between 0.5 and 0.8 is preferred. Further analysis on the influence of the switching probability p has been performed by even considering a time-variable switching probability [27]. Optimal switching probability was found to be constant p = 0.8. Therefore, all optimizations were performed by using n = 10 and p = 0.8 for this study. The main steps for classic FPA are summarized in Algorithm 1.

Remark 2
The choice of FPA as an optimization tool was mainly due to its simple structure and low parameter number to be tuned. However, it should be noted that any other optimization algorithm could also be used to approximate the parameters, such as genetic algorithm or particle swarm optimization.

Algorithm 1: Flower Pollination Algorithm
Initialize a population of n flowers with random solutions Find the best solution g best in the initial population Define a switch probability p while t < Max I terations do for i = 1 : n do if rand < p then Global Pollination with equation (17)  3.1 Approximation of the series impedances Z 1 and Z 2

Asymptotic approximation
As it was seen in the previous section, the series impedance Z 1 and Z 2 are equal and behave like a resistance for low frequencies and like a constant-phase element at high frequencies. For simplicity sake, only Z 1 is used. By combining both asymptotic behaviors, a simple approximation of Z 1 is now proposed as the parallel combination of a resistance and a fractance, which leads to: This first proposition will be called asymptotic approximation, as it is obtained by only considering both extremes of the frequency domain. The cutoff frequency obtained at −3 d B under the static gain corresponds to (see Appendix B): Figures 3 and 4, respectively, show the gain and phase for the approximation Z 1−asymp compared with the exact expression of Z 1 over a frequency range [0.1, 100] rad/s. As expected, both curves get closer toward the limits in low and high frequencies, however, high gain and phase errors occur in the middle range frequencies. Another important error by using the approximation Z 1−asymp , is the cutoff frequency. The true −3 d B frequency under the static gain is almost one decade greater than the one obtained with the asymptotic approximation, which means this approximation leads to a slower dynamic.

Contribution with zero-pole approximation
In order to improve the approximation Z 1−asymp , a correction term is proposed to be added to the asymptotic approximation. This correction term will adjust the impedance for a mid-band frequency range, as the asymptotic term can properly handle both low and high frequency limits. This means that the additional term should not contribute to the gain or to the phase in the low-or the high-frequency range. A first proposition for the corrected impedance is to use zero-pole cells that will be placed inside the frequency range of interest: The parameter vector θ is given by: where p and z, respectively, are vectors containing all the poles and zeros. The number of cells N cells , used for the correction, will depend on each specific case. Each cell adds a single pole and zero to the transfer function Z 1 . A large number of cells may give pretty accurate approximations, but the thermal model may become too complex and identifiability problems may occur (large number of parameters in the transfer functions).
In order to determine an optimal number of cells to add for the zero-pole approximation, the error criteria J (θ ) were optimized for a number of cells N cells in order to estimate an optimal number of poles and zeros. However, the choice of the optimal number of cells may not be trivial as is well known in optimization that the higher the number of parameters, the lower the quadratic criterion J (see Fig. 5). In order to have a more objective criterion, the novel PIC criterion proposed in [29] is used and is defined as: where SS E is the sum of squared errors, k the number of parameters and N the number of datapoints. The PIC criterion evolution is plotted in Fig. 6.
The PIC criterion provides an optimal number of zero-pole cells for N cells = 3. Even though this technique might succeed in fitting the true frequency response, the dimension of the parameter vector of the correction term should be taken into account. In this case, dim(θ ) = 6. If this thermal impedance is only a section of a global thermal model, it is possible that an approximation of this type will lead to extremely complex expressions for global transfer functions in thermal systems.

Contribution with fractional Butterworth approximation
By observing the gain plot in Fig. 3, another approximation can be made for the thermal impedance Z 1 . The true frequency response (blue curve) is considerably flat before the cutoff frequency. This suggests an impedance that behaves like a Butterworth filter [4]. However, the high-frequency gain slope is not an integer multiple of −20 dB/dec, but −10 dB/dec. Therefore, inspired by [1], a fractional Butterworth filter can be used as an alternative approximation model: This new type of model is not a correction term multiplied by the asymptotic approximation, as it can include both the low-and high-frequency behaviors within its expression: A term-by-term identification gives: For this type of approximation, the parameter vector will always be limited to 3 parameters: The fractional Butterworth filter may be a simpler solution with respect to the parameter vector dimension. However, the stability of a system with such a transfer function should be discussed a priori. The BIBO stability for fractional commensurate transfer functions was established by Matignon [23]. An extended criterion for incommensurate systems was later developed [31]. In [1], it is stated that parameters a and b need to be equal or inferior to zero in order to guarantee the fractional stability of the system. On the other hand, the differentiation orders α and β are positive and their sum is 0.5. They are therefore both bounded between 0 and 0.5. The parameters obtained by using the fractional Butterworth filter are given in Table 2.
It is interesting to observe that parameter a is zero, which simplifies even more the transfer function of the impedance. For this particular case, the obtained thermal impedance is of the form: or even, its thermal admittance is: From this last expression, impedance Z 1 may actually be reinterpreted as being the parallel combination of a half-capacitance, a resistance and a fractional-order new element which can be expressed as: Further studies are necessary in order to analyze and give a possible physical meaning to this capacitancelike expression, particularly its negative sign.

Comparison of the approximations of the series impedance Z 1
In order to compare the previous approximations, the absolute error is defined as:  Figure 7 illustrates the gains for the exact case and the proposed approximations of Z 1 , whereas the error signals are plotted in Fig. 8.
The minimum, maximum and mean values of the error obtained for the three approximations are provided in Table 3.
Moreover, the mean error for the asymptotic approximation is almost 3 dB and the lowest error is around 1 dB, which translates an inaccurate approximation for all mid-range frequencies. When zero-pole cells are added, the mean error is reduced by a factor of 35, and the error signal can get really close to zero, as shown by its minimum error value (which is around 10 −4 ). The Butterworth approximation offers intermediate results, as there is a rough factor of 4 between the asymptotic and the Butterworth mean error. Its peak is still acceptable, but 4 times higher than the one obtained with the zero-pole approximation. If one needs an extremely precise result, the complexity of the zero-pole cells might be adequate, but the simpler Butterworth structure may offer accurate enough results for many applications without adding too much parameters.
3.2 Approximation of the shunt impedance Z 3

Capacitance approximation
The shunt impedance Z 3 behaves like a pure capacitance for low frequencies, but its high-frequency behavior is more complex than that of the series impedance, namely: with C t = k L S w a = 1 J · K −1 . The gain frequency response of Z 3 and its capacitance approximation are shown in Fig. 9.
The −20 dB/dec slope, seen in low frequencies, confirms the capacitance behavior of Z 3 . However, the slope of the gain curve increases as the frequency increases. Even though a capacitance tends to behave like a short-circuit on high frequencies, the slope of the gain is not increased and the error curve increases indefinitely in high frequency (see Fig. 10).

Contribution with a fractional slope approximation
If the true thermal system is excited in the mid-range frequencies, it is suggested to modify the classic capacitance model to include a slope increment.
The difference between the initial and the final slopes is of almost −60 dB/dec. In order to correct this error, a simple fractional filter is proposed as a multiplicative term to the low-frequency capacitance: This correction term by using a fractional transfer function increases the high-frequency slope magnitude by the factor φ, namely −20φ dB/dec. 1/τ indicates the transitional frequency from where the correction term begins to influence. An initial estimate for φ would be 2.90 in order to be close to the final slope of −78.28 dB/dec around the higher frequency. The parameter vector for this case is: (39) The correction filter term in Eq. (38) has the following form: a filter that reminds a commensurate transfer function with order α and Matignon's stability theorem may be applied [23], namely: where for λ = s α , λ i are the roots of the characteristic polynomial in λ.
For the general structure presented in Eq. (40), we have: (42) Therefore, for this type of filter, the stability condition reads α < 2 for b > 0.

Contribution with multiple fractional slope approximation
According to Matignon's stability theorem, the fractional order α being limited to 2, it is not possible to go beyond this limit and to get a more suited highfrequency slope with such a model (40). Therefore, a multiple fractional slope filter may also be proposed: The number of cells N m depends on the required additional slope. As each cell may add up to −40 dB/dec, in this case N m = 2.

Contribution with reduced pole approximation
Looking at the gain curve and its slope in high frequency, an alternative correction term could be the addition of poles to take into account the slope increment.
The difference of slopes in the frequency band of interest, namely: · denoting the nearest integer function, provides the number of poles to be added. A new approximation of the shunt impedance can be expressed as: with parameter vector: If the required number of poles is low enough, one could simply use this approximation and take each parameter τ i as an additional variable to approximate in order to get the correction term. However, there may be cases in which the number of additional poles might be too high and its determination will be too complex in terms of direct optimization (local minima convergence problem). In this case, an arbitrary relation such as a recurrent relation may be established between poles as a way to reduce the number of parameters: where all the poles are related to the first pole τ 1 .
The main advantage of this correction is that the parameter vector will always be limited to the following: The parameter space will always be 2-D in this case. On the other hand, the main drawback is the loss of flexibility with respect to the general pole addition as detailed in paragraph Sect. 3.2.3.

Comparison of the approximations of the shunt impedance Z 3
The approximations obtained by using the fractional slope, the multiple fractional poles and the reduced poles are compared to the pure capacitance model. Note that to help the convergence of the FPA algorithm for a better fit in high frequency, the weighting factor λ was fixed to 0.95 in the criterion (16) for the fractional and multiple fractional slope models. Figure 11 illustrates the gains for the exact case and the proposed approximations of Z 3 , whereas the error signals are plotted in Fig. 12.
The minimum, maximum and mean values of the error obtained for the four approximations are provided in Table 4.
Even though there is an error that may go up to almost 4.16 dB by using the fractional slope model, the error is clearly reduced when compared to the pure capacitance model. For this example, the fractional order was optimized to φ = 1.45, which is a rather surprising result. The optimization did not lead to a sat-  Fig. 12).
Concerning the multiple fractional slope model, the optimization have led to better results. There is roughly a factor of 3 between the mean errors of this model and the simple fractional slope one. The commensurate order was found to be ν = 1.13, which means a highfrequency order of 2ν = 2.26. The additional flexibility provided by the multiple cell fractional filter permits to go beyond the limit of 2 such as imposed by Matignon's stability criterion in order to get a better fit.
The reduced pole approximation provides even more accurate results and the dimension of the research space for the optimization stays the same. There is a dramatic difference at the point in which the approximation starts to modify the curve, as the recursive pole gave a first pole located at τ 1 = 0.033 rad/s, which is two decades below the optimal value obtained by the fractional slope. This approximation may allow a more predictive approach to error reduction than the fractional slope and multiple-cell fractional filter.

Comparison of the reduced pole approximation versus the all-pole approximation
In view of the previous results, let us consider the latter approximation given in Eq. (49). Its optimized model can be used as an initial starting point to perform a further optimization: Instead of using reduced poles (which indeed has the advantage to reduce the number of parameters in the optimization procedure), this constraint defined in Eq. (47) can be released so that all the parameters can be optimized independently. Such an approximation of the shunt impedance can be expressed as: with the more general parameter vector: where all τ i are independent and gathered in τ . The FPA algorithm was performed with the reduced and all poles by using a gradient descent. As expected, both approximations give really good results. Figure 13 illustrates the comparison of the associated errors for the reduced and the all-pole approximation cases.
It is indeed difficult to distinguish a further improvement for this new case by allowing the poles to be located at any position. Table 5 shows the pole locations as well as the criterion J for both cases.
Even if there is a slight improvement regarding the criterion J , the non-recursive approximation increased the vector parameter size (three times more than in the recursive case) and the criterion improvement proved to be very slight. On the other hand, it can be seen that all pole approximation locations are almost the same as the ones with the reduced pole approximation. This means that, for this case, the optimal pole location for the approximation was really close to the results provided by the reduced pole case. It may also justify the use of a reduced pole model in more complex cases.

Application example: heat transfer in a lung bronchus
A potential application for these approximations is to model heat transfers in the human body. During surgery, heat transfer and temperature control of a tissue may help in reducing necrosis and preserving a greater amount of a given organ. Patients can suffer from lung injuries during cardiopulmonary bypass, such as ischemia [3]. Perfusion techniques involving some type of indirect temperature control exist [21]. Some lung injuries may be treated by applying mild hypothermia [33]. The precision required to keep a fairly precise temperature in the tissue justifies our interest in enhancing the knowledge regarding dynamic of thermal models for physiological scenarios.
For the application, the heat transfer is considered in the right middle bronchus. Under normal conditions, the nasal cavity and larynx play a significant role regarding air heating and humidification (see [10,11]) and procedures impairing this conditioning may change thermal conditions inside the lungs. Human bronchi may be exposed to thermal stress in such conditions. For this analysis, an intermediate length of L = 0.0236 m was chosen as it is a mean value for right middle bronchus length [25] and radius r ≈ 1 mm (see [24]).
In terms of frequency, even though human respiration is not a naturally fast mechanism, it may have nonnegligible dynamics for heat transfers. Human breathing can go up to 20 breaths per minute in normal conditions, this breathing being around 12 breaths per minute for an adult in average. Slow breathing is considered in the range [0.07, 0.16] Hz (see [32]), and above 14 breaths per minute (or 0.23 Hz), such breathing can already be considered for some as an abnormal value (see [7]). Therefore, the optimization will be carried out in the frequency range [0.01, 5] rad/s in order to go slightly over a normal breathing condition. For the series impedance Z 1 , the optimal number of cells for the zero-pole correction was found to be 6, as shown in Fig. 14. The gain diagrams of the different approximations proposed in paragraph Sect. 3.1 (asymptotic, zero-pole and fractional Butterworth approximations) and their associated errors are, respectively, shown in Figs. 15 and 16. Table 6 summarizes a comparison between the obtained errors. As expected, the asymptotic approximation offers a significant error for the whole frequency range, which may indicate that a human bronchus cannot be considered in the low-frequency limit as human breathing cycles are in the mid-band frequency range. Once  again, it can be seen that the zero-pole approximation offers the most accurate results. However, for this scenario, the improvement obtained by using the Butterworth is more significant than in the academic case. For this case, the zero-pole approximation has 8 parameters, which includes 5 additional parameters with respect to the Butterworth case. However, the Butterworth impedance still provides an improved estimation than the asymptotic one. It can be seen in Appendix B that for this case, the Butterworth impedance was once again reduced to the form shown in equation (31), which offers further advantages in terms of parameter space (2D).

Shunt impedance Z 3
The gain diagrams of the different approximations proposed in paragraph Sect.3.2 (capacitance, fractional slope, multiple fractional slope and pole approximations) and their associated errors are, respectively, shown in Figs. 17 and 18. Table 7 summarizes a comparison between the obtained errors.  The pure capacitance model provides the worst approximation: The mean error is of 7.58 dB in the considered frequency range, and it can go up to 42.05 dB as we get closer to the higher frequencies. For this particular geometry, three integer poles were required to provide a good approximation. It should be noted that for a frequency of 1 rad/s (close to normal breathing), the slope of the Z 3 gain curve is already −38 dB/dec, which means that a normal breathing is already in its mid-band frequency range. The best results are obtained with the reduced pole approximation with a mean error of 0.14 dB.

Time simulations
Consider a simple scenario for the bronchus model with its end insulated at x = L, which means φ out = 0W. A sine temperature variation is applied at the input at x = 0: with T 0 = 1K and ω 0 = 0.45 rad/s. Such a frequency ω 0 is coherent with a normal human breathing frequency. The transfer function which relates input and output temperatures is given by: By combining approximations of Z 1 and Z 3 , five different models were tested and the different combinations are shown in Table 8. The results of the simulation are shown in Fig. 19. Note that for the simulation of the fractional models, the Grünwald-Letnikov definition is used [30]: where p = d dt denotes the differential operator, and T s = 0.33s is the sampling time. The error terms are proportional to the sampling time [30,Section 7.4]. Consequently, the smaller the sampling time, the lesser the approximation errors.  Steady-state peak-to-peak values for each model were compared to the exact response. Errors are summarized in Table 9.
Model A exhibits the largest gain difference with respect to the exact solution due to lack of accuracy for both used impedances. The other four models show improved gain accuracy. As expected, models B and D are, respectively, outperformed by models C and D thanks to the zero-pole correction term in Z 1 . On the other hand, one can also see that multiple fractional slope leads to more accurate results than using a single one by comparing models C and E. The main drawback seen is that the approximations do not take into account the thermal lag present in the exact solution. This can be explained as only gain was optimized and lag phenomena are only visible in the phase. Optimizing both gain and phase could be studied for future works.

Conclusion and final remarks
Different propositions for approximating thermal impedances were proposed and compared in terms of accuracy in the frequency band of interest. The series impedance of a plane wall proved to have an interesting behavior in the frequency domain, as its usual simplification as a simple resistance is inaccurate beyond the low-frequency spectrum. The presence of a halfcapacitance indicates a possible additional thermal accumulation that cannot be taken into account by a simple capacitance model. Even if the resistance-half capacitance parallel model succeeds in approximating the frequency limits, there is a significative error in the mid-band frequency that can be reduced by using zero-pole cells or a fractional Butterworth expression. The zero-pole approximation proved to be better in terms of accuracy, but the impedance expression is complex and the number of parameters required for an optimal approximation may be too large. The zeropole cell approximation is an exclusively mathematical construction and lacks physical meaning. The fractional Butterworth approximation may not offer such accurate results, but it still greatly improves the error when compared to the resistance-half-capacitance parallel model. The parameter vector is limited to a maximum dimension of 3, and it is suggested that this type of expression may actually be reinterpreted in a circuit model as the parallel addition of elements to the resistance-half capacitance model. Additional physical meaning may be obtained from this model in future studies.
The shunt impedance has a more complex frequency behavior and only a low-frequency analysis was done. In low-frequency, this impedance behaves like a capacitance. However, the quick slope decreases (both in gain and phase) in high frequency, a circuit element nor a simple fractional model is sufficient. This is the reason why the approximations presented for Z 3 are limited to a specific frequency range. The fractional slope helped in reducing the errors and has little parameters to be determined. However, the approximation with recursive poles proved to be more accurate without adding more parameters to the approximation. For the academic example, the accuracy was not significantly improved by allowing the poles to be independent. A simple recursive pole approximation gives accurate enough results and provides a simple expression. However, it should be noted that Z 1 and Z 3 are inseparable elements of a same system. Therefore, equivalent transfer functions obtained from the T circuit will involve fractional-order derivatives as it cannot be avoided for the series impedances Z 1 and Z 2 .
All the proposed approximations for modeling the series and shunt impedance were applied in a simulation example of a human bronchus. As previously stated, the simulation results well illustrate that the zero-pole approximation well fits the series impedance Z 1 and Z 2 . For the shunt impedance Z 3 , it could be fitted by using a single integer-order pole.
Research perspectives of this study include analyzing different types of scenarios with the thermal twoport networks and its transfer function properties when applying the proposed approximation techniques. The possibility of getting transfer function models with physically meaningful parameters may improve the analysis of system identification for thermal systems. On the other hand, time simulations showed that the inclusion of a thermal lag may be required to fully emulate the exact solution of the heat equation. Another perspective could be to analyze impedance approximations by means of a two-step optimization in order to also take into account the phase.
Author contributions In relation to your submission: J.F. Duhé has developed the results with S. Victor. He also made the simulations. S. Victor, P. Melchior, Y. Abdelmoumen and F. Roubertie have contributed to the biological and physiological results.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.

A Impedance expressions from the change of a matrix formulation to a T circuit
For a T network as shown in Fig. 2, the input-output relationship may be expressed by means of a transmission matrix: The transmission matrix can be expressed in function of the impedances Z 1 , Z 2 and Z 3 : Therefore, impedances: which lead to expressions in Eqs. (8) and (9). It should be noted that Z 1 = Z 2 for this specific geometry as A = D, but it should not be taken as a general assumption.

B Proof of fractional first-order system bandwidth
The frequency response of the asymptotic approximation for the series impedance is: from which one gets the gain: ; it can be seen that the highest gain in d B for this system is: The −3 dB frequency is given by: By replacing Eqs. (57) in (59), one gets the following polynomial: By taking the single positive and real-valued solution for ω b : which is the expression for the bandwidth.

C Estimated transfer functions for the academic example
See Table 10.

D Estimated transfer functions for human bronchus
See Table 11.