Full-state feedback LQR with integral gain for control of induction heating of steel billet

Engineering Science and Technology, an International Journal 55 (2024) 101721 Available online 4 June 2024 2215-0986/© 2024 THE AUTHORS. Published by Elsevier BV on behalf of Karabuk University. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Full Length Article Full-state feedback LQR with integral gain for control of induction heating of steel billet Sohaibullah Zarghoon a , Samuel Emebu b , c , e , * , Radek Matuˇ sů b , Cyril Belavý a , Luk´ aˇ s Bartalský a , Stanislav ˇ Duriˇ s a , Sabir Husnain f , Clara Mendoza Martinez d a Institute of Automation, Informatization, and Measurement, Faculty of Mechanical Engineering, Slovak University of Technology, N´ am. slobody 17, 812 31 Bratislava, Slovakia b Department of Automation and Control Engineering, Faculty of Applied Informatics, Tomas Bata University in Zlín, Nad Str´ anˇ emi 4511, 760 05 Zlín, Czech Republic c Department of Chemical Engineering, Faculty of Engineering, University of Benin, PO Box 1154, Benin City, Nigeria d Department of Energy, Lappeenranta-Lahti University of Technology LUT, Yliopistonkatu 34, FI-53850 Lappeenranta, Finland e Department of Separation Science, Lappeenranta-Lahti University of Technology LUT, Yliopistonkatu 34, FI-53850 Lappeenranta, Finland f Department of Mechanical, Mechatronics and Manufacturing Engineering, University of Engineering and Technology Lahore, Lahore, Punjab, Pakistan ARTICLE INFO Keywords: Induction heating Skin effects Carbon steel billet Semi-numerical method Full-state feedback Linear Quadratic Regulator (LQR) ABSTRACT Induction heating is widely used in industrial furnaces due to its rapid response and energy efficiency. Computeraided modelling and simulation are necessary for the design and optimisation of these furnaces. This article focuses on modelling and simulating induction heating of a one-dimensional spatial partial differential model for a rectangular carbon steel billet. The simulation considered temperature dynamics, including the skin effect and heat loss via radiation. The partial differential model was converted into a state-space model for designing a fullstate feedback Linear Quadratic Regulator (LQR) with integral action. This controller was effectively applied to regulate the billet’s core temperature from 1000 ◦C to 1200 ◦C, under different input, R and state, Q, weighing matrices for the LQR, as well as in the presence of disturbances. 1. Introduction Induction heating is a widely used heating technique in various industries due to its environmental benefits. Traditional heating methods, such as burner heating furnaces, consume more energy when compared to induction heating equipment. As a result, incorporating induction heating systems in hot-forging plants helps reduce energy consumption, as well as air pollution [1]. This method efficiently generates heat directly within the molecular level of metallic inductive materials such as aluminium, brass, copper, steel, and semiconducting materials like silicon carbide [2], through electromagnetic induction. Therefore, induction heat treatment has become a globally recognised technique for repairing, welding both microscopic components and larger parts, offering economic as well as practical advantages [3]. In practice, the design of induction heating furnaces for various applications often involves a lengthy and expensive trial-and-error process. Measuring and evaluating key factors such as current frequency, temperature, and heat distribution can pose challenges. And computer modelling is a suitable approach for such evaluation, as well as the design of induction heating furnaces [4]. The modelling of the induction heating process is based on the Maxwell equations, which constitute a set of Partial Differential Equations (PDEs) as such a nonlinear model that can be difficult to solve. However, the solutions can be obtained or simplified using modelling techniques such as analytical, numerical, and semi-analytical approaches [5–8]. Developing and simulating models of induction heating furnaces are essential for optimising and controlling the process to achieve the desired temperature for the heated billet. This primarily aims to prevent hazards and ensure efficient energy utilisation [9]. The induction heating process can be controlled using either nonlinear or linear controllers. Linear controllers are commonly employed in industrial control systems since many nonlinear processes can be accurately described by linear models under specific nominal conditions through approximation methods like transfer functions and state-space models [10]. Additionally, linear controllers have a wellestablished foundation, and identifying linear models based on process data is relatively straightforward. They provide a good approximation in the vicinity of the designed operating conditions, offer fast response times, and are cost-effective to implement [11]. * Corresponding author. E-mail address: [email protected] (S. Emebu). Contents lists available at ScienceDirect Engineering Science and Technology, an International Journal journal homepage: www.elsevier.com/locate/jestch https://doi.org/10.1016/j.jestch.2024.101721 Received 20 December 2023; Received in revised form 20 April 2024; Accepted 21 May 2024 Engineering Science and Technology, an International Journal 55 (2024) 101721 2 Although the importance of modelling, simulation, and control has been discussed, it is worth noting that most literature studies have primarily focused on the modelling and simulation of induction heating furnaces. These studies have often compared the results obtained from numerical models, particularly the Finite Element Method (FEM), with experimental data, and have generally found good agreement [8,12–17]. However, Demidovich et al. [18] not only considered FEM but also explored and compared other numerical methods such as finitedifference methods (FDM), boundary element method (BEM), infinite element method (IEM), as well as their combinations, without the use of experimental data. Furthermore, some reports have compared analytical models with experimental data [5,15]. For instance, Li et al. [5] developed an analytical model to predict the temperature profile in a planar moving induction heating process. The study demonstrated that the temperature profile and current intensity calculated by the analytical model were similar to the results obtained from the numerical model. Importantly, the analytical model exhibited higher computational efficiency compared to the numerical model. The advantage of reduced experimental and computational complexity offered by the analytical model provides support for studies such as Patidar et al.’s[6] report on the combination of analytical and numerical methods, as well as Areitioaurtena et al.’s[7] semi-analytical method for modelling the induction heating process. Based on the discussion thus far, there are limited reports on the control of the induction heating process using the developed non-linear PDEs models. Camber et al. [19] employed the distributed parameter theory for the control of a continuous steel billet induction heater system as lumped-inputs (i.e., alternating current) and distributed-outputs (temperature profile of billet). Goodwin et al. [20] modelled the induction heating process via a non-linear PDE, controlled the process using nonlinear model predictive control (MPC), as well as compared experimental and simulated results. Furthermore, Roetzer et al. [21] also modelled the induction heating process via a non-linear PDE, controlled and compared the process with different types of controllers (i.e., Proportional integral (PI), and linear-quadratic Gaussian (LQG) controllers). However, as highlighted by Goodwin et al. [20], the incorporation of models to account for skin effect phenomena, and the temperature dependency of parameters, are still limited in literature. Particularly the control of material with the incorporation of skin effect [22] and the temperature dependency of material properties rather than constant values (such as density, emissivity, electrical conductivity, thermal conductivity, and specific heat capacity) [23], can result in a very complex nonlinear system that is difficult to solve as well as control. In practice, this kind of model can be computationally demanding, slow, and expensive to control using nonlinear controllers [24]. Therefore, to reduce cost and increase the controller response time, arises the need to explore linear controllers such as the LQR, to control the induction heating process using simplified approximations of non-linear PDEs models, with both consideration of skin effect, and temperature dependency of material properties. Although attempts have been made to cover this research gas [25,26], however results were still inadequate. For instance, Asadzadeh et al. [26] applied an analytical approximation, specifically the inverse model to simulate the heating process with temperature dependency of the material properties and attempted to consider skin effect. However, the inverse model negated the spatial distribution of temperature along the material, hence being unable to effectively portray the skin effect on the material. Furthermore, Kapusta et al. [25] applied a numerical model for the simulation of induction heating of steel billet and control via distributed-input and distributedparameter-output systems. Although Kapusta et al. [25] considered the temperature dependency of the material on the model, however, the model may not have adequately considered skin effect on the spatial temperature distribution within the depth of the material. Therefore, it is this research gap that this work attempts to cover. Hence the novelty of this work entails the control via linear controller to facilitate the rapid response of the highly nonlinear model of induction heating of steel billet with adequate consideration of skin effect, along with temperature dependency of its material properties. Specifically, the LQR is a reliable linear controller that can be applied, because of its ability to deliver optimal as well as rapid control response to the system per the designer’s specifications (e.g. actuators characteristics or limitations) [27]. LQR is however more suitable for systems with limited steady-state error on application of the LQR feedback gain matrix to the controlled system [28]. Therefore, for systems with high steady-state error the LQR is usually incorporated with integral gain, proportional integral gain, proportional derivative gain, etc [27,29,30]. This additional gain to the LQR gain matrix serves to eliminate the high steady-state error, as well as increase the controller response. However, this may result in an undesirable overshoot on the controller response [31]. Therefore, to limit the effect of controller overshoot, it is necessary to optimise the LQR weighing matrices [32]. The induction heating of steel billet is characterised by high steady-state error, hence in this study, the LQR with an integral gain is a potential linear controller that can be applied. Which to the best of the authors’knowledge, have not been reported in literature. Therefore, this study aims to model, simulate, and control the induction heating process of a rectangular steel billet. This will be achieved through the following objectives: Source data on the properties of steel as a function of temperature; Develop a model and simulate the induction heating process using a semi-numerical method, incorporating skin effect phenomena and the properties of steel as a function of temperature; Linearise the model using the Jacobian method; Apply a linear control scheme (i.e., LQR with integral action in conjunction with optimised weighing matrices) to regulate the core temperature of the steel billet. By accomplishing these objectives, the study intends to provide a comprehensive understanding of the induction heating process and propose an effective control strategy for maintaining the desired temperature in the steel billet. Furthermore, the innovative contribution of this study would entail a predictive estimation and control of the nonlinear spatial temperature variation along the depth of steel billet because of skin effect, as well as the non-uniform energy or current distribution, especially regarding power losses associated with this phenomenon. 2. Methodologies 2.1. Process description and model development The Fig. 1a illustrates an induction heating furnace, in which the coil carries high-frequency alternating current (AC) and generates timevarying electromagnetic fields. These varying fields, in turn, induce eddy currents in the billet, leading to heat generation through the Joule effect. The governing model to describe the induction heating can be deduced from Ampere’s law, Equations (1) a subset of Maxwell’s law [33] in relation to Equations (2), (3), and (4) for magnetic flux density, B (N.A -1 m −1 ), electric flux density, D (C.m −2 ), and conduction current density, J(A.m −2 ) respectively. ∇ × H=J+ ∂ D ∂ t(1) D= ε E (2) B= μ H (3) J= σ E (4) In solving for electric field intensity, E (N.C -1 ), and magnetic field intensity, H (A.m −1 ), it is more convenient to introduce the magnetic vector potential, A(N.A -1 ), and electric scalar potential, φ(V), as given by Equation (5) –(6). The parameters ε , μ ,and σ are respectively emissivity, magnetic permeability (H.m −1 ), and electrical conductivity (S.m −1 ) of the billet. S. Zarghoon et al. Engineering Science and Technology, an International Journal 55 (2024) 101721 3 B= ∇ × A(5) E= ∇φ+ ∂ A ∂ t(6) Substituting Equations (2), (3), and (4) into (1) yields Equation (7d). Noting that ∂ε / ∂ t=0, since ε is a constant, and ( σ E+ ε∂ E/ ∂ t) ≈ σ E. 1 μ ∇B= σ E+ ∂ ( ε E) ∂ t(7a) 1 μ ∇B= σ E+ ε∂ E ∂ t+E ∂ε ∂ t(7b) 1 μ ∇B= σ E+ ε∂ E ∂ t(7c) 1 μ ∇B= σ E (7d) Furthermore, substituting Equations (5) and (6) into (7d), yields Equation (8b) upon rearranging. Noting that J = − σ ∇φ. 1 μ ∇2A= σ (∇φ+ ∂ A ∂ t)= σ ∇φ+ σ∂ A ∂ t(8a) 1 μ ∇2A= − J+ σ∂ A ∂ t(8b) Assuming heating of the steel billet is implemented by alternating current, typically a sinusoidal harmonic, so that the magnetic vector potential, A, and current density, J, can be described as A=Aoej ω tand J=Jsej ω trespectively. Noting that Aois the maximum magnetic vector potential and Jsis the maximum current density on billet surfaces [8]. Therefore, Equation (8b) can be expressed in terms of the sinusoidal harmonic, Equation (9c). In this equation, j implies an imaginary state, ω =2 π f, is the angular frequency and, t is the time of observation. In solving this Equation (9c), the Dirichlet (Ao=0) and Neumann boundary conditions are respectively used in the core and outer surfaces of the billet. 1 μ ∇2Aoej ω t= − Jsej ω t+ σ∂ (Aoej ω t) ∂ t(9a) 1 μ ∇2Aoej ω t= − Jsej ω t+ σ Aoj ω ej ω t(9b) 1 μ ∇2Ao= − Js+j ωσ Ao(9c) The Equation (9c) can be modified to incorporate the skin effect, i.e., the tendency of the source current to be greatest at the surface. This implies that a higher current is induced at the surface of the billet but decreases towards its core as described by Equation (10) [4,7] and Fig. 2. In Equation (10),dis the distance of the core of the billet to its surface, δis the skin depth, Equation (11) [4,7], i.e., the depth below the surface of the billet at which the current density has fallen by the inverse exponent of the maximum source current at the billet outer surface, Jso. In Equation (11), f (Hz) is the frequency of the source current, and μ ris the relative magnetic permeability of the billet to vacuum. Js=Jsoe−d/δ(10) δ= 1/( πμσ f) √≅503  1/( σμ rf) √(11) The generic magnetic vector potential on the billet surface, Ao, can be deduced from Equation (9c) –(11), and subsequently the source heat Fig. 1. Illustration of induction heating of a rectangular steel billet. Fig. 2. Illustration of the skin effect for a rectangular billet. S. Zarghoon et al. Engineering Science and Technology, an International Journal 55 (2024) 101721 4 flux, q, can be calculated as given by Equation (12) in terms of the conduction current density, J, which can also be expressed in terms of Aoas given by Equation (13) [8]. In this equation, Rsignifies the real state, and (AoAo)*indicates a product of complex conjugate. J= σ E= − j ωσ Ao(12) q=[R(J)]2 σ =1 2 ω 2 σ (AoAo)*(13) The heat diffusion equation of induction heating on the billet can be expressed by Equation (14). In the equation, ρ (Kg.m −3 ) is the density, Cp(J. Kg −1 . K −1 ) is the specific heat capacity, and k (W.m −1 . k −1 ) is the thermal conductivity of the billet. ρ Cp ∂ T ∂ t=k∇2T+q (14) The Equations (9c) and (14) are solved using the ordinary differential equation solver in MATLAB (specifically, ode15s) in collaboration with Equations (10), (11), and (13). In solving these equations the spatial temperature distribution or partial differential, ∇, was approximated using the Central scheme of the Finite Difference Method (FDM). This results in N th number of equations that represents each space point, i = 1, 2, 3…, Nxas shown in Fig. 1b. Furthermore, in Equation (14) the Neumann boundary condition is applied to the core, and for the outer surface of the billet, radiation heat loss, i.e., Equation (15), is applied. The Equation (15) relates the environment or surrounding air temperature, TE, to the outer temperature of the billet, T, via the StefanBoltzmann constant, σ *, thermal conductivity, and emissivity. k∇T= σ * ε (T4−T4 E)(15) These boundary conditions applied are based on the operating condition of industrial induction heating furnaces as well as popularly reported in literature, for which both convection and radiative heat losses are considered [34,35]. However, considering that the radiative heat loss is ≫convective heat loss[36], only radiative heat loss was considered in this work. 3. Process control theory To effectively control the system, it is necessary to initially investigate the system properties such as the system’s dynamics, stability, and performance (i.e., the difference in steady-state error). Prior observations of induction heating systems from literature [7,15] indicate that the system is generally stable but may exhibit significant steady-state error. To mitigate this steady-state error, a full-state feedback Linear Quadratic Regulator (LQR) can be employed. In general, the LQR control aims to control linear system state variables using optimal control cost function in finite horizon or infinite horizon. It aims to deduce gain matrices that minimises a cost function to ensure adequate system stability and control performance. Therefore, the choice of weighing matrices (i.e., state weighing and input weighing matrices) is crucial to adequately design the LQR [37]. The full-state feedback LQR methodology is based on the concept of a cost function, denoted as, j, wherein the objective is to determine an optimal input or actuator effort, u, that changes the state vector, x, through minimisation of Equation (16) by adjusting the state weighing matrix, Q. The input weighing matrix, R, is utilised to penalise the actuator effort and achieve the desired control performance. j=∫(xTQx +uTRu)dt (16) Q and R are diagonal matrices, which must be positive-definite to ensure that when multiplied by the state, x(i.e., temperature), and input, u(i.e., source heat flux), the column vectors yield xTQx ≥0 and uTRu >0. The dimension Q and R matrices are defined by the number of states, s, and inputs, w, as given by Equation (17). Q=⎡ ⎢ ⎢ ⎣ Q1,10⋯0 0 Q2,2⋯0 0 0 ⋱ ⋮ 0 0 ⋯Qs,s ⎤ ⎥ ⎥ ⎦ &R=⎡ ⎢ ⎢ ⎣ R1,10⋯0 0 R2,2⋯0 0 0 ⋱ ⋮ 0 0 ⋯Rw,w ⎤ ⎥ ⎥ ⎦(17) Typically, the LQR using the infinite horizon cost function, Equation (16), can be structured as given by Equation (18). This follows the traditional zero-terminal quadratic cost function, applied in conjunction with decretised state-space model, Equation (22) –(23) as constraints. min uj=1 2∑ N k=1 (x(kt)*−x(kt))TQ(x(kt)*−x(kt))+u(kt)TRu(kt)(18) Subject to: x(kt+1) = Ax(kt) + Bu(kt) y(kt) = Cx(kt) In Equation (18), the time scalar of the process model is divided into discrete points, kt∈ [1 Nt], and x(k)*is the reference value of outputs. The solution to Equation (18) is equivalent to solving the Riccati Fig. 3. Illustration of a full-state feedback system with integral action. Table 1 Simulation parameters utilised for the inductive heating furnace. Parameter, symbol Value, unit Duration of simulation,t 5000 s Thickness of billet,ϑ0.12 m Number of discretised spatial points,Nx5 Number of discretised time points,Nt100 Differential distance between discretised spatial points, Δxϑ/2Nx Stefan-Boltzmann constant, σ *5.6704 ×10 8 Curie temperature of steel,Tc750 ℃ Permeability factor,n 0.5 Magnetic permeability of vacuum, μ o4 π ×10 -7 H.m −1 Relative permeability of carbon steel at room temperature, μ ro 10H.m −1 Frequency of alternating current, f [17] 100 Hz Maximum current density on billet surface,Jso 50 ×10 5 A.m −2 Surrounding air temperature,TE35 ℃ Initial temperature of billet,Tin 1000 ℃ Initial maximum magnetic vector potential,Ao0 N.A -1 S. Zarghoon et al. Engineering Science and Technology, an International Journal 55 (2024) 101721 5 equation [38]. In solving the Riccati equation, to deduce the input, u, Equation (21), the so-called “state transition matrix”, P, must be deduced via Equation (19), to estimate the controller gain, K, Equation (20). P=ATP+PA−PBR−1BTP+Q (19) K=R−1(BTP)(20) u= − R−1(BTPx)(21) The full-state feedback LQR with integral action within the feedforward pathway connecting the error comparator and the system is depicted in Fig. 3. In summary of the diagram shown in Fig. 3, r is the desired setpoint, u is the predicted controller input, A, B,and C are matrices representing the process model, Kxis the LQR feedback gain matrix and Kiis Fig. 4. Properties of carbon steel as a function of temperature. S. Zarghoon et al. Engineering Science and Technology, an International Journal 55 (2024) 101721 6 the integral gain. The addition of integral action to the control system allows the controller to track the setpoint (instead of relying on the classical null stabilisation) and, as such ensure a better perturbation robustness [39]. The aim is to stabilise the controlled system as well as integrate the steady-state error. However, the integral action can lead to controller overshoots, therefore to resolve this issue[40], adequate measures must be taken to optimise the LQR through its weighting matrices, i.e. Q and R. Adequate design or tuning of these two matrices are the most important parameters of the LQR which can be tedious and challenging[41]. To apply the process model to the LQR, its state and output dynamics are usually expressed as state-space models, Equation (22) and (23). The state-space model defined by matrices, A, B,and C can be deduced from the linearisation of the dynamic models, Equation (9c) and (14). Note that C is a matrix of zeros and ones that defines the outputs from the states. ˙ x=Ax +Bu (22) y=Cx (23) ˙ I=r−y=r−Cx (24) The resulting error from the system, ˙ I, can be deduced from the difference between the desired setpoint, r, and output, y, of the system. Equations (22) and (24), can be used to design a new state-space model with integral action as given by the state-space matrices, Equation (25), Fig. 5. Illustration of temperature dynamics for the discretised billet surfaces. Table 2 Evaluation metric for the comparison of developed model and data for carbon steel properties. Carbon steel properties R 2 Density, ρ 0.9985 Specific heat capacity,Cp0.9997 Thermal conductivity, k 0.9998 Emissivity, ε 1.0000 Electrical conductivity, σ 0.9950 Fig. 6. Illustration of the resulting magnetic vector potential from the current density on billet surfaces. S. Zarghoon et al. Engineering Science and Technology, an International Journal 55 (2024) 101721 7 which can be further simplified into Equation (26). [˙ x ˙ I]=[A 0 −C 0][x I]+[B 0]u+[0 1]r (25) ˙ x=AIxI+BIu+Brr (26) The controller input, u, of the full-state feedback LQR can be expressed as given in Equation (27). Substituting Equation (27) into (26) yields Equation (28). The gain matrices of the controller are given as K = [Kx⋮Ki]. Typically Kxis a single-row matrix whose column size is the same as the number of state variables, x, i.e. the five temperatures of each discretised point, Kx=[Kx1, Kx2, Kx3, Kx4, Kx5] and Kia single value. Therefore, K =[Kx1, Kx2, Kx3, Kx4, Kx5, Ki]. u= − KxI(27) ˙ x= (AI−KBI)xI+Brr (28) In the aforementioned equations, the gain matrix, K, can be determined through the utilisation of the “lqr”function in MATLAB R2021b by solving the Riccati equation, Equation (19), and the control system is implemented with SIMULINK R2021b in the same design framework given by Fig. 3. 4. Results and discussion 4.1. Simulation of induction heating To simulate the induction heating process, Equations (9c) and (14) are considered as one-dimensional spatial differential, ∇, along the xaxis, and discretised into five segments (i.e., Nx=5), Fig. 1b. The central difference approximations are used to derive five equations to represent Equations (9c) and (14) each. These equations are then solved simultaneously with Equations (10), (11), and (13). The simulation results, obtained using the parameters specified in Table 1, and Equations (29) – (34), are presented in Figs. 4 –5. The Equation (29) –(34) shows the temperature dependence (i.e., T in kelvins, K) of the density, ρ [42,43], specific heat capacity, Cp [42,44–48], thermal conductivity, k [42–48], emissivity, ε [49], electrical conductivity, σ [50], and relative magnetic permeability, μ r [34,51,52], of carbon steel. The Equations (29) –(33) were curve-fitted as a function of temperature, using data from the respective highlighted literature. However, Equation (34) is as reported in literature. The resulting comparison of the developed models and data is shown in Fig. 4, with their respective R-squared (R 2 ) evaluation metrics given in Table 2, and the result showed adequate fit (R 2 >0.95). ε =⎧ ⎨ ⎩ T≤663→0.28 653 ≤T≤793→ T≥793→0.69 0.6999exp(−(T−848.7278 203.6169 )2)(32) μ r= μ μ o =⎧ ⎪ ⎨ ⎪ ⎩ T≤1023→1+ ( μ ro −1)(1−(T−273 TC(◦C))n) T≥1023→1 (34) In Fig. 5a, upon the application of source current to the billet surface, Js, the outermost surface (i.e., Nx5) initially experiences a decrease in temperature due to radiation heat loss, Equation (15) to the surrounding air film. This heat loss occurs because of the sharp temperature contrast ρ (kg.m3)=⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 273 ≤T≤1673K→7849exp(−(T−146.6 4858 )2)+367.2exp(−(T−1914 468 )2) 1673 ≤T≤1803K→12420exp( − 0.0003T) 1803 ≤T≤1973K→7120exp(1.635E−8T) (29) Cp(J.kg−1K−1)=⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 273 ≤T≤993K→4780exp(−(T−996.9 13.39 )2)+710.2exp(−(T−2046 2627 )2) 993 ≤T≤1973K→4412exp(−(T−990.6 42.22 )2)+603.7exp(−(T−21320 356500 )2)(30) k(W.m−1K−1)=⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 273 ≤T≤1673K→12.13exp(−(T−349.3 384.4)2)+3.292E+15exp(−(T−641800 113100 )2) 1673 ≤T≤1773K→23.32exp( − 0.0002T) + 1.706E−13exp(0.0197T) 1773 ≤T≤1973K→259.7exp(7.473E−7T) (31) σ (S.m−1) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 273 ≤T≤1773K→5364015.153exp(−(T−4977.4071 4746.0594 )2)+239071.5586exp(−(T−781.0628 317.3960 )2) 1773 ≤T≤1803K→1.4860exp(0.0074T) 1803 ≤T≤1973K→2.2269T2−8614.6218T+9202966.681 (33) S. Zarghoon et al. Engineering Science and Technology, an International Journal 55 (2024) 101721 8 between the billet surface (1000 ℃) and the surrounding air temperature (35 ℃). This heat loss takes place even though the skin effect theorem predicts that Jsshould be maximum at Nx5as such should experience an instantaneous increase in temperature (i.e., in the absence of heat loss). Therefore, due to this counteracting effect, after this initial heat loss, the surface begins to steadily increase in temperature. In contrast to the observation at the outermost surface, the inner surfaces (i.e., Nx4to Nx1) experience an instantaneous increase in temperature. During the initial phase of heating (≤205 s), the middle surface, Nx3, experiences the highest heating rate, followed by Nx2and Nx1. This is because the skin effect has a greater influence on the Nx3surface compared to radiation heat loss. On the other hand, the Nx4surface experiences the lowest heating rate due to its proximity to the outer surface, resulting in more radiation heat loss to the surrounding surface, Nx5. However, as the heating process continues (>205 s), there is more heat accumulation toward the core, Nx1, of the billet due to negligible radiation heat loss to the surroundings. This continues even though, according to the skin effect theorem, Jsdecreases towards the core, Nx1. Consequently, as the heating progresses, when all surfaces reach their steady-state temperatures at approximately 4400 s, the core temperature finally reaches the highest temperature. The preceding surface temperatures follow the order: Nx1(≅1259 ◦C) >Nx2(≅1258 ◦C) > Nx3(≅1246 ◦C) >Nx4(≅1217 ◦C) >Nx5(≅1164 ◦C). Considering that the spatial discretisation of steel billet was linear, i.e., Δx=ϑ/2Nx. It would have been expected that the spatial distribution of temperature along the depth of the steel billet should be correspondingly linear. The results, (Tx1−Tx2≅1◦C, Tx2−Tx3≅12 ◦C, Tx3−Tx4≅29 ◦C, Tx4−Tx5≅53 ◦C) and Fig. 5a, however shows otherwise. This highly nonlinear spatial temperature distribution can be attributed to the skin effect, as well as the nonlinear boundary conditions (i.e., the Dirichlet boundary condition via radiative heat loss on Nx5, and the Neumann boundary on Nx1). Furthermore, Fig. 5b shows a surface plot illustrating the dynamic changes in temperature across the discretised surfaces as time progresses. It visually represents how the temperature evolves throughout the billet. In addition, Fig. 6 showcases the resulting magnetic vector potential, A, on the billet surfaces. As anticipated, A, decreases towards the core of the billet, following the influence of the skin effect on Js. This demonstrates the distribution and behavior of the magnetic field within the induction heating process. The results of the simulation of this model thus far conform to those reported in literature. Luozzo et al.[34] for carbon steel tube, Hansson and Fisk [35] for stainless steel tube, Jelicic [53] and Baldan [54] for solid steel billet, all reported similar results that verify the result of this work. Considering that both radiation and convection heat losses on the surface were considered, the inner surface also experienced higher temperatures than the outer surface, as found in this study. However, in this work, the temperature difference was higher compared to those reported in literature. The difference can be attributed to the obvious difference in initial condition (i.e. 1000 ◦C in this work as compared to about 20 to 35 ◦C for other work), emissivity, as well as the assumption of convective heat transfer considered. Furthermore Jankowski et al. [55] although considered a cylindrical steel billet, the trend of the result also showed similarity with that of this work, even though the spatial Fig. 8. LQR temperature control for billet core surface based on changes in input weighing matrix. Fig. 7. LQR temperature control for billet core surface based on changes in state weighing matrix. S. Zarghoon et al. Engineering Science and Technology, an International Journal 55 (2024) 101721 9 temperature distribution wasn’t considered in their work. Janutien˙ e and Maˇ zeika [56] however reported higher temperatures on the billet surface than at the core of a solid cylindrical carbon steel billet, contrary to the result of this work. This discrepancy is a result of negligible heat loss assumed at the surface, as such due to the prevailing skin effect, the result is expected for such an assumption. 5. Process control of induction heating To control the induction heating process, the state-space model described by Equation (22) –(28) is utilised. To deduce the A and B matrices required for this purpose, the dynamic models, Equations (9c) and (14) in conjunction with Equations (10), (11), and (13) are linearised using Taylor’s series or Jacobian linearisation theorem, as extensively described in literature [11]. The value of Jso given in Table 1 was used as the steady-state value of the maximum current density, Jso, for the linearisation process, along with the preliminary steady-state temperature values of surfaces Nx1to Nx5i.e., T1=1259 ◦C,T2= 1258 ◦C, T3=1246 ◦C, T4=1217 ◦C andT5=1164 ◦C. These temperature values are used for the preliminary calculation of carbon steel properties as highlighted in Equation (29) –(34). The resulting matrices obtained from the linearisation process are given by Equation (35), with the corresponding output matrix, C = [ 10000], and potentially the disturbance or direct transmission matrix, D =0. The typical application of the matrices A, B, and C to the state-space model, Equation (22) – (23), is given in Equation (30b). To verify the earlier assumption of system stability, the open-loop stability of the induction heating system was evaluated based on the eigenvalues, λ, obtained from the determinant of the matrix A and its identity matrix, I, i.e, |A−λI|= 0. The system stability is confirmed if all eigenvalues have negative real parts, i.e., λ<0. The resulting eigenvalues (−0.0017, −0.0308, −0.0954, −0.1850, −0.1599) affirm the system is stable. With the matrices A and B determined, and stability confirmed, the LQR control system can be designed. The resulting AI, and BImatrices are given by Equation (38), with the previously specified matrix, Br=[ 0 1 ]T. The LQR, irrespective of the specific values assigned to the Q and R parameters, serves as design elements to penalise both the state variables and the control signals [57,58]. In this study, the Q and R matrices were tested under two case studies. In casestudy-1, R was kept constant while Q was varied, and in case-study-2, Q was kept constant while R was assigned different values. The objective of these case studies was to evaluate the system’s output when attempting to reach a setpoint of 1200 ℃at the core of the billet. In case-study-1, the following values of Q and R matrices were respectively applied, i.e., Q1=diag (1,1,1,1,1,100), Q2=diag (1,1,1,1,1,150), Q3=diag (1,1,1,1,1,200), Q4=diag (1,1,1,1,1,300), and Q5=diag (1,1,1,1,1,350), Q6=diag (1,1,1,1,1,600) as well as R = 2, i.e., the R matrix is kept constant. Fig. 7a illustrates the response of the system upon the implementation of a full-state feedback LQR, with the resulting gain matrices given by Equation (39). The observed behaviour indicates that as the value of Q increases, the system approaches the setpoint with minimal overshoot and an impressive settling time both before and after the occurrence of disturbance. Generally, in a control system, a faster controller response or settling time to reach the setpoint is desirable, while overshoot is considered undesirable. Among the tested values of Q, Q3demonstrates the best performance as it results in minimal overshoot and an adequate settling time. Increasing the value of Q generally leads to faster settling time but at the expense of higher overshooting due to the initial overreaction of the controller [59], which may have resulted in or amplified by the integral action added to the LQR. This can be observed in the response of Q6(i.e., the highest Q value), which achieves the fastest settling time but with the highest overshoot, while Q1(i.e., the lowest Q value) exhibits the opposite behaviour. Overall, the controller proves to be effective even in the presence of disturbance at 4000 s, as the system quickly returns to its setpoint. The resulting maximum current density, Jso, for achieving temperature control of the billet core surface is depicted in Fig. 7b. A=⎡ ⎢ ⎢ ⎢ ⎢ ⎣ −0.0908 0.0908 0.0000 0.0000 0.0000 0.0454 −0.0908 0.0454 0.0000 0.0000 0.0000 0.0455 −0.0909 0.0455 0.0000 0.0000 0.0000 0.0456 −0.0911 0.0456 0.0000 0.0000 0.0000 0.0915 −0.1091 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ &B =⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 5.1363E −8 1.7614E −7 3.2884E −7 4.6366E −7 5.3544E −7 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (35) ˙ x=⎡ ⎢ ⎢ ⎢ ⎢ ⎣ dT1/dt dT2/dt dT3/dt dT4/dt dT5/dt ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ =⎡ ⎢ ⎢ ⎢ ⎢ ⎣ −0.0908 0.0908 0.0000 0.0000 0.0000 0.0454 −0.0908 0.0454 0.0000 0.0000 0.0000 0.0455 −0.0909 0.0455 0.0000 0.0000 0.0000 0.0456 −0.0911 0.0456 0.0000 0.0000 0.0000 0.0915 −0.1091 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⏟⏞⏞⏟ A ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ T1 T2 T3 T4 T5 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⏟⏞⏞⏟ x +⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 5.1363E −8 1.7614E −7 3.2884E −7 4.6366E −7 5.3544E −7 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⏟⏞⏞⏟ B q ⏟⏞⏞⏟ u (36) Fig. 9. LQR temperature control for billet core surface based on specific state and input weighing matrix. S. Zarghoon et al.