# KPCA-ESN Soft-Sensor Model of Polymerization Process Optimized by Biogeography-Based Optimization Algorithm

Authors: Wen-hua Cui; Jie-sheng Wang; Shu-xia Li
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

For solving the problem that the conversion rate of vinyl chloride monomer (VCM) is
hard for real-time online measurement in the polyvinyl chloride (PVC) polymerization production
process, a soft-sensor modeling method based on echo state network (ESN) is put forward. By
analyzing PVC polymerization process ten secondary variables are selected as input variables of the
soft-sensor model, and the kernel principal component analysis (KPCA) method is carried out on
the data preprocessing of input variables, which reduces the dimensions of the high-dimensional
data. Thek-means clustering method is used to divide data samples into several clusters as inputs of
each submodel. Then for each submodel the biogeography-based optimization algorithm (BBOA)
is used to optimize the structure parameters of the ESN to realize the nonlinear mapping between
input and output variables of the soft-sensor model. Finally, the weighted summation of outputs of
each submodel is selected as the final output. The simulation results show that the proposed
soft-sensor model can significantly improve the prediction precision of conversion rate and
conversion velocity in the process of PVC polymerization and can satisfy the real-time control
requirement of the PVC polymerization process.


---

## Body

## 1. Introduction

Polyvinyl chloride (PVC) is one of the largest plastic products in the world. Because PVC has characteristics of high strength and good stability, it has become a widely used synthetic material in the world today. According to different purposes, by adding the different additives or plasticizers, all kinds of PVC plastic products can be produced, such as plates, doors and windows profiles, pipe fittings, foam materials, and electrical parts. These products have widespread applications in industry, agriculture, health care and daily necessities, and other fields. As the embodiment of the superiority of PVC, the PVC production and improvement of the technology get more and more attention of people.PVC is mainly produced by vinyl chloride monomer (VCM), so the quality of VCM directly affects the quality of PVC resin, production costs, and economic benefits [1]. With the large-scale development trend of polymerization kettle and the in-depth research of vinyl chloride polymerization, further improving the conversion rate of polymerization kettle for productivity improvement and reducing the production cost have important significance. The different VCM conversion has a certain impact on the molecular weight of PVC resin, thermal stability, porosity, the residues of VCM, the absorptivity of plasticizers, and processing liquidity. As a result of the immature detection device, the complexity of the suspension polymerization reaction, and the strong nonlinear and strong coupling, in the actual production process, vinyl chloride conversion rate and conversion velocity are hard to acquire real-time, so it is difficult to achieve direct closed-loop control [2].Echo state network (ESN) is a new type of recursion neural network [3, 4]. In terms of network structure and learning mechanism, it is different from the previous cycle networks. ESN has better nonlinearity approximation capability, which makes it get good performance in the nonlinear prediction fields. A kind of wavelet decomposition echo state network predictive model was proposed, which adopts the wavelet decomposition method to match different ESN models at each scale with different properties. Then the least square method realizes the optimal integration of the predictive components through weight coefficients so as to reach accurate prediction of each scale and integration [5]. ESN was adopted to predict the conditions of flue gas turbine, and the singular decomposition method was used to carry on the modification of the linear regression training algorithm of ESN [6]. The autocorrelation function method was used to construct the input sequence of ESN for setting up the time series forecast method [7], which is used in the field of the mobile communication traffic prediction. The ESN prediction model based on the principal component analysis method was established in order to reduce the training time and improve the forecast speed [8]. However, the calculation of the ESN output weights based on the standard linear regression algorithm may easily lead to the pathological solutions when dealing with the practical problems. On the other hand, the output weights are often with larger amplitudes. In order to conquer the ill-conditioned matrix problems of traditional echo state network model, the evolutionary algorithms, such as genetic algorithms (GAs) [9, 10], particle swarm optimization (PSO) algorithm [11, 12], memetic algorithm (MA) [13–15], and ant colony algorithm (ACA) [16, 17], are applied to train the output weights of ESN.In this paper, a kind of echo state network (ESN) soft-sensor model for the VCM conversion rate and conversion velocity in the PVC production process based on the biogeographic algorithm is put forward and the simulation results verify the effectiveness of the proposed method. The paper is organized as follows. In Section2, the technique flowchart of the PVC polymerization process is introduced. The data dimension reduction based on KPCA method is presented in Section 3. In Section 4, the ESN soft-sensor model based on BBO algorithm is introduced. The simulation experiments and results analysis are introduced in detail in Section 5. Finally, the conclusion illustrates the last part.


## 2. PVC Polymerization Process

