Modeling and Estimation of Thermal Flows Based on Transport and Balance Equations

Heat transfer in counter ﬂ ow heat exchangers is modeled by using transport and balance equations with the temperatures of cold ﬂ uid, hot ﬂ uid, and metal pipe as state variables distributed along the entire pipe length. Using such models, boundary value problems can be solved to estimate the temperatures over all the length by means of measurements taken only at the boundaries. Conditions for the stability of the estimation error given by the di ﬀ erence between the temperatures and their estimates are established by using a Lyapunov approach. Toward this end, a method to construct nonlinear Lyapunov functionals is addressed by relying on a polynomial diagonal structure. This stability analysis is extended in case of the presence of bounded modeling uncertainty. The theoretical ﬁ ndings are illustrated with numerical results, which show the e ﬀ ectiveness of the proposed approach.


Introduction
The mathematical modeling of heat exchangers is really important since these devices are employed in chemical processes, petroleum refineries, power systems, and heating/refrigeration of air-conditioning plants [1,2]. In this paper, we address the modeling of a counterflow heat exchanger based on transport and balance PDEs and use such models for the purpose of monitoring the thermodynamic process by solving boundary value problems that exploit only few temperature measurements at the boundaries. A rigorous stability analysis is addressed to prove the effectiveness of the resulting temperature estimates from the theoretical and numerical points of view in line with the widespread methods based on the Lyapunov approach in a number of different applications [3][4][5][6].
Heat exchangers for industrial processes allow recovering lost energy from hot fluid streams or to heat a cold fluid for the purpose of air conditioning. Finite-dimensional models of such processes are used for fouling detection by using Kalman filtering methods [1,2]. As compared with these models, the proposed modeling framework is in an infinite dimension since it relies on hyperbolic PDEs. Based on such equations, first we will focus on a model with the dynamics of the temperatures of cold fluid, hot fluid, and metal pipe. In addition, a reduced model is considered with only the temperatures of cold and hot fluid as state variables in line with [7]. Since such temperatures cannot be determined by direct, distributed measurements, it is fundamental to estimate them by using a few numbers of measures. Toward this end, using such models, we address the solution of boundary value problems to estimate the distributed state, which provides such estimate by having at our disposal measurements of temperature only at the boundaries. We will refer to them also as boundary-output Luenberger observers or simply state observers (see [8] and the references therein). Stability conditions on the estimation error (i.e., the difference between the temperatures and their estimates) will be presented to construct such estimators in case of perfect modeling and assuming bounded uncertainty. In such a case, we will rely on the notion of quadratic boundedness (QB, for short) [9]. Finally, numerical results will be given to showcase the theoretical findings.
The stability of the estimation error is analyzed in two different settings. First, we will consider perfect modeling and hence asymptotic stability results will be provided. Next, we will assume the presence of additive noises, which account for bounded modeling uncertainties and stability is studied by using QB. Generally speaking, QB allows dealing with positively invariant sets and derive upper bounds on the trajectories of the state of a dynamic system subject to bounded disturbances [9]. Thus, QB has been successfully used for both output feedback control [10] and state estimation [11].
Here, we will extend the use of QB for the purpose of estimation in the considered infinite dimensional context. It is worth noting that, in lieu of searching for exact solutions of the flow equations (see, e.g., [12,13]), here we deal with the construction of boundary-output observers with guaranteed stability properties on the estimation error.
All the stability conditions that will be presented demand the selection of nonlinear Lyapunov functionals. In this paper, we will address this task by using the SOS (sum-ofsquares) approach [14,15], which was originally developed for the analysis and control of polynomial systems described by ODEs (see [16] and the references therein) and recently applied also to systems described by PDEs [17][18][19]. SOS methods are quite computationally efficient since they are based on semidefinite programming (SDP) [20,21]. In this respect, the SOS approach allows finding Lyapunov functionals just like the methods based on LMIs (linear matrix inequalities) enable computing Lyapunov functions for linear systems based on ODE models [22]. As compared with the search of classical Lyapunov functions to investigate local stability [23], we do not rely on linearized models but provide conditions of global stability.
The paper is organized as follows: Section 2 reports the basic definitions that will be used in the following discussions. Modeling of thermal flows in one-dimensional heat exchangers is given in Section 3. The proposed estimation approach is presented in Section 4, where stability is analyzed in a noise-free modeling framework. This analysis is extended to models affected by bounded disturbances in Section 5. A brief account on the use of the SOS method to find the Lyapunov functionals and thus ensure stability is given in Section 6. Section 7 illustrates the numerical results, while conclusions are drawn in Section 8.

