# System Identification of Heat-Transfer Process of Frequency Induction Furnace for Melting Copper Based on Particle Swarm Algorithm

Authors: Zhi-gang Jia; Xing-xuan Wang
Journal: Journal of Control Science and Engineering
 (2015)

---

## Abstract

An adaptive evolutionary strategy in standard particle swarm optimization is introduced. Adaptive evolution particle swarm optimization is constructed to improve the capacity of global search. A method based on adaptive evolution particle swarm optimization for identification of continuous system with time delay is proposed. The basic idea is that the identification of continuous system with time delay is converted to an optimization of continuous nonlinear function. The adaptive evolution particle swarm optimization is utilized to find an optimal solution of continuous nonlinear function. Convergence conditions are given by the convergence analysis based on discrete time linear dynamic system theory. Numerical simulation results show that the proposed method is effective for a general continuous system with time delay and the system of heat-transfer process of frequency induction furnace for melting copper.


---

## Body

## 1. Introduction

Frequency induction furnace has been developed into a kind of smelting equipment widely used. It is composed of four main parts: furnace body, electrical system, water cooling system, and hydraulic system [1]. Its power supply is industrial frequency (50 or 60 Hz) current. The research object in this paper is a vertical frequency induction furnace in a copper company. After adding raw material, workers regulate transformer output-voltage so that material melts. Copper smelting and heat preservation need lower transformer output-voltage. The molten copper is kept in proper temperature. Finally pour the molten copper [2].The heat-transfer process in frequency induction furnace system consists of copper melting process and heat-transfer process of molten copper to furnace wall. The approximation of most industrial processes is low-order continuous system. For the sake of convenience, the approximations of those two processes are first-order continuous systems, respectively [2]. The following mathematical model can be set up by frequency induction furnace input voltage and measured data of furnace wall temperature
(1)
G
s
=
K
T
0
s
+
1
T
1
s
+
1
,
where K, T
0, and T
1 denote heat-transfer gain, time constant of the process of copper melting, and time constant of heat-transfer process of molten copper to furnace wall [2]. But this industrial process has time delay. Therefore time delay is introduced to model (1) in literature [1]. The mathematical model of frequency induction furnace heat-transfer process is proposed as follows:
(2)
G
s
=
K
T
0
s
+
1
T
1
s
+
1
e
-
L
s
.
The controller design of frequency induction furnace heat-transfer process needs identification of system (2).Identification methods of second-order system with time delay in literature [3–6] set up identification model by integral computation of step response and measured data of input signal in time domain. Then identify parameters though Newton-Raphson method and least square method. Methods in literature [3–6] can only be used for second-order system with time delay.The continuous system identification methods based on particle swarm optimization (PSO) proposed in this paper are effective for general continuous system with time delay. The input signal has no limitation.


## 2. Adaptive Evolution PSO (AEPSO)

### 2.1. Standard PSO