In PVC polymerization process, various raw materials and additives are added to the reaction kettle, which are full evenly dispersed under the mixing action. Then the suitable amounts of the initiators are added to the kettle and start to react. The cooling water is constantly poured into the jacket and baffle of reaction kettle to remove the reaction heat. The reaction will be terminated and the final products are obtained when the conversion ratio of the vinyl chloride (VCM) reaches a certain value and a proper pressure drop appears. Finally, after the reaction completed and VCM contained in slurry separated by the stripping technique, the remaining slurry is fed into the drying process for dewatering and drying. A typical PVC polymerization kettle technological process is shown in Figure1 [2].Figure 1
Technique flowchart of polymerization kettle.According to characteristics of polymerization process, 10 process variables related to VCM conversion rate and conversion velocity are selected as secondary variables of soft-sensor modeling. They, respectively, are kettle temperature (TIC-P101), kettle pressure (PIC-P102), baffle water flow (FIC-P101), jacket water flow (FIC-P102), injection water flow (FIC-P104), seal water flow (FIC-P105), inlet temperature of cooling water (for jacket water and baffles water sharing, TI-P107), outlet temperature of jacket water (TI-P109), outlet temperature of baffle water (TI-P110), and inlet temperature of injection water and seal water (namely, the outlet temperature of the cold water tank, TIC-WA01).


## 3. Data Dimension Reduction Based on KPCA Method

According to the above process characteristics, 10 process variables as auxiliary variables are determined. But for neural network soft-sensor model, too large dimensions of input vector will cause dimension disaster which can lead to the fact that network topology becomes tedious and training is complicated. In this paper, the kernel principal component analysis (KPCA) method is adopted to reduce the dimensions of the input vector. KPCA, introducing the concept of kernel function into principal component analysis (PCA), is a kind of nonlinear extension of PCA. KPCA as a nonlinear form of PCA has the ability to deal with nonlinear problem better than PCA [18, 19]. Data sets composed by input variables of the model are given a kernel principal component analysis and the analysis results are shown in Table 1. The variance contribution rate of the first five principal components has reached more than 95%. Original data primary variables disposed by KPCA are selected as inputs of the neural network model, which not only retains the characteristic information of original variables but also simplifies the structure of neural network.Table 1
Cumulative contribution ratio of different principal component number.
Principal component number
Percentage of variance (%)
Cumulative percentage of variance (%)
1
59.34
59.34
2
16.49
75.83
3
8.35
84.18
4
6.27
90.45
5
4.91
95.36
6
2.52
97.88
7
1.01
98.89
8
0.59
99.48
9
0.52
100.00
10
0
100.00


## 4. ESN Soft-Sensor Model Based on BBO Algorithm

### 4.1. Structure of Soft-Sensor Model

Soft-sensor technology is mainly composed of four parts: data acquisition and processing, choice of auxiliary variables, soft-sensor modeling, and online correction. The framework of the proposed multiple model soft-sensor modeling based on clustering is shown in Figure2.Figure 2
Structure of soft-sensor model.The dimension of the process parameters space is reduced based on KPCA method and the five input variables (kettle pressurex1, baffle water flow x2, jacket water flow x3, outlet temperature of baffle water x4, and outlet temperature of jacket water x5) are selected. The K-means clustering method is adopted to divide the sample data into k classes and each class will be as inputs of a sub-soft-sensor model. Thus the multiple models soft-sensor modeling method based on K-means clustering method is established. The conversion rate and velocity of VCM are the output variables. The BBO algorithm is used to optimize the ESN parameters to realize the nonlinear relationship f{·} between them. Therefore, the soft-sensor model of VCM conversion rate is set up. Experiment proves that this method can effectively improve the prediction accuracy of models.


### 4.2.K-Means Clustering Method

K-means clustering method is a kind of widely used clustering algorithm. Its main thought is that the data set is divided into different clusters through the iteration process, making the criterion function of evaluation clustering performance achieve optimization. Objects within the cluster have a high similarity, and the objects between clusters have a low similarity. The general steps of the algorithm are described as follows.Step 1.
For the sample data setX=xm∣m=1,2,…, randomly select k data samples {c1,c2,…,ck} as the initial cluster centers.Step 2.
Calculate the distanced of each residual sample and the cluster center and assign it to the nearest cluster. The distance between the samples expresses the similarity of the two samples. The smaller the distance, the more similar the two samples, and the smaller the difference degree. The distance between the samples is calculated by Euclidean distance, and its formula is shown as follows:(1)dxi,xj=∑k=1nxik-xjk2,where xi and xj are two samples and n is the sample dimension.Step 3.
Through Step 2, update and get the new clusters. Calculate the square error criterion until the overall average error function is met. Consider(2)E=∑i=1k∑p∈Xip-ci2,where ci is the clustering center and p is sample data.Step 4.
The average value of all objects of clusters is chosen as a new cluster centerci. Repeat Steps 2 and 3 and stop the iteration until the center does not change.According to the chosen number of clusters and cluster center, the data objects are divided intok different clusters. For different clustering data, respectively, establish submodels shown in Figure 3. Through weighted summation of predicted value of each submodel, the final output forecast value is got.Figure 3
Structure of submodel.


### 4.3. Echo State Network

