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

Authors: Jianshuo Li; Qingyu Xiong; Kai Wang; Xin Shi; Shan Liang; Min Gao
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

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.


---

## Body

## 1. 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–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–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 Section2, 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.


## 2. Model Building

The power distribution model applied in this analysis is shown in Figure1. The heating material is placed in free space, with permittivity ε
n and permeability μ
n, 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 normally upon an isotropic material which fills the region 0
<
z
<
L, the transverse electric (TE) wave propagates along the z-direction, with perpendicular electric (E
x) along the x-direction and magnetic (H
y) along the y-direction. The wave satisfies the first order Mur absorbing boundary condition. The electromagnetic field satisfies the following formulas:(1)
∂
E
→
x
∂
z
=
μ
r
μ
0
∂
H
y
∂
t
∂
H
→
y
∂
z
=
ε
r
ε
0
∂
E
x
∂
t
.Figure 1
Wave propagation model.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 haveN layers, so each layer can be represented by(2)
d
2
E
x
,
0
d
z
2
+
k
0
2
E
x
,
0
=
0
z
<
0
(
surrounding
material
)
d
2
E
x
,
n
d
z
2
+
k
n
2
E
x
,
n
=
0
0
<
z
<
L
(
heating
material
)
d
2
E
x
,
N
+
1
d
z
2
+
k
N
+
1
2
E
x
,
N
+
1
=
0
z
>
L
surrounding
material
,where z
n
-
1
≤
z
≤
z
n and n
=
1
⋯
N. Here z
0
,
z
1
⋯
z
N denote the interface position. k
n
=
ω
μ
n
μ
0
(
ε
n
′
(
T
)
+
i
ε
n
′
′
(
T
)
)
ε
0 is the propagation constant which depends on dielectric constant ε
n
′
(
T
) and dielectric loss ε
n
′′
(
T
).  μ
n 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(3)
E
x
,
0
=
A
0
e
i
k
0
z
+
B
0
e
-
i
k
0
z
z
<
0
E
x
,
n
=
A
n
e
i
k
n
z
+
B
n
e
-
i
k
n
z
0
<
z
<
L
,
n
=
1
⋯
N
E
x
,
N
+
1
=
A
N
+
1
e
i
k
N
+
1
z
+
B
N
+
1
e
-
i
k
N
+
1
z
z
>
L
.The electric and magnetic field are continuous at the interface between different layers. That is as follows:(4)
E
x
,
n
=
E
x
,
n
+
1
n
=
0
⋯
N
,
d
E
x
,
n
d
z
=
d
E
x
,
n
+
1
d
z
z
=
z
0
⋯
z
N
.By using the interface conditions, the general solution coefficients can be obtained via solving the set of algebraic equations:(5)
A
n
e
i
×
k
n
z
n
+
B
n
e
-
i
×
k
n
z
n
-
A
n
+
1
e
i
×
k
n
+
1
z
n
-
B
n
+
1
e
-
i
×
k
n
+
1
z
n
=
0
k
n
A
n
e
i
×
k
n
z
n
-
k
n
B
n
e
-
i
×
k
n
z
n
-
k
n
+
1
A
n
+
1
e
i
×
k
n
+
1
z
n
+
k
n
+
1
B
n
+
1
e
-
i
×
k
n
+
1
z
n
=
0
h
h
h
h
h
h
h
h
h
h
h
h
h
h
n
=
0
⋯
N
.Let us denote byI
0 the incident power, which is the control variable in actual applications. Then the electric amplitude at z
<
0 can be calculated by(6)
A
0
=
E
0
=
2
I
0
c
ε
0
,where c is the velocity of light.Based on (5) and (6), the absorbed power in the nth layer, calculated by Poynting vector theorem, is(7)
P
z
=
1
2
ω
ε
0
ε
n
′′
E
x
,
n
z
E
x
,
n
∗
z
,where the superscript “∗” denotes the complex conjugate.