The particle swarm optimization algorithm is a new evolutionary computation method proposed by Kenndy and Eberhart in 1995. It is a simulation of the natural phenomenon of a flock of birds searching for food. As a key optimization technique, PSO has been used extensively in many fields, including function optimization, neutral network, and fuzzy system control [7]. PSO algorithm belongs to the evolutionary algorithm similar to genetic algorithm. It starts from the random solutions and searches for the optimal solution through iteration [8].InD dimensional solution space, the position of ith particle can be expressed as a D dimensional vector
(3)
X
i
=
x
i
1
x
i
2
x
i
3
⋯
x
i
D
.
The velocity vector of ith particle is expressed as follows:
(4)
V
i
=
v
i
1
v
i
2
v
i
3
⋯
v
i
D
.
The position of the optimal solution of ith particle can be expressed as a D dimensional vector as follows:
(5)
P
i
=
p
i
1
p
i
2
p
i
3
⋯
p
i
D
.
The position of the optimal solution of the swarm which has the maximum fitness value in the swarm can be expressed as follows:
(6)
P
g
=
p
g
1
p
g
2
p
g
3
⋯
p
g
D
.
The dynamic range of the particles is as follows:
(7)
X
min
⁡
=
x
min
⁡
1
x
min
⁡
2
x
min
⁡
3
⋯
x
min
⁡
D
X
max
⁡
=
x
max
⁡
1
x
max
⁡
2
x
max
⁡
3
⋯
x
max
⁡
D
.
The dynamic range of the velocity of the particles is as follows:
(8)
V
max
⁡
=
v
max
⁡
1
v
max
⁡
2
v
max
⁡
3
⋯
v
max
⁡
D
.
f is a real function which is defined in D dimensional space. The fitness value of ith particle can be expressed as follows:
(9)
fitnes
s
i
=
f
X
i
.
Particles’ velocity and position are updated from generation t to t
+
1 according to the following equations:
(10)
V
i
t
+
1
=
ω
*
V
i
t
+
c
1
*
rand
1
*
P
i
t
-
X
i
t
+
c
2
*
rand
2
P
g
t
-
X
i
t
(11)
X
i
t
+
1
=
X
i
t
+
V
i
t
+
1
,
where rand
1 and rand
2 are random number from 0 to 1. c
1 and c
2 are local learning factor and global learning factor, respectively. c
2 is larger than c
1 generally. ω is inertia weight.Studies show that the standard PSO algorithm can find the local optimal solution quickly, but it also can produce premature convergence phenomenon. Each particle adjusts its search direction based on individual optimal position and global optimal position in standard PSO algorithm. But particles are more and more close to the global optimal particle along with iteration. Thus particles’ velocity decreases and tends to zero gradually. So this particle swarm loses ability to evolve further. The standard PSO is considered to be convergent at this time. In many cases, the convergence results do not converge to the global optimal solution, not even converge to a local optimal solution sometimes [9].In order to avoid premature convergence phenomenon, many improved PSO algorithms have been proposed. These algorithms tend to focus on improving inertia weightω or introduce special mutation operation to particles [10]. This paper will introduce adaptive inertia weight ω
(
t
) and special mutation operation to the swarm to avoid premature convergence phenomenon.


### 2.2. The Introduction of Adaptive Inertia Weight

PSO algorithm which has introduced adaptive inertia weight needs to do the following modifications. Equation (10) can be modified as follows:
(12)
V
i
t
+
1
=
w
(
t
)
*
V
i
t
+
c
1
*
rand
1
*
P
i
t
-
X
i
t
+
c
2
*
rand
2
*
P
g
t
-
X
i
t
.Time-varying inertia weightω
(
t
) has been used in (12). If a linear decreasing inertia weight ω
(
t
) is used in iterative process, the PSO algorithm has good global search capability in the beginning and has a good local search ability in the later stage [11]. A good global search capability makes the PSO algorithm locate the area around the global optimal solution quickly in the beginning. The good local search ability makes the PSO algorithm obtain the global optimal solution accurately in the later stage. The adaptive inertia weight can be determined by the following equation in iteration:
(13)
ω
t
=
ω
max
⁡
-
ω
max
⁡
-
ω
min
⁡
G
*
t
,
where G is the number of evolution generation. ω
max
⁡ and ω
min
⁡ are constants. Let adaptive inertia weight ω
(
t
) vary in the range of [
ω
min
⁡
,
ω
max
⁡
]. ω
max
⁡ and ω
min
⁡ are usually 0.9 and 0.1, respectively, according to experience.


### 2.3. Mutation Operation

Mutation operation is a natural biological model of gene mutation induced by various accidental factors in genetic environment. It takes very small probability to change the position of the particle randomly. In order to obtain high quality optimal solution as far as possible, mutation operation must be taken. Mutation operation plays an auxiliary role in the genetic algorithm but is a necessary operation step of genetic algorithm. In this paper, the basic mutation strategy of genetic algorithm is introduced to the standard PSO algorithm [12]. The particles which will do mutation operation are selected according to the probability. Then the positions of the selected particles do mutation operation. The position density distribution model of each particle is considered to be uniform distribution density function in this paper [13]
(14)
p
X
i
j
=
1
X
max
⁡
j
-
X
min
⁡
j
,
where X
i
j denotes jth parameter of ith mutation particle. Then a new particle can have the following representation:
(15)
X
i
=
X
min
⁡
+
X
max
⁡
-
X
min
⁡
*
r
,
where r is a random variable which satisfies uniform distribution in the range of [
0,1
].AEPSO consists of the following steps.(1)
Initialization: generateN particles in D dimensional solution space which have random position and velocity vectors.
(2)
Fitness evaluation: calculate the fitness value ofith particle and find out the optimal particle in the swarm which has the maximum fitness value, i
=
1,2
,
3
,
…
,
N.
(3)
Updating of particles’ position and velocity: update the position and velocity of each particle with (12) and (11), respectively. Check whether the position and velocity of each particle are out of bounds.
(4)
Mutation operation: the positions of particles do mutation operation with (15).
(5)
Comparison of current fitness value ofith particle f
(
X
i
) and individual optimal fitness value of ith particle f
(
P
i
): if f
(
X
i
)
>
f
(
P
i
), then P
i
=
X
i.
(6)
Comparison of current fitness value ofith particle f
(
X
i
) and global optimal fitness value of the swarm f
(
P
g
): if f
(
X
i
)
>
f
(
P
g
), then P
g
=
X
i.
(7)
Checking criterion for stopping iteration of algorithm: if the criterion is satisfied, the iteration of the algorithm is stopped. Otherwise, go to Step 3 and continue searching for the optimal solution until the criterion is satisfied.


