The Application of MPC in Microwave Heating Process Based on Model Constructed by Lambert’s Law Combined with Temperature

Microwave heating has been the research hot spot for many years. It has the capability of volumetric heating, which makes it save energy and time compared to traditional heating by conduction. A lot of work has been done to easily and exactly describe power distribution in the heating material. Maxwell’s equations and Lambert’s law are the most common ways. Maxwell’s equation is complicated and hard to apply, while Lambert’s law ignores the temperature influence. For material thickness less than penetration depth, only Maxwell’s equation can accurately solve power distribution. For large thickness material, Lambert’s law combined with regional temperature proposed in this paper can be more precise than only Lambert’s law. But there also exist some differences. To precisely control the heating process and make the whole process safe, this paper proposes the use of model predictive control (MPC) algorithm to make the maximum temperature follow a preset reference trajectory. The simulation results demonstrate that the algorithm can well control the heating process with little difference between the reference trajectory and the practical output.


Introduction
Microwave heating has been widely used in our lives, from domestic utilized microwave oven to industrial high-power applications.Compared with traditional heating method, microwave heating saves energy and time and is easier to control.But there exist two main problems holding back its wide applications.Firstly, during the heating process, due to the microwave's inhomogeneous heating, temperature distribution nonuniformity will make the heating results unsatisfied.Secondly, at some parts of the heating material, sharp temperature rising may happen.Then it may lead to thermal runaway, which is very dangerous in practical applications.It may cause reactants burning or even blow them up.
Maxwell's equation and Lambert's exponential law form the basis for modeling microwave assisted heating process.According to Lambert's law, microwave power decays exponentially by the penetration depth into the material.
The correctness to calculate power distribution for large thickness material has been proved.While Maxwell's equation is based on space and time, its calculation characterizes the exact microwave propagation behavior.But the calculation process is complex and needs exact temperature distribution.
Maxwell's equation is commonly used to theoretically explain some phenomenon or validate proposed models, for example, the simulation analysis of the heating process in microwave ovens [1,2], the relationship between the emerging of resonance and material thickness [3,4], the temperature field distribution during the thawing process [5], and the verification of impulse microwave on making heating process temperature distribution uniform [6]; those researches help people know microwave heating characteristic.But, in practical applications, the exact microwave distribution is hard to obtain.Microwave field distribution is influenced by material permittivity, which directly relates to temperature.The usual way to obtain temperature is by optical fiber thermometer, infrared thermometer, or infrared camera [7].Those methods are suitable for detecting surface or point temperature.Another way to detect temperature is by ultrasonic sensors, which can obtain the temperature field in the heating material.But there are some problems unsolved to be applied in practical applications [8].Because of the incomplete temperature field, the exact microwave distribution is basically unknown during the heating process.So, Maxwell's equation cannot be used to solve power distribution in real time practical applications.
Lambert's law is the common way to calculate power distribution in actual applications, such as food heating, thawing, and drying process [9][10][11][12][13].But it disregards the influence of temperature which directly relates to permittivity and makes the power distribution quite different with reality.This paper uses Lambert's law combined with material temperature to analyze power distribution.It makes up the shortcoming of ignoring material regional permittivity.Although Lambert's law combined with temperature can have a better solution to solve power distribution, there still exist some differences.So, to make the heating process safe and avoid thermal runaway, some control methods have been proposed by former researchers, for example, the global linearizing control of multiple-output-multiple-input microwave assisted thawing process [14] and the automatic control of microwave heating process [15].But because of the lacking of power distribution information, some problems remain unsolved to exactly control the heating process.
Model predictive control (MPC), an advanced model based control strategy, is well dedicated to solve this constrained problem of ordinary differential equation systems [16][17][18].MPC or receding horizon control is used to predict and optimize process performance.The main idea is to solve the manipulated variable input value over a finite prediction horizon at each sampling time.The procedure is reiterated at the next sampling time with the updated process measurements and model parameters.Today, MPC has become an advanced control strategy widely used in industry.
According to the above, to analyze microwave assisted heating process and guarantee its safety, MPC algorithm based on Lambert's law combined with temperature is proposed in this paper.This paper is organized as follows.In Section 2, it mainly deals with the difference among the three ways (Maxwell's equation, Lambert's law, and Lambert's law combined with local temperature) to calculate power distribution.Section 3 introduces the way how control algorithm can be applied in the proposed model to control the heating process.At last, the simulation results of this model and the effect of the control algorithm are revealed.

