Burden Surface Decision Using MODE With TOPSIS in Blast Furnace Ironmkaing

Burden surface distribution plays a key role in achieving an energy-efficient status of blast furnace (BF). However, actual adjustment of burden surface usually depends on the operator’s experience when the production status changes. Meanwhile, due to the characteristics of high dimension, strong coupling, and distributed parameters, it is difficult to establish the accurate mechanism model for BF ironmaking process. Considering the aforementioned issues, this paper proposes an integrated multi-objective optimization framework for optimizing burden surface distribution based on the analysis of BF operation characteristics. Firstly, data-driven models are constructed for two objectives, i.e., gas utilization ratio (GUR) and coke ratio (CR), and two constraints using adaptive particle swarm optimization (APSO) based extreme learning machine (ELM), named APSO-ELM. Multi-objective optimization is subsequently carried out between GUR and CR using the multi-objective differential evolution algorithm (MODE) to generate the Pareto optimal solutions. Finally, TOPSIS is applied to select a best compromise solution among the Pareto optimal solutions for this optimization problem. Comprehensive experiments are presented to illustrate the performance of the proposed integrated multi-objective optimization framework. The experimental results demonstrate that the proposed framework can give a reasonable burden surface profile according to the production status changes to guarantee the BF operation more efficient and stable.

Hence, the optimal setting of burden surface is of vital to stable operation and improvement of energy efficiency, which is still an open problem both in academic research and industrial production.
In the past decades, many research works paid much attention to the burden distribution, especially in the mathematical model for the burden distribution process. Park et al. [8], [9] developed a burden distribution analysis model based on burden trajectory, descent model, and stock model to calculate the burden surface. Fu et al. [10] demonstrated a mathematical model with the combination of falling curve model, stock line profile formation model and burden descending model, and discussed the effects of the non-uniform descending velocity on burden distribution. Zhao et al. [11] established a comprehensive model from flow control gate to stock surface in detail and analyzed the non-uniformity phenomenon. Xu et al. [12] proposed the circumferential burden distribution behaviors at bell-less top blast furnace with parallel type hoppers to reduce the uneven degree of burden distribution. The aforementioned contributions make us understand the ''black box'' system more clearly, and provide theoretical support for decision-making. In terms of the optimal burden surface decision, it is basically given by the results of experience and empirical knowledge. Generally speaking, platform plus funnel mode is recognized as the ideal burden surface profile, which can satisfy the demand of gas flow distribution to achieve the improvement of gas utilization efficiency [13]- [15]. To our best knowledge, researchers rarely attempt to find the optimal burden surface matching the production status to achieve optimal key production indicators. Uniquely, Li et al. [16], [17] constructed multiple models set of burden surface by k-means clustering algorithm and established a selection mechanism, but they did not consider the effects of production status changes on burden surface and this approach cannot guarantee the optimization of key production indicators.
The purpose of BF ironmaking operation is high efficiency, high quality and low consumption, which are interrelated and differentiated [18]. Actually, the burden surface optimal setting can be treated as a multi-objective optimization problem (MOP) under the premise that the burden surface is discretized. Meanwhile, the methods to solve the MOPs have made considerable progresses in recent years. Many multi-objective evolutionary algorithms (MOEAs) have been introduced to solve MOPs, which can be divided into two categories. The first category, including the non-dominated sorting genetic algorithm (NSGA) [19], does not provide an elitism mechanism. The second category, known as NSGA-II [20], gains much attention from the researchers due to their effectiveness and easy implementation. In addition, there are several multi-objective metaheuristic algorithms, such as multi-objective ant colony optimization (MOACO) [21] and multi-objective particle swarm optimization (MOPSO) [22]. A review of multi-objective metaheuristic algorithms has been presented by Jones et al. [23]. Since the development of MOEAs, there has been a growing interest in obtaining the Pareto optimal solutions using different evolutionary algorithms. Among these attended MOEAs, multi-objective differential evolution algorithm (MODE) [24] preforms better in solving the MOP. In the past decade, MODE has been improved and widely applied in optimization of various problems [25], [26]. These algorithms lay a solid theoretical foundation for the burden surface optimization problem.
Accordingly, in this paper, we focus on the optimal setting of burden surface. Due to the complexity of BF and the installation of radars in the top of BF, a multi-objective optimization framework for burden surface is established based on data-driven technique. The original burden surface data collected by radars are analyzed, and the characteristics of the burden surface are carefully considered to be discretized according to the experience of operator. After data processing, considering the constrains of the actual situation, and upper and lower bounds of the variables, data-driven multiobjective optimization models are constructed for optimizing burden surface with the aim to minimize the difference between actual values and target values of gas utilization ratio (GUR) and coke ratio (CR). We combine extreme learning machine (ELM), which is a new kind of machine learning approach with fast learning speed and excellent generalization performance, and adaptive particle swarm optimization (APSO) to establish the accurate data-driven models for objectives and constraints, named APSO-ELM, in which APSO is implemented to help determine appropriate hidden layer parameters in ELM. On this basis, MODE algorithm is used to optimize the goals to obtain the Pareto optimal solutions. Then, a popular multiple criteria decision making method called TOPSIS is adopted to rank the Pareto optimal solutions to select a best compromise solution according to production needs.
The main contributions of this paper can be summarized as follows: (1) At present, the optimal setting of burden surface has been seldom investigated. In terms of the production indicator optimization, the modeling for optimal burden surface is considered as a MOP during ironmaking process.
(2) An integrated optimization framework is proposed to solve the MOP for burden surface. The modified APSO-ELM is applied to establish the process models between production indicators and decision variables. The inputs of these process models are the characteristic parameters of burden surface and the operating status parameters. Next, MODE is utilized to generate the Pareto optimal solutions, and then a best compromise solution can be obtained using TOPSIS for this optimization problem.
(3) Comprehensive experiments have been carried out to validate the effectiveness of the proposed optimization framework using actual data collected from a BF.
The remaining parts of this paper are arranged as follows: Section II introduces the multi-objective problem for burden surface based on the analysis of BF operation characteristics. Multi-objective optimization strategy of burden surface is detailed in Section III, including data-driven models, MODE solved Pareto optimal solutions, and the application of TOPSIS for ranking the solutions. Section IV presents the simulation results to verify the applicability and effectiveness of the proposed multi-objective optimization framework for burden surface. Conclusions are given in Section V.