## 2.1. Standard PSO

The particle swarm optimization algorithm is a new evolutionary computation method proposed by Kenndy and Eberhart in 1995. It is a simulation of the natural phenomenon of a flock of birds searching for food. As a key optimization technique, PSO has been used extensively in many fields, including function optimization, neutral network, and fuzzy system control [7]. PSO algorithm belongs to the evolutionary algorithm similar to genetic algorithm. It starts from the random solutions and searches for the optimal solution through iteration [8].InD dimensional solution space, the position of ith particle can be expressed as a D dimensional vector
(3)
X
i
=
x
i
1
x
i
2
x
i
3
⋯
x
i
D
.
The velocity vector of ith particle is expressed as follows:
(4)
V
i
=
v
i
1
v
i
2
v
i
3
⋯
v
i
D
.
The position of the optimal solution of ith particle can be expressed as a D dimensional vector as follows:
(5)
P
i
=
p
i
1
p
i
2
p
i
3
⋯
p
i
D
.
The position of the optimal solution of the swarm which has the maximum fitness value in the swarm can be expressed as follows:
(6)
P
g
=
p
g
1
p
g
2
p
g
3
⋯
p
g
D
.
The dynamic range of the particles is as follows:
(7)
X
min
⁡
=
x
min
⁡
1
x
min
⁡
2
x
min
⁡
3
⋯
x
min
⁡
D
X
max
⁡
=
x
max
⁡
1
x
max
⁡
2
x
max
⁡
3
⋯
x
max
⁡
D
.
The dynamic range of the velocity of the particles is as follows:
(8)
V
max
⁡
=
v
max
⁡
1
v
max
⁡
2
v
max
⁡
3
⋯
v
max
⁡
D
.
f is a real function which is defined in D dimensional space. The fitness value of ith particle can be expressed as follows:
(9)
fitnes
s
i
=
f
X
i
.
Particles’ velocity and position are updated from generation t to t
+
1 according to the following equations:
(10)
V
i
t
+
1
=
ω
*
V
i
t
+
c
1
*
rand
1
*
P
i
t
-
X
i
t
+
c
2
*
rand
2
P
g
t
-
X
i
t
(11)
X
i
t
+
1
=
X
i
t
+
V
i
t
+
1
,
where rand
1 and rand
2 are random number from 0 to 1. c
1 and c
2 are local learning factor and global learning factor, respectively. c
2 is larger than c
1 generally. ω is inertia weight.Studies show that the standard PSO algorithm can find the local optimal solution quickly, but it also can produce premature convergence phenomenon. Each particle adjusts its search direction based on individual optimal position and global optimal position in standard PSO algorithm. But particles are more and more close to the global optimal particle along with iteration. Thus particles’ velocity decreases and tends to zero gradually. So this particle swarm loses ability to evolve further. The standard PSO is considered to be convergent at this time. In many cases, the convergence results do not converge to the global optimal solution, not even converge to a local optimal solution sometimes [9].In order to avoid premature convergence phenomenon, many improved PSO algorithms have been proposed. These algorithms tend to focus on improving inertia weightω or introduce special mutation operation to particles [10]. This paper will introduce adaptive inertia weight ω
(
t
) and special mutation operation to the swarm to avoid premature convergence phenomenon.


## 2.2. The Introduction of Adaptive Inertia Weight

