Fractional-Order Identification and Control of Heating Processes with Non-Continuous Materials

Abstract

The paper presents a fractional order model of a heating process and a comparison of fractional and standard PI controllers in its closed loop system. Preliminarily, an enhanced fractional order model for the heating process on non-continuous materials has been identified through a fitting algorithm on experimental data. Experimentation has been carried out on a finite length beam filled with three non-continuous materials (air, styrofoam, metal buckshots) in order to identify a model in the frequency domain and to obtain a relationship between the fractional order of the heating process and the different materials’ properties. A comparison between the experimental model and the theoretical one has been performed, proving a significant enhancement of the fitting performances. Moreover the obtained modelling results confirm the fractional nature of the heating processes when diffusion occurs in non-continuous composite materials, and they show how the model’s fractional order can be used as a characteristic parameter for non-continuous materials with different composition and structure. Finally, three different kinds of controllers have been applied and compared in order to keep constant the beam temperature constant at a fixed length.

1. Introduction

The association between non-continuous materials and Fractional Order Systems [1,2,3,4] is an actual research issue involving investigation at the micro- and macroscopic scale [5,6]. This research follows the recent decades’ stream, according to which an increasing interest has arisen for non-continuous, heterogeneous, and composite materials. They are in fact becoming more common in a wide range of functional devices—particularly those involving energy transport and conversion [7]. Their physical properties can be designed as a function of their composition in terms of materials and structure. Modeling and controlling such properties therefore becomes crucial for the optimal materials design.

In particular, the heat transfer model optimization in the frequency domain represents a topic of great interest, considering the common situations of periodic heat flux in engineering problems such as cyclic heating of the cylinder surface of internal combustion engines, diurnal heating and nocturnal cooling of building structures, and periodic (pulse) laser heating of solid surfaces in materials processing [8].

The fractional order calculus was especially useful for modelling thermal processes, as shown in [9,10,11,12,13]. The physics phenomena of such a distributed system model have been approached in literature through different methods: finite difference method by using Grunwald–Letnikov (see [1]), definition for the fractional time derivative [14], and other numerical methods [15,16], such as Podlubny’s matrix approach for discretization of integrals and derivatives of non-integer order [17] and implicit difference approximation to solve the time fractional diffusion equation have been used in [18]. Furthermore, non-Markovian stochastic process—related to a phenomenon of slow anomalous diffusion—have been studied with discrete models of random walk [19].

The aim of this paper is twofold. The first goal is to investigate the fractional order nature of such phenomena in non-continuous, heterogeneous, porous, or composite materials. In particular, the order variation of the heat transfer model has been investigated on non-homogeneous beam heating process by exploiting optimization methods for model identification and validation in the frequency and time domains. The second goal consists of the comparison of classical Proportional Integrative (PI) controllers with the Fractional Order one (FO-PI), keeping constant the beam temperature at a fixed point.

2. Fractional Theoretical Model

Non-continuous material beam heating is associated with the effect of anomalous diffusion, due to the discontinuous interfaces within the material and to part of the heat flux dispersed into the neighbouring environment [8]. Theoretical and experimental models are here presented and compared. The diffusion heat equation applied to a semi-infinite bar is described in Equation (1) in the time domain. Considering the initial condition in Equation (2), the boundary conditions in Equations (3) and (4), and applying the Laplace transform, it is possible to obtain the transfer function $G(x,s)$ in Equation (5), as in [8]. H is the heat flux, T is the temperature, and λ is a normalized distance from the heat source.

$$\frac{\partial T(x,t)}{\partial t}=a^2\frac{\partial T^2(x,t)}{\partial x^2}$$

(1)

$$T(x,t=0)=0$$

(2)

temperature along the tube at $t=0$

$$H(x=0,t)=-K\frac{\partial T}{\partial x}(x=0)$$

(3)

heat flux at $x=0$

$$T(x\to \infty ,t=0)=0$$

(4)

temperature at $x\to \infty$

$$G(x,s)=\frac{T(x,s)}{H(x,s)}=\frac{T_1}{{\left(T_{3}s\right)}^{0.5}+1}{e}^{-\lambda {\left(T_2\right)}^{0.5}}$$

(5)

3. Fractional Experimental Model Identification

3.1. Experimental Setup