### 2.1. Case of Lambert’s Law

Lambert’s law is valid for semi-infinite material, and the power absorbed by the material per unit volume can be represented by [19](8)
P
(
z
)
=
1
2
ω
ε
0
ε
1
′′
T
r
,
1
2
E
0
2
e
-
2
β
z
,where T
r
,
1 represents the penetration coefficient calculated by the sampling temperature on material surface and β represents attenuation factor:(9)
T
r
,
1
=
2
μ
1
k
r
μ
1
k
r
+
μ
r
k
1
(10)
β
=
π
2
ε
1
′
1
+
ε
1
′′
/
ε
1
′
2
-
1
λ
0
.


### 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: (11)
E
1
(
z
)
=
T
r
,
1
E
0
e
-
β
1
z
(12)
E
n
(
z
)
=
E
n
-
1
(
z
n
-
1
)
e
-
β
n
(
z
-
z
n
-
1
)
n
=
2
⋯
N
(13)
P
n
z
=
1
2
ω
ε
0
ε
n
′′
E
n
2
z
n
=
1
⋯
N
.


### 2.3. The Comparison among the Three Methods

#### 2.3.1. Material with Constant Permittivity

The permittivity and thermophysical properties of the heating material are referred to in Table1. The microwave frequency used in this analysis is chosen to be 2.45 Ghz. The initial material temperature is set to be 280 K. The surrounding material is filled with air. The penetration depth can be calculated by 1
/
β, equal to 3.64 cm. For the entire material with the same temperature, the calculation of Lambert’s law and Lambert’s law combined with temperature has the same value. So Maxwell’s equation and Lambert’s law are utilized to calculate the difference in power distribution. As Figure 2(a) shows, for small thickness material, power distribution calculated by Lambert’s law has a huge difference with Maxwell’s equation. That means Lambert’s law cannot be utilized to material thickness less than penetration depth. For thickness a little larger than penetration depth, as shown in Figure 2(b), Lambert’s law can approximately describe power distribution. But it cannot reveal the phenomenon of power fluctuation. For thickness much larger than penetration depth, seeing Figure 2(c), Lambert’s law can exactly describe power distribution.Table 1
Dielectric and thermophysical properties.
Property
Variables (T in Kelvin)
Dielectric constant(ε
′)
-
0.0017
T
2
+
0.5299
T
+
73.3099
Dielectric loss(ε
′
′)
-
0.0021
T
2
+
1.2297
T
-
169.6084
Thermal conductivity
k
t (W·m−1·K−1
)
-
6.6964
×
1
0
-
6
T
2
+
0.0054
T
-
0.4095
Specific heat
C
p (J·K−1·Kg−1)
4080
Density
ρ (Kg·m−3)
990Figure 2
Power distribution with the input powerI
0
=
1.5
W
·
c
m
-
2, (a) thickness of 2 cm; (b) thickness of 5 cm; (c) thickness of 10 cm.
(a)
(b)
(c)


#### 2.3.2. Material with Varying Permittivity