II. PROBLEM FORMULATIONS
Creating a reasonable burden distribution in a BF is a highly complicated task in reality, even though it is of utmost importance and absolutely essential for the smooth and efficient operation of the furnace [27]. In this section, the BF ironmaking process and burden distribution mechanism are firstly described. In addition, the actual adjustment status of burden surface is analyzed. Then, according to the characteristics of BF ironmaking and smelting experience, optimization objectives and constraint conditions are determined, and the mathematical model of optimization problem for burden surface is constructed, respectively.

A. BF IRONMAKING PROCESS
As shown in Fig. 1, BF is a giant shaftlike countercurrent reactor used for smelting to produce molten iron. The solid reactants, including iron ore and coke, are fed into the top layer by layer iteratively, while the preheated air and some auxiliary fuels are blown into the bottom through the tuyeres. During the ironmaking process, complex gas-solid, gasliquid, and solid-solid chemical reactions occur in different zones under different temperatures, which are shown in the right low corner of Fig. 1. The radial distribution of the charged solid raw materials, i.e., burden surface (see the right top corner of Fig. 1), influences the pressure loss and the local mass flows of solid and gas inside the furnace, and further affects the indirect reduction degree of the ore [15]. In addition, the burden surface is closely related to operating status.
The actual adjustment status of burden surface in BF ironmaking process is presented in Fig. 2. According to the current production status, the upper limit y u , lower limit y l and detection value y (t) of main production indicators and production boundary conditions B, the operator gives the burden surface X sp with experience. In practice, the operator cannot accurately adjust the burden surface when the production status changes in most situations. Therefore, an appropriate adjustment of the burden surface is required for smooth BF operation, energy efficiency and quality improvement of molten iron according to the production status changes.