The thermal system approximating a semi-infinite beam consists of a steel beam (length $L=40$ cm, outer section diameter $S=4$ cm, inner section diameter $s=3$ cm, thickness $t=5$ mm) filled with different non-continuous materials. The entire experimental setup, schematic and devices, are shown in Figure 1a,b. One end of the metal beam is connected to a Peltier cell (representing the heat source), which is integrated in the heat pump assembly $THP51B$. Such assembly contains a heat sink and the fan: they are used to control the temperature of the cell cold side and to maintain it approximately constant at the ambient temperature. The heat flux out of the Peltier module is controlled by the input current signal coming from a signal generator, and is amplified by a power amplifier $OPA549$. Three temperature sensors $LM35DH$ are positioned at a fixed distance from the heat source ($sensor1$ at $d_{{\lambda }_1}$, $sensor2$ at $d_{{\lambda }_2}$, and $sensor3$ at $d_{{\lambda }_3}$), as detailed in the picture in Figure 1c. Three different experimental setups have been created by varying the beam filling material. The idea was to create composite materials through the composition of granular-shaped materials with air. Such a choice would provide an easy solution for inserting and removing the heterogeneous material from the tube without changing the base of the experimental setup. The maintenance of the basic apparatus, consisting of the metal tube and the Peltier module, ensures the repeatability of the experiments when using different materials. In particular: air, Styrofoam (spheres of 2 mm diameter) and lead buckshot (spheres of 1 mm diameter) have been used in order to investigate how the fractional order of the thermal phenomena varies for the different materials with different thermal conductivity. The lead buckshots filling the tube represent a composite material made of a metal, with good thermal conductivity (35 W/(m·K)), and air, having thermal conductivity of $0.025$ W/(m·K). The Styrofoam spheres, which represent a thermal insulator with thermal conductivity of $0.04$ W/(m·K) with air, represent the second composite materials; while the air filling is used for the control experiment.

3.2. Experimental Model

In this paper, an enhanced experimental model for thermal diffusion has been introduced. The theoretical ideal model in Equation (5) has been thus modified, leading to Equation (6) in order to better fit the experimental data. In particular, the fractional order of the system (α) has been left as a parameter to be identified.

The ${\alpha }_1$ parameter has been included in order to discriminate within the mathematical model the component of the energy transport related to the metal tube. This is important in the parameter identification procedure, in order to decouple the fractional component related to the composite sample from the one related to the experimental setup artifact.

Moreover, the spatial variable λ contribution has been considered in the pole position. Such parameters to be identified are at the basis for the investigation of a possible relationship between the thermal diffusion fractional order and the non-continuous material of the beam.

$$T(s,\lambda )=\frac{T_1}{{\left(T_3{\lambda }^{0.5}s\right)}^{\alpha }+1}{e}^{-\lambda {\left(T_{2}s\right)}^{\alpha {\alpha }_1}}H(s,0)$$

(6)

3.3. Experimental Campaign

Sinusoidal signals at different frequencies have been used for the model identification in the frequency domain. In particular, input current sinusoids with frequencies $0.1\mathrm{mHz}$, $0.3\mathrm{mHz}$, $0.5\mathrm{mHz}$, $0.8\mathrm{mHz}$, $1\mathrm{mHz}$, $2\mathrm{mHz}$, $3\mathrm{mHz}$, $4\mathrm{mHz}$, $5\mathrm{mHz}$, $6\mathrm{mHz}$, $7\mathrm{mHz}$, $8\mathrm{mHz}$, $9\mathrm{mHz}$, $10\mathrm{mHz}$, and $20\mathrm{mHz}$ were used, causing a heat flux peak-to-peak amplitude of about 45 W. For each frequency, the sinusoidal signal has been applied as a current input, and the temperature measurement was collected on three sinusoids when the steady state condition was reached for the output temperature sinusoid off-set. At the end of each sinusoid measurement, the system temperature was allowed to return to the environmental temperature.

Three sets of experimental data acquisition have been conducted by varying the beam filling material: air, Styrofoam, and lead buckshot.

3.4. Identification and Validation Method

Given the spatial dependence of the model, the identification process has been performed in two steps. First, a multi-objective Nelder–Mead simplex optimization algorithm has been used in order to identify the model parameters ($T1$, $T2$, $T3$, α, and ${\alpha }_1$) by fitting the model response (module and phase) to the experimental data for $sensor1$ at the normalized distance ${\lambda }_1$. The objective function $J_{err}$ combines the error functions on the module and phase. Then, in order to validate the identified model at different distances from the heat source, the same algorithm was used in order to validate the ${\lambda }_2$ and ${\lambda }_3$ parameters by comparing—in the frequency domain—the model output T to the experimental temperature measured, respectively, at distances $d_{{\lambda }_2}$ and $d_{{\lambda }_3}$. The quality of the identified model on the spatial domain has been measured by comparing estimated distances $d_{{\lambda }_2}$ and $d_{{\lambda }_3}$ (corresponding to the identified ${\lambda }_2$ and ${\lambda }_3$) with the real ones.