In recent years, neural network has been widely applied in nonlinear systems. The most common types of neural network are the feed-forward neural network and recurrent neural network. Most feed-forward neural networks are static neural network, without the ability of dealing with information dynamically. The recurrent neural network is joined by the dynamic mechanism of processing information on the basis of feed-forward neural network, which makes the whole network with dynamic characteristics and approximates the target value better. Common recurrent neural networks are Elman network and Hopfield network. Echo state network (ESN) is a new type of recurrent neural network, whose typical structure is shown in Figure4.Figure 4
Structure of echo state neural network.It can be seen from Figure4 that the structure of ESN is similar to that of most neural network, which is composed of three parts: input layer, hidden layer, and output layer. Unlike other neural networks, the hidden layer of the ESN is a larger dynamic reservoir (DR) and the number of neurons in the ESN is much more than that of other neural networks. The dynamic reservoirs (DR) can unceasingly store a large number of teachers signals and have short-term memory ability. Although there is no input signal after the training, ESN still can predict for a short period of time; thus this ability can make the network reach the approximation effect for learning system [20].Suppose the input layer containsK neurons, DR contains N neurons, and the output layer contains L neurons. Input sample of network is a K-dimensional vector u(n)=[u1(n),…,uK(n)]T, the state vector of DR is N-dimensional vector x(n)=[x1(n),…,xN(n)]T, and the output sample is L-dimensional vector y(n)=[y1(n),…,yL(n)]T. Between the input neurons and DR a link weight matrix Win exists, whose dimension is N×K; N neurons of DR are connected to a sparse network, and the number of connections is N2 (including the self-connection neurons). The matric W expresses the link weight between DR neurons, and it usually keeps the sparse connection of 1%~5%, so W is the N×N sparse matrix and the element Wij of W expresses the link weight between the i-neuron and the j-neuron in the DR. Between DR and output layer there is an output weight matrix Wout, whose dimension is (K+N+L)×L. In addition, between the output layer and DR there is a feedback connection weight Wfb, and its dimension is N×L.Suppose now we haveM samples (u(i),d(i)) (i=1,2,…,M), where u(i), d(i), respectively, are the i-input and output sample, and its dimensions are, respectively, K and L. The basic equation of echo state network can be represented as follows:(3)xn+1=fWinun+1+Wxn+Wfbdn,(4)yn+1=Woutxn+1;un+1,where f(·) is the activation function of reservoir neuron, generally taking sigmoid function tanh⁡(). It can be seen from (3) x(n+1) is associated with u(n+1), x(n), and d(n). When n=0, there is no definition for d(0), so use d(0)=0 as the output sample. In the training process, Win, Wfb, and W always remain the same and Wout is not involved in network training process, so the value Wout is calculated at the end of the network training. The error of the network forecast output y(n) and actual output of test samples d(n) is smaller; the prediction performance of the network is better.


### 4.4. BBO Algorithm

In recent years, due to the nature’s inspiration, many scholars have proposed some optimization algorithms based on swarm intelligence to solve complex optimization problems, such as genetic algorithm, simulated annealing algorithm, and artificial fish swarm algorithm. Although these intelligent algorithms have appeared for a short time, their good ability to solve complex optimization problems makes them extensively applied into many actual production processes. Biogeography-based optimization algorithm was put forward by an American scholar Simon inspired by the biogeography [21].Biogeography is a kind of natural science that researches species distribution, migration and extinction, and so forth. In nature, the distribution of biology population is different. The place where population lives is named as the “habitat.” And every habitat living environment is not the same, such as rainfall, temperature, humidity, and geology (appropriate index variables, SIV). The habitat suitability index (HSI) is adopted to describe whether the habitat living environment is good or bad. The HSI is higher in which the environment is suitable for the survival of species; on the contrary, it is not suitable for the survival of species. The corresponding relation between biogeography mathematical model and BBO algorithm variables is listed in Table2.Table 2
Corresponding relation between biogeography mathematical model and BBO algorithm variables.
Biogeography mathematical model
Biogeography-based optimization (BBO) algorithm
Habitat
Individual
Habitat suitability index (HSI)
Fitness of the individual
Suitability index variables (SIVs)
Individual variables
High HSI habitats
Excellent individualEvery habitat space is limited and the number of accommodated species is also limited. When the high HSI habitat accommodates more than the largest number of species and the habitat resources are not enough for allocation, competition between species becomes anabatic, making HSI low. At this time, some species will choose to leave the habitat and migrate to a place whose resources are relatively rich, thus improving the HSI lower habitat. So, due to species of the space with high HSI which is multifarious and relatively stable, it has a low immigration rate and high emigration rate. That is to say, low HSI habitats have higher immigration rate and low emigration rate.


#### 4.4.1. Migration Operation

BBO algorithm adopts the migration operation to share information between habitats, making the achievement of the global optimal objective faster. In the high HSI habitats, because of the diversity of species and large number of species, the emigration rate is higher. Due to its high emigration rate, the characteristics information of better habitat is shared to lower HSI habitats. This does not mean that the HSI of better habitat reduces, but its features are copied to the habitats with lower HSI. A species migration model of a single island is shown in Figure5. In order to illustrate the basic principle of biogeography, this model is adopted to describe BBO algorithm.Figure 5
Migration ratio of single island species.In Figure5, Smax⁡ is the largest number of species. Suppose Smax⁡=n; for the kth island the immigration rate is λk, the emigration rate is μk, I is the biggest immigration rate, and E is the biggest emigration rate. When the number of species of an island is zero (S=0), λk=I. As the number of species in the island increases, the number of species in the habitat gradually tends to saturation value; therefor, the immigration rate λk decreases linearly. When the species number reaches maximum Smax⁡, λk=0; that is to say, no species move in. For the emigration rate μk, when S=0, no species move out; namely, μk=0. With the increase of species number, in order to find more suitable survival habitats, the emigration rate increases. When the species number reaches maximum, μk=E. When λk=μk, the species number of the habitat reaches an equilibrium state. According to Figure 5, the immigration rate and emigration rate can be calculated as follows:(5)λk=I×1-kn,μk=E×kn.