B. MAIN PRODUCTION INDICATORS OF BF IRONMAKING PROCESS
The production indicators can reflect the operation level of production comprehensively. In terms of the optimization problem of burden surface, only the suitable indicators can be optimized to achieve a satisfactory effect.

1) ENERGY UTILIZATION RATIO
In the actual production, GUR is the main indicator for energy utilization ratio. It represents the ratio of CO to CO 2 . The improvement of GUR is an important embodiment of technical progress in BF operation [28], [29]. The corresponding mathematical description is as follows: where, η CO represents the GUR, CO 2 (%) and CO(%) are the content of CO 2 and CO, respectively. During the ironmaking process, CO is the main reactant to reduce iron ore, the utilization degree directly affects the process of chemical reactions inside the furnace as described in Fig. 1, leading to the rise of product yield. The burden surface and the operating status have impacts on GUR. If the burden surface profile is reasonable and the chemical reactions are sufficient, the gas utilization degree will become high. In addition, the operating status parameters, such as blast temperature, blast volume, blast pressure, top temperature, top pressure, and oxygen enrichment, are also closely related to GUR.

2) ECONOMIC COST
Besides energy utilization ratio, economic cost is another important target in BF ironmaking process. CR is the main indicator for economic cost, which represents the amount of coke consumed by smelting one ton of qualified pig iron. In order to reduce the production cost, it is urgent to reduce the use of coke. Similarly, the burden surface and the operating status also have effects on CR. In the practical ironmaking process, energy utilization ratio and economic cost are two coupling and contradictory objectives.

3) PERMEABILITY INDEX
Permeability index (PI) is a significant symbol to measure the smooth operation of BF, which can be defined as the ratio of blast volume and pressure difference. The corresponding mathematical representation can be shown as where BV is the blast volume, P is the pressure difference between the blast pressure and the top pressure. In normal production, PI is stable, which means that it achieves a dynamic balance between ascending hot gas and descending charge material. If the gas permeability inside the furnace becomes worse, it will cause the charge material fall to be difficult. In severe case, it will lead to hanging. On the other hand, increased PI indicates that the reaction between hot gas and charge material is not enough. Generally speaking, in order to ensure the BF running smoothly, PI should be kept within 30 ∼ 34.

4) QUALITY INDEX
As mentioned previously, the basic task of BF production is to smelt iron ore into qualified hot metal. Hot metal silicon content (HMSC) is a main parameter by which product quality of pig iron is measured [30]. It actually refers to the percentage of [Si] elements in hot metal and is generally controlled between 0.5% ∼ 0.7% to satisfy production demand.