PSO algorithm which has introduced adaptive inertia weight needs to do the following modifications. Equation (10) can be modified as follows:
(12)
V
i
t
+
1
=
w
(
t
)
*
V
i
t
+
c
1
*
rand
1
*
P
i
t
-
X
i
t
+
c
2
*
rand
2
*
P
g
t
-
X
i
t
.Time-varying inertia weightω
(
t
) has been used in (12). If a linear decreasing inertia weight ω
(
t
) is used in iterative process, the PSO algorithm has good global search capability in the beginning and has a good local search ability in the later stage [11]. A good global search capability makes the PSO algorithm locate the area around the global optimal solution quickly in the beginning. The good local search ability makes the PSO algorithm obtain the global optimal solution accurately in the later stage. The adaptive inertia weight can be determined by the following equation in iteration:
(13)
ω
t
=
ω
max
⁡
-
ω
max
⁡
-
ω
min
⁡
G
*
t
,
where G is the number of evolution generation. ω
max
⁡ and ω
min
⁡ are constants. Let adaptive inertia weight ω
(
t
) vary in the range of [
ω
min
⁡
,
ω
max
⁡
]. ω
max
⁡ and ω
min
⁡ are usually 0.9 and 0.1, respectively, according to experience.


## 2.3. Mutation Operation

Mutation operation is a natural biological model of gene mutation induced by various accidental factors in genetic environment. It takes very small probability to change the position of the particle randomly. In order to obtain high quality optimal solution as far as possible, mutation operation must be taken. Mutation operation plays an auxiliary role in the genetic algorithm but is a necessary operation step of genetic algorithm. In this paper, the basic mutation strategy of genetic algorithm is introduced to the standard PSO algorithm [12]. The particles which will do mutation operation are selected according to the probability. Then the positions of the selected particles do mutation operation. The position density distribution model of each particle is considered to be uniform distribution density function in this paper [13]
(14)
p
X
i
j
=
1
X
max
⁡
j
-
X
min
⁡
j
,
where X
i
j denotes jth parameter of ith mutation particle. Then a new particle can have the following representation:
(15)
X
i
=
X
min
⁡
+
X
max
⁡
-
X
min
⁡
*
r
,
where r is a random variable which satisfies uniform distribution in the range of [
0,1
].AEPSO consists of the following steps.(1)
Initialization: generateN particles in D dimensional solution space which have random position and velocity vectors.
(2)
Fitness evaluation: calculate the fitness value ofith particle and find out the optimal particle in the swarm which has the maximum fitness value, i
=
1,2
,
3
,
…
,
N.
(3)
Updating of particles’ position and velocity: update the position and velocity of each particle with (12) and (11), respectively. Check whether the position and velocity of each particle are out of bounds.
(4)
Mutation operation: the positions of particles do mutation operation with (15).
(5)
Comparison of current fitness value ofith particle f
(
X
i
) and individual optimal fitness value of ith particle f
(
P
i
): if f
(
X
i
)
>
f
(
P
i
), then P
i
=
X
i.
(6)
Comparison of current fitness value ofith particle f
(
X
i
) and global optimal fitness value of the swarm f
(
P
g
): if f
(
X
i
)
>
f
(
P
g
), then P
g
=
X
i.
(7)
Checking criterion for stopping iteration of algorithm: if the criterion is satisfied, the iteration of the algorithm is stopped. Otherwise, go to Step 3 and continue searching for the optimal solution until the criterion is satisfied.


## 3. System Identification of General Continuous System with Time Delay

Consider the general continuous system as follows:(16)
g
s
=
b
m
s
m
+
b
m
-
1
s
m
-
1
+
⋯
+
b
1
s
+
b
0
a
n
s
n
+
a
n
-
1
s
n
-
1
+
⋯
+
a
1
s
+
a
0
e
-
L
s
,
where m
⩽
n. The time delay of such a general continuous system needs approximate processing according to the identification method of second-order continuous system with time delay in literature [3]. The time delay can be represented by Taylor series as follows:
(17)
e
-
L
s
=
1
-
L
s
+
L
s
2
2
+
⋯
.Now the parameter identification ofL becomes very difficult. After approximate processing, L has high math power. The form of the general continuous system will become complicated after multiplication. Parameters of second-order continuous system with time delay can be identified well by Newton-Raphson method. But it does not work well for high-order system. The identification of the general continuous system can be converted to optimization of nonlinear function. Then the AEPSO algorithm mentioned in this paper can be used to solve this optimization problem. The parameters which need to be identified can be expressed as the following vector:
(18)
θ
=
b
m
b
m
-
1
⋯
b
1
a
n
a
n
-
1
⋯
a
1
L
.
The nonlinear function which needs to be optimized can be expressed as follows:
(19)
J
θ
^
=
Y
-
Y
^
T
Y
-
Y
^
,
where Y is the measured output vector of real system. Y
^ is calculated output vector of estimation model. The fitness function f can be expressed as follows:
(20)
f
=
1
J
θ
^
.
Therefore the parameter identification problem of general continuous system can be regarded as a standard optimization problem, minimization of (19) according to the input and output sequences. Search for the parameter vector which can minimize (19) in n
+
m
+
1 dimensional space.