#### 4.4.2. Mutation Operation

As species migrate between habitats, the number of species in each island is constantly changing. Suppose the probability of habitats including species numberS is Ps; the change value formula of the probability Ps at time t to t+Δt [21] is described as follows:(6)pst+Δt=1-λs-μsPsΔt+λs-1Ps-1Δt+μs+1Ps+1Δt,where λs and μs express immigration rate and emigration rate when the habitat contains species S.Assuming thatΔt is small enough, the calculation of probability Ps is shown as the following formula:(7)Ps=-λs+μsPs+μs+1Ps+1,S=0,-λs+μsPs+λs-1Ps-1+μs+1Ps+1,1≤S<Smax⁡-1,-λs+μsPs+λs-1Ps-1,S=Smax⁡.The above formula can be abbreviated toP=AP. When the largest species number of the habitat is equal to Smax⁡, the probability Ps function corresponding to different habitats is a symmetric function about equilibrium point. Individuals with larger species number and less species number all have low stable probability; that is to say, the probability is small, the number of species near the equilibrium point is relatively stable, and the existence probability is higher. On this basis, the variation rate mi is designed as the following equation:(8)A=-λ0+μ0μ1λ0-λ1+μ1μ2λ1⋯⋯⋯-λn-1+μn-1μn-1λn-1-λn+μn,mi=mmax⁡1-PSiPmax⁡,where mmax⁡ is the biggest mutation rate and Ps is the probability of habitat accommodating species S.For the mutation operation, whether the habitat needs to mutate first should be determined. If the random number is less than the mutation probabilitymi, it means the habitat needs to mutate. Then a group of randomly generated vectors are used to replace the original vector. In nature incidents (such as volcanic eruptions, tsunamis, and disease) are often unavoidable; the occurrence of these events will affect the number of species, makes the ecological environment unstable, and reduces the habitat suitability index. If islands of low suitability index are given a variation, the chance of getting a better solution will increase; if islands of higher suitability index are given a variation, it may not get a better solution, so retain islands with high suitability index and make a mutation for islands with low suitability index.


### 4.5. Algorithm Procedure

First, main parameters of ESN are the input matrixWin, the reservoir weight matrix W, the output feedback weight matrix Wfb, and the output weight matrix Wout. So optimizing the ESN is equivalent to optimizing the four matrixes. In the network learning process, Win, Wfb, and W always remain the same, and Wout is not involved in the network training process and its value is calculated after the end of the network training. So in this paper, the habitat of biogeographic optimization algorithm is in correspondence with the output connection weight of the ESN. Through biogeography-based optimization algorithm, the weight of ESN is optimized to realize the ideal predicted values. The algorithm flowchart of ESN prediction model optimized by BBO algorithm is shown in Figure 6.Figure 6
Algorithm flowchart of ESN prediction model optimized by BBO algorithm.Step 1 (initialization parameters).
Initialize the following algorithm parameters: the largest species number of islandSmax⁡, the number of island N, the emigration rate λk and immigration rate μk, the maximum variation rate mmax⁡ and the number of iterations itermax⁡. Initialize a group of islands, each island, namely, each habitat, all represents an individual, which is the solution of the problem.Step 2 (calculate the fitness value).
Use suitability index HSI of the island as the fitness value. Calculate the fitness value of each islandfitness function. Judge whether the termination condition is met or not. If satisfied, output the optimal solution, otherwise continue Step 3.Step 3.
According to the fitness values, the individuals are arranged in a descending order and the highest HSI individual is stored. Calculate species numberS corresponded by each island, the emigration rate λk, the immigration rate μk, and the mutation probability mmax⁡. Then the optimal individual is noted as p.Step 4 (migration operation).
According to the emigration rateλk and the immigration rate μk, judge whether the islands k need to perform the migration operation or not. If needed, perform the migration operations, generate new individuals, and recalculate the fitness value of island k. Compare it with the optimal individual p and retain the best individual, noted as p′.Step 5 (mutation operation).
Make a mutation operation for islands with lower HSI and recalculate the fitness value of each island. Record the best individual now. Compare with the optimal individualp′ and keep the optimal individual noted as p′′.Step 6.
Check the same individual and use random vectors instead of the same individuals. Reorder all individuals and maintain the individual corresponded by the highest HSI. Record the best individual at this time. If the terminating condition is satisfied, exit iteration. Otherwise, repeat from Step 3.Step 7 (End).
Return the vector corresponding to the highest HSI individual.


## 4.1. Structure of Soft-Sensor Model