Model Building
The power distribution model applied in this analysis is shown in Figure 1.The heating material is placed in free space, with permittivity   and permeability   , which are varying with temperature. 0 is the free space dielectric constant, and  0 is the free space permeability.Assuming that a plane, time harmonic electromagnetic wave of frequency  impinges Layer N

Absorbing boundary
Absorbing boundary normally upon an isotropic material which fills the region 0 <  < , the transverse electric (TE) wave propagates along the -direction, with perpendicular electric (  ) along the direction and magnetic (  ) along the -direction.The wave satisfies the first order Mur absorbing boundary condition.The electromagnetic field satisfies the following formulas: To solve the electric field within the sample, the solutions of following coupled Maxwell's equations are required.The heating sample may be made by material with different permittivity.Let us set it to have  layers, so each layer can be represented by where  −1 ≤  ≤   and  = 1 ⋅ ⋅ ⋅ .Here  0 ,  1 ⋅ ⋅ ⋅   denote the interface position.  = √   0 (   () +    ()) 0 is the propagation constant which depends on dielectric constant    () and dielectric loss    ().  is the material permeability, which equals one for nonmagnetic heating material.Let us assume each layer has a constant dielectric property; therefore, the solutions of (2) can be represented by propagation and reflection wave as The electric and magnetic field are continuous at the interface between different layers.That is as follows: By using the interface conditions, the general solution coefficients can be obtained via solving the set of algebraic equations: Let us denote by  0 the incident power, which is the control variable in actual applications.Then the electric amplitude at  < 0 can be calculated by where  is the velocity of light.
Based on ( 5) and ( 6), the absorbed power in the th layer, calculated by Poynting vector theorem, is where the superscript " * " denotes the complex conjugate.

Case of Lambert's Law.
Lambert's law is valid for semiinfinite material, and the power absorbed by the material per unit volume can be represented by [19] where  ,1 represents the penetration coefficient calculated by the sampling temperature on material surface and  represents attenuation factor: 2.2.Case of Lambert's Law Combined with Regional Temperature.The main characteristic of Lambert's law combined with regional temperature is that it considers the regional permittivity and attenuation factor.The algorithm uses electric amplitude attenuation instead of power attenuation.The eclectic amplitude in each layer can be calculated by ( 10) and (11).The penetration depth of each layer is calculated by the corresponding permittivity, and the power dissipation in the material can be solved by (12).The use of this equation needs to sample material temperature.Space distribution of temperature sensors can be set as demanded: Lambert's law cannot exactly describe power distribution even for thickness much larger than penetration depth.For 2-centimeter-thick material, with the temperature distribution of Figure 3(a), the power distribution among the three methods has a colossal difference.The reflected wave between material and right surrounding material has relatively high amplitude compared with penetrated wave from input.So, Lambert's law and Lambert's law combined with temperature cannot be applied to exactly calculate power distribution.In the method of Lambert's law combined with temperature, the length of each layer is set to be 1 mm.
For 5-centimeter-thick material, with the temperature distribution of Figure 4(a), the power distribution among the three methods also has many differences.The exact power value solved by Maxwell's equation has a large fluctuation, which cannot be seen by Lambert's law.Lambert's law combined with temperature can approximately describe the fluctuation.Compared with Lambert's law, it has a more precise solution.In the method of Lambert's law combined with temperature, the length of each layer is also set to be 1 mm.
For 10-centimeter-thick material, with the temperature distribution of Figure 5(a), the exact power distribution solved by Maxwell's equation has a maximum value around 1.3 cm, which cannot be seen by Lambert's law.Lambert's law combined with temperature can approximately describe the inner maximum phenomena.Knowing maximum power can help the heating process, making the whole process safe.In the method of Lambert's law combined with temperature, the length of each layer is set to be 2.5 mm.