## 4. Analysis of Convergence

Analyze one element of the parameter vector without loss of generality:(21)
V
i
j
t
+
1
=
ω
*
V
i
j
t
+
c
1
*
rand
1
*
P
i
j
t
-
X
i
j
t
+
c
2
*
rand
2
P
g
j
t
-
X
i
j
t
X
i
j
t
+
1
=
X
i
j
t
+
V
i
j
t
+
1
.
Equations (21) can be simplified as follows:
(22)
v
t
+
1
=
ω
v
t
+
c
1
*
rand
1
*
p
i
-
x
t
+
c
2
*
rand
2
*
p
g
-
x
t
x
t
+
1
=
x
t
+
v
t
+
1
.
Equations (22) can be converted into the following discrete dynamic system:
(23)
v
t
+
1
=
ω
v
t
-
c
1
*
rand
1
+
c
2
*
rand
2
·
x
t
-
c
1
*
rand
1
*
p
i
+
c
2
*
rand
2
*
p
g
c
1
*
rand
1
+
c
2
*
rand
2
x
t
+
1
-
c
1
*
rand
1
*
p
i
+
c
2
*
rand
2
*
p
g
c
1
*
rand
1
+
c
2
*
rand
2
=
ω
v
t
+
1
-
c
1
*
rand
1
+
c
2
*
rand
2
·
x
t
-
c
1
*
rand
1
*
p
i
+
c
2
*
rand
2
*
p
g
c
1
*
rand
1
+
c
2
*
rand
2
.
The individual optimal position of a particle and the global optimal position of the swarm are usually unchangeable. The above two equations can be simplified as follows:
(24)
v
t
+
1
=
ω
v
t
-
c
z
t
z
t
+
1
=
ω
v
t
+
1
-
c
z
t
.
So the above system can be rewritten as a discrete autonomous system:
(25)
Y
t
+
1
=
A
Y
t
,
where Y
t
=
v
t
z
t
T, A
=
w
-
c
w
1
-
c.The equilibrium in Lyapunov sense can be obtained by the following equations:(26)
I
-
A
Y
e
=
0
.
If det
⁡
I
-
A
=
c
≠
0, there is only one equilibrium. The characteristic equation of autonomous system can be expressed as follows:
(27)
λ
2
-
1
-
c
+
ω
λ
+
ω
=
0
.
Characteristic root is the two following equations:
(28)
λ
1
=
ω
+
1
-
c
+
ω
+
1
-
c
2
-
4
ω
2
λ
2
=
ω
+
1
-
c
-
ω
+
1
-
c
2
-
4
ω
2
.Theorem 1.
The sufficient and necessary condition of autonomous system (25) converging to the equilibrium state is as follows:
(29)
max
⁡
λ
1
λ
2
<
1
.
Therefore the sufficient and necessary condition can be obtained as follows:
(30)
ω
<
1
0
<
c
<
4
.


## 5. Numerical Simulation

### 5.1. Experiment One