Soft-sensor technology is mainly composed of four parts: data acquisition and processing, choice of auxiliary variables, soft-sensor modeling, and online correction. The framework of the proposed multiple model soft-sensor modeling based on clustering is shown in Figure2.Figure 2
Structure of soft-sensor model.The dimension of the process parameters space is reduced based on KPCA method and the five input variables (kettle pressurex1, baffle water flow x2, jacket water flow x3, outlet temperature of baffle water x4, and outlet temperature of jacket water x5) are selected. The K-means clustering method is adopted to divide the sample data into k classes and each class will be as inputs of a sub-soft-sensor model. Thus the multiple models soft-sensor modeling method based on K-means clustering method is established. The conversion rate and velocity of VCM are the output variables. The BBO algorithm is used to optimize the ESN parameters to realize the nonlinear relationship f{·} between them. Therefore, the soft-sensor model of VCM conversion rate is set up. Experiment proves that this method can effectively improve the prediction accuracy of models.


## 4.2.K-Means Clustering Method

K-means clustering method is a kind of widely used clustering algorithm. Its main thought is that the data set is divided into different clusters through the iteration process, making the criterion function of evaluation clustering performance achieve optimization. Objects within the cluster have a high similarity, and the objects between clusters have a low similarity. The general steps of the algorithm are described as follows.Step 1.
For the sample data setX=xm∣m=1,2,…, randomly select k data samples {c1,c2,…,ck} as the initial cluster centers.Step 2.
Calculate the distanced of each residual sample and the cluster center and assign it to the nearest cluster. The distance between the samples expresses the similarity of the two samples. The smaller the distance, the more similar the two samples, and the smaller the difference degree. The distance between the samples is calculated by Euclidean distance, and its formula is shown as follows:(1)dxi,xj=∑k=1nxik-xjk2,where xi and xj are two samples and n is the sample dimension.Step 3.
Through Step 2, update and get the new clusters. Calculate the square error criterion until the overall average error function is met. Consider(2)E=∑i=1k∑p∈Xip-ci2,where ci is the clustering center and p is sample data.Step 4.
The average value of all objects of clusters is chosen as a new cluster centerci. Repeat Steps 2 and 3 and stop the iteration until the center does not change.According to the chosen number of clusters and cluster center, the data objects are divided intok different clusters. For different clustering data, respectively, establish submodels shown in Figure 3. Through weighted summation of predicted value of each submodel, the final output forecast value is got.Figure 3
Structure of submodel.


## 4.3. Echo State Network

In recent years, neural network has been widely applied in nonlinear systems. The most common types of neural network are the feed-forward neural network and recurrent neural network. Most feed-forward neural networks are static neural network, without the ability of dealing with information dynamically. The recurrent neural network is joined by the dynamic mechanism of processing information on the basis of feed-forward neural network, which makes the whole network with dynamic characteristics and approximates the target value better. Common recurrent neural networks are Elman network and Hopfield network. Echo state network (ESN) is a new type of recurrent neural network, whose typical structure is shown in Figure4.Figure 4
Structure of echo state neural network.It can be seen from Figure4 that the structure of ESN is similar to that of most neural network, which is composed of three parts: input layer, hidden layer, and output layer. Unlike other neural networks, the hidden layer of the ESN is a larger dynamic reservoir (DR) and the number of neurons in the ESN is much more than that of other neural networks. The dynamic reservoirs (DR) can unceasingly store a large number of teachers signals and have short-term memory ability. Although there is no input signal after the training, ESN still can predict for a short period of time; thus this ability can make the network reach the approximation effect for learning system [20].Suppose the input layer containsK neurons, DR contains N neurons, and the output layer contains L neurons. Input sample of network is a K-dimensional vector u(n)=[u1(n),…,uK(n)]T, the state vector of DR is N-dimensional vector x(n)=[x1(n),…,xN(n)]T, and the output sample is L-dimensional vector y(n)=[y1(n),…,yL(n)]T. Between the input neurons and DR a link weight matrix Win exists, whose dimension is N×K; N neurons of DR are connected to a sparse network, and the number of connections is N2 (including the self-connection neurons). The matric W expresses the link weight between DR neurons, and it usually keeps the sparse connection of 1%~5%, so W is the N×N sparse matrix and the element Wij of W expresses the link weight between the i-neuron and the j-neuron in the DR. Between DR and output layer there is an output weight matrix Wout, whose dimension is (K+N+L)×L. In addition, between the output layer and DR there is a feedback connection weight Wfb, and its dimension is N×L.Suppose now we haveM samples (u(i),d(i)) (i=1,2,…,M), where u(i), d(i), respectively, are the i-input and output sample, and its dimensions are, respectively, K and L. The basic equation of echo state network can be represented as follows:(3)xn+1=fWinun+1+Wxn+Wfbdn,(4)yn+1=Woutxn+1;un+1,where f(·) is the activation function of reservoir neuron, generally taking sigmoid function tanh⁡(). It can be seen from (3) x(n+1) is associated with u(n+1), x(n), and d(n). When n=0, there is no definition for d(0), so use d(0)=0 as the output sample. In the training process, Win, Wfb, and W always remain the same and Wout is not involved in network training process, so the value Wout is calculated at the end of the network training. The error of the network forecast output y(n) and actual output of test samples d(n) is smaller; the prediction performance of the network is better.


