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

Riccardo Caponetto 1,*, Francesca Sapuppo 1, Vincenzo Tomasello 1, Guido Maione 2 and Paolo Lino 2 1 Department of Electrical, Electronic and Information Engineering, University of Catania, Viale A. Doria 6, 95125 Catania, Italy; francesca.sapuppo@dieei.unict.it (F.S.); vincenzo.tomasello@dieei.unict.it (V.T.) 2 Department of Electrical and Information Engineering, Politecnico di Bari, Via Amendola 126/b, 70126 Bari, Italy; guido.maione@poliba.it (G.M.); paolo.lino@poliba.it (P.L.) * Correspondence: riccardo.caponetto@dieei.unict.it; Tel.: +39-095-7382336


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.

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.
temperature along the tube at t = 0

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 λ 1 , sensor2 at d λ 2 , and sensor3 at d λ 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.

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 α 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.

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 mHz, 0.3 mHz, 0.5 mHz, 0.8 mHz, 1 mHz, 2 mHz, 3 mHz, 4 mHz, 5 mHz, 6 mHz, 7 mHz, 8 mHz, 9 mHz, 10 mHz, and 20 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.

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 α 1 ) by fitting the model response (module and phase) to the experimental data for sensor1 at the normalized distance λ 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 λ 2 and λ 3 parameters by comparing-in the frequency domain-the model output T to the experimental temperature measured, respectively, at distances d λ 2 and d λ 3 .The quality of the identified model on the spatial domain has been measured by comparing estimated distances d λ 2 and d λ 3 (corresponding to the identified λ 2 and λ 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. with 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.

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 Tables 1 and 2, respectively, for the theoretical model (T1, T2, T3) and for the experimental model (T1, T2, T3, α, and α 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 Figures 3-5.The experimental model's enhanced performances in fitting experimental data at the sensor1 location are visible by inspection of the Bode diagrams (see Figures 3-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.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 α = 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.

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 ∆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.

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 λ 2 , d λ 3 with the real sensor distances d λ 2 = 4.2 cm, d λ 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 λ 2 and d λ 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.

Fractional Order and Non-Continuous Materials
The identified values of the fractional order of the experimental model (see Table 4) 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.

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 λ 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 α = 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.

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.

Figure 1 .
Figure 1.Experimental setup.(a) Schematic, where SM (V/K) is the Seebeck coefficient , RM (Ohm) is the Electrical Resistance, and KM the Thermal Conductance (W/K) of the Peltier cell.Such parameters where calculated using the Ferrotec method (http://www.ferrotec.com)through measurements of the module input current I, the hot side temperature th and the cold side temperature tc.NI DAQ: National Instrument c Data Acquisition Board; (b) Entire acquisition setup; (c) Beam photo with details on the sensors position.

Figure 2 .
Figure 2. Nelder-Mead multi-objective optimization algorithm for heat transfer model identification.(a) Identification of model parameters in the frequency domain; (b) Validation of the identified model in the spatial domain.

Figure 3 .
Figure 3. Metal beam filled with air.Theoretical and experimental model comparison in the frequency domain.(a) Module diagram; (b) Phase diagram.

Figure 4 .
Figure 4. Metal beam filled with Styrofoam.Theoretical and experimental model comparison in the frequency domain.(a) Module diagram; (b) Phase diagram.

Figure 5 .
Figure 5. Metal beam filled with metal buckshots.Theoretical and experimental model comparison in the frequency domain.(a) Module diagram; (b) Phase diagram.

Figure 7 .
Figure 7. Step response for the three different controllers: PI, anti-windup PI and FO-PI.

Table 1 .
Theoretical model identified parameters at sensor1.

Table 2 .
Experimental model identified parameters at sensor1.

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

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