Consider the following third-order continuous system with time delay. The input signal is random signal sequence with amplitude of 4:(31)
g
s
=
4
s
2
+
6
s
4
s
3
+
15
s
2
+
17.5
s
+
1
e
-
0.8
s
.
There are six parameters needed to be identified. They can be expressed as the following vector:
(32)
θ
=
b
2
b
1
a
3
a
2
a
1
L
.The settings of identification method based on AEPSO algorithm are as follows. The identification windows is set to 1000. The number of particles N
=
200. The number of iteration generationsTimes = 1000. The local learning factor c
1
=
1.3. The global learning factor c
2
=
1.7. The inertia weight ω decreases form 0.9 to 0.1 with the iteration. The search range of parameters mentioned above is [
0,20
].The results of the nonlinear function needing to be optimized are shown in Figure1.Figure 1
The results of nonlinear function needing to be optimized versusTimes.The errorδ of parameters estimation can be calculated by the following equation:
(33)
δ
=
θ
^
-
θ
θ
×
100
%
.
The error δ of parameters estimation is shown in Figure 2.Figure 2
The errorδ of parameters estimation versusTimes.The parameters estimation is shown in Table1. From Figures 1-2, we can see that the results of the nonlinear optimized function and the error of parameters estimation δ are becoming smaller asTimes increases. From Table 1, the parameters estimation is becoming closer to the true values with iteration. The identification method based on AEPSO algorithm is effective for continuous system with time delay according to the results of experiment one.Table 1
Estimation of parameters.
Times
b
2
b
1
a
3
a
2
a
1
L
δ (%)
1
4.3001
3.6775
4.7041
15.8040
4.4133
0.8062
54.46
100
3.9306
3.0649
3.8284
12.0910
9.8059
0.8043
35.67
200
3.6863
3.2186
3.6900
11.7114
10.7010
0.8015
32.91
300
3.6091
3.3414
3.6158
11.6022
10.7973
0.8000
32.62
400
2.9136
4.1961
2.9262
10.7083
12.8079
0.8003
27.70
500
2.5226
5.4856
2.5382
10.9166
15.6515
0.7986
20.28
600
2.6645
5.7444
2.6757
11.5040
16.3083
0.7994
16.95
700
3.7783
6.0790
3.7766
14.5121
17.5467
0.7996
2.40
800
3.9280
6.0538
3.9294
14.8759
17.6169
0.8000
0.84
900
3.9867
6.0446
3.9872
15.0075
17.6055
0.8000
0.47
1000
3.9997
6.0106
4.0000
15.0090
17.5290
0.8000
0.14
True value
4
6
4
15
17.5
0.8Now we will use MATLAB system identification toolbox to identify that system, because the MATLAB identification toolbox identifies process system by fitting the transfer function of process system. But due to the limitation of MATLAB identification toolbox, it can identify the process system model as follows:(34)
K
1
+
T
z
*
s
1
+
2
ζ
*
T
w
s
+
T
w
2
s
2
1
+
T
p
*
s
e
-
T
d
*
s
.The process system transfer function obtained by MATLAB system identification toolbox can be shown as follows:(35)
5.78
×
1
0
4
s
+
0.005793
2.949
s
3
+
9.687
s
2
+
16.66
s
+
1
e
-
0.00182
s
.It can be obviously seen that the transfer function obtained by MATLAB system identification toolbox is completely wrong. The parameters which we need cannot be identified efficiently by the identification toolbox. So we can draw the conclusion that identification algorithm mentioned in this paper is more effective than traditional identification toolbox on identification of nonlinear system.


### 5.2. Experiment Two

Identify a second-order system which is similar to the heat-transfer system of frequency induction furnace. The second-order continuous system with time delay is as follows. The input signal is unit step signal(36)
G
s
=
6
6
s
+
1
20
s
+
1
e
-
4
s
.There are four parameters needed to be identified. They can be expressed as the following vector:(37)
θ
=
K
T
1
T
2
L
.The settings of identification method based on AEPSO algorithm are as follows. The identification windows is set to 1000. The number of particles N
=
100. The number of iteration generations T
i
m
e
s
=
200. The local learning factor c
1
=
1.3. The global learning factor c
2
=
1.7. The inertia weight ω decreases form 0.9 to 0.1 with the iteration. The search range of parameters mentioned above is [
0
,
30
].The results of the nonlinear optimized function and the errorδ of parameters estimation are shown in Figures 3 and 4, respectively.Figure 3
The results of nonlinear function needing to be optimized versusTimes.Figure 4
The errorδ of parameters estimation versusTimes.The parameters estimation of the second-order system with iteration is shown in Table2.Table 2
Estimation of parameters.
Times
K
T
1
T
2
L
δ (%)
1
6.0794
5.1238
17.8664
6.3655
14.96
20
5.9288
6.1143
17.8270
5.1525
11.15
40
6.0121
6.5824
20.1386
3.5830
3.30
60
6.0044
6.2609
19.8547
3.8584
1.50
80
5.9985
6.2603
19.8246
3.8509
1.57
100
5.9962
6.2567
19.8341
3.8470
1.55
120
6.0000
6.0510
19.9795
3.9592
0.31
140
6.0000
6.0381
19.9858
3.9694
0.23
160
5.9998
5.9995
19.9996
3.9979
0.01
180
5.9999
6.0057
19.9969
3.9968
0.03
200
6.0000
6.0040
19.9973
3.9978
0.02
True value
6
6
20
4From Table2, we can see that the step response of continuous system can be used to identify the system. The proposed identification method in this paper has no restrictions on input signal.Now we will identify that system by MATLAB system identification toolbox. Fit the following transfer function model:(38)
K
1
+
T
p
1
*
s
1
+
T
p
2
*
s
e
-
T
d
*
s
.The transfer function obtained by MATLAB system identification toolbox can be shown as follows:(39)
5.9902
1
+
19.9651
s
1
+
5.9793
s
e
-
4.9955
s
.The errorδ of parameters estimated by the identification method mentioned by this paper is 0.02%. The error δ
1 of parameters obtained by MATLAB identification toolbox is 4.51%. Because of δ
<
δ
1, the method mentioned in this paper is obviously superior to the system identification toolbox. In order to observe this kind of superiority more directly, the step response of the system transfer function identified by the mentioned method in this paper and MATLAB identification toolbox can be painted in Figure 5.Figure 5
The step response of the system transfer function identified.The step responses of transfer function obtained by different methods are painted in Figure5. The blue circles represent step response of system identified by MATLAB identification toolbox. The green solid line represents step response of system identified by method mentioned in this paper. The red rectangles represent step response of the original system. The system transfer function obtained by system identification toolbox in MATLAB is very accurate and bestfits has reached 95.22%. But by comparison of the error of parameter estimated and step response in Figure 5, we can draw the conclusion that the identification method mentioned in this paper is more effective than traditional methods.