C. BURDEN SURFACE ADJUSTMENT AS A MULTI-OBJECTIVE OPTIMIZATION PROBLEM
The purpose of this paper is to determine the reasonable burden surface according to the production status with highest yield and lowest consumption among those satisfying all constraints. Based on the aforementioned analysis, GUR and CR are the two most important production indicators for BF ironmaking process. Meanwhile, through analyzing the monthly statistics of indexes, the relationship between GUR and CR are depicted in Fig. 3. As observed from Fig. 3, the GUR-CR relationship can be illustrated clearly and they cannot be optimal at the same time. In practice, they need to be maintained at their desired ranges, and their deviation between real and target values should be less than a required value.
With the development of automatic detection technique for burden surface, radar-type instrument has been widely used to obtain the burden surface data to establish burden surface distribution model [31], [32]. In order to better characterize the burden surface, feature extraction is performed for burden surface discretization. In the charging operation, the width  of platform, the depth of funnel and the width of funnel are the main features of burden surface. Hence, 7 features are extracted to represent the burden surface according to operator's cognition and charging operation, including the width of funnel l 1 , the width of platform l 2 , the distance between zero position and burden surface h 1 , the depth of funnel h 2 , the inclination angle of funnel α, the central angle of funnel β and the inclination angle of edge γ , as depicted in Fig. 4. Therefore, 7 characteristic parameters are used as decision variables in the following optimization problem.
According to the above analysis, more specifically, the burden surface optimization task can be summarized as • Satisfaction of constraint on PI and HMSC, keeping it within prescribed bounds.
Due to the high complexity of the furnace interior, the mechanism model of the iromaking process is difficult to be established. Fortunately, data-driven intelligent modeling approach can effectively deal with modeling problem for complex industrial processes. In this paper, we try to construct the data-driven process models for the two objectives as well as the two constraints as where y GUR , y CR , y PI and y HMSC are GUR, CR, PI and HMSC, respectively. X = (l 1 , l 2 , h 1 , h 2 , α, β, γ ) is the burden surface features, and S is the operating status parameters. VOLUME 8, 2020 Accordingly, burden surface optimization can be regarded as a typical MOP. The corresponding mathematical description is as follows: where i = 1, 2, . . . , 7, x b i and x u i are the boundary constraint and upper bounds of x i , y l and y u are the lower and upper bounds of y . In the constraints of Eq.(4), the first inequality is to constrain PI to ensure the smooth and stable production, the second inequality is to constrain HMSC to ensure the production of qualified hot metal, and the third constraint guarantees the 7 burden surface decision variables in the given ranges. It should be noted that x b i is the burden surface features of the last time charging operation, which is to ensure that the optimized burden surface is higher than that before optimization.
In order to solve the MOP for burden surface, an integrated framework is proposed, which is illustrated in Fig. 5. Firstly, this MOP is descripted and its mathematical model is presented (Phase 1). Next, the data-driven models between control targets and variables are established using APSO-ELM algorithm (Phase 2). Finally, in Phase 3, MODE is utilized to generate the Pareto optimal solutions, and then TOPSIS is applied to select a best compromise solution among the Pareto optimal solutions for this optimization problem. Different from traditional manual decision, it can avoid the subjectivity and randomness. In addition, the proposed optimization framework can search a trade-off solution that is consistent with production needs and provides proper and feasible solutions rapidly for decision makers.

III. MULTI-OBJECTIVE OPTIMIZATION STRATEGY OF BURDEN SURFACE
As mentioned above, burden surface optimization can be treated as a MOP. In addition, the process is difficult to be modeled due to its complex nonlinear characteristic in nature. In order to tackle these issues, we apply the proposed multi-objective optimization strategy to determine the optimal burden surface according to production status based on the established data-driven process models.

