Fractional Order Model of the Two Dimensional Heat Transfer Process

: In this paper, a new, state space, fractional order model of a heat transfer in two dimensional plate is addressed. The proposed model derives directly from a two dimensional heat transfer equation. It employes the Caputo operator to express the fractional order differences along time. The spectrum decomposition and stability of the model are analysed. The formulae of impluse and step responses of the model are proved. Theoretical results are veriﬁed using experimental data from thermal camera. Comparison model vs experiment shows that the proposed fractional model is more accurate in the sense of MSE cost function than integer order model.


Introduction
It is known that the non integer order calculus can be applied in modeling of processes and phenomena hard to analyse with the use of other tools. Non integer order (NIO) or fractional order (FO) models of different physical phenomena have been presented by many authors. The amount of FO models of various processes is collected in the book [1]. The book [2] presents fractional order models of chaotic systems and Ionic Polymer Metal Composites (IPMC). Fractional models of ultracapacitors are "classics" of FO modeling. They are given for example by [3]. Distributed parameter systems can be also described using FO approach. As an example diffusion processes discussed in [4][5][6], can be given. A collection of recent results employing a new Atangana-Baleanu operator can be found in [7]. In this book, i.e., the FO blood alcohol model, the Christov diffusion equation and fractional advection-dispersion equation for groundwater transport process are presented.
Different thermal processes can also be described using FO approach. For example a temperature-heat flux relationship for heat flow in semi-infinite conductor is presented in [1], the beam heating problem is given in [3]. Theoretical analysis of fractional thermal equation is given, e.g., in [8][9][10][11].
Models of temperature fields obtained using thermal cameras are presented by [12,13], temperature models in three dimensional solid body are given in [14]. The use of fractional order approach to the modeling and control of heat systems is also presented in [15]. The analytical solution of the two dimensional heat equation is proposed in the paper [16]. Numerical methods of solution partial differential equation can be found, e.g., in [17]. Fractional Fourier integral operators are analyzed, u.a., in [18].
This paper is devoted to present a new, discrete, FO model of heat transfer process in two dimensional plate. Such a process is described by a partial differential equation (PDE) of parabolic type. The partial derivative along time is expressed by Caputo operator, both partial, spatial derivatives are integer order. Such a model of thermal process has not been previously published.
The paper is organized as follows: Preliminaries recall some elementary ideas and definitions from fractional calculus. Next, the state space model of two dimensional heat transfer process is proposed and discussed. Finally the experimental validation of theoretical results is given.