Notations and Definitions
The set of real numbers is denoted by ℝ and hence ℝ + will be used for the set of strictly positive real numbers. The set of the integer numbers equal to or greater than zero is denoted by ℕ.
The minimum and maximum eigenvalues of a symmetric matrix P ∈ ℝ n×n are denoted by λ min ðPÞ and λ max ðPÞ, respectively. Moreover, P > 0 (P < 0) means that it is also positive (negative) definite. Given a generic matrix M, jMj ≔ ðλ max ðM ⊤ MÞÞ 1/2 = ðλ max ðMM ⊤ ÞÞ 1/2 , and hence, for a vector x = ðx 1 , x 2 , ⋯, x n Þ ∈ ℝ n , |x| ≔ ðx ⊤ xÞ 1/2 is its Euclidean norm. ℝ½x denotes the ring of real polynomials in x ∈ ℝ n ; ℝ½x d is the set of real polynomials of a degree equal or less than d ∈ ℕ.

Modeling Heat Transfer in Heat Exchangers
Consider the one-dimensional model of counterflow heat exchangers with three distributed state variables, i.e., T c ðx, tÞ, T h ðx, tÞ, and T m ðx, tÞ for x ∈ ½0, L and t ≥ 0, as depicted in Figure 1; such variables represent the temperatures of cold fluid, hot fluid, and pipe, respectively (the subscript indicating that it is made of some metal in such a way to maximize the heat transfer). Based on the balance of enthalpy for the hot fluid, we get where ρ h is the density (in kg/m 3 ), c h is the specific heat (in is the internal surface of the pipe, and p h is its perimeter (in m). Concerning the cold fluid, the balance of enthalpy is given by where ρ c is density (in kg/m 3 ), c c is the specific heat (in J/(kg K)), and U c is the transfer coefficient (in W/(m 2 )). Moreover, S c (in m 2 ) is the external surface of the pipe and p c is its perimeter (in m). Finally, for the metal pipe it follows that where ρ m is the density (in kg/m 3 ) of the metal, c m is the specific heat (in J/(kg K)) of the metal, and S m ≔ S c − S h (in m 2 ) is the surface of the pipe section.
After dividing each of the three previous equations for the first coefficient in each respective l.h.s., we get where or, more concisely, as follows: where Tðx, tÞ ≔ colðT h ðx, tÞ, T c ðx, tÞ, T m ðx, tÞÞ ∈ ℝ 3 and The system (11) is solved with initial conditions and Dirichlet boundary conditions for all t ≥ 0. The boundary conditions on T m are not necessary since equation (9) does not contain spatial derivatives of T m . Likewise, in [7], a reduced model can be easily derived from (11) by neglecting the dynamics of the temperature in the pipe, i.e., assuming that such a temperature is fixed by the values of the temperatures of cold and hot fluids, i.e., imposing that the r.h.s. of (6) is null. Therefore, we obtain or, more compactly, where Tðx, tÞ ≔ colðT h ðx, tÞ, T c ðx, tÞÞ ∈ ℝ 2 and Advances in Mathematical Physics with initial conditions and Dirichlet boundary conditions for all t ≥ 0.
Model (11) is of a hyperbolic type but not strictly hyperbolic since matrix A 1 has a zero eigenvalue. By contrast, (16) is strictly hyperbolic since the eigenvalues of A 2 have a nonnull-real part. In Section 4, we will consider the problem to reconstruct all the temperatures by using a generic hyperbolic model.

Estimation for Heat Exchangers with Stable Estimation Error
Consider the hyperbolic equation where Tðx, tÞ ∈ ℝ n , A ∈ ℝ n×n diagonal, and B ∈ ℝ n×n . Consider also the observer for (20) given by whereTðx, tÞ ∈ ℝ n is the estimate of Tðx, tÞ. Such PDEs need some suitable boundary conditions, which will be given later. The stability of the estimation error can be guaranteed under suitable conditions as follows. and for all x ∈ ½0, L.
Proof. Consider the Lyapunov functional and note that the asymptotic stability of the estimation error is ensured if for all eðx, tÞ ∈ ℝ n since PðxÞ > 0 for all x ∈ ½0, L. Since P does not depend on time, we get and, owing to the diagonal structure of A and P, Thus, it follows and hence (25) holds owing to (22) and (23).
Note that if, instead of (23), we consider with α > 0 for all x ∈ ½0, L, the estimation error is L 2 exponentially stable with a rate of decrease equal to α. Thus, a design goal may consist in maximizing α, which can be regarded as a sort of generalized eigenvalue problem for Lyapunov functionals instead of Lyapunov functions [22]. In case of the reduced model (16) (i.e., with A = A 2 and B = B 2 ), the stability conditions of Theorem 1 can be greatly simplified according to [7] by using and boundary conditionŝ for all t ≥ 0 with ℓ 1 , ℓ 2 , μ 1 , μ 2 ∈ ℝ to be suitably chosen. More specifically, in [7] it is proposed to satisfy (23) by using an LMI-based discretized approach to get μ 1 , μ 2 and then select ℓ 1 , ℓ 2 such that Advances in Mathematical Physics by choosing, for example, The construction of Lyapunov functionals is not an easy task to solve in general and will be addressed later in Section 6.

Estimation under Modeling Uncertainties
Let us focus on estimation in the presence of bounded disturbances. More specifically, instead of (20) we consider where D ∈ ℝ n×q and wðx, tÞ ∈ ℝ q such that jw i ðx, tÞj ≤ 1, i = 1, ⋯, q for all x ∈ ½0, L and t ≥ 0 without loss of generality (different upper bounds can be taken into account by suitably scaling the coefficients of D), where q is an integer number no larger than n. Estimation will be performed by using the same observer adopted in the noise-free case, i.e., (21). The presence of the system noises prevents from ensuring asymptotic stability, but since such disturbances are bounded, an invariant set for the estimation error exists and can be studied by using QB [9].
Toward this end, first of all we need to define QB in the L 2 sense. More specifically, the estimation error is said to be L 2 quadratically bounded with the Lyapunov functional VðtÞ as defined in (24) As a matter of fact, any positive constant can be chosen instead of L, but with such a choice we reduce the notational burden. This does not entail loss of generality since the conditions ensuring QB are homogeneous in the design parameters, which scales with L and thus allows for this simplification.
Owing to (36), the set turns out to be positively invariant and it is attractive (i.e., if the error is out of E, it approaches E asymptotically).
Proof. Likewise in the proof of Theorem 1, we easily compute the time derivative of the Lyapunov functional (24) and get that _ VðtÞ < 0 subject to (22) is equivalent to Since −e ⊤ Pe + 1 > 0 implies VðtÞ > L and −w 2 i + 1 > 0 holds, using the well-known S-procedure (see [22], p. 23) we obtain that QB holds if for some α ∈ ℝ q >0 and β > 0 or equivalently if (38) and (39) hold.

Clearly, (38) implies
for all x ∈ ½0, L, which ensures that the estimation error in the absence of noises is L 2 exponentially stable with a rate of decrease equal to β. Thus, a convenient design goal may consist in maximizing β, as pointed out in [7]. Another, quite popular objective of design is related to the L 2 gain between disturbance and estimation error. Theorem 2 provides conditions ensuring that the set E is attractive. Thus, one can design the estimator in such a way to keep E as small as possible. Toward this end, note that for all e ∈ E and hence, as a design objective, we may maximize the minimum eigenvalue of PðxÞ since for all t ≥ 0, x ∈ ½0, L. The inequality above allows proving the boundedness in the L 2 sense as for all t ≥ 0.
To compare the effectiveness of a given observer in terms of rejection of the noise, it is quite popular to rely on the notion of the L 2 gain by assuming finite-energy noises, i.e., wð·, tÞ ∈ L 2 ðΩÞ. More specifically, we say that the proposed 5 Advances in Mathematical Physics observer admits and L 2 gain between disturbance and estimation error is equal to γ > 0 if for T > 0. It is straightforward to show that the bound of the L 2 gain holds if Note also that this condition is quite similar to (38). Moreover, minimizing γ is equivalent to maximizing 1/γ 2 in accordance with the goal of maximizing β.

Search of Lyapunov Functionals Using the SOS Approach
Stability conditions such as (23), (29), (38), or (47) (each together with PðxÞ > 0) demand the satisfaction of LMIs in PðxÞ for all x ∈ ½0, L, i.e., the solution of semi-infinite programming problems. Such problems arise in a variety of applications and are usually approximately solved in discretized form on a sufficiently fine mesh of points (see, e.g., [24]). The main difficulty to address semi-infinite programming problems concerns both the choice of a sufficiently large number of points and especially the local minima trapping, by which nonlinear programming solvers may be affected in trying to find the solution. A more appealing way to solve such problems consists in resorting to the SOS paradigm, which enables turning the construction of a Lyapunov functional into a convex problem that can be efficiently solved by using SDP without encountering such issues due to local minima.
The idea behind such an approach is the SOS decomposition of a candidate Lyapunov functional as well as of the opposite of its time derivative by using a positivity certification, which does not depend on the characteristics of the chosen polynomial [21,25]. More specifically, the following result holds.

Theorem 3.
A polynomial pðxÞ ∈ ℝ½x 2d in x = ðx 1 , ⋯, x n Þ ∈ ℝ n has sum-of-squares decomposition (or is said to be SOS) if and only if there exists a real symmetric and positive semide- where v d ðxÞ is the vector of all the monomials in the components of x ∈ ℝ n of degree equal to or less than d ∈ ℕ, i.e., of dimension Proof. See Proposition 2.1 in [26] (p. 17).
Consider, for example, the estimation of the state variables of (7), (8), and (9) by using the results of Theorem 1. Based on Theorem 3, a procedure can be adopted to find suitable SOS polynomials as diagonal elements of PðxÞ = diag ðP 11 ðxÞ, P 22 ðxÞ, P 33 ðxÞÞ, all to be taken in ℝ½x 2d with x ∈ ℝ for increasing values of d. Toward this end, we say that pðxÞ ∈ ℝ½x 2d is ε-SOS polynomial if pðxÞ − ε∑ n i=1 x 2 i ∈ ℝ½x 2d with d ∈ ℕ, d ≥ 1 is SOS for some "small" tolerance ε > 0 [27]. In addition, a polynomial, square matrix MðxÞ with M ij ðxÞ ∈ ℝ½x 2d for i, j = 1, ⋯, m is said to be an ε-SOS matrix if y Τ MðxÞy is ε-SOS in ℝ½x, y with x ∈ ℝ and y ∈ ℝ m .
Given d ∈ ℕ with d ≥ 1 and ε > 0, the problem to solve when dealing, for example, with the stability conditions of Theorem 1 is the following: Problem 4. Find P 11 ðxÞ, P 22 ðxÞ, P 33 ðxÞ ∈ ℝ½x 2d such that for all x ∈ ½0, L.
If Problem 4 admits a solution, we can follow the same reasoning that leads to (34) to ensure the satisfaction of (22) by choosing Similar problems may be formulated by replacing the SOS description of (52) with those corresponding to (29), (38), and (47).

Numerical Results
For the sake of brevity, in the following discussion, we will describe only the numerical solution of the boundary value problem that results from the discretization of the model (11) since this includes that of the reduced model (16) as a special case. We have discretized the domain ½0, L with M equally spaced points x j for j = 1, ⋯, M, and both (16) and The system of ODEs (52) is solved by using a fourthorder Runge-Kutta method (implemented in the ode45 function of MATLAB) with variable time step. In the first instants of the simulation, the fast dynamics of the convective part prevails, while later the reaction term predominates. To accurately simulate the model, it is therefore necessary to select a very small time step in the first time instants, while it may be chosen larger afterwards, but of course, always respecting the stability condition of the numerical scheme. Implicit methods (e.g., the scheme implemented in the routine ode15s of MATLAB) are tested, since they generally allow using a larger time step. This resulted in a longer running time, due to the higher complexity to complete each iteration. Then, explicit methods seem to be the most appropriate choice for this problem.
In the numerical case study, we have fixed L = 5 m, a 1 = 1 m/s, a 2 = 1 m/s, b 1 = 0:01 1/s, b 2 = 0:01 1/s, b 3 = 0:01 1/s, and b 4 = 0:01 1/s. Concerning the proposed estimation approach with (11) and (16), we have solved Problem 4 by using the SOS toolbox [27] and the numerical solver Yalmip [31] with ε = 10 −6 . In the first case (i.e., A = A 1 and B = B 1 ), we have got a solution to such a problem with d = 1 given by and taking ℓ 1 , ℓ 2 again according to (52). We have chosen M = 500, T h ð0, 0Þ = 55, T c ðL, 0Þ = −15, and T m = 50. Figure 2 shows the behavior of the temperatures T h , T c , and T m and their estimates given byT h ,T c , andT m at different time instants, all based on (11). Figure 3 shows the behav-ior of the temperatures T h , T c and the corresponding estimatesT h ,T c , based on the reduced model (16).

Conclusions
Two models have been presented to account for the temperature dynamics in heat exchangers. Such models are based on hyperbolic PDEs with the most complete having the temperatures of cold fluid, hot fluid, and pipe as state variables, whereas the reduced one relies only on the temperatures of cold and hot fluids. Estimators resulting from the solution of boundary value problems based on such models have been proposed and provided with a rigorous stability analysis with the Lyapunov approach. To this end, the existence of nonlinear polynomial Lyapunov functionals has been addressed by using the SOS approach.  Advances in Mathematical Physics Future work will concern the investigation on the stability of systems described by nonlinear hyperbolic equations (see, e.g., [32]). Another direction of research is the study of stability of cascaded hyperbolic and other nonlinear PDEs such as the normal flow equation to deal with multiphase problems, which turn out to be increasingly difficult but with a wide range of potential applications [33].

Data Availability
The simulation code is available upon request to Angelo Alessandri, email: alessandri@dime.unige.it.

Conflicts of Interest
The authors declare that they have no conflicts of interest.