In order to determine the model parameters, the following $J_{err}$ index error—consisting of the sum of two terms—has been considered. The first term, $J_{module}$, is related to the normalized error between the module of the experimental ratio $G=T/H$ and the simulated one of Equation (6), while the $J_{phase}$ represents the corresponding error for the phase.

$$J_{err}=J_{module}+J_{phase}$$

(7)

with

$$J_{module}=\frac{\sqrt{\sum \left(\right|\left|G_{exp}\right|-\left|G_{sim}\right|{|}^2)}}{\sqrt{\sum \left(\right|G_{exp}{|}^2)}}$$

(8)

$$J_{phase}=\frac{\sqrt{\sum \left(\right|\angle G_{exp})-\angle G_{sim}{|}^2)}}{\sqrt{\sum \left(\right|\angle G_{exp}{|}^2)}}$$

(9)

The minimization of the cost function has been performed sequentially using the Nelder–Mead simplex method [20] implemented via the Matlab function fminsearch [21], see also Figure 2.

4. Modelling Results

4.1. Model Identification at sensor1

Three couples of identified models, theoretical and experimental, have been obtained by applying the identification and the validation algorithm on the experimental data acquired at the $sensor1$ location on the metal tube filled with the three different materials. The identified parameters are presented in

Table 1. Theoretical model identified parameters at $sensor1$.

Material	T1	T2	T3	α	${\mathit{\alpha }}_1$
Air	13.8610	289.7543	1222344.6	$N/D$	$N/D$
Styrofoam	13.7534	293.0218	1248787.9	$N/D$	$N/D$
Metal Buckshots	8.1853	282.1817	481148.7	$N/D$	$N/D$

and

Table 2. Experimental model identified parameters at $sensor1$.

Material	T1	T2	T3	α	${\mathit{\alpha }}_1$
Air	2.0589	117.6201	9260.2	0.7408	0.7077
Styrofoam	2.1589	109.8189	11252.4	0.7382	0.7649
Metal Buckshots	2.5013	456.3768	13193.5	0.5771	0.6869

, respectively, for the theoretical model ($T1$, $T2$, $T3$) and for the experimental model ($T1$, $T2$, $T3$, α, and ${\alpha }_1$). Comparison between the theoretical model, the experimental model, and the experimental data at the $sensor1$ location is performed in the frequency domain through the Bode diagrams—both module and phase—as presented in Figure 3, Figure 4 and Figure 5.

The experimental model’s enhanced performances in fitting experimental data at the $sensor1$ location are visible by inspection of the Bode diagrams (see Figure 3, Figure 4 and Figure 5). More evidence of such enhancement is given by the mean square error ($J_{err}$) calculated between the identified models module and the experimental data at $sensor1$, as shown in

Table 3. Experimental and theoretical model validation in the spatial domain.

Material	Model	${\mathit{J}}_{\mathit{err}}$	${\mathit{\lambda }}_1$	${\mathit{d}}_{{\mathit{\lambda }}_1}$	${\mathit{\lambda }}_2$	${\mathit{d}}_{{\mathit{\lambda }}_2}$	${\mathit{\lambda }}_3$	${\mathit{d}}_{{\mathit{\lambda }}_3}$
Air	Theo.	10.5569	1	3.2	1.1542	3.69	1.3496	4.32
	Exp.	3.5595	1	3.2	1.2031	3.85	1.50	4.8
Styrofoam	Theo.	9.073	1	3.2	1.1622	3.72	1.3607	4.35
	Exp.	3.087	1	3.2	1.2104	3.87	1.51	4.83
Metal Buckshots	Theo.	6.9831	1	3.2	1.2001	3.84	1.4096	4.51
	Exp.	2.3032	1	3.2	1.2929	4.14	1.6062	5.14

. Moreover, as it is shown in Table 2, the fractional order of the system for the metal buckshots is close to the theoretical value of $\alpha =0.5$. The module of the Bode diagram in Figure 5 confirms such evidence, since the theoretical model fits the experimental data with a smaller $J_{err}$ with respect to the other pipe filling material.

4.2. Time Domain Model Validation at sensor1