Temperature Field Calculation.
After power distribution is deduced, heat generation and transfer in the material can be calculated by the following equation: with the initial condition (,  = 0) =  0 (), where  is material density (Kg ⋅ m −3 );   is specific heat of the material (J ⋅ K −1 ⋅ Kg −1 );   is thermal conductivity of the material (W ⋅ m −1 ⋅ K −1 ); (, ) is material temperature (K); () is electromagnetic power dissipated per unit volume (W ⋅ m −3 ).
The temperature satisfies the surface heat balance; boundary condition can be set by where ℎ is the convection heat transfer parameter (W⋅m −2 ⋅K) and  0 is the temperature of the surround external material.
During the microwave heating process, the temperature rising can be represented by the following algorithms by discretizing ( 14)-( 15) as where   and   ( = 1, . . ., ) represent the temperature and power distribution at the th layer.The power distribution   at each layer is calculated by solving (1)- (7).To make the simulation approximate reality, the length of each layer is set to be 0.1 mm and the time step length is chosen to be 0.5 ms.
The simulation process can refer to paper [20].However, during the actual heating process, the sampling time step usually takes a few seconds and the number of temperature sensors is limited.In this analysis, the sampling time step is chosen to be one second; the temperature rising does not make a huge difference in   and ().The number of temperature sensors is set to be   , and those sensors distributed uniformly along the heating material.So, by using the assumption that, at each prediction horizon, heat and electromagnetic parameters keep the same, the following equations can be deduced: where   () is the sampling temperature at the th sampling point,   = /  ,  represents the length of the heating material,  represents the sampling moment,   represents the prediction time, and   is the power distribution at the th sampling point and is solved by using ( 11)-( 13).

Temperature Field Calculation by the Three Methods.
The mean square error is used to estimate the difference between Maxwell's equation and the other two ways.This calculation will clearly show the correctness between the two ways: where  Max  presents the temperature calculated by Maxwell's equation,    presents the temperature calculated by Lambert's law, and    presents the temperature calculated by Lambert's law combined with temperature.
The initial temperature of the heating material and environment is set to be 280 K.After heating time of 100 seconds with the power of 20 W ⋅ cm −2 , the temperature field in the material is shown in Figure 6.The difference between Lambert's law and Lambert's law combined with temperature can be referred to in Table 2. Lambert's law combined with temperature has a better ability to simulate microwave heating process than Lambert's law.The difference between Lambert's law combined with temperature and Maxwell's equation can be quite small.

Model Predictive Control
In actual applications, due to the constraint of temperature sensors, Maxwell's equation cannot be used.The most common way to solve power distribution is Lambert's law.By combining temperature sampled by sensors, Lambert's law combined with temperature can be used to improve calculation precision, but there still exists some difference.To make the whole heating process under control by using the inaccurately calculated power distribution, MPC has been proved to be suitable to such situation.The mean thought of MPC is to minimize the following equation, under some constraints, the cost function : where  is the actual discrete time index,  is the discrete time index, and   is the receding horizon. ref contains the description of reference behavior;   is the sequence of the future input.Consider the following: (20) where   is the control horizon and the other inputs are the same in the future step: For the problem of this paper, the input power has the constraint as follows: In this analysis, the control model can be observed in Figure 6.Maxwell's equation is used to simulate the actual