A. PROCESS MODEL BASED ON APSO-ELM
The premise of burden surface optimization is to establish the accurate process models. Due to the high complexity of the furnace interior, data-driven process models are constructed for the two objectives and the two constraints, i.e., GUR, CR, PI, and HMSC.
ELM is a competitive machine learning method for training single hidden layer feedforward neural networks (SLFNs) with fast learning speed and good generalization performance [33]- [35]. Different from other conventional machine learning methods, e.g., back propagation neural network (BPNN), or support vector machine (SVM), its hidden layer parameters can be generated randomly, and the output weights are analytically determined by Moore-Penrose generalized inverse. The universal approximation capability of ELM has been proved theoretically [36], [37]. Due to the above advantages, ELM has been widely used in many fields [38]- [40].
In order to establish the data-driven models between the variables and the output shown in Eq.(3) and Fig. 6, we employ the following modified APSO-ELM algorithm.
For N given samples ∈ n represents n-dimensional input attributes of the ith sample, and y i = [y i1 , y i2 , . . . , y im ] T ∈ m represents m-dimensional output variables. As four data-driven models are expected to be established, there are four output variables, i.e., m = 4. The output function of ELM isŷ 35716 VOLUME 8, 2020 where L is the number of hidden nodes, a i = [a i1 , a i2 , . . . , a in ] and b i are the learning parameters of the ith hidden node, β = [β 1 , β 2 , . . . β L ] T is the output weight vector between the hidden layer and the output layer, and h(x) is the hidden-layer output. h(x) actually maps the data from the n-dimensional input space to L-dimensional hidden-layer feature space and is decided beforehand.
The learning process of ELM aims to minimize the training error as well as the norm of the output weight. Thus, it can be represented as a constrained optimization problem to identify β: where e i is the training error of the ith sample. According to Karush-Kuhn-Tucker (KKT) theorem, a Lagrangian function can be constructed as where α i denotes the Lagrangian multiplier. The optimality conditions of Eq.(7) are given in the following equation: After solving Eq.(8), the estimated output weight of β is obtained as follows:β After obtainingβ, the predicted output of ELM can be represented by Eq.(5).
However, due to the random determination of the hidden layer learning parameters, some un-optimal hidden layer parameters may be generated. It should be noted that the parameters may impose negative impacts on the performance of ELM [41]. Adaptive particle swarm optimization (APSO) can perform a global search over the entire search space, which has the advantages of less parameters, low computational complexity and fast convergence speed. In addition, the inertia weight is dynamically adapted for every particle according to the fitness [42], [43]. Therefore, we adopt APSO algorithm to search the optimal hidden layer parameters of ELM.
The performance of ELM is evaluated using the following root mean squared error (RMSE): In normal case, RMSE on the whole training dataset is used as the fitness. However, it may cause overfitting problem to ELM [44]. Hence, we set the fitness to RMSE on the validation dataset (randomly selected from the training dataset) instead of using the whole training dataset [41], [44]. Firstly, all the hidden layer parameters a i = [a i1 , a i2 , . . . , a in ] and b i are randomly initialized within the range of [-1, 1]. Then, we establish the ELM model, and meanwhile, APSO is used to search the permitted area to decrease the fitness value η gradually. Finally, η converges to the minimum, and the corresponding a i and b i are the optimal parameters.
As mentioned above, we can summarize the proposed APSO-ELM in the following steps: Inputs: a training dataset, a validation dataset, a testing dataset, number of hidden nodes L, the maximum iteration times T max , the maximum and minimum of inertia weight ω (ω max and ω min ), population size .
Step1: Randomly generate particles for population. Initialize iteration counter k = 0, generate initial velocities ν k iτ of each particle. Initialize particle position ϕ k iτ as the best position P best , and the position with the best fitness η τ (calculated by Eq.(10)) of all particles as the best position of the entire swarm G best .
Step2: k = k +1, update the velocity and position for each particle: Step3: Evaluate fitness η τ of new position and update the best position of the entire swarm. If η τ < η τ , then η τ = η τ and the best position of position P k best = ϕ k iτ . Then denoting η e = min τ =1 η τ , and the corresponding position G k best is the new best position of the entire swarm. If η e < η e , then η e = η e and G best = G k best . Step4: If k < T max , turn back to Step 2. Otherwise, the global best position G best is the optimal parameter (a i , b i ) of ELM.
Step5: Obtain the updated ELM model with the optimal parameters by Eq. (5).
Then, we employ APSO-ELM to establish the data-driven models in regarding with GUR, CR, PI and HMSC, respectively. In this sense, the burden surface optimization problem VOLUME 8, 2020 shown in Eq.(4) is formulated. Next, the selected objectives are optimized according to problem formulation by a MOEA.