In practical applications, different temperature leads to different permittivity. 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 Figure3(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.Figure 3
Material with thickness of 2 cm and input powerI
0
=
1.5
W
·
c
m
-
2, (a) temperature distribution; (b) power distribution.
(a)
(b)For 5-centimeter-thick material, with the temperature distribution of Figure4(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.Figure 4
Material with thickness of 5 cm and input powerI
0
=
1.5
W
·
c
m
-
2, (a) temperature distribution; (b) power distribution.
(a)
(b)For 10-centimeter-thick material, with the temperature distribution of Figure5(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.Figure 5
Material with thickness of 10 cm and input powerI
0
=
1.5
W
·
c
m
-
2, (a) temperature distribution; (b) power distribution.
(a)
(b)


### 2.4. Temperature Field Calculation

After power distribution is deduced, heat generation and transfer in the material can be calculated by the following equation:(14)
ρ
C
p
∂
T
(
z
,
t
)
∂
t
=
k
t
∂
2
T
(
z
,
t
)
∂
z
2
+
P
(
z
)with the initial condition T
(
z
,
t
=
0
)
=
T
0
(
z
), where ρ is material density (Kg
·
m
-
3); C
p is specific heat of the material (J
·
K
-
1
·
K
g
-
1); k
t is thermal conductivity of the material (W
·
m
-
1
·
K
-
1); T
(
z
,
t
) is material temperature (K); P
(
z
) is electromagnetic power dissipated per unit volume (W
·
m
-
3).The temperature satisfies the surface heat balance; boundary condition can be set by(15)
k
t
∂
T
(
t
)
∂
z
=
h
T
t
-
T
0
,
z
=
0
-
k
t
∂
T
t
∂
z
=
h
T
t
-
T
0
,
z
=
L
,where h is the convection heat transfer parameter (W
·
m
-
2
·
K) and T
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(16)
d
T
1
d
t
=
T
0
-
T
1
p
C
p
Δ
d
h
+
k
t
×
T
2
-
T
1
p
C
p
Δ
d
2
+
P
1
p
C
p
d
T
i
d
t
=
k
t
×
T
i
+
1
-
2
T
i
+
T
i
-
1
p
C
p
Δ
d
2
+
P
i
p
C
p
i
=
2
,
…
,
N
-
1
d
T
N
d
t
=
T
0
-
T
N
p
C
p
Δ
d
h
+
k
t
×
T
N
-
1
-
T
N
p
C
p
Δ
d
2
+
P
N
p
C
p
,where T
i and P
i (i
=
1
,
…
,
N) represent the temperature and power distribution at the ith layer. The power distribution P
i 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 ink
t and P
(
z
). The number of temperature sensors is set to be N
s, 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:(17)
T
1
t
+
T
s
=
T
1
t
+
T
0
-
T
1
ρ
C
p
d
s
h
T
s
+
k
t
t
×
T
2
-
T
1
ρ
C
p
d
s
2
T
s
+
P
1
ρ
C
p
T
s
T
i
t
+
T
s
=
T
i
t
+
k
t
t
×
T
i
+
1
-
2
T
i
+
T
i
-
1
ρ
C
p
d
s
2
T
s
+
P
i
ρ
C
p
T
s
i
=
2
,
…
,
N
s
-
1
T
N
s
t
+
T
s
=
T
N
s
t
+
T
0
-
T
N
s
ρ
C
p
d
s
h
T
s
+
k
t
t
×
T
N
s
-
1
-
T
N
s
ρ
C
p
d
s
2
T
s
+
P
N
s
ρ
C
p
T
s
,where T
i
(
t
) is the sampling temperature at the ith sampling point, d
s
=
L
/
N
s, L represents the length of the heating material, t represents the sampling moment, T
s represents the prediction time, and P
i is the power distribution at the ith sampling point and is solved by using (11)–(13).


### 2.5. 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:(18)
E
L
,
T
=
∑
i
=
1
N
T
i
Max
⁡
-
T
i
L
,
T
2
N
,where T
i
Max
⁡ presents the temperature calculated by Maxwell’s equation, T
i
L presents the temperature calculated by Lambert’s law, and T
i
T 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 of20
W
·
c
m
-
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.Table 2
The mean square error of different thickness.
Thickness (cm)
E
L
E
T
2
1.2102
0.5011
5
0.2744
0.1345
10
0.2128
0.1278Figure 6
Temperature field distribution byI
0
=
20
W
·
c
m
-
2 after 100 seconds, (a) thickness of 2 cm; (b) thickness of 5 cm; (c) thickness of 10 cm.
(a)
(b)
(c)


## 2.1. Case of Lambert’s Law

Lambert’s law is valid for semi-infinite material, and the power absorbed by the material per unit volume can be represented by [19](8)
P
(
z
)
=
1
2
ω
ε
0
ε
1
′′
T
r
,
1
2
E
0
2
e
-
2
β
z
,where T
r
,
1 represents the penetration coefficient calculated by the sampling temperature on material surface and β represents attenuation factor:(9)
T
r
,
1
=
2
μ
1
k
r
μ
1
k
r
+
μ
r
k
1
(10)
β
=
π
2
ε
1
′
1
+
ε
1
′′
/
ε
1
′
2
-
1
λ
0
.


## 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: (11)
E
1
(
z
)
=
T
r
,
1
E
0
e
-
β
1
z
(12)
E
n
(
z
)
=
E
n
-
1
(
z
n
-
1
)
e
-
β
n
(
z
-
z
n
-
1
)
n
=
2
⋯
N
(13)
P
n
z
=
1
2
ω
ε
0
ε
n
′′
E
n
2
z
n
=
1
⋯
N
.


## 2.3. The Comparison among the Three Methods

### 2.3.1. Material with Constant Permittivity

The permittivity and thermophysical properties of the heating material are referred to in Table1. The microwave frequency used in this analysis is chosen to be 2.45 Ghz. The initial material temperature is set to be 280 K. The surrounding material is filled with air. The penetration depth can be calculated by 1
/
β, equal to 3.64 cm. For the entire material with the same temperature, the calculation of Lambert’s law and Lambert’s law combined with temperature has the same value. So Maxwell’s equation and Lambert’s law are utilized to calculate the difference in power distribution. As Figure 2(a) shows, for small thickness material, power distribution calculated by Lambert’s law has a huge difference with Maxwell’s equation. That means Lambert’s law cannot be utilized to material thickness less than penetration depth. For thickness a little larger than penetration depth, as shown in Figure 2(b), Lambert’s law can approximately describe power distribution. But it cannot reveal the phenomenon of power fluctuation. For thickness much larger than penetration depth, seeing Figure 2(c), Lambert’s law can exactly describe power distribution.Table 1
Dielectric and thermophysical properties.
Property
Variables (T in Kelvin)
Dielectric constant(ε
′)
-
0.0017
T
2
+
0.5299
T
+
73.3099
Dielectric loss(ε
′
′)
-
0.0021
T
2
+
1.2297
T
-
169.6084
Thermal conductivity
k
t (W·m−1·K−1
)
-
6.6964
×
1
0
-
6
T
2
+
0.0054
T
-
0.4095
Specific heat
C
p (J·K−1·Kg−1)
4080
Density
ρ (Kg·m−3)
990Figure 2
Power distribution with the input powerI
0
=
1.5
W
·
c
m
-
2, (a) thickness of 2 cm; (b) thickness of 5 cm; (c) thickness of 10 cm.
(a)
(b)
(c)


### 2.3.2. Material with Varying Permittivity

In practical applications, different temperature leads to different permittivity. 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 Figure3(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.Figure 3
Material with thickness of 2 cm and input powerI
0
=
1.5
W
·
c
m
-
2, (a) temperature distribution; (b) power distribution.
(a)
(b)For 5-centimeter-thick material, with the temperature distribution of Figure4(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.Figure 4
Material with thickness of 5 cm and input powerI
0
=
1.5
W
·
c
m
-
2, (a) temperature distribution; (b) power distribution.
(a)
(b)For 10-centimeter-thick material, with the temperature distribution of Figure5(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.Figure 5
Material with thickness of 10 cm and input powerI
0
=
1.5
W
·
c
m
-
2, (a) temperature distribution; (b) power distribution.
(a)
(b)


## 2.3.1. Material with Constant Permittivity

The permittivity and thermophysical properties of the heating material are referred to in Table1. The microwave frequency used in this analysis is chosen to be 2.45 Ghz. The initial material temperature is set to be 280 K. The surrounding material is filled with air. The penetration depth can be calculated by 1
/
β, equal to 3.64 cm. For the entire material with the same temperature, the calculation of Lambert’s law and Lambert’s law combined with temperature has the same value. So Maxwell’s equation and Lambert’s law are utilized to calculate the difference in power distribution. As Figure 2(a) shows, for small thickness material, power distribution calculated by Lambert’s law has a huge difference with Maxwell’s equation. That means Lambert’s law cannot be utilized to material thickness less than penetration depth. For thickness a little larger than penetration depth, as shown in Figure 2(b), Lambert’s law can approximately describe power distribution. But it cannot reveal the phenomenon of power fluctuation. For thickness much larger than penetration depth, seeing Figure 2(c), Lambert’s law can exactly describe power distribution.Table 1
Dielectric and thermophysical properties.
Property
Variables (T in Kelvin)
Dielectric constant(ε
′)
-
0.0017
T
2
+
0.5299
T
+
73.3099
Dielectric loss(ε
′
′)
-
0.0021
T
2
+
1.2297
T
-
169.6084
Thermal conductivity
k
t (W·m−1·K−1
)
-
6.6964
×
1
0
-
6
T
2
+
0.0054
T
-
0.4095
Specific heat
C
p (J·K−1·Kg−1)
4080
Density
ρ (Kg·m−3)
990Figure 2
Power distribution with the input powerI
0
=
1.5
W
·
c
m
-
2, (a) thickness of 2 cm; (b) thickness of 5 cm; (c) thickness of 10 cm.
(a)
(b)
(c)


## 2.3.2. Material with Varying Permittivity

In practical applications, different temperature leads to different permittivity. 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 Figure3(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.Figure 3
Material with thickness of 2 cm and input powerI
0
=
1.5
W
·
c
m
-
2, (a) temperature distribution; (b) power distribution.
(a)
(b)For 5-centimeter-thick material, with the temperature distribution of Figure4(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.Figure 4
Material with thickness of 5 cm and input powerI
0
=
1.5
W
·
c
m
-
2, (a) temperature distribution; (b) power distribution.
(a)
(b)For 10-centimeter-thick material, with the temperature distribution of Figure5(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.Figure 5
Material with thickness of 10 cm and input powerI
0
=
1.5
W
·
c
m
-
2, (a) temperature distribution; (b) power distribution.
(a)
(b)


## 2.4. Temperature Field Calculation

After power distribution is deduced, heat generation and transfer in the material can be calculated by the following equation:(14)
ρ
C
p
∂
T
(
z
,
t
)
∂
t
=
k
t
∂
2
T
(
z
,
t
)
∂
z
2
+
P
(
z
)with the initial condition T
(
z
,
t
=
0
)
=
T
0
(
z
), where ρ is material density (Kg
·
m
-
3); C
p is specific heat of the material (J
·
K
-
1
·
K
g
-
1); k
t is thermal conductivity of the material (W
·
m
-
1
·
K
-
1); T
(
z
,
t
) is material temperature (K); P
(
z
) is electromagnetic power dissipated per unit volume (W
·
m
-
3).The temperature satisfies the surface heat balance; boundary condition can be set by(15)
k
t
∂
T
(
t
)
∂
z
=
h
T
t
-
T
0
,
z
=
0
-
k
t
∂
T
t
∂
z
=
h
T
t
-
T
0
,
z
=
L
,where h is the convection heat transfer parameter (W
·
m
-
2
·
K) and T
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(16)
d
T
1
d
t
=
T
0
-
T
1
p
C
p
Δ
d
h
+
k
t
×
T
2
-
T
1
p
C
p
Δ
d
2
+
P
1
p
C
p
d
T
i
d
t
=
k
t
×
T
i
+
1
-
2
T
i
+
T
i
-
1
p
C
p
Δ
d
2
+
P
i
p
C
p
i
=
2
,
…
,
N
-
1
d
T
N
d
t
=
T
0
-
T
N
p
C
p
Δ
d
h
+
k
t
×
T
N
-
1
-
T
N
p
C
p
Δ
d
2
+
P
N
p
C
p
,where T
i and P
i (i
=
1
,
…
,
N) represent the temperature and power distribution at the ith layer. The power distribution P
i 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 ink
t and P
(
z
). The number of temperature sensors is set to be N
s, 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:(17)
T
1
t
+
T
s
=
T
1
t
+
T
0
-
T
1
ρ
C
p
d
s
h
T
s
+
k
t
t
×
T
2
-
T
1
ρ
C
p
d
s
2
T
s
+
P
1
ρ
C
p
T
s
T
i
t
+
T
s
=
T
i
t
+
k
t
t
×
T
i
+
1
-
2
T
i
+
T
i
-
1
ρ
C
p
d
s
2
T
s
+
P
i
ρ
C
p
T
s
i
=
2
,
…
,
N
s
-
1
T
N
s
t
+
T
s
=
T
N
s
t
+
T
0
-
T
N
s
ρ
C
p
d
s
h
T
s
+
k
t
t
×
T
N
s
-
1
-
T
N
s
ρ
C
p
d
s
2
T
s
+
P
N
s
ρ
C
p
T
s
,where T
i
(
t
) is the sampling temperature at the ith sampling point, d
s
=
L
/
N
s, L represents the length of the heating material, t represents the sampling moment, T
s represents the prediction time, and P
i is the power distribution at the ith sampling point and is solved by using (11)–(13).


## 2.5. 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:(18)
E
L
,
T
=
∑
i
=
1
N
T
i
Max
⁡
-
T
i
L
,
T
2
N
,where T
i
Max
⁡ presents the temperature calculated by Maxwell’s equation, T
i
L presents the temperature calculated by Lambert’s law, and T
i
T 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 of20
W
·
c
m
-
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.Table 2
The mean square error of different thickness.
Thickness (cm)
E
L
E
T
2
1.2102
0.5011
5
0.2744
0.1345
10
0.2128
0.1278Figure 6
Temperature field distribution byI
0
=
20
W
·
c
m
-
2 after 100 seconds, (a) thickness of 2 cm; (b) thickness of 5 cm; (c) thickness of 10 cm.
(a)
(b)
(c)


## 3. 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 functionJ:(19)
min
⁡
J
U
p
=
∑
j
f
y
ref
j
,
y
p
j
,
u
j
-
1
∀
j
∈
k
+
1
,
…
,
k
+
N
p
,where k is the actual discrete time index, j is the discrete time index, and N
p is the receding horizon. y
ref contains the description of reference behavior; U
p is the sequence of the future input. Consider the following:(20)
U
p
=
⋯
u
j
⋯
T
∀
j
∈
k
,
…
,
k
+
N
c
-
1
,where N
c is the control horizon and the other inputs are the same in the future step:(21)
u
j
=
u
k
+
N
c
-
1
∀
j
∈
k
+
N
c
,
…
,
k
+
N
p
-
1
.For the problem of this paper, the input power has the constraint as follows:(22)
u
min
⁡
≤
u
j
≤
u
max
⁡
∀
j
∈
k
,
…
,
k
+
N
p
-
1
.In this analysis, the control model can be observed in Figure6. Maxwell’s equation is used to simulate the actual 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(23)
e
k
=
y
ref
k
-
y
p
k
.At each sampling timek, the error e
(
k
) is updated and remains the same over the prediction horizon N
p (Figure 7).Figure 7
Control model.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 functionJ can take the following form:(24)
min
⁡
J
=
∑
j
=
1
N
p
q
j
y
p
k
+
j
-
y
ref
k
+
j
2
+
∑
j
=
1
N
p
r
j
u
2
k
+
j
-
1
.Let us setQ
=
diag
[
q
1
,
q
2
,
…
,
q
N
p
] and R
=
diag
[
r
1
,
r
2
,
…
,
r
N
p
] as the weight coefficient. Y
ref
(
k
+
1
)
=
[
y
ref
(
k
+
1
)
,
…
,
y
ref
(
k
+
N
p
)
]
T and Y
p
(
k
)
=
[
y
p
(
k
+
1
)
,
…
,
y
p
(
k
+
N
p
)
]
T represent the reference temperature and the prediction temperature. B
=
[
P
I
/
ρ
C
ρ
,
…
,
P
I
/
ρ
C
ρ
]
T is the coefficient matrix between input power and temperature variation, where P
I represents the power value under per unit input power at the location of the maximum temperature. H
=
[
h
1
,
h
2
,
…
,
h
N
p
]
T is the error coefficient matrix. Let us set the control horizon equal to the prediction horizon; then the optimal input at each sampling time can be calculated by the following algorithms:(25)
U
p
=
B
T
Q
B
+
R
-
1
·
[
Y
ref
(
k
+
1
)
-
Y
p
(
k
)
+
H
×
e
(
k
)
-
W
(
k
)
]
u
k
=
U
p
1
B
1
,where (26)
W
k
=
T
0
-
T
1
k
ρ
C
ρ
d
s
h
+
k
t
t
×
T
2
k
-
T
1
k
ρ
C
p
d
s
2
,
…
,
h
l
×
T
2
k
-
T
1
k
ρ
C
p
d
s
2
,
…
,
N
p
h
l
·
T
0
-
T
1
k
ρ
C
ρ
d
s
h
+
k
t
t
×
T
2
k
-
T
1
k
ρ
C
p
d
s
2
Tis the coefficient matrix for the maximum temperature located at material left surface. And(27)
W
k
=
k
t
t
×
T
i
+
1
k
-
2
T
i
k
+
T
i
-
1
k
ρ
C
p
d
s
2
,
h
l
h
…
,
N
p
k
t
t
×
T
i
+
1
k
-
2
T
i
k
+
T
i
-
1
k
ρ
C
p
d
s
2
h
h
h
l
h
h
h
h
h
l
h
×
T
i
+
1
k
-
2
T
i
k
+
T
i
-
1
k
ρ
C
p
d
s
2
Tis the coefficient matrix for the maximum temperature located at the ith layer in the internal of the heating material. The coefficient matrix for the maximum temperature located at material right surface is as follows:(28)
W
k
=
T
0
-
T
N
s
k
ρ
C
ρ
d
s
h
+
k
t
t
×
T
N
s
-
1
k
-
T
N
s
k
ρ
C
p
d
s
2
,
h
h
…
,
N
p
T
0
-
T
N
s
k
ρ
C
ρ
d
s
h
+
k
t
t
×
T
N
s
-
1
k
-
T
N
s
k
ρ
C
p
d
s
2
h
h
h
h
h
h
h
h
l
h
×
T
N
s
-
1
k
-
T
N
s
k
ρ
C
p
d
s
2
T
.
U
p
(
1
) and B
(
1
) represent the first element in the columns U
p and B.


## 4. 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 Section2, 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, N
p
=
N
c
=
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
≤
u
≤
20
W
·
c
m
-
2. The distance between two adjacent sensors is selected to be d
s
=
1
mm.For the material with the thickness of 5 cm, the matrix parameters are set asQ
=
diag
[
5
×
1
0
9
,
5
×
1
0
9
,
5
×
1
0
9
], R
=
diag
[
0.6,0.6,0.6
], and H
=
[
0.1,0.1,0.1
]
T. The reference trajectory is set as y
ref
(
k
)
=
280
+
0.6
×
k
-
0.001
×
k
2, for 0
≤
k
≤
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.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.
(a)
(b)
(c)
(d)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 asQ
=
diag
[
1
0
9
,
1
0
9
,
1
0
9
], R
=
diag
[
0.08,0.08,0.08
], and H
=
[
-
0.9
,
-
0.9
,
-
0.9
]
T. 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).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.
(a)
(b)
(c)
(d)


## 5. 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.



---

*Source: 610490-2015-10-25.xml*