## 4.4. BBO Algorithm

In recent years, due to the nature’s inspiration, many scholars have proposed some optimization algorithms based on swarm intelligence to solve complex optimization problems, such as genetic algorithm, simulated annealing algorithm, and artificial fish swarm algorithm. Although these intelligent algorithms have appeared for a short time, their good ability to solve complex optimization problems makes them extensively applied into many actual production processes. Biogeography-based optimization algorithm was put forward by an American scholar Simon inspired by the biogeography [21].Biogeography is a kind of natural science that researches species distribution, migration and extinction, and so forth. In nature, the distribution of biology population is different. The place where population lives is named as the “habitat.” And every habitat living environment is not the same, such as rainfall, temperature, humidity, and geology (appropriate index variables, SIV). The habitat suitability index (HSI) is adopted to describe whether the habitat living environment is good or bad. The HSI is higher in which the environment is suitable for the survival of species; on the contrary, it is not suitable for the survival of species. The corresponding relation between biogeography mathematical model and BBO algorithm variables is listed in Table2.Table 2
Corresponding relation between biogeography mathematical model and BBO algorithm variables.
Biogeography mathematical model
Biogeography-based optimization (BBO) algorithm
Habitat
Individual
Habitat suitability index (HSI)
Fitness of the individual
Suitability index variables (SIVs)
Individual variables
High HSI habitats
Excellent individualEvery habitat space is limited and the number of accommodated species is also limited. When the high HSI habitat accommodates more than the largest number of species and the habitat resources are not enough for allocation, competition between species becomes anabatic, making HSI low. At this time, some species will choose to leave the habitat and migrate to a place whose resources are relatively rich, thus improving the HSI lower habitat. So, due to species of the space with high HSI which is multifarious and relatively stable, it has a low immigration rate and high emigration rate. That is to say, low HSI habitats have higher immigration rate and low emigration rate.


### 4.4.1. Migration Operation

BBO algorithm adopts the migration operation to share information between habitats, making the achievement of the global optimal objective faster. In the high HSI habitats, because of the diversity of species and large number of species, the emigration rate is higher. Due to its high emigration rate, the characteristics information of better habitat is shared to lower HSI habitats. This does not mean that the HSI of better habitat reduces, but its features are copied to the habitats with lower HSI. A species migration model of a single island is shown in Figure5. In order to illustrate the basic principle of biogeography, this model is adopted to describe BBO algorithm.Figure 5
Migration ratio of single island species.In Figure5, Smax⁡ is the largest number of species. Suppose Smax⁡=n; for the kth island the immigration rate is λk, the emigration rate is μk, I is the biggest immigration rate, and E is the biggest emigration rate. When the number of species of an island is zero (S=0), λk=I. As the number of species in the island increases, the number of species in the habitat gradually tends to saturation value; therefor, the immigration rate λk decreases linearly. When the species number reaches maximum Smax⁡, λk=0; that is to say, no species move in. For the emigration rate μk, when S=0, no species move out; namely, μk=0. With the increase of species number, in order to find more suitable survival habitats, the emigration rate increases. When the species number reaches maximum, μk=E. When λk=μk, the species number of the habitat reaches an equilibrium state. According to Figure 5, the immigration rate and emigration rate can be calculated as follows:(5)λk=I×1-kn,μk=E×kn.


### 4.4.2. Mutation Operation

As species migrate between habitats, the number of species in each island is constantly changing. Suppose the probability of habitats including species numberS is Ps; the change value formula of the probability Ps at time t to t+Δt [21] is described as follows:(6)pst+Δt=1-λs-μsPsΔt+λs-1Ps-1Δt+μs+1Ps+1Δt,where λs and μs express immigration rate and emigration rate when the habitat contains species S.Assuming thatΔt is small enough, the calculation of probability Ps is shown as the following formula:(7)Ps=-λs+μsPs+μs+1Ps+1,S=0,-λs+μsPs+λs-1Ps-1+μs+1Ps+1,1≤S<Smax⁡-1,-λs+μsPs+λs-1Ps-1,S=Smax⁡.The above formula can be abbreviated toP=AP. When the largest species number of the habitat is equal to Smax⁡, the probability Ps function corresponding to different habitats is a symmetric function about equilibrium point. Individuals with larger species number and less species number all have low stable probability; that is to say, the probability is small, the number of species near the equilibrium point is relatively stable, and the existence probability is higher. On this basis, the variation rate mi is designed as the following equation:(8)A=-λ0+μ0μ1λ0-λ1+μ1μ2λ1⋯⋯⋯-λn-1+μn-1μn-1λn-1-λn+μn,mi=mmax⁡1-PSiPmax⁡,where mmax⁡ is the biggest mutation rate and Ps is the probability of habitat accommodating species S.For the mutation operation, whether the habitat needs to mutate first should be determined. If the random number is less than the mutation probabilitymi, it means the habitat needs to mutate. Then a group of randomly generated vectors are used to replace the original vector. In nature incidents (such as volcanic eruptions, tsunamis, and disease) are often unavoidable; the occurrence of these events will affect the number of species, makes the ecological environment unstable, and reduces the habitat suitability index. If islands of low suitability index are given a variation, the chance of getting a better solution will increase; if islands of higher suitability index are given a variation, it may not get a better solution, so retain islands with high suitability index and make a mutation for islands with low suitability index.