## 5.1. Experiment One

Consider the following third-order continuous system with time delay. The input signal is random signal sequence with amplitude of 4:(31)
g
s
=
4
s
2
+
6
s
4
s
3
+
15
s
2
+
17.5
s
+
1
e
-
0.8
s
.
There are six parameters needed to be identified. They can be expressed as the following vector:
(32)
θ
=
b
2
b
1
a
3
a
2
a
1
L
.The settings of identification method based on AEPSO algorithm are as follows. The identification windows is set to 1000. The number of particles N
=
200. The number of iteration generationsTimes = 1000. The local learning factor c
1
=
1.3. The global learning factor c
2
=
1.7. The inertia weight ω decreases form 0.9 to 0.1 with the iteration. The search range of parameters mentioned above is [
0,20
].The results of the nonlinear function needing to be optimized are shown in Figure1.Figure 1
The results of nonlinear function needing to be optimized versusTimes.The errorδ of parameters estimation can be calculated by the following equation:
(33)
δ
=
θ
^
-
θ
θ
×
100
%
.
The error δ of parameters estimation is shown in Figure 2.Figure 2
The errorδ of parameters estimation versusTimes.The parameters estimation is shown in Table1. From Figures 1-2, we can see that the results of the nonlinear optimized function and the error of parameters estimation δ are becoming smaller asTimes increases. From Table 1, the parameters estimation is becoming closer to the true values with iteration. The identification method based on AEPSO algorithm is effective for continuous system with time delay according to the results of experiment one.Table 1
Estimation of parameters.
Times
b
2
b
1
a
3
a
2
a
1
L
δ (%)
1
4.3001
3.6775
4.7041
15.8040
4.4133
0.8062
54.46
100
3.9306
3.0649
3.8284
12.0910
9.8059
0.8043
35.67
200
3.6863
3.2186
3.6900
11.7114
10.7010
0.8015
32.91
300
3.6091
3.3414
3.6158
11.6022
10.7973
0.8000
32.62
400
2.9136
4.1961
2.9262
10.7083
12.8079
0.8003
27.70
500
2.5226
5.4856
2.5382
10.9166
15.6515
0.7986
20.28
600
2.6645
5.7444
2.6757
11.5040
16.3083
0.7994
16.95
700
3.7783
6.0790
3.7766
14.5121
17.5467
0.7996
2.40
800
3.9280
6.0538
3.9294
14.8759
17.6169
0.8000
0.84
900
3.9867
6.0446
3.9872
15.0075
17.6055
0.8000
0.47
1000
3.9997
6.0106
4.0000
15.0090
17.5290
0.8000
0.14
True value
4
6
4
15
17.5
0.8Now we will use MATLAB system identification toolbox to identify that system, because the MATLAB identification toolbox identifies process system by fitting the transfer function of process system. But due to the limitation of MATLAB identification toolbox, it can identify the process system model as follows:(34)
K
1
+
T
z
*
s
1
+
2
ζ
*
T
w
s
+
T
w
2
s
2
1
+
T
p
*
s
e
-
T
d
*
s
.The process system transfer function obtained by MATLAB system identification toolbox can be shown as follows:(35)
5.78
×
1
0
4
s
+
0.005793
2.949
s
3
+
9.687
s
2
+
16.66
s
+
1
e
-
0.00182
s
.It can be obviously seen that the transfer function obtained by MATLAB system identification toolbox is completely wrong. The parameters which we need cannot be identified efficiently by the identification toolbox. So we can draw the conclusion that identification algorithm mentioned in this paper is more effective than traditional identification toolbox on identification of nonlinear system.