B. OPTIMIZATION STRATEGY FOR BURDEN SURFACE BASED ON MODE AND TOPSIS
At present, the most popular multi-objective algorithms based on meta-heuristic methods are the MOEAs [45]. Many MOEAs have emerged and obtained satisfactory achievement [46]. In recent years, MODE is overwhelmingly applied for several engineering problems [27], [47], which uses DE [48] as the underlying optimization technique. The salient features of MODE are: (a) it has a low computational complexity, (b) it is easy to implement with simple structure, (c) it has fast convergence speed and strong robustness, and (d) it has efficient constraint processing method. Therefore, MODE is employed to search the Pareto optimal solutions. MODE can be outlined by Algorithm 1 (referred to [24] for more details of MODE).
Pareto optimal solutions can be obtained through MODE algorithm. However, even though the results are informative, the number of solutions may still be prohibitive for a decision maker to make suitable choices. At this point, it is very important to select a representative solution from the Pareto optimal solutions [49], [50]. TOPSIS is a multi-criteria decision analysis method, which was original developed by Hwang and Masud [51]. The basic principle of this method is that choosing the best alternative should have the shortest vector distance from the positive ideal solution (PIS) and the longest vector distance from the negative ideal solution (NIS) [52]. Accordingly, TOPSIS is used to rank the Pareto optimal solutions to select a best compromise solution.
An evaluation matrix R = (x ij ) m×n is firstly formulated, which consists of m alternatives and n criteria. Then, after initialization and weighting, the weighted normalized evaluation matrix is constructed as follows: µ ij = r ij × ω j , i = 1, 2, . . . , m; j = 1, 2, . . . , n (14) where , i = 1, 2, . . . , m; j = 1, 2, . . . , n and n j=1 ω j = 1. Subsequently, the evaluation distance between each alternative, and PIS and NIS, respectively, are calculated, which denote as d + i and d − i . Finally, the closeness coefficient C * i is calculated using Eq. (15). Ranking all the alternatives according to C * i and the optimal solution is obtained.
During the optimization process of burden surface, the operating status parameters are taken as the current detection value and not involved in the optimization link. The burden surface features are decision variables. Then, MODE enhanced with TOPSIS is employed to search for burden surface features satisfying the optimization conditions under Algorithm 1 MODE Algorithm 1: Iteration times t = 0 2: Create a random initial population P i,t , ∀i, i = 1, . . . , N pop 3: for t = 0 to I t max do 4: for i = 1 to N pop do 5: U i,t+1 =P i,t 6: end for 7: for i = 1 to N pop do 8: Select randomly three different chromosomes r 1 , r 2 , r 3 9: Generate a random integer value i rand from 1 to D 10: for i = 1 to N pop do 11: Generate a random real value rand j belongs to [0, 1] 12: if rand j < CR or j = i rand then 13:  if n + k < N pop then 25: for i = n + 1 to n + k do 26:  29: Apply crowding distance sorting to V p,t+1 30: for i = n + 1 to N pop do 31:

IV. SIMULATION RESULTS
In this section, we present the simulation results based on the actual production data from a BF to verify the effectiveness of the proposed multi-objective optimization strategy for burden surface. Data-driven process models are firstly established based on APSO-ELM algorithm. After that, the MODE enhanced with TOPSIS is applied to obtain the optimal solutions of the MOP.