## 4.4.1. Migration Operation

BBO algorithm adopts the migration operation to share information between habitats, making the achievement of the global optimal objective faster. In the high HSI habitats, because of the diversity of species and large number of species, the emigration rate is higher. Due to its high emigration rate, the characteristics information of better habitat is shared to lower HSI habitats. This does not mean that the HSI of better habitat reduces, but its features are copied to the habitats with lower HSI. A species migration model of a single island is shown in Figure5. In order to illustrate the basic principle of biogeography, this model is adopted to describe BBO algorithm.Figure 5
Migration ratio of single island species.In Figure5, Smax⁡ is the largest number of species. Suppose Smax⁡=n; for the kth island the immigration rate is λk, the emigration rate is μk, I is the biggest immigration rate, and E is the biggest emigration rate. When the number of species of an island is zero (S=0), λk=I. As the number of species in the island increases, the number of species in the habitat gradually tends to saturation value; therefor, the immigration rate λk decreases linearly. When the species number reaches maximum Smax⁡, λk=0; that is to say, no species move in. For the emigration rate μk, when S=0, no species move out; namely, μk=0. With the increase of species number, in order to find more suitable survival habitats, the emigration rate increases. When the species number reaches maximum, μk=E. When λk=μk, the species number of the habitat reaches an equilibrium state. According to Figure 5, the immigration rate and emigration rate can be calculated as follows:(5)λk=I×1-kn,μk=E×kn.


## 4.4.2. Mutation Operation

As species migrate between habitats, the number of species in each island is constantly changing. Suppose the probability of habitats including species numberS is Ps; the change value formula of the probability Ps at time t to t+Δt [21] is described as follows:(6)pst+Δt=1-λs-μsPsΔt+λs-1Ps-1Δt+μs+1Ps+1Δt,where λs and μs express immigration rate and emigration rate when the habitat contains species S.Assuming thatΔt is small enough, the calculation of probability Ps is shown as the following formula:(7)Ps=-λs+μsPs+μs+1Ps+1,S=0,-λs+μsPs+λs-1Ps-1+μs+1Ps+1,1≤S<Smax⁡-1,-λs+μsPs+λs-1Ps-1,S=Smax⁡.The above formula can be abbreviated toP=AP. When the largest species number of the habitat is equal to Smax⁡, the probability Ps function corresponding to different habitats is a symmetric function about equilibrium point. Individuals with larger species number and less species number all have low stable probability; that is to say, the probability is small, the number of species near the equilibrium point is relatively stable, and the existence probability is higher. On this basis, the variation rate mi is designed as the following equation:(8)A=-λ0+μ0μ1λ0-λ1+μ1μ2λ1⋯⋯⋯-λn-1+μn-1μn-1λn-1-λn+μn,mi=mmax⁡1-PSiPmax⁡,where mmax⁡ is the biggest mutation rate and Ps is the probability of habitat accommodating species S.For the mutation operation, whether the habitat needs to mutate first should be determined. If the random number is less than the mutation probabilitymi, it means the habitat needs to mutate. Then a group of randomly generated vectors are used to replace the original vector. In nature incidents (such as volcanic eruptions, tsunamis, and disease) are often unavoidable; the occurrence of these events will affect the number of species, makes the ecological environment unstable, and reduces the habitat suitability index. If islands of low suitability index are given a variation, the chance of getting a better solution will increase; if islands of higher suitability index are given a variation, it may not get a better solution, so retain islands with high suitability index and make a mutation for islands with low suitability index.


## 4.5. Algorithm Procedure

First, main parameters of ESN are the input matrixWin, the reservoir weight matrix W, the output feedback weight matrix Wfb, and the output weight matrix Wout. So optimizing the ESN is equivalent to optimizing the four matrixes. In the network learning process, Win, Wfb, and W always remain the same, and Wout is not involved in the network training process and its value is calculated after the end of the network training. So in this paper, the habitat of biogeographic optimization algorithm is in correspondence with the output connection weight of the ESN. Through biogeography-based optimization algorithm, the weight of ESN is optimized to realize the ideal predicted values. The algorithm flowchart of ESN prediction model optimized by BBO algorithm is shown in Figure 6.Figure 6
Algorithm flowchart of ESN prediction model optimized by BBO algorithm.Step 1 (initialization parameters).
Initialize the following algorithm parameters: the largest species number of islandSmax⁡, the number of island N, the emigration rate λk and immigration rate μk, the maximum variation rate mmax⁡ and the number of iterations itermax⁡. Initialize a group of islands, each island, namely, each habitat, all represents an individual, which is the solution of the problem.Step 2 (calculate the fitness value).
Use suitability index HSI of the island as the fitness value. Calculate the fitness value of each islandfitness function. Judge whether the termination condition is met or not. If satisfied, output the optimal solution, otherwise continue Step 3.Step 3.
According to the fitness values, the individuals are arranged in a descending order and the highest HSI individual is stored. Calculate species numberS corresponded by each island, the emigration rate λk, the immigration rate μk, and the mutation probability mmax⁡. Then the optimal individual is noted as p.Step 4 (migration operation).
According to the emigration rateλk and the immigration rate μk, judge whether the islands k need to perform the migration operation or not. If needed, perform the migration operations, generate new individuals, and recalculate the fitness value of island k. Compare it with the optimal individual p and retain the best individual, noted as p′.Step 5 (mutation operation).
Make a mutation operation for islands with lower HSI and recalculate the fitness value of each island. Record the best individual now. Compare with the optimal individualp′ and keep the optimal individual noted as p′′.Step 6.
Check the same individual and use random vectors instead of the same individuals. Reorder all individuals and maintain the individual corresponded by the highest HSI. Record the best individual at this time. If the terminating condition is satisfied, exit iteration. Otherwise, repeat from Step 3.Step 7 (End).
Return the vector corresponding to the highest HSI individual.