## 5.2. Experiment Two

Identify a second-order system which is similar to the heat-transfer system of frequency induction furnace. The second-order continuous system with time delay is as follows. The input signal is unit step signal(36)
G
s
=
6
6
s
+
1
20
s
+
1
e
-
4
s
.There are four parameters needed to be identified. They can be expressed as the following vector:(37)
θ
=
K
T
1
T
2
L
.The settings of identification method based on AEPSO algorithm are as follows. The identification windows is set to 1000. The number of particles N
=
100. The number of iteration generations T
i
m
e
s
=
200. The local learning factor c
1
=
1.3. The global learning factor c
2
=
1.7. The inertia weight ω decreases form 0.9 to 0.1 with the iteration. The search range of parameters mentioned above is [
0
,
30
].The results of the nonlinear optimized function and the errorδ of parameters estimation are shown in Figures 3 and 4, respectively.Figure 3
The results of nonlinear function needing to be optimized versusTimes.Figure 4
The errorδ of parameters estimation versusTimes.The parameters estimation of the second-order system with iteration is shown in Table2.Table 2
Estimation of parameters.
Times
K
T
1
T
2
L
δ (%)
1
6.0794
5.1238
17.8664
6.3655
14.96
20
5.9288
6.1143
17.8270
5.1525
11.15
40
6.0121
6.5824
20.1386
3.5830
3.30
60
6.0044
6.2609
19.8547
3.8584
1.50
80
5.9985
6.2603
19.8246
3.8509
1.57
100
5.9962
6.2567
19.8341
3.8470
1.55
120
6.0000
6.0510
19.9795
3.9592
0.31
140
6.0000
6.0381
19.9858
3.9694
0.23
160
5.9998
5.9995
19.9996
3.9979
0.01
180
5.9999
6.0057
19.9969
3.9968
0.03
200
6.0000
6.0040
19.9973
3.9978
0.02
True value
6
6
20
4From Table2, we can see that the step response of continuous system can be used to identify the system. The proposed identification method in this paper has no restrictions on input signal.Now we will identify that system by MATLAB system identification toolbox. Fit the following transfer function model:(38)
K
1
+
T
p
1
*
s
1
+
T
p
2
*
s
e
-
T
d
*
s
.The transfer function obtained by MATLAB system identification toolbox can be shown as follows:(39)
5.9902
1
+
19.9651
s
1
+
5.9793
s
e
-
4.9955
s
.The errorδ of parameters estimated by the identification method mentioned by this paper is 0.02%. The error δ
1 of parameters obtained by MATLAB identification toolbox is 4.51%. Because of δ
<
δ
1, the method mentioned in this paper is obviously superior to the system identification toolbox. In order to observe this kind of superiority more directly, the step response of the system transfer function identified by the mentioned method in this paper and MATLAB identification toolbox can be painted in Figure 5.Figure 5
The step response of the system transfer function identified.The step responses of transfer function obtained by different methods are painted in Figure5. The blue circles represent step response of system identified by MATLAB identification toolbox. The green solid line represents step response of system identified by method mentioned in this paper. The red rectangles represent step response of the original system. The system transfer function obtained by system identification toolbox in MATLAB is very accurate and bestfits has reached 95.22%. But by comparison of the error of parameter estimated and step response in Figure 5, we can draw the conclusion that the identification method mentioned in this paper is more effective than traditional methods.


## 6. Conclusion

The paper discusses how to identify the heat-transfer process system of frequency induction furnace. An identification method of general continuous system with time delay based on AEPSO algorithm has been proposed. Experiment one demonstrates the proposed identification method’s effectiveness for general continuous system. Experiment two shows that the heat-transfer process system of frequency induction system can be identified by the proposed method. The proposed method has advantage over the method in literature [3–6]. There are no special restrictions on the input signal. Through the comparison experiment with traditional MATLAB identification toolbox, we can see that the proposed method is more effective.



---

*Source: 502828-2015-02-25.xml*