A. DATA DESCRIPTION AND PREPROCESSING
The data are collected from a medium-size BF with an inner volume of about 2500 m 3 . Some important and concerned variables are considered as the operating status parameters, including blast temperature ( • C), blast volume (m 3 /min), blast pressure (kPa), top pressure (kPa), top temperature (including four-point temperature) ( • C), differential pressure (kPa), and oxygen enrichment (%). The daily production data from February 2014 to June 2014 are selected from the historical database for analysis. Fig. 7 demonstrates the series of GUR, blast volume and PI. According to Fig. 7, there is a great difference in the magnitude of the variables clearly. Considering the impacts of convergence and complexity on modeling, all the samples are normalized into [0, 1] to eliminate the influence of magnitude before applying in the experiments. The method is as follows: wherex k and x k are the kth variable before and after change respectively, min(x k ) and max(x k ) are the minimum and maximum values of the kth variable before normalization. In addition, the raw data contain outliers, which may affect the performance of the created models. These may be due to irregular behaviors of the furnace interior in a certain period of time, furnace shutdown or some wrong readings [53]. In order to ensure the reliability of data-driven models and decision making, the elimination of outliers is performed. The histogram of blast volume and blast temperature is depicted in Fig. 8. From this figure, it can be seen that the data distribution does not obey normal distribution. In this point, it is unreasonable to eliminate the outliers by the criterion of 3σ . Therefore, the box-plot method is adopted to eliminate outliers, which has no requirements for data distribution. Fig. 9 represents the box-plot diagram of variables. The red crossing mark is the extreme abnormal points, which are choose as the outliers. The data can be used for modeling and decision making after the aforementioned preprocessing.   In this section, data-driven process models are constructed for the two objectives and the two constraints based on APSO-ELM algorithm. We randomly select 800 samples from the preprocessed data. The first 600 samples are chosen as training samples ℵ tr = {(x i , y i )} 600 i=1 , and the remaining 200 samples are used as testing samples ℵ te = ,where x i ∈ 17 is the input vector and y i ∈ 4 is the output for GUR, CR, PI and HMSC, respectively. Because CR is a statistic once a day, 100 samples are used as training samples and 50 samples are used as testing samples for CR. The sigmoidal function G(a, b, x) = 1/(1+exp(−(a·x +b))) is adopted as the activation function. In order to construct the accurate data-driven models using APSO-ELM algorithm, the number of hidden nodes L needs to be determined. We gradually increase the hidden nodes number, and select the one with the minimum testing error as the final one. Taking the data-driven model for GUR as an example, Fig. 10 presents the changing trends of training error, testing error and training time with the increase of L. In Fig. 10, the green and blue curves correspond to training and testing errors, and the red dotted curve is the training time. As observed from Fig. 10, when L is less than 500, both errors are relatively large. With the increase of L, the errors (i.e., RMSE) of the model become smaller gradually and the training time is increased. The minimum testing error is achieved when L is within the range [500, 1000]. Considering the computational complexity and testing error, L is set as 600. The same method is performed for CR, PI and HMSC model to select the corresponding optimal L.
The parameters of APSO used in APSO-ELM are set as follows: population size is 25; the maximum iteration times T max is set to 150 to get enough iterations to search for the optimal solution; the maximum and minimum of inertia weight ω (ω max and ω min ) are set as 0.9 and 0.4, respectively; the acceleration coefficient c 1 and c 2 are both set as 2.
As a result, four data-driven process models in regarding with GUR, CR, PI and HMSC with ELM and APSO-ELM, respectively, are depicted in Fig. 11. The number of hidden nodes L are set to 600 for both ELM and APSO-ELM in this simulation to make it more reasonable and fair. Fig. 12 shows the scatter diagram of actual value and predicted value. According to Fig. 11, the predicted values of APSO-ELM are closer to the actual values than ELM. As observed from Fig. 12, it is apparent that the scattered points obtained by APSO-ELM are more concentrated in the vicinity of the diagonal, i.e., y = x. For GUR model, the most predicted values are between y 1 = x + 0.025 and y 1 = x − 0.025. It shows that APSO-ELM provides more actuate predicted results and lays a solid foundation for subsequent optimization.
Furthermore, the training RMSE, testing RMSE and testing mean absolute percentage error (MAPE) are summarized in Table 1. Twenty trials are carried out for each approach. The average results are adopted as the comparison results. According to Table 1, APSO-ELM outperforms ELM and SVM, which again shows that the four data-driven process models can obtain satisfactory performance. It proves that the effectiveness of hidden layer parameters determination based on APSO and implies that optimal parameters are searched. In this way, the created model can be used to provide the basis for burden surface optimization.

