Application of nonlinear model predictive control to an industrial ...
Application of Nonlinear Model Predictive Control to an Industrial Induction Heating Furnace G.C. Goodwin, R.H. Middleton, M.M. Seron and B. Campos
A BSTRACT Induction Heating Furnaces are used extensively in industry. The basic principle is that induced eddy currents are used to heat a ferromagnetic material as it passes through a series of coils. Because of the importance of such systems, there has been on-going interest in their design and operation. Past work includes model development from physical principles and optimal design of operational practices. However, previous work has invariably been based on open-loop strategies. Our work is aimed at the design of a closed-loop control strategy incorporating feedback from the available measurements. This paper reports initial work including model development and calibration together with preliminary control system design. Proposed future work includes full scale industrial implementation. I. I NTRODUCTION Induction Heating Furnaces are frequently used in industry. These systems consume extremely large amounts of energy. Indeed, the annual cost of electrical energy for a typical furnace can exceed US$1M. Hence, there is substantial financial motivation to reduce energy consumption. Beyond the cost of energy, there also exist strong environmental drivers since electrical energy is a major contributor to green house gas emissions and related issues. Simulation and control of induction heating furnaces are topics which have been extensively studied. There are presently two main ways in which these systems are designed. The first is where a computationally expensive simulation is run over a long period of time in order to find a theoretically optimal physical setup ([1], [2], [3], [4], [5], [6], [7], [8]). The models typically involve finite element analysis with a very large number of elements ([9], [10], [6], [11], [12]). The major shortcoming of these methods is that there is no flexibility in the design. Thus minor miscalculations can result in large errors appearing when physically implemented. The other main method used is closed loop temperature control. In this, the surface temperature of the bar is measured via pyrometers and any deviations from desired values are fed into a controller ([13], [14], [15]). While this offers some compensation for model error and disturbances, current implementations are of an empirical nature. Also, as we discuss later, surface temperature may be a poor guide to internal states of the process. The models in current use are unsuitable for on-line control studies since they are either excessively complex (i.e. based on a full partial differential equation description) or too simplistic (i.e. aimed at standard operational scenarios and/or design of PID controllers). Our proposal is to use nonlinear model predictive control based on a simplified model with feedback. This is done by creating a low complexity, “control relevant” model and using its predictions as part of a closed-loop control strategy. To date, the following steps have been carried out: (i) model development (ii) model calibration using industrial data (iii) simulated control studies using nonlinear MPC (iv) simulated observer design using nonlinear optimization (v) comparison of initial proposed control solutions with current industrial practice This paper is based on the plenary paper ”Opportunities and Challenges in the Application of Nonlinear MPC to Industrial Problems” presented at NMPC 2012. G.C. Goodwin, R.H. Middleton and M.M. Seron are with School of Electrical Engineering and Computer Science, University of Newcastle, 2308, Australia Graham.Goodwin, Richard.Middleton, [email protected] B. Campos is with Donhad Pty Ltd, Newcastle. [email protected]
Future work is aimed at full scale industrial implementation of these ideas. This phase is planned over a time horizon of several years since substantial practical issues are involved. This problem has a range of associated challenges including: (i) Whilst one can reasonably assume radial symmetry of temperatures within the rod, the temperatures are not constant in the longitudinal direction. This complicates model development and can lead to excessive model complexity. For example, if one were to divide the radial and longitudinal spatial dimension into 1000 segments each, then one would have a model containing a million state variables! (ii) Due to skin effect phenomena ([16]), much of the heating takes place in a narrow surface region. This region can be of the order of 1/2 mm or less in thickness and is associated with rapid thermal heating. By way of contrast, the core heats slowly. Hence the model is inherently “stiff” with a wide range of associated time constants. Moreover, the skin effect region progressively moves into the core as the outer temperature passes the Curie point. (iii) The available measurements are very limited being restricted, in essence, to surface temperature measurements (made by pyrometers) and the external currents and voltages. It may also be possible to use the forces measured in subsequent forming processes. (iv) The system has many degrees-of-freedom in the input including line speed, voltages, currents and frequencies in the induction coils (of which there may be 10 or 20 in a given set-up). (v) The system needs to satisfy a large number of constraints. These include, but are not restricted to: 1) A lower bound on minimum exit temperature. 2) An upper bound on maximum exit temperature. 3) Constraints on temperature gradients in both time and space to avoid metallurgical stress on the material. 4) Bounds on operating voltages, currents, frequencies and switching times dictated by the driving power electronics. (vi) All of the parameters appearing in the model are highly temperature dependent. Moreover, the temperatures within the rod range from ambient to above the Curie point ([17]). Some parameters, e.g. relative permeability, change by several orders of magnitude as the rod traverses a given coil. Hence the system is highly nonlinear. The layout of the remainder of the paper is as follows: In Section II, we give a brief overview of the physical apparatus. In Section III, we describe a control relevant model. Section IV discuss model parameter variations. Section V describes model validation against industrial data. Section VI describes the proposed control law and gives simulated performance results. Section VII describes the proposed observer design. Section VIII describes challenges associated with the problem. Conclusions are presented in Section IX. II. T HE P HYSICAL A PPARATUS A schematic diagram of a typical furnace is given in Figure 1. Velocity v u1 Fig. 1.
u2
u3
uNc
Schematic Diagram of Induction Heating Furnace {ui ≡ (Voltage, Current, Frequency); i = 1, . . . , Nc }
As shown in Figure 1, the key features of the apparatus are that a cylindrical rod is passed through a succession of coils. These coils are capable of being operated at different voltages, currents and frequencies. Moreover, the rod can be moved at different speeds through the furnace. The coils induce currents in the rod which, due to hysteresis effects, heat the rod. Most of the current in the rod is initially confined to the surface due to skin effect. Once the temperature passes the Curie point, then the relative permeability drops from about 200 to near 1 (i.e. that of air). The skin depth then increases. This, together with radial heat conduction, continues until the rod temperatures are all above the Curie temperature. Heating still continues but the skin depth will be now an order of magnitude larger; being 5 to 10 mms if the frequency is unchanged. For this reason, and also to improve the efficiency of the process with reduced permeability, it is usual to increase the frequency in this region to aid heating. The main goal in operating the furnace is to satisfy the constraints described earlier in Section I whilst minimizing energy consumption.
III. A C ONTROL R ELEVANT M ODEL The system is described by a set of partial differential equations which, inter alia, include Fourier’s Law (for heat conduction), Newton’s Law (for Thermal heating) and Maxwell’s equations (for electromagnetic phenomena). Many authors have begun modelling at this level. However, this inevitably leads to high order models with many thousands of state variables (see observation (i) in Section I). Our approach is to develop a control-relevant model. After much initial work with a larger model, we were led to the following: (i) We assume the electromagnetic aspects are in steady state and are thus described by algebraic relationships. (ii) We assume that the majority of heat conduction occurs in the radial direction. (iii) We focus on a single longitudinal zone and track it as it passes through the furnace. (iv) We assume that downstream zones have a minimal impact on electromagnetic effects in the current zone. (This is justified, in part, by the skin effect phenomena.) This approach leads to a nonlinear-state-space model in algebraic-difference form. Details of the model are outlined in the sequel. We note that all physical parameters are highly temperature dependent. This aspect will be discussed further in Section V. Our model focuses on one coil. However, the model applies mutatis mutandis to multiple coils and the region between coils (where the external coil current is zero). We consider a slice through the rod as shown in Figure 2. We assume this slice is in a longitudinal “zone” located between (j − 1) ∗ d and j ∗ d from the entry to the coil where d is the width of each slice. radiated heat loss
θN j
θ2j θ1j J1j r1
r2
rN
r
J2j JN j
Fig. 2.
Radial distribution of temperature and current density in the jth zone through the bar
We assume Nc coils. The concentric rings are assumed to have temperature, radius, and current density θ ij , ri , Jij , ; i = 1, . . . , N, j = 1, . . . , M respectively, where r N is the radius of the rod and M d is the length of the coil. A. Heat Transfer Model We denote the temperature in the i th ring in zone j at time step k as θ ij (k). There are three mechanisms for heat transfer: (i) heat conductivity due to temperature gradients (ii) heat generated due to circulating current (iii) radiation from the surface
The heat, Q1 , due to temperature gradients satisfies Fourier’s Law ([18], page 23), i.e., dθ dQ1 = −KA dt dx
where
(1)
dQ1 dt
is heat flux per unit time is heat conductivity (which is a function of temperature) A is the area Since the area is a function of radius, we have K
dθ dQ1 = −2Kπd r dt dr
(2)
where d is the width of the zone. Dividing (2) by r and integrating we obtain the heat flux into a region per unit time dQ1ij dt
=−
2Kπd (θij − θi−1,j ) n (ri /ri−1 )
(3)
Next we note that the heat generated due to circulating current is ρ|J| 2 (W/cm3 ). To use this, we need the volume of the ith ring in zone j , ρij ρij 2πri = Aij (ri − ri−1 )d 2 2 Vij = π(ri − ri−1 )d
Rij =
(4) (5)
Let ρij = ρ(θij ) be the resistivity in the ith ring in zone j . Hence heat per unit time due to circulating current, Q 2ij , satisfies dQ2ij = ρij |Jij |2 Vij dt Finally the temperature rise satisfies Newton’s Law ([18]), namely (for all but the outer ring): dQ2ij dQ1ij dθij 1 = + dt c(θij )δVij dt dt
(6)
(7)
where c is the temperature dependent specific heat capacity [nominally 0.45 J/(gK) for iron], δ is the metal density [7g/cm3 for iron]. If we use a time increment of Δ seconds, then (7) can be approximated by θij (k + 1) = θij (k) dQ2ij dQ1ij 1 + +Δ c(θij (k))δVij dt dt
(8)
θij (k + 1) = θij (k) d +Δ c(θij (k))δVij θi+1,j (k) − θij (k) · K(θij (k))2π d n (ri+1 /ri ) θij (k) − θi−1,j (k) d − K(θij (k))2π n (ri /ri−1 ) + ρij |Jij |2 Vij
(9)
Combining (3), (7), (8) we have
Equation (9) holds for i = 2, . . . , N − 1 ; j = 1, . . . , M .
For i = 1, we have no heat conducted out. For i = N , we replace the heat conducted into the cell by heat 4 − θ4 radiation loss. We model heat radiation by; a “grey body” law: σ2πr N d θN j amb , where σ is the Stefan Boltzman constant ([19]) and θ amb is the ambient temperature. We note that the above equations, in particular (6), depend upon the current density J ij (k). Equations for the latter variables are developed in the sequel. B. Electromagnetic Modelling We model the electromagnetics within the bar via the complex sinusoidal steady state version of Maxwell’s equations, in cylindrical coordinates. These equations are: d E = jωo μH (10) dr d H = −J (11) dr where j is the imaginary unit, E is the circular component of the electric field, H is the axial component of the field intensity, ωo is the excitation angular frequency, μ = μ o μr is the magnetic permeability, and J is the current density given by 1 E (12) ρ Note that E , H and J in (10) - (12) are complex valued functions of radius, r . To solve (10), (11), (12) we do a spatial discretization of E and H , as vectors representing the fields in each zone. For a fixed zone, this is J=
¯ [E(0), E(r1 ), . . . , E(rN )]T E ¯ [H(0), H(r1 ), . . . , H(rN )]T H d The spatial derivative operator, dr , can then be represented by a matrix multiplication with ⎡ ⎤ 1 1 − r(1) 0 · · · · · · 0 r(1) ⎢ 0 ⎥ 1 1 0 ... 0 − r(2)−r(1) ⎢ ⎥ r(2)−r(1) ⎢ ⎥ Γ⎢ .. .. .. .. .. .. ⎥ . . . . . . ⎣ ⎦ 1 1 0 ... ... . . . − r(N )−r(N −1) r(N )−r(N −1)
Then (10), (11), (12) can be approximated by the implicit, sparse, linear system of equations: ¯ ¯ 0 diag{jω0 pi } Γ 0 E E ¯ = diag{− 1 } ¯ 0 0 Γ H H pi
(13) (14)
(15)
(16)
together with equality constraints E(0) = 0 and H(r N ) = Nc Ic /L (the field intensity due to the Nc turn coil current, Ic , over a length, L). Solution of (16) in each zone allows computation of the current densities for (12). C. Modeling the Motion of the Bar The final step in the model is to describe the motion of the bar. We do this by noting that, if we choose the width of the rings as d = vΔ where v is the bar velocity and Δ is the time step, then θi,j+1 (k + 1) = θi,j (k) for i = 1, . . . , N ; j = 1, . . . , M − 1
(17)
The boundary conditions for this equation are ambient (for the first coil) or the exit temperatures from the previous coil. The above strategy allows us to combine the subscripts j and k, i.e. we use θ i (k) to denote the temperature in ring i after it has passed k steps into the coil.
D. Final Model Format We note that the model as described above has two key features: (i) Progressively heat is generated via induction as the rod transits the furnace. (ii) There is a mixture of difference equations (for temperature) and algebraic equations (for electromagnetic variables). j So as to express the model in a familiar notation, we use x jk to denote the vector of temperatures θ 1j (k) . . . θN (k) in the j th zone. Also, let ejk denote the corresponding vector of electromagnetic variables. Then the model described above takes the following form:
(18) xjk+1 = f 1 xjk , ejk where ejk is a nonlinear function of xjk . Next we note that the electromagnetic model is linear in e jk but highly nonlinear in the temperature. Thus these equations can be expressed in the form of a linear matrix equation: Dk ejk = bjk
(19)
where Dk depends nonlinearly on temperature i.e. x jk . The right hand side of (19), i.e. b jk , depends linearly on the external current. For future use, we denote the “control” variables i.e. the external current and frequency by u j . Then, for each fixed temperature xjk , equation (19) can be inverted to express e jk as a function of uj and xjk , i.e. ejk = Dk−1 bjk
(20)
Finally, substituting (20) into (18) we obtain a model having the familiar form:
xjk+1 = f xjk , uj ; j = 1, . . . , Nc
(21)
with boundary conditions: x11 = ambient temperature j+1 x1 = xjM ; j = 1, . . . , Nc − 1 IV. M ODEL PARAMETER VARIATIONS One of the key challenges in this problem is the fact that all physical parameters are highly temperature dependent. As an illustration, typical temperature dependencies for heat conductivity, resistivity and relative permeability are shown in Figure 3, Figure 4 and Figure 5, respectively. Note that the dependence of permeability on temperature is particularly significant. k W/(mK) 50
25
800 100 Fig. 3.
temp o C
Heat conductivity as a function of temperature
In particular, the relative permeability changes by 200:1 over a very narrow temperature region as the material approaches and passes the Curie temperature. This phenomena is a major factor in the control problem.
ρ
1.2*10−4
0.2*10−4 800 100 Fig. 4.
temp o C
Resistivity as a function of temperature
μr
200
1 800 100 Fig. 5.
temp o C
Relative permeability as a function of temperature
V. M ODEL VALIDATION AGAINST E XPERIMENTAL DATA The model parameters were chosen based on published empirical data. The model was then run using typical values used in industry for the control variables. Figure 6 compares the predicted surface temperature and real surface temperature data obtained from a process of similar characteristics in a local industry. We see from Figure 6 that the model performs qualitatively the same as the real process. However, some important discrepancies are evident. These discrepancies will be addressed in Section VII. VI. P ROPOSED O N - LINE C ONTROL L AW The basic principle behind the development of the control law is to use finite horizon optimization in a nonlinear MPC control law. However, there are several departures for the usual MPC type controllers described in the literature (see for example [20]). In particular, once a control sequence u 1 , . . . , uNc has been determined, then it can all be applied since the “time index” actually relates to the spatial dimension of distance through the furnace. Details are given below. We use a prediction model of the form:
x ¯jk+1 = f x ¯jk , uj ; j − 1, . . . , Nc (22) We use a cost function which captures energy utilization in the N c coils. Nc J= uj j=1
(23)
Real Data and Uncorrected Model 1200 Real Data Uncorrected Model 1000
Temperature (C)
800
600
400
200
0
0
200
400
600
800
1000
1200
1400
1600
Axial Distance (cm)
Fig. 6.
Real surface temperature data (black stars) and output of the uncorrected model (solid line).
We optimize (23) subject to (22) and the constraints listed in Section I, (v). For example, the constraint on the minimum exit temperature would be (assuming the minimum temperature is at the core) c [ 1 0 . . . 0 ]¯ xN M θmin
(24)
The available control signals are the current into each coil, the frequency used in each coil and the rod velocity. Figure 7, top plot, shows an initial guess for the currents and frequencies in 10 coils at a fixed rod velocity of 13 cm/s. The resulting temperatures within the rod as a function of distance into the furnace are shown in the bottom plot. The normalized energy corresponding to this design is 1.9643. Figure 8 shows the optimized solution obtained using currents and frequencies as decision variables and fixing the rod velocity at 13 cm/s. Note that the minimum temperature in the rod, on exit from the furnace, is above the desired minimum value (θ min =1000 C) for both the initial and optimized designs. The normalized energy utilization for the optimized solution is 1.7479. This represents an improvement of 11% over the initial guess. It is worth noting that control input settings currently utilized in industry achieve energy utilization values within 3% of this optimal cost for a furnace of similar characteristics. Initial Design H (A/cm), f0 (Hz)
3000 2500 2000 1500
500 0
Temperature (C)
H (A/cm) Frequency (Hz)
1000
0
200
400
600
800
1000
1400
1600
1200 1000 800 Core half R Surface
600 400 200 0
200
400
600
800
1000
Axial Distance (cm)
Fig. 7.
1200
Temperature profile for an initial control sequence
1200
1400
1600
Optimal Design H (A/cm), f0 (Hz)
3000 2500 2000 1500
500 0
Temperature (C)
H (A/cm) Frequency (Hz)
1000
0
200
400
600
800
1000
1200
1400
1600
1200 1000 800 Core half R Surface
600 400 200 0
200
400
600
800
1000
1200
1400
1600
Axial Distance (cm)
Fig. 8.
Temperature profiles for optimal current and frequency distribution in multiple coils
VII. N ONLINEAR O BSERVER The model described above depends on several simplifying assumptions. Hence, it cannot exactly describe the system behaviour. (Recall the comparison made in Figure 6.) Consequently it is desirable that the model be “calibrated” on-line by comparing the model response with real data. This will be done via a nonlinear observer driven by the available measurements. The most common measurement is that of surface temperature (typically measured by one or multiple pyrometers). As can be seen from Figures 7 and 8, this may give some indication of the internal states of the system. However, it may be useful to seek additional measurements. We believe that an important measurement is that of energy supplied to the system. This can be readily measured in the external circuitry. It can also be found as an output of the model. (It is the sum of the energy going into heating the rod plus energy radiated from the surface). Both these quantities are functions of the model states. Our proposal is to use optimization to estimate the state with a cost function of the form Jo =
N
(yk − yˆk )r (yk − yˆk ) + vkT Qvk
(25)
k−1
where yk is the sequence of available measurements along the furnace, Q ≥ 0 is a weighting matrix and v k is the driving sequence into a model of the form:
x ˆjk+1 = f x ˆjk , uj + E vkj (26) for some appropriate matrix E . The basic idea of the observer, is to adjust {v kj } to minimize Jo . However, depending on the choice of the matrix E , there may be a major complexity issue since the dimension of v kj may be of the order of a thousand! Hence, some “parametrization” of the process noise is necessary to regularize the state estimation problem. In preliminary studies, we have reduced the complexity of the process noise by restricting it to errors in the heat radiated from the surface. Our first test considers data from a hypothetical “real” process simulated by changing the coefficient in the heat radiation loss from 4 to 3.5 (see the “grey body” law after equation (9)). For this process we “measure” two different outputs, as discussed above: surface temperature and energy supplied into the system. Figure 9 shows the results for the observer that employs measurements of energy supplied into the system only. In this case, there is good matching between the estimated and “real” states up to when the rod exits the last coil (at 1200 cm); this is to be expected since no more more energy is supplied into the system and, hence no more measurements are available, along the last section of the process. Figure 10 shows the results for the observer that employs surface temperature measurements. The “real” process surface and core temperatures are plotted in black dotted lines whilst the corresponding temperatures achieved by
Observer using Energy Measurements 1400 Observer Real Process 1200
Temperature (C)
1000
800
600
400
200
0
0
200
400
600
800
1000
1200
1400
1600
Axial Distance (cm)
Fig. 9.
Surface and core temperature profiles for the simulated “real” process and the observer employing energy measurements
the observer are shown by red solid lines. We see from this figure that there is a very good matching between the estimated and “real” states. Observer using Surface Temperature Measurements 1400 Real Process Observer 1200
Temperature (C)
1000
800
600
400
200
0
0
200
400
600
800
1000
1200
1400
1600
Axial Distance (cm)
Fig. 10. Surface and core temperature profiles for the simulated “real” process and the observer employing surface temperature measurements
Our next test uses real surface temperature data obtained from a process of similar characteristics in a local industry as shown earlier in Figure 6. The results obtained when the model states are adjusted using the observer are shown in Figure 11. We can see from this figure that the observer interpolates the real data with great accuracy, whereas the uncorrected model fails to interpolate most of the measurement points. This experiment highlights the importance of using the available measurements to correct the model and compensate for uncertainties, varying parameters, disturbances, etc. In practice, we envisage combining measurements of both energy and surface temperature in the observer designs. VIII. C HALLENGES There are many challenges with the control and state estimation problem described above, including:
Observer using Real Measurements 1200 Real Data Observer Uncorrected Model 1000
Temperature (C)
800
600
400
200
0
0
200
400
600
800
1000
1200
1400
1600
Axial Distance (cm)
Fig. 11. Real surface temperature data (black stars) and resulting surface and core temperature profiles for the observer employing these measurements (red solid curves) and the uncorrected model (Eq. (26) with vk = 0, blue dotted curves)
(i) The optimal solution is prone to local minimum effects. (The only remedy we currently have for this problem is to begin with many different initial guesses for the control.) (ii) The state dimension is very high. (We have used 1000 states.) (iii) The available measurements are scarce in number. (iv) The model equations are very stiff, i.e. there is a wide range in “time constraints”. Hence numerical issues are prevalent. (v) A key step is the development of a control relevant model. This can be viewed as a special form of model order reduction for a nonlinear system. Currently this step requires a trial and error procedure guided by physical intuition. (vi) Full scale implementation of the idea requires attention to a large number of practical details which fall outside the scope of the current paper. This aspect is still under development. IX. C ONCLUSIONS This paper has described initial work (spanning a 2 year period) to apply nonlinear model predictive control to a complex industrial process. The system is challenging for many reasons including: • high state dimension (thousands) • wide range in time constraints (100:1) • state dependent parameter variations (permeability changes by 200:1 when Curie temperature is passed) • many constraints • local minima problems in the optimization • few available measurements Not withstanding these challenges, the incentives to achieve better control are extremely high since the annual energy costs for just one furnace are of the order of a million dollars. On going work includes improvements to the numerical aspects of the problem, further work on the control strategy and full scale implementation. More generally, we believe that the induction heating furnace problem is a quintessential example of a complex system where the potential rewards achievable from advanced control warrant the effort needed to carry out the design and implement the final solution. We also believe that some of the challenges associated with this problem may be relevant to other problems. We thus hope some of the discussion presented here may be of broader interest in the systems and control community.
R EFERENCES [1] F. Bay, V. Labbe, Y. Favennec, and J. Chenot, “A numerical model for induction heating processes coupling electromagnetism and thermomechanics,” International Journal for Numerical Methods in Engineering, vol. 58, pp. 839–867, 2003. [2] T. Skoczkowski and M. Kalus, “The mathematical model of induction heating of ferromagnetic pipes,” IEEE Transactions on Magnetics, vol. 25, no. 3, pp. 2745–2750, 1989. [3] R. Pascal, P. Conraux, and J.-M. Bergheau, “Coupling between finite elements and boundary elements for the numerical simulation of induction heating processes using a harmonic balance method,” IEEE Transactions on Magnetics, vol. 39, no. 3, pp. 1535–1538, 2003. [4] X. Yang, Y. Wang, and W. Yan, “Simulation of induction heating device with double inductors for continuously heating up steel bars,” in World Automation Congress. WAC0’8., Hawaii, Sept. 28-Oct. 2 2008. [5] Z. Qianzhe, L. Yibing, L. Yanping, and Z. Weisong, “Numerical simulation of pc steel bar induction heating process,” Advanced Materials Research, vol. 228-229, pp. 270–275, 2011. [6] C. Chaboudez, S. Clain, R. Glardon, D. Mari, J. Rappaz, and M. Swierkosz, “Numerical modeling in induction heating for axisymmetric geometries,” IEEE Transactions on Magnetics, vol. 33, no. 1, pp. 739–745, 1997. [7] M. Fisk, “Simulation of induction heating in manufacturing,” Licentiate Thesis, Department of Applied Physics and Mechanical Engineering, Lule˚a University of Technology, 2008. [8] E. ter Maten and J. Melissen, “Simulation of inductive heating,” IEEE Transactions on Magnetics, vol. 28, no. 2, pp. 1287–1290, 1992. [9] L. Egan and E. Furlani, “A computer simulation of an induction heating system,” IEEE Transactions on Magnetics, vol. 27, no. 5, pp. 4343–4354, 1991. [10] S.-M. Jang, S. Cho, S.-H. Lee, H. Cho, and H. Park, “Thermal analysis of induction heating roll with heat pipes,” IEEE Transactions on Magnetics, vol. 39, no. 5, pp. 3244–3246, 2003. [11] V. Rudnev, “Intricacies of computer simulation of induction heating processes,” Science and Technology, Inductoheat Inc., http://www.inductoheat.com/pdf/Intricacies-of-Computer-Simulation-of-Induction-HeatingProcesses Valery Rudnev Chicago 2011.pdf, Tech. Rep., 2011. [12] C. Chaboudez, S. Clain, R. Glardon, J. Rappaz, M. Swierkosz, and R. Touzani, “Numerical modelling of induction heating of long workpieces,” IEEE Transactions on Magnetics, vol. 30, no. 6, pp. 5028–5037, 1994. ¨ [13] H. Unver and M. Aydemir, “Power and frequency control in a 60kw induction steel heating furnaces through plc,” Department of Electrical and Electronics Engineering, Kirikkale University, Turkey, http://www.kku.edu.tr/ unver/ISHF.pdf, Tech. Rep., 2005. [14] V. Rudnex, D. Loveless, R. Cook, and M. Black, Handbook of Induction Heating. Marcel Dekker AG, Switzerland, 2003. [15] S. Zinn and S. Semiatin, Elements of Induction Heating: Design, Control and Applications. Electric Power Research Institute, Inc., 1987. [16] W. Hayt and J. Buck, Engineering Electromagnetics. McGraw-Hill, 2006. [17] C. Kittel, Introduction to Solid State Physics. Wiley and Sons, New York, 1986. [18] M. Bailyn, A Survey of Thermodynamics. American Institute of Physics, New York, 1994. [19] Y. Cengel, Heat and Mass Transfer: A Practical Approach. McGraw-Hill Series in Mechanical Engineering, 2007. [20] J. Rawlings and D. Mayne, Model Predictive Control: Theory and Design. Nob Hill Publishing, 2009.
A BSTRACT Induction Heating Furnaces are used extensively in industry. The basic principle is that induced eddy currents are used to heat a ferromagnetic material as it passes through a series of coils. Because of the importance of such systems, there has been on-going interest in their design and operation. Past work includes model development from physical principles and optimal design of operational practices. However, previous work has invariably been based on open-loop strategies. Our work is aimed at the design of a closed-loop control strategy incorporating feedback from the available measurements. This paper reports initial work including model development and calibration together with preliminary control system design. Proposed future work includes full scale industrial implementation. I. I NTRODUCTION Induction Heating Furnaces are frequently used in industry. These systems consume extremely large amounts of energy. Indeed, the annual cost of electrical energy for a typical furnace can exceed US$1M. Hence, there is substantial financial motivation to reduce energy consumption. Beyond the cost of energy, there also exist strong environmental drivers since electrical energy is a major contributor to green house gas emissions and related issues. Simulation and control of induction heating furnaces are topics which have been extensively studied. There are presently two main ways in which these systems are designed. The first is where a computationally expensive simulation is run over a long period of time in order to find a theoretically optimal physical setup ([1], [2], [3], [4], [5], [6], [7], [8]). The models typically involve finite element analysis with a very large number of elements ([9], [10], [6], [11], [12]). The major shortcoming of these methods is that there is no flexibility in the design. Thus minor miscalculations can result in large errors appearing when physically implemented. The other main method used is closed loop temperature control. In this, the surface temperature of the bar is measured via pyrometers and any deviations from desired values are fed into a controller ([13], [14], [15]). While this offers some compensation for model error and disturbances, current implementations are of an empirical nature. Also, as we discuss later, surface temperature may be a poor guide to internal states of the process. The models in current use are unsuitable for on-line control studies since they are either excessively complex (i.e. based on a full partial differential equation description) or too simplistic (i.e. aimed at standard operational scenarios and/or design of PID controllers). Our proposal is to use nonlinear model predictive control based on a simplified model with feedback. This is done by creating a low complexity, “control relevant” model and using its predictions as part of a closed-loop control strategy. To date, the following steps have been carried out: (i) model development (ii) model calibration using industrial data (iii) simulated control studies using nonlinear MPC (iv) simulated observer design using nonlinear optimization (v) comparison of initial proposed control solutions with current industrial practice This paper is based on the plenary paper ”Opportunities and Challenges in the Application of Nonlinear MPC to Industrial Problems” presented at NMPC 2012. G.C. Goodwin, R.H. Middleton and M.M. Seron are with School of Electrical Engineering and Computer Science, University of Newcastle, 2308, Australia Graham.Goodwin, Richard.Middleton, [email protected] B. Campos is with Donhad Pty Ltd, Newcastle. [email protected]
Future work is aimed at full scale industrial implementation of these ideas. This phase is planned over a time horizon of several years since substantial practical issues are involved. This problem has a range of associated challenges including: (i) Whilst one can reasonably assume radial symmetry of temperatures within the rod, the temperatures are not constant in the longitudinal direction. This complicates model development and can lead to excessive model complexity. For example, if one were to divide the radial and longitudinal spatial dimension into 1000 segments each, then one would have a model containing a million state variables! (ii) Due to skin effect phenomena ([16]), much of the heating takes place in a narrow surface region. This region can be of the order of 1/2 mm or less in thickness and is associated with rapid thermal heating. By way of contrast, the core heats slowly. Hence the model is inherently “stiff” with a wide range of associated time constants. Moreover, the skin effect region progressively moves into the core as the outer temperature passes the Curie point. (iii) The available measurements are very limited being restricted, in essence, to surface temperature measurements (made by pyrometers) and the external currents and voltages. It may also be possible to use the forces measured in subsequent forming processes. (iv) The system has many degrees-of-freedom in the input including line speed, voltages, currents and frequencies in the induction coils (of which there may be 10 or 20 in a given set-up). (v) The system needs to satisfy a large number of constraints. These include, but are not restricted to: 1) A lower bound on minimum exit temperature. 2) An upper bound on maximum exit temperature. 3) Constraints on temperature gradients in both time and space to avoid metallurgical stress on the material. 4) Bounds on operating voltages, currents, frequencies and switching times dictated by the driving power electronics. (vi) All of the parameters appearing in the model are highly temperature dependent. Moreover, the temperatures within the rod range from ambient to above the Curie point ([17]). Some parameters, e.g. relative permeability, change by several orders of magnitude as the rod traverses a given coil. Hence the system is highly nonlinear. The layout of the remainder of the paper is as follows: In Section II, we give a brief overview of the physical apparatus. In Section III, we describe a control relevant model. Section IV discuss model parameter variations. Section V describes model validation against industrial data. Section VI describes the proposed control law and gives simulated performance results. Section VII describes the proposed observer design. Section VIII describes challenges associated with the problem. Conclusions are presented in Section IX. II. T HE P HYSICAL A PPARATUS A schematic diagram of a typical furnace is given in Figure 1. Velocity v u1 Fig. 1.
u2
u3
uNc
Schematic Diagram of Induction Heating Furnace {ui ≡ (Voltage, Current, Frequency); i = 1, . . . , Nc }
As shown in Figure 1, the key features of the apparatus are that a cylindrical rod is passed through a succession of coils. These coils are capable of being operated at different voltages, currents and frequencies. Moreover, the rod can be moved at different speeds through the furnace. The coils induce currents in the rod which, due to hysteresis effects, heat the rod. Most of the current in the rod is initially confined to the surface due to skin effect. Once the temperature passes the Curie point, then the relative permeability drops from about 200 to near 1 (i.e. that of air). The skin depth then increases. This, together with radial heat conduction, continues until the rod temperatures are all above the Curie temperature. Heating still continues but the skin depth will be now an order of magnitude larger; being 5 to 10 mms if the frequency is unchanged. For this reason, and also to improve the efficiency of the process with reduced permeability, it is usual to increase the frequency in this region to aid heating. The main goal in operating the furnace is to satisfy the constraints described earlier in Section I whilst minimizing energy consumption.
III. A C ONTROL R ELEVANT M ODEL The system is described by a set of partial differential equations which, inter alia, include Fourier’s Law (for heat conduction), Newton’s Law (for Thermal heating) and Maxwell’s equations (for electromagnetic phenomena). Many authors have begun modelling at this level. However, this inevitably leads to high order models with many thousands of state variables (see observation (i) in Section I). Our approach is to develop a control-relevant model. After much initial work with a larger model, we were led to the following: (i) We assume the electromagnetic aspects are in steady state and are thus described by algebraic relationships. (ii) We assume that the majority of heat conduction occurs in the radial direction. (iii) We focus on a single longitudinal zone and track it as it passes through the furnace. (iv) We assume that downstream zones have a minimal impact on electromagnetic effects in the current zone. (This is justified, in part, by the skin effect phenomena.) This approach leads to a nonlinear-state-space model in algebraic-difference form. Details of the model are outlined in the sequel. We note that all physical parameters are highly temperature dependent. This aspect will be discussed further in Section V. Our model focuses on one coil. However, the model applies mutatis mutandis to multiple coils and the region between coils (where the external coil current is zero). We consider a slice through the rod as shown in Figure 2. We assume this slice is in a longitudinal “zone” located between (j − 1) ∗ d and j ∗ d from the entry to the coil where d is the width of each slice. radiated heat loss
θN j
θ2j θ1j J1j r1
r2
rN
r
J2j JN j
Fig. 2.
Radial distribution of temperature and current density in the jth zone through the bar
We assume Nc coils. The concentric rings are assumed to have temperature, radius, and current density θ ij , ri , Jij , ; i = 1, . . . , N, j = 1, . . . , M respectively, where r N is the radius of the rod and M d is the length of the coil. A. Heat Transfer Model We denote the temperature in the i th ring in zone j at time step k as θ ij (k). There are three mechanisms for heat transfer: (i) heat conductivity due to temperature gradients (ii) heat generated due to circulating current (iii) radiation from the surface
The heat, Q1 , due to temperature gradients satisfies Fourier’s Law ([18], page 23), i.e., dθ dQ1 = −KA dt dx
where
(1)
dQ1 dt
is heat flux per unit time is heat conductivity (which is a function of temperature) A is the area Since the area is a function of radius, we have K
dθ dQ1 = −2Kπd r dt dr
(2)
where d is the width of the zone. Dividing (2) by r and integrating we obtain the heat flux into a region per unit time dQ1ij dt
=−
2Kπd (θij − θi−1,j ) n (ri /ri−1 )
(3)
Next we note that the heat generated due to circulating current is ρ|J| 2 (W/cm3 ). To use this, we need the volume of the ith ring in zone j , ρij ρij 2πri = Aij (ri − ri−1 )d 2 2 Vij = π(ri − ri−1 )d
Rij =
(4) (5)
Let ρij = ρ(θij ) be the resistivity in the ith ring in zone j . Hence heat per unit time due to circulating current, Q 2ij , satisfies dQ2ij = ρij |Jij |2 Vij dt Finally the temperature rise satisfies Newton’s Law ([18]), namely (for all but the outer ring): dQ2ij dQ1ij dθij 1 = + dt c(θij )δVij dt dt
(6)
(7)
where c is the temperature dependent specific heat capacity [nominally 0.45 J/(gK) for iron], δ is the metal density [7g/cm3 for iron]. If we use a time increment of Δ seconds, then (7) can be approximated by θij (k + 1) = θij (k) dQ2ij dQ1ij 1 + +Δ c(θij (k))δVij dt dt
(8)
θij (k + 1) = θij (k) d +Δ c(θij (k))δVij θi+1,j (k) − θij (k) · K(θij (k))2π d n (ri+1 /ri ) θij (k) − θi−1,j (k) d − K(θij (k))2π n (ri /ri−1 ) + ρij |Jij |2 Vij
(9)
Combining (3), (7), (8) we have
Equation (9) holds for i = 2, . . . , N − 1 ; j = 1, . . . , M .
For i = 1, we have no heat conducted out. For i = N , we replace the heat conducted into the cell by heat 4 − θ4 radiation loss. We model heat radiation by; a “grey body” law: σ2πr N d θN j amb , where σ is the Stefan Boltzman constant ([19]) and θ amb is the ambient temperature. We note that the above equations, in particular (6), depend upon the current density J ij (k). Equations for the latter variables are developed in the sequel. B. Electromagnetic Modelling We model the electromagnetics within the bar via the complex sinusoidal steady state version of Maxwell’s equations, in cylindrical coordinates. These equations are: d E = jωo μH (10) dr d H = −J (11) dr where j is the imaginary unit, E is the circular component of the electric field, H is the axial component of the field intensity, ωo is the excitation angular frequency, μ = μ o μr is the magnetic permeability, and J is the current density given by 1 E (12) ρ Note that E , H and J in (10) - (12) are complex valued functions of radius, r . To solve (10), (11), (12) we do a spatial discretization of E and H , as vectors representing the fields in each zone. For a fixed zone, this is J=
¯ [E(0), E(r1 ), . . . , E(rN )]T E ¯ [H(0), H(r1 ), . . . , H(rN )]T H d The spatial derivative operator, dr , can then be represented by a matrix multiplication with ⎡ ⎤ 1 1 − r(1) 0 · · · · · · 0 r(1) ⎢ 0 ⎥ 1 1 0 ... 0 − r(2)−r(1) ⎢ ⎥ r(2)−r(1) ⎢ ⎥ Γ⎢ .. .. .. .. .. .. ⎥ . . . . . . ⎣ ⎦ 1 1 0 ... ... . . . − r(N )−r(N −1) r(N )−r(N −1)
Then (10), (11), (12) can be approximated by the implicit, sparse, linear system of equations: ¯ ¯ 0 diag{jω0 pi } Γ 0 E E ¯ = diag{− 1 } ¯ 0 0 Γ H H pi
(13) (14)
(15)
(16)
together with equality constraints E(0) = 0 and H(r N ) = Nc Ic /L (the field intensity due to the Nc turn coil current, Ic , over a length, L). Solution of (16) in each zone allows computation of the current densities for (12). C. Modeling the Motion of the Bar The final step in the model is to describe the motion of the bar. We do this by noting that, if we choose the width of the rings as d = vΔ where v is the bar velocity and Δ is the time step, then θi,j+1 (k + 1) = θi,j (k) for i = 1, . . . , N ; j = 1, . . . , M − 1
(17)
The boundary conditions for this equation are ambient (for the first coil) or the exit temperatures from the previous coil. The above strategy allows us to combine the subscripts j and k, i.e. we use θ i (k) to denote the temperature in ring i after it has passed k steps into the coil.
D. Final Model Format We note that the model as described above has two key features: (i) Progressively heat is generated via induction as the rod transits the furnace. (ii) There is a mixture of difference equations (for temperature) and algebraic equations (for electromagnetic variables). j So as to express the model in a familiar notation, we use x jk to denote the vector of temperatures θ 1j (k) . . . θN (k) in the j th zone. Also, let ejk denote the corresponding vector of electromagnetic variables. Then the model described above takes the following form:
(18) xjk+1 = f 1 xjk , ejk where ejk is a nonlinear function of xjk . Next we note that the electromagnetic model is linear in e jk but highly nonlinear in the temperature. Thus these equations can be expressed in the form of a linear matrix equation: Dk ejk = bjk
(19)
where Dk depends nonlinearly on temperature i.e. x jk . The right hand side of (19), i.e. b jk , depends linearly on the external current. For future use, we denote the “control” variables i.e. the external current and frequency by u j . Then, for each fixed temperature xjk , equation (19) can be inverted to express e jk as a function of uj and xjk , i.e. ejk = Dk−1 bjk
(20)
Finally, substituting (20) into (18) we obtain a model having the familiar form:
xjk+1 = f xjk , uj ; j = 1, . . . , Nc
(21)
with boundary conditions: x11 = ambient temperature j+1 x1 = xjM ; j = 1, . . . , Nc − 1 IV. M ODEL PARAMETER VARIATIONS One of the key challenges in this problem is the fact that all physical parameters are highly temperature dependent. As an illustration, typical temperature dependencies for heat conductivity, resistivity and relative permeability are shown in Figure 3, Figure 4 and Figure 5, respectively. Note that the dependence of permeability on temperature is particularly significant. k W/(mK) 50
25
800 100 Fig. 3.
temp o C
Heat conductivity as a function of temperature
In particular, the relative permeability changes by 200:1 over a very narrow temperature region as the material approaches and passes the Curie temperature. This phenomena is a major factor in the control problem.
ρ
1.2*10−4
0.2*10−4 800 100 Fig. 4.
temp o C
Resistivity as a function of temperature
μr
200
1 800 100 Fig. 5.
temp o C
Relative permeability as a function of temperature
V. M ODEL VALIDATION AGAINST E XPERIMENTAL DATA The model parameters were chosen based on published empirical data. The model was then run using typical values used in industry for the control variables. Figure 6 compares the predicted surface temperature and real surface temperature data obtained from a process of similar characteristics in a local industry. We see from Figure 6 that the model performs qualitatively the same as the real process. However, some important discrepancies are evident. These discrepancies will be addressed in Section VII. VI. P ROPOSED O N - LINE C ONTROL L AW The basic principle behind the development of the control law is to use finite horizon optimization in a nonlinear MPC control law. However, there are several departures for the usual MPC type controllers described in the literature (see for example [20]). In particular, once a control sequence u 1 , . . . , uNc has been determined, then it can all be applied since the “time index” actually relates to the spatial dimension of distance through the furnace. Details are given below. We use a prediction model of the form:
x ¯jk+1 = f x ¯jk , uj ; j − 1, . . . , Nc (22) We use a cost function which captures energy utilization in the N c coils. Nc J= uj j=1
(23)
Real Data and Uncorrected Model 1200 Real Data Uncorrected Model 1000
Temperature (C)
800
600
400
200
0
0
200
400
600
800
1000
1200
1400
1600
Axial Distance (cm)
Fig. 6.
Real surface temperature data (black stars) and output of the uncorrected model (solid line).
We optimize (23) subject to (22) and the constraints listed in Section I, (v). For example, the constraint on the minimum exit temperature would be (assuming the minimum temperature is at the core) c [ 1 0 . . . 0 ]¯ xN M θmin
(24)
The available control signals are the current into each coil, the frequency used in each coil and the rod velocity. Figure 7, top plot, shows an initial guess for the currents and frequencies in 10 coils at a fixed rod velocity of 13 cm/s. The resulting temperatures within the rod as a function of distance into the furnace are shown in the bottom plot. The normalized energy corresponding to this design is 1.9643. Figure 8 shows the optimized solution obtained using currents and frequencies as decision variables and fixing the rod velocity at 13 cm/s. Note that the minimum temperature in the rod, on exit from the furnace, is above the desired minimum value (θ min =1000 C) for both the initial and optimized designs. The normalized energy utilization for the optimized solution is 1.7479. This represents an improvement of 11% over the initial guess. It is worth noting that control input settings currently utilized in industry achieve energy utilization values within 3% of this optimal cost for a furnace of similar characteristics. Initial Design H (A/cm), f0 (Hz)
3000 2500 2000 1500
500 0
Temperature (C)
H (A/cm) Frequency (Hz)
1000
0
200
400
600
800
1000
1400
1600
1200 1000 800 Core half R Surface
600 400 200 0
200
400
600
800
1000
Axial Distance (cm)
Fig. 7.
1200
Temperature profile for an initial control sequence
1200
1400
1600
Optimal Design H (A/cm), f0 (Hz)
3000 2500 2000 1500
500 0
Temperature (C)
H (A/cm) Frequency (Hz)
1000
0
200
400
600
800
1000
1200
1400
1600
1200 1000 800 Core half R Surface
600 400 200 0
200
400
600
800
1000
1200
1400
1600
Axial Distance (cm)
Fig. 8.
Temperature profiles for optimal current and frequency distribution in multiple coils
VII. N ONLINEAR O BSERVER The model described above depends on several simplifying assumptions. Hence, it cannot exactly describe the system behaviour. (Recall the comparison made in Figure 6.) Consequently it is desirable that the model be “calibrated” on-line by comparing the model response with real data. This will be done via a nonlinear observer driven by the available measurements. The most common measurement is that of surface temperature (typically measured by one or multiple pyrometers). As can be seen from Figures 7 and 8, this may give some indication of the internal states of the system. However, it may be useful to seek additional measurements. We believe that an important measurement is that of energy supplied to the system. This can be readily measured in the external circuitry. It can also be found as an output of the model. (It is the sum of the energy going into heating the rod plus energy radiated from the surface). Both these quantities are functions of the model states. Our proposal is to use optimization to estimate the state with a cost function of the form Jo =
N
(yk − yˆk )r (yk − yˆk ) + vkT Qvk
(25)
k−1
where yk is the sequence of available measurements along the furnace, Q ≥ 0 is a weighting matrix and v k is the driving sequence into a model of the form:
x ˆjk+1 = f x ˆjk , uj + E vkj (26) for some appropriate matrix E . The basic idea of the observer, is to adjust {v kj } to minimize Jo . However, depending on the choice of the matrix E , there may be a major complexity issue since the dimension of v kj may be of the order of a thousand! Hence, some “parametrization” of the process noise is necessary to regularize the state estimation problem. In preliminary studies, we have reduced the complexity of the process noise by restricting it to errors in the heat radiated from the surface. Our first test considers data from a hypothetical “real” process simulated by changing the coefficient in the heat radiation loss from 4 to 3.5 (see the “grey body” law after equation (9)). For this process we “measure” two different outputs, as discussed above: surface temperature and energy supplied into the system. Figure 9 shows the results for the observer that employs measurements of energy supplied into the system only. In this case, there is good matching between the estimated and “real” states up to when the rod exits the last coil (at 1200 cm); this is to be expected since no more more energy is supplied into the system and, hence no more measurements are available, along the last section of the process. Figure 10 shows the results for the observer that employs surface temperature measurements. The “real” process surface and core temperatures are plotted in black dotted lines whilst the corresponding temperatures achieved by
Observer using Energy Measurements 1400 Observer Real Process 1200
Temperature (C)
1000
800
600
400
200
0
0
200
400
600
800
1000
1200
1400
1600
Axial Distance (cm)
Fig. 9.
Surface and core temperature profiles for the simulated “real” process and the observer employing energy measurements
the observer are shown by red solid lines. We see from this figure that there is a very good matching between the estimated and “real” states. Observer using Surface Temperature Measurements 1400 Real Process Observer 1200
Temperature (C)
1000
800
600
400
200
0
0
200
400
600
800
1000
1200
1400
1600
Axial Distance (cm)
Fig. 10. Surface and core temperature profiles for the simulated “real” process and the observer employing surface temperature measurements
Our next test uses real surface temperature data obtained from a process of similar characteristics in a local industry as shown earlier in Figure 6. The results obtained when the model states are adjusted using the observer are shown in Figure 11. We can see from this figure that the observer interpolates the real data with great accuracy, whereas the uncorrected model fails to interpolate most of the measurement points. This experiment highlights the importance of using the available measurements to correct the model and compensate for uncertainties, varying parameters, disturbances, etc. In practice, we envisage combining measurements of both energy and surface temperature in the observer designs. VIII. C HALLENGES There are many challenges with the control and state estimation problem described above, including:
Observer using Real Measurements 1200 Real Data Observer Uncorrected Model 1000
Temperature (C)
800
600
400
200
0
0
200
400
600
800
1000
1200
1400
1600
Axial Distance (cm)
Fig. 11. Real surface temperature data (black stars) and resulting surface and core temperature profiles for the observer employing these measurements (red solid curves) and the uncorrected model (Eq. (26) with vk = 0, blue dotted curves)
(i) The optimal solution is prone to local minimum effects. (The only remedy we currently have for this problem is to begin with many different initial guesses for the control.) (ii) The state dimension is very high. (We have used 1000 states.) (iii) The available measurements are scarce in number. (iv) The model equations are very stiff, i.e. there is a wide range in “time constraints”. Hence numerical issues are prevalent. (v) A key step is the development of a control relevant model. This can be viewed as a special form of model order reduction for a nonlinear system. Currently this step requires a trial and error procedure guided by physical intuition. (vi) Full scale implementation of the idea requires attention to a large number of practical details which fall outside the scope of the current paper. This aspect is still under development. IX. C ONCLUSIONS This paper has described initial work (spanning a 2 year period) to apply nonlinear model predictive control to a complex industrial process. The system is challenging for many reasons including: • high state dimension (thousands) • wide range in time constraints (100:1) • state dependent parameter variations (permeability changes by 200:1 when Curie temperature is passed) • many constraints • local minima problems in the optimization • few available measurements Not withstanding these challenges, the incentives to achieve better control are extremely high since the annual energy costs for just one furnace are of the order of a million dollars. On going work includes improvements to the numerical aspects of the problem, further work on the control strategy and full scale implementation. More generally, we believe that the induction heating furnace problem is a quintessential example of a complex system where the potential rewards achievable from advanced control warrant the effort needed to carry out the design and implement the final solution. We also believe that some of the challenges associated with this problem may be relevant to other problems. We thus hope some of the discussion presented here may be of broader interest in the systems and control community.
R EFERENCES [1] F. Bay, V. Labbe, Y. Favennec, and J. Chenot, “A numerical model for induction heating processes coupling electromagnetism and thermomechanics,” International Journal for Numerical Methods in Engineering, vol. 58, pp. 839–867, 2003. [2] T. Skoczkowski and M. Kalus, “The mathematical model of induction heating of ferromagnetic pipes,” IEEE Transactions on Magnetics, vol. 25, no. 3, pp. 2745–2750, 1989. [3] R. Pascal, P. Conraux, and J.-M. Bergheau, “Coupling between finite elements and boundary elements for the numerical simulation of induction heating processes using a harmonic balance method,” IEEE Transactions on Magnetics, vol. 39, no. 3, pp. 1535–1538, 2003. [4] X. Yang, Y. Wang, and W. Yan, “Simulation of induction heating device with double inductors for continuously heating up steel bars,” in World Automation Congress. WAC0’8., Hawaii, Sept. 28-Oct. 2 2008. [5] Z. Qianzhe, L. Yibing, L. Yanping, and Z. Weisong, “Numerical simulation of pc steel bar induction heating process,” Advanced Materials Research, vol. 228-229, pp. 270–275, 2011. [6] C. Chaboudez, S. Clain, R. Glardon, D. Mari, J. Rappaz, and M. Swierkosz, “Numerical modeling in induction heating for axisymmetric geometries,” IEEE Transactions on Magnetics, vol. 33, no. 1, pp. 739–745, 1997. [7] M. Fisk, “Simulation of induction heating in manufacturing,” Licentiate Thesis, Department of Applied Physics and Mechanical Engineering, Lule˚a University of Technology, 2008. [8] E. ter Maten and J. Melissen, “Simulation of inductive heating,” IEEE Transactions on Magnetics, vol. 28, no. 2, pp. 1287–1290, 1992. [9] L. Egan and E. Furlani, “A computer simulation of an induction heating system,” IEEE Transactions on Magnetics, vol. 27, no. 5, pp. 4343–4354, 1991. [10] S.-M. Jang, S. Cho, S.-H. Lee, H. Cho, and H. Park, “Thermal analysis of induction heating roll with heat pipes,” IEEE Transactions on Magnetics, vol. 39, no. 5, pp. 3244–3246, 2003. [11] V. Rudnev, “Intricacies of computer simulation of induction heating processes,” Science and Technology, Inductoheat Inc., http://www.inductoheat.com/pdf/Intricacies-of-Computer-Simulation-of-Induction-HeatingProcesses Valery Rudnev Chicago 2011.pdf, Tech. Rep., 2011. [12] C. Chaboudez, S. Clain, R. Glardon, J. Rappaz, M. Swierkosz, and R. Touzani, “Numerical modelling of induction heating of long workpieces,” IEEE Transactions on Magnetics, vol. 30, no. 6, pp. 5028–5037, 1994. ¨ [13] H. Unver and M. Aydemir, “Power and frequency control in a 60kw induction steel heating furnaces through plc,” Department of Electrical and Electronics Engineering, Kirikkale University, Turkey, http://www.kku.edu.tr/ unver/ISHF.pdf, Tech. Rep., 2005. [14] V. Rudnex, D. Loveless, R. Cook, and M. Black, Handbook of Induction Heating. Marcel Dekker AG, Switzerland, 2003. [15] S. Zinn and S. Semiatin, Elements of Induction Heating: Design, Control and Applications. Electric Power Research Institute, Inc., 1987. [16] W. Hayt and J. Buck, Engineering Electromagnetics. McGraw-Hill, 2006. [17] C. Kittel, Introduction to Solid State Physics. Wiley and Sons, New York, 1986. [18] M. Bailyn, A Survey of Thermodynamics. American Institute of Physics, New York, 1994. [19] Y. Cengel, Heat and Mass Transfer: A Practical Approach. McGraw-Hill Series in Mechanical Engineering, 2007. [20] J. Rawlings and D. Mayne, Model Predictive Control: Theory and Design. Nob Hill Publishing, 2009.