Preliminaries
At the beginning the non integer order, integro-differential operator is presented (see, e.g., [1,[19][20][21]). Definition 1 (The elementary non integer order operator). The non integer order integrodifferential operator is defined as follows: (1) where a and t denote time limits for operator calculation, α ∈ R denotes the non integer order of the operation.
Next, an idea of complete Gamma Euler function is recalled (see for example [20]): Definition 2 (The complete Gamma function).
An idea of Mittag-Leffler function needs to be given next. It is a non integer order generalization of exponential function e λt and it plays crucial role in the solution of fractional order (FO) state equation. The two parameter Mittag-Leffler function is defined as follows: Definition 3 (The two parameter Mittag-Leffler function).
For β = 1 we obtain the one parameter Mittag-Leffler function: Definition 4 (The one parameter Mittag-Leffler function).
The fractional order, integro-differential operator (1) is described by different definitions, given by Grünwald and Letnikov (GL definition), Riemann and Liouville (RL definition) and Caputo (C definition). Relations between Caputo and Riemann-Liouville, between Riemann-Liouville and Grünwald-Letnikov operators are given, e.g., in [1,22]. Discrete versions of these operators are analysed with details in [23]. The C definition has a simple interpretation of an initial condition (it is analogical as in integer order case) and intuitive Laplace transform. Additionally, its value from a constant equals to zero, in contrast to, e.g., RL definition. That is why in the further consideration the C definition along time will be used. It is recalled beneath.

Definition 5 (The Caputo definition of the FO operator).
In (5) M is a limiter of the non integer order: M − 1 ≤ α < M. If M = 1 then, consequently, 0 ≤ α < 1 is considered and the definition (5) takes the form: Finally a fractional linear state equation using Caputo definition should be recalled. It is as follows: where α ∈ (0, 1) is the fractional order of the state equation, x(t) ∈ R N , u(t) ∈ R L , y(t) ∈ R P are the state, control and output vectors respectively, A, B, C are the state, control, and output matrices, respectively. Figure 1 shows the simplified scheme of the considered heat system. It has the form of thin rectangular metallic surface (PCB plate) heated by flat heater, denoted by H. Coordinates of heater are denoted by x h1 , x h2 , y h1 and y h2 , respectively. The temperature is read using thermal camera; the area of measurement is configurable. The area of measurement is denoted by S and its coordinates are equal x s1 , x s2 , y s1 and y s2 , respectively. More details about the construction of this laboratory system are given in the section "Experimental Results".

Measuring Area
Heater Y X x s1 ,y s1 x s2 ,y s2 x h1 ,y h1 x h2 ,y h2 (0,0) The fundamental mathematical model describing the heat transfer in the surface is the Partial Differential Equation (PDE) of the parabolic type. All the side surfaces of plate are much smaller than its frontal surface. This allows to assume the homogeneous Neumann boundary conditions at all edges of the plate as well as the heat exchange on the surface needs to be considered too. It is expressed by coefficient R a . The control and observation are distributed due to the size of heater and size of temperature field read by camera. The heat conduction along both directions x and y is the same and described by coefficient a w . Two dimensional, IO heat transfer equation has been considered in many papers (e.g., [24][25][26]). Its fractional version proposed by authors is as beneath: In (8), α is non integer order of the system, a w > 0, R a ≥ 0 are coefficients of heat conduction and heat exchange, k 0 is a steady-state gain of the model, b(x, y) and c(x, y) are heater and sensor functions. Denote the area of heater by H and the area of measurement by S. Then the heater and sensor functions take the following form: Let Ω ⊂ R N it be an appropriate restricted area. The Laplace operator ∆ = ∂ 2 (..) ∂x 2 + ∂ 2 (..) ∂y 2 in L 2 (Ω) with Dirichlet or Neumann boundary conditions is a discrete operator. The discrete operator has only a point spectrum (see, e.g., [27], pp. 204, 460). Without going into detail, it is generally known from the spectral theorem for compact and self-adjoint operators that all eigenvalues λ m,n of the Laplace operator ∆ in L 2 (Ω) (with Dirichlet or Neumann boundary conditions) are non-negative, with finite multiplicities and λ n → ∞ for n → ∞. Additionally, there is in L 2 (Ω) an orthonormal basis (complete system) composed of eigenfunctions of the appropriate operator ∆ . In special cases of the area Ω ⊂ R N (e.g., for a rectangle on a plane) analytical formulae for eigenvalues and eigenfunctions of the appropriate Laplace operator ∆ can be given (see, e.g., [28], pp. 21, 26 for the Dirichlet problem or [29], pp. 133, 138, 301, 305).
However in the considered case the construction of the experimental system requires to assume the homogenous Neumann boundary conditions. Then eigenfunctions and eigenvalues are as follows: Consequently, the 2D heat Equation (8) can be expressed as an infinite dimensional state equation: . (13) where: The state vector Q(t) is defined as beneath: The state operator A takes the following form: The control operator B is as beneath: where: Taking into account (11) we obtain: where: The output operator C is defined analogically: where: Each element c m,n is expressed analogically, as (19): In (23) h xm,yn are expressed by (20).

The Decomposition of the Model
The form of state operator A expressed by (16) allows to decompose the system (13)-(23) to infinite number of independent scalar subsystems, associated to particular eigenvalues. This can be expressed as follows: q m,n + Bu D α q m,n = a w ∂ 2 q m,n ∂x 2 + ∂ 2 q m,n ∂y 2 − R a q m,n + b m,n u.
The form of Equation (24) implies the decomposition of the system (7) into systems related to single eigenvalues λ m,n , m, n = 0, 1, 2, . . . This decomoposition allows to easily compute the step and impulse responses of the system as a sum of responses of particular modes. This is presented below. Remark 1. Consider the system described by (13)- (22). Its step response takes the following form: where m, n-th mode of response is as follows: In (26) E α (..) is the one parameter Mittag-Leffler function (4), λ m,n , b m,n and c m,n are expressed by (12), (18) and (22) respectively.
Proof. The m, n-th mode of the model is expressed as follows: C 0 D α t q m,n (t) = λ m,n q m,n (t) + b m,n 1(t).
The steady-state response to the Heaviside function 1(t) is as follows: Remark 2. Consider the system described by (13)- (22). Its steady-state response to the Heaviside function 1(t) takes the following form: Proof. The Laplace transform from m, n-th mode of step response is given by (28). Its steady-state value is as beneath: The sum of elements y ss m,n gives directly (30).
The impulse response of the proposed model is described by the following remark: Remark 3. Consider the system described by (13)- (22). Its impulse response is as beneath: where: g m,n (t) = t α−1 E α,α (λ m,n t α )b m,n c m,n .
The above remark can be proved with the use of Laplace transform, analogically as (25).
The crucial difference to one-dimensional heat system presented in the paper [31] is that the eigenvalues can be multiple.
The existence of multiple eigenvalues is determined by the dimensions X and Y of the plate and it can be described by the following remark: Remark 4 (The existence of multiple eigenvalues in the model). Consider the model of two dimensional heat transfer described by state equation with state operator (16). Assume that the eigenvalues of this state operator are expressed by (12) and the size of plate is equal X × Y. Then two different eigenvalues λ m 1 ,n 1 and λ m 2 ,n 2 (where m 1 = m 2 , n 1 = n 2 ) are equal iff: Proof. Taking into consideration (12) we obtain: Comparing: λ m 1 ,n 1 = λ m 2 ,n 2 and making some elementary transformations we obtain directly (34).
Analogically, as in one dimensional case discussed, u.a., in [31] the model (13)-(33) is infinite dimensional. Its practical implementation requires to use finite dimensional approximation obtained via truncation of further modes of step or impulse responses. Then the infinite dimensional operators A, B and C can be interpreted as matrices. Consequently, step and impulse responses are the following finite sums: In (35) and (36) M and N denote the size of finite dimensional approximation. It can be estimated numerically or analytically.

The Stability
The spectrum of the proposed model is defined as the set of all eigenvalues of the system: where λ m,n are defined by (12). The spectrum contains negative, single or multiple, pure real eigenvalues. This implies that the proposed FO model is asymptotically stable for each value of plant parameters α, a w and R a . Another interesting problem is the Mittag-Leffler stability of the model. It is determined by location of eigenvalues in the left semiplane. In the one dimensional case all eigenvalues in the spectrum are ordered in the decreasing order. The most poorly damped eigenvalue is always most close to imaginary axis, next eigenvalues follow it in the decreasing order.
In the two dimensional case, we deal with the situation is rather more complicated due to the value of λ m,n is determined by two indices m and n as well as dimensions of the plate. The most poorly damped eigenvalue is still the eigenvalue λ 0,0 , but next eigenvalues are not ordered in the decreasing order. Determing the damping rate requires to order all others eigenvalues.

The Convergence
The convergence of the proposed model will be tested using approach close to presented in papers [32,33]. This can be performed by estimating the orders M and N assuring a predefined value of Rate Of Convergence (ROC). In the considered case the ROC is defined as the increment of steady-state response y ss m,n as a function of orders M and N. This increment is equal to the absolute value of M, N-th mode of the steady-state response (31): The ROC as a function of orders M and N is expressed by the following remark: Remark 5. Consider the system (13)- (22). The Rate of Convergence (ROC) as a function of both orders M and N is expressed as beneath: where λ M,N is expressed by (12), the coefficients H R , P h and P s are as follows: In (40) h xm and h yn are described by (20).
The form of ROC(M, N), expressed by (39), is rather complicated. To simplify it assume the maximum values of P h and P s equal 4. This allows to give the upper estimation of (39) in the following form: The values of orders M and N assuring the keeping a predefined value ∆ of ROC can be computed via numerical solution of the following inequality:

Experimental Validation of Results
Experiments were conducted using laboratory system shown in the Figure 2. The dimensions of the PCB plate in pixels are following: X = 380, Y = 290. The PCB is heated by the heater 170 × 20 pixels with maximum power 10 W located in points: x h1 = 100, y h1 = 40 The temperature of the plate is measured using thermal camera OPTRIS PI 450, connected to computer via USB. The measured range of temperature is 0-250 • C, the frequency of sampling is 80 Hz. Data from camera are collected using dedicated software Optris PI Connect. The signal powering the heater is given from computer using NI LabView, NI MyRIO and amplifier. The maximum current from amplifier equal 400 [mA] at a voltage of 12 V gives a power of 4.8 W. The tested PCB plate not isolated from the environment. This implies that measurements strongly depend on ambient temperature. The presented experiment was conducted during the heat of summer. Estimation the quality of the models only by comparing of diagrams is not enough accurate. The use of cost function is necessary, because only it allows to precisely estimate the accuracy of a model. Additionally such a cost function is employed to parameters identification. In this paper the Mean Square Error (MSE) cost function is used. It describes the mean difference between step response of plant and model at the same time-spatial mesh: In (43), K is the number of all collected samples, y(k) is the step response of the model, computed using (35), y e (k) is the experimental response measured in the same place and at the same time instants k with the use of thermal camera. The sample time during a step response measurement was equal 1 [s]. In each case the mean temperature of the whole area is measured. During experiments, the step response of the system was investigated. The "zero" level denotes the heater switched off, the "one" level is the full power of the heater. The temperature fields for both states are shown in the Figure 3. This figure shows also the points of measurement of the step response, marked as "Area 1-4". Areas 1-3 are located in different points of plate, area 4 covers the heater and it describes its mean temperature. Coordinates of all measuring areas are described by the Table 1. During calculations these coordinates x .. and y .. were given relative to X and Y. For example, x s1 = 75 during calculating elements of C matrix with respect to (23) was equal: x s1 = 75 380 = 0.1974.   The parameters of the fractional model were identified via minimization the cost function (43) as a function of parameters α, a w and R a . Optimization was conducted with the use of MATLAB function fminsearch, the step response was calculated using Formula (35). Results of identification are given in the Table 2       Next, the accuracy of the proposed FO model was compared to accuracy of analogical Integer Order (IO) model, obtained for α = 1.0. The step response of IO model was computed using MATLAB function step. Results are given in the Table 3 and illustrated by Figures 9-12.

Discussion of Results and Final Conclusions
At the beginning limitations of the proposed model need to be discussed. Firstly, the proposed model uses fractional operator only to express derivative along the time and the derivative along the space coordinates is being still integer order. Results for one dimensional heat transfer given in [34] point that the most accurate model uses fractional derivatives along time and space coordinates. On the other hand, the analysis of the proposed model, using only FO operator along time is simplier.
Next, the surface which the temperature is measured, is non homogenous in the sense of values of parameters α, a w and R a . This was observed during identification, when we obtained different values of parameters for different points of measurement. This is visible well in area 3. In addition, this point is the most distant from heater and the impact of disturbances is more significant than in other tested places. A more precise model should employ parameters dependent on spatial coordinates x and y. Such a a model is planned to construct during further investigations.
An another limitation are dimensions M and N of the model. They can not be too high due to the summarized size of the model equal MN. Increasing of these dimensions causes increasing numerical complexity of the model. Fortunately, the proposed model assures relatively good accuracy for M, N < 10. The numerical complexity is also worth researching during future work.
Next, the shape of analyzed surface is relatively simple (rectangular). In reality, a temperature of more complicated elements needs to be modeled. This is not a trivial issue due to it requires to consider much more complicated form of border conditions. Furthermore, some properties, e.g., positivity of this model need o be more deeply analyzed.For one dimensional case the positivity was analyzed with details in the paper [35] and results presented here need to be generalized to the two-dimensional plant considered here. This is planned to do during future work.
The main final conclusion from this paper is that the proposed FO model is more accurate in the sense of the considered MSE cost function than the analogical IO model. This is visible by comparing Table 2 to Table 3.