C. OPTIMIZATION USING MODE ENHANCED WITH TOPSIS
In this section, the proposed MODE enhanced with TOPSIS, presented in Section III-B, are applied to solve the MOP of burden surface formulated as Eq.(4). There are 7 decision variables in this optimization problem, represented by X = (l 1 , l 2 , h 1 , h 2 , α, β, γ ). Before running the MODE algorithm, its parameters need to be specified to control its convergence rate. The parameters are set as follows: the number of population N pop is set as 150, the maximum iteration times I t max are assumed to be 80,100 and 150, crossover probability CR is 0.5, mutation factor F is 0.6.
The Pareto optimal solutions for GUR and CR in BF ironmaking process obtained by MODE with different I t max    are given in Fig. 13. From Fig. 13(a), the Pareto optimal solutions are not ideal. According to Fig. 13(b) and Fig. 13(c), with the increase of the maximum iterations, the solutions trend to convergence gradually. Therefore, I t max is set as 150. Table 2 presents a summary of the results obtained by MODE. As observed from Fig. 13(c) and Table 2, it is possible to generate optimized operating conditions from the Pareto optimal solutions to improve production performance. Furthermore, the Pareto optimal solutions are prioritized using TOPIS technique. For this purpose, objective functions are considered as criteria. The weights of the criterion are assumed to be equal. The rank of each solution is obtained after applying the TOPSIS. The top 10 optimal efficient solutions obtained by TOPSIS are presented in Table 3. In addition, the top 10 ranked optimal solutions are represented in Fig. 14. Actually, considering the specific requirements, operators can select appropriate burden surface features to maximize efficiency. In this paper, we select the first ranked solution as the final compromise solution. Furthermore, considering the production scheduling, the weights of the criterion can be set to different values at a certain time period by operators to give different attention. Thus, it will give a solution that is more suitable for production planning. Finally, according to the obtained burden surface features, the burden surface profile can be drawn out.
In order to evaluate the performance of the proposed integrated multi-objective optimization framework in dealing with burden surface optimization, we develop it in three  different production situations. Table 4 shows the comparison results between the actual production values and the multi-objective optimization results. In order to more intuitively illustrate the burden surface optimization results, the burden surface profile are depicted in Fig. 15. From Table 4, we can find that the optimal burden features are adjusted with the different production situations. Compared with pre-optimization, GUR has increased, while CR has decreased. It is consistent with production requirements for energy efficiency. The phenomenon indicates that the multi-objective optimization strategy can provide an effective and feasible solution for the operators. In addition, as observed from Fig. 15, the optimized burden surface profile is ''platform + funnel'' type, which accords with the charging operation and again shows the reliability of the optimization results.
In addition, in order to better verify the effectiveness and feasibility of the proposed integrated optimization framework for this optimization problem in BF ironmaking process, we choose some relevant state-of-the-art MOEAs as comparison, including NSGA-II [20], PESA-II [54], and NSGA-III [55]. The usual Hypervolume (HV) [56], [57] is used as the performance metric in the following experiments, which calculates the volume of the objective space between the obtained solution set and a reference point. The larger the HV, the better performance of the solution set in terms of convergence, diversity and uniformity can be obtained. For the fairness, the maximum iteration times are assumed to be the same value in all comparison algorithms. Twenty trials are carried out in each approach. The comparison results of statistical properties including mean and standard deviation (SD) are summarized in Table 5. As observed from Table 5, MODE has smaller SD than other three methods, which shows that MODE runs more stable. In addition, for the measurements of mean of HV, compared with PESA-II and NSGA-II, MODE can achieve higher HV value. Meanwhile, the HV value of MODE is little lower than that of NSGA-III. Therefore, MODE has satisfactory performance on convergence, diversity and uniformity. The main reason is that differential evolution generates new parameter vector by adding the weighted difference between two population vectors to a third vector to ensure the diversity and convergence in mutation phase.
Furthermore, If K , D, N pop and I t max represent the number of objectives, the number of optimal parameters, the number of population and the number of iterations, respectively, the running complexity of MODE is O K · D · N pop · I t max [24]. In addition to the calculated   performance measure, the running time of the above optimization algorithms over ten independent runs are also calculated to evaluate the computational complexity. The average running time with 150 iterations are obtained as 607.416s, 873.365s, 890.652s, and 875.500s, respectively. It can be found that MODE takes the least time, which may be mainly due to its parallel search mechanism. Overall, it is obvious that MODE achieves better results as compared to other algorithms, further indicating that MODE is more suitable and effective for dealing with this kind of burden surface optimization problem.

V. CONCLUSION
In this paper, the optimal setting of burden surface is converted to a MOP, and an integrated multi-objective optimization framework is proposed for solving the corresponding MOP. Due to the high complexity of the furnace interior, data-driven process models between objectives and variables are established using APSO-ELM. Then, MODE enhanced with TOPSIS method is used to optimize the goals to obtain a best compromise solution for GUR and CR, meanwhile, PI is guaranteed to vary in prescribed bounds to ensure the smooth and stable production and HSMC is limited to a given range to ensure the production of qualified hot metal. The actual production data are employed to validate the effectiveness of the proposed integrated multi-objective optimization framework. The experimental results indicate that the established data-driven models are more accurate and the obtained optimal solution provides important guidance to optimize BF ironmaking process.
In the present work, the optimization framework mainly focuses on analysis of deterministic model and only considers two objectives. However, in practice, variations in raw material properties and production process lead to frequent changes of model parameters including the constraints. In addition, there may be more than two objectives that need to be optimized for further improving the optimal performance. The future works worth exploring would be establishing online model, introducing efficient constraint handling mechanism and improving the multi-objective optimization algorithm.