Optimization algorithm
Power distribution calculated by Lambert's law combined with temperature heating process.Lambert's law combined with temperature is used to calculate power distribution.The error between the actual output and the reference can be represented as At each sampling time , the error () is updated and remains the same over the prediction horizon   (Figure 7).
The heating process is nonuniform, which means the temperature is different at different location.To make the heating process safe, the maximum temperature must be under control.So the maximum temperature is selected as the control variable.The cost function  can take the following form: where ))] . ( (1) and (1) represent the first element in the columns   and .

Simulation Results
To prove the correctness of the proposed algorithm, the materials with thickness of 5 cm and 2 cm are chosen as the control objects.As expounded in Section 2, for the material with the thickness of 5 cm, Lambert's law combined with temperature can be used to solve approximate power distribution.For material thickness larger than 5 cm, more precise solution can be obtained.So the thickness of 5 cm is suitable to verify the effectiveness of the proposed algorithm under approximate power calculation.For the thickness of 2 cm, Lambert's law combined with temperature cannot solve power distribution precisely.Under such situation, the model has a large difference with reality.So, it is selected to ascertain the commonality of this algorithm.In this paper,   =   = 3 is chosen as the prediction and control horizon.The initial temperature of the material and environment is set to be 280 K.The convection heat transfer parameter is chosen to be 10 W ⋅ m −2 ⋅ K.The input power has the constraint of 0 ≤  ≤ 20 W ⋅ cm −2 .The distance between two adjacent sensors is selected to be   = 1 mm.For the material with the thickness of 5 cm, the matrix parameters are set as  = diag[5 × 10 9 , 5 × 10 9 , 5 × 10 9 ],  = diag[0.6,0.6, 0.6], and  = [0.1,0.1, 0.1]  .The reference trajectory is set as  ref () = 280 + 0.6 ×  − 0.001 ×  2 , for 0 ≤  ≤ 100. Figure 8(a) indicates the temperature field calculated by Maxwell's equation; Figure 8(b) shows temperature field calculated by Lambert's law combined with temperature.There is little temperature difference between the two algorithms.The maximum temperature calculated by Lambert's law combined with temperature is just a little higher than Maxwell's equation.So as Figure 8(c) shows, the maximum temperature can well follow the reference trajectory.Figure 8(d) shows the variation of input power, and an optimal input with a slight fluctuation can be obtained.The emerging of input fluctuation is caused by the variation of the maximum temperature location, which slowly moves from the material surface to the interior.
For the material thickness larger than penetration depth, the model based on Lambert's law combined with temperature can approximately solve power distribution, and a satisfied result can be obtained.But for small thickness material, the power distribution cannot be exactly solved by Lambert's law combined with temperature.MPC algorithm can still have a satisfied result.The matrix parameters are set as  = diag[10 9 , 10 9 , 10 9 ],  = diag[0.08,0.08, 0.08], and  = [−0.9,−0.9, −0.9]  .The reference trajectory is the same with 5-centimeter-thick material.There is a large difference between temperature distribution calculated by Maxwell's equation and Lambert's law combined with temperature, as shown in Figures 9(a) and 9(b).Temperature field calculated by Maxwell's equation has a considerable fluctuation.The maximum temperature is located inside the material.However, Lambert's law combined with temperature only has a gradually declining temperature field from the left side to the right side. Figure 9(c) shows the maximum temperature can approximately track the reference trajectory.But for the inaccurately calculated power distribution, the input power cannot always be optimal.Large fluctuations appear and the time to get steady is large, as shown in Figure 9(d).

Conclusion
This paper compares Lambert's law with Maxwell's equation.Lambert's law can have a good power approximation for material with thickness larger than penetration depth and constant permittivity.But for material in which temperature has a great influence on permittivity, Lambert's law is not suitable.The simulation results show Lambert's law combined with regional temperature is more precise than Lambert's law.
To exactly control the whole heating process, MPC algorithm is used to make up the difference between power distribution calculated by Lambert's law combined with temperature and reality.The simulation results show that this algorithm can have a good control effect under model inaccuracy.But for material thickness less than penetration depth, input cannot always be optimal and the time needed to get steady is quite large.So how to quickly get the optimal input under model inaccuracy is the future research emphasis.

Figure 8 :
Figure 8: Algorithm verification for the material with the thickness of 5 cm, (a) temperature field calculated by Maxwell's equation; (b) temperature field calculated by Lambert's law combined with temperature; (c) maximum temperature variation versus the reference trajectory; (d) input power at each sampling step.

Figure 9 :
Figure 9: Algorithm verification for the material with the thickness of 2 cm, (a) temperature field calculated by Maxwell's equation; (b) temperature field calculated by Lambert's law combined with temperature; (c) maximum temperature variation versus reference trajectory; (d) input power at each sampling step.

Table 2 :
The mean square error of different thickness.