The experimental model has been validated in the time domain: the step response analysis has been performed, the steady state value of the heat flux has been measured providing the value of 40 W over a $\mathrm{\Delta }T$ = 100 K of the input step amplitude (see Figure 6). This leads to a static gain of 2.5 K/W for all the materials, which is comparable to the $T1$ parameter found through experimental model identification, ranging from 2 and 2.5 K/W. Such analysis suggests that the theoretical model in this case fails in fitting the experimental data having a static gainranging from 8 to 13 K/W. It is also worth noticing that the position of the pole for the theoretical model is at a very low frequency, and this affects such a mismatch between the static gains of the two considered models.

4.3. Model Validation at sensor2 and sensor3

The model validation results in the spatial domain has been presented for the three materials by comparing the estimated distances $d_{{\lambda }_2},d_{{\lambda }_3}$ with the real sensor distances $d_{{\lambda }_2}=4.2$ cm, $d_{{\lambda }_3}=5.2$ cm. Also in the validation step, the experimental model performs better than the theoretical model.

It is clear how the estimated distances $d_{{\lambda }_2}$ and $d_{{\lambda }_3}$ in Table 3 for the experimental model are closer to the $sensor2$ and $sensor3$ actual distances with respect to the theoretical model estimated distances.

4.4. Fractional Order and Non-Continuous Materials

The identified values of the fractional order of the experimental model (see

Table 4. Relationship between thermal conductivity of the non-continuous beam material and the experimental fractional order model α.

Material	Thermal Conductivity	α
Air	0.025	0.7408
Styrofoam	0.04	0.7382
Metal Buckshots	35	0.5771

) for the different non-continuous materials suggests that there is a relationship between the material composing the beam and the fractional behavior. This might be due to the multiple interfaces within the material leading to anomalous diffusion. The morphological composition and the relative thermal properties of the non-continuous composite materials can therefore be considered the basis for designing complex materials with desired thermal properties.

4.5. Controllers Implementation and Comparison

In this section, the comparison of three different approaches for the control of the temperature of the beam at the distance $d_{{\lambda }_1}=3.2$ cm is reported. The control systems have been implemented by using the HIL system DSpace ACE-kit DS1103. Its main features are the PowerPC $604e@400MHz$, 2 MByte local SRAM, 128 MByte global DRAM, 20 ADC channels (16/12 bit), and 8 DAC channels (14 bit).

All the considered controllers have been designed in Matlab/Simulink and then downloaded and tested in real time using the $ControlDesk$ tool as part of the HIL system.

The compared controllers are a conventional PI, an anti-windup PI, and a fractional order (FO) PI controller. The PI has been designed according to the Ziegler–Nichols tuning rules, and the following parameters have been used during the simulation: $K_p=81.08$, $K_i=0.19$; for the anti-windup controller, see [22]. A feedback gain equal to $0.02$ has been added, while the FO-PI has been designed according to the tuning table given in [23], with $K_p=57.62$, $K_i=0.27$, and the fractional order $\alpha =0.7$. In particular, the fractional order controller has been implemented using the Grunwald–Letnikov discrete approximation, as in [1].

As it can be seen in Figure 7, the best performances have been obtained with the fractional controller. This result can be supported by the evidence that a fractional order system can be better controlled with a fractional controller instead of a classical one. The optimization of all the controllers is foreseen as a next step, as well as the implementation and comparison of the control laws on a low-cost microcontroller.

5. Conclusions

In this paper, as a first step, the fractional order of the heat transfer phenomenon on a finite beam made of non-continuous composite materials has been investigated.

An optimized heat transfer model has been obtained by identifying the fractional order by means of a Nelder–Mead optimization algorithm. The model has been identified for three different non-continuous materials: air, Styrofoam, and metal buckshot.

The proposed model is a first step towards a complete enhanced experimental model. The convection analysis is to be included in future models, considering that convective flow can occur in the air and Styrofoam samples. Such an additive effect could affect the fractional order of the model, especially in the case of air and Styrofoam filling.

With this in view, the experimental set-up could be improved taking into account the orientation of the tube to evaluate the effect of the convection and the effect of the tube insulation on the energy transport between the composite sample and the tube itself.

The model optimization paves the way for future work in the design of heating process control for engineering applications. Such results allow us to investigate a simplified thermal model of complex composite materials with a view towards defining the relationship between fabrication parameters, the thermal capacitance and possible mixtures, and the fractional exponent value. This concept opens the way for investigation of the technological control of material thermal properties by composite material structure, concentration, and mixtures, making such design more simple and flexible.

The second step shows the comparison among fractional order and standard PI controllers, outlining a better performance of the fractional one. With respect to the control strategies, future steps consist of the optimization of all the controllers, as well as the implementation and comparison of the control laws on microcontrollers.