## 5. Simulation Results

In this paper, the polymerization industrial process of a chemical factory with 40000 tons/year polyvinyl chloride (PVC) production device is taken as background, whose technology is introduced by America B·F·G company, taking vinyl chloride monomer (VCM) as raw material and using suspension polymerization technology to produce polyvinyl chloride (PVC) resin. A soft-sensor model of the VCM conversion rate and velocity in the polyvinyl chloride (PVC) production process based on BBOA-ESN algorithm is put forward. After reducing ESN dimensions, the number of input variables is 5 and the output dimension is 2. In addition, suppose the reservoir size is 100, the sparse connection rate of reservoir weight matrix is 5%, the activation function of reservoir istanh⁡(), and the output unit adopts linear activation function. Suppose the initial values of ESN parameters: Win=0.3, W=0.2, and Win=0.03.The initialized parameters of adopted BBO algorithm are described as follows: the habitat sizeN=100, the largest species number Smax⁡=100, the largest immigration rate I=1, the largest emigration rate E=1, the largest mutation probability mmax⁡=0.005, and the maximum iterations number is 200.Before setting up the soft-sensor model of VCM conversion rate and velocity in the PVC polymerization process, in order to measure the performances of prediction models, several performance indicators are defined in Table3, where y(t) is predicted value and yd(t) is actual value.Table 3
Definition of model performance index.
Mean absolute percentage error
MAPE
=
100
T
∑
t
=
1
T
y
t
-
y
d
(
t
)
y
d
(
t
)
Mean-square error
MSE
=
1
T
∑
t
=
1
T
y
t
-
y
d
t
2
Root mean square error
RMSE
=
1
T
∑
t
=
1
n
y
d
t
-
y
t
2
1
/
2
Normalized root mean square error
NRMSE
=
1
T
y
d
2
∑
t
=
1
T
y
t
-
y
d
t
2
Sum of squared error
SSE
=
∑
t
=
1
n
y
d
t
-
y
(
t
)
2The production historical data of PVC polymerization process are collected and 2 kettles including 1600 sets historical data with the uniformity and representativeness are chosen. Then after data preprocessing the data is divided into two parts, in which the front 1350 data are the training data, and the rest 250 data are used to validate the performance of soft-sensor models. The simulation results are shown in Figures7–10. Figure 7 shows the output comparison curves of conversion velocity of VCM, respectively, predicted by ESN soft-sensor model, BBO-ESN soft-sensor model, and BBO-MESN soft-sensor model. Figure 8 shows the predicted error curves of VCM conversion velocity. Figure 9 shows the output comparison curves of VCM conversion rate, respectively, predicted by ESN soft-sensor model, BBO-ESN soft-sensor model, and BBO-MESN soft-sensor model. Figure 10 shows the predicted error curves of VCM conversion rate.Figure 7
Predicted output curves of VCM conversion velocity.Figure 8
Predicted output curves of VCM conversion velocity.Figure 9
Predicted output curves of VCM conversion ratio.Figure 10
Predicted error curves of VCM conversion ratio.The model performances comparison is shown in Table4. From the simulation results, the prediction accuracy of BBO-MESN soft-sensor model proposed in this paper is higher than that of ESN and BBO-ESN soft-sensor model. Its application of predicting the VCM conversion rate and velocity in PVC polymerization process has great significance in improving the capacity of equipment and reducing the production cost.Table 4
Performances comparison of different soft-sensor models.
MAPE
MSE
RMSE
NMSE
SSE
ESN
5.3509
0.0078
0.0880
0.0201
1.9381
BBO-ESN
5.0907
0.0073
0.0856
0.0044
1.8322
BBO-MESN
4.1694
0.0051
0.0713
0.0036
1.2692


## 6. Conclusions

Based on that ESN has good capability of nonlinear approximation and biogeography-based optimization algorithm (BBOA) is simple and easy-to-implement, which can obtain the global extreme and avoid falling into local extreme, a BBO-MESN soft-sensor model is proposed to predict VCM conversion rate and conversion velocity. BBOA is used to optimize the output weights of the ESN network. The simulation results show that the neural network soft-sensor model based on BBO-MESN has higher prediction accuracy.



---

*Source: 493248-2015-04-20.xml*