Assessing the Predictability for Blast Furnace System through

Mar 28, 2008 - This paper is concerned with assessing the predictability for a blast furnace (BF) system through ... of silicon content in blast furna...
0 downloads 0 Views 362KB Size
Download PDF
Ind. Eng. Chem. Res. 2008, 47, 3037-3045

3037

Assessing the Predictability for Blast Furnace System through Nonlinear Time Series Analysis Chuanhou Gao Institute of Systems Optimization Techniques, Department of Mathematics, Zhejiang UniVersity, Hangzhou 310027, People’s Republic of China

Zhimin Zhou College of Mathematics and Statistics, Zhejiang UniVersity of Finance and Economics, Hangzhou 310018, People’s Republic of China

Jiming Chen* Institute of Industrial Process Control, Department of Control, Zhejiang UniVersity, Hangzhou 310027, People’s Republic of China

This paper is concerned with assessing the predictability for a blast furnace (BF) system through nonlinear time series analysis of silicon content in blast furnace hot metal, in three representative blast furnaces with different volumes. The results indicate that the predictability of silicon content in hot metal of the selected BFs approaches that of a totally predictable system while it departs from that of a totally unpredictable system greatly. The predictability of silicon sequences conversely renders a strong indication of the presence of a deterministic mechanism ruling the dynamics of the ironmaking blast furnace process. In the meantime, the results also provide important information on selecting the minimum number of macroscopic variables and the delay time needed for reconstructing the dynamics of these blast furnace systems. Moreover, the above information is found to be independent of the size of these blast furnaces. However, there is a slight difference in the complexity of the silicon sequence collected from different types of blast furnaces. As one might expect, larger difficulty will be encountered in predicting the silicon content of a pint-sized blast furnace because of the presence of the larger complexity. 1. Introduction Ironmaking blast furnace process (IBFP) is a crucial operation unit in the manufacture of iron and steel. For a long time, the modeling and control of IBFP have been always the important subjects in the metallurgical field since a blast furnace (BF) is a sophisticated reactor and is also a significant item in the economy of any country. Among these subjects, exactly predicting the evolution of silicon content in blast furnace hot metal constitutes the important issues for the economic operation of BFs. The silicon content in BF hot metal usually acts as an indicator to represent the inner thermal state of BF and further determines the quality of BF hot metal. After decades of extensive research on the thermodynamic and kinetic behaviors occurring inside BFs,1-9 a fundamental understanding of silicon transfer has been progressively established. However, a fully reliable predictive model based on ironmaking mechanism has not emerged so far. Thus, data-driven models have been investigated quite intensively in recent years10-20 and exhibit good performance in determining the input-output relationships. Despite being superior to analytical models in this regard, a distinct shortcoming for data-driven models lies in the difficulty of choosing an appropriate tool and algorithm to establish the input-output relationships, because there are many tools that can approximate arbitrary function in arbitrary accuracy. Practical applied tools for silicon prediction include neural net,10-15 nonlinear time series analysis,16 partial least-squares,17,18 fuzzy theory,19 and recursive algorithm (see ref 20 and references * To whom correspondence should be addressed. Tel.: +86-57187953762. Fax: +86-571-87951879. E-mail: [email protected].

therein). Notwithstanding the progress in the development of predictive techniques of silicon content in BF hot metal achieved using the above methods, the understanding of the dynamical behavior of the evolution of silicon content is still far from complete. Since determining whether dynamical behavior of systems exhibit regular, chaotic, or random features is usually a goal in many problems,21 it is of prime importance to understand the inherent mechanism hidden behind the time series of silicon content in BF hot metal for constructing a suitable predictive model. In the early work, we have proven quantitatively that there is a chaotic mechanism ruling the evolution of silicon content in BF hot metal by calculating the maximum Lyapunov exponent,22 which frequently acts as an indicator to represent the chaotic feature of systems, power spectrum23 and fractal dimension.24 More details about how to determine the chaotic mechanism ruling the evolution of silicon content in hot metal may be found in refs 20-24. For this reason, chaos seems, at least in principle, to be the optimal tool to predict silicon content in BF hot metal. However, this paper is not concerned with how to establish a chaotic model for predicting silicon content in BF hot metal but rather with measuring how difficult it is to predict the next point, i.e., measuring the predictability of silicon content in BF hot metal, and it also provides some preliminary information for constructing a chaosbased data-driven model. The ironmaking BF is a giant complex reactor in which a large number of gas-solid, solid-solid, and solid-liquid chemical reactions and countercurrent heat exchange coupled with complicated heat and mass transfer phenomena take place.25 To quantify the complexity of the dynamical behaviors of such

10.1021/ie070879s CCC: $40.75 © 2008 American Chemical Society Published on Web 03/28/2008

3038

Ind. Eng. Chem. Res., Vol. 47, No. 9, 2008

The ironmaking BF process is a complex industrial semibatch process. During the operation of a BF, the lump solid raw materials consisting of iron ore, coke, and fluxes are alternatively in the shape of layers charged into the top of the furnace, while the preheated compressed air and auxiliary fuels including pulverized coal, oil, or natural gas are introduced through tuyeres equipped just above the hearth. As the various materials move down through the furnace, they are progressively heated and reduced by the hot ascending gases that contain rich carbon monoxide, then form the final products, molten hot metal and slag. In the meantime, the generated carbon dioxide, along with the unused carbon monoxide and the part of nitrogen originally introduced, forms the BF gas as a byproduct. Since this BF gas contains a great deal of chemical energy in the form of carbon oxide, it is usually introduced into the hot stoves as a heat resource to heat the fresh blown air after being released from the top of the BF. The molten hot metal and slag continuously trickle down into the hearth and form two separate layers due to the difference of density. The upper layer is slag, while the lower layer is hot metal. Finally, the liquid hot metal and slag are tapped at regular intervals through tapholes and sent for further processing. The whole ironmaking blast furnace process will last ∼6-8 h. Figure 1. Schematic diagram of a blast furnace.

3. Method for Measuring the Predictability a system is one of the most challenging endeavors. Since the direct experimental measurement for BF is exceedingly difficult in the case of high temperature, pressurized operation, existence of packed particles, and distributed parameter characteristics, it is not possible to monitor all relevant physical variables pertaining to the ironmaking process, to say nothing of the exact relationships among them. Thus, data analysis of time series of some physical variables of system becomes a useful tool for characterizing the dynamical behavior of the system, partly based on the fact that every physical variable contains enough information on the complex dynamics of system.26 For instance, the time series of silicon content in BF hot metal can be used to measure the complexity of the dynamical behavior of the IBFP. There are various approaches to measure the complexity of a time series data set.27 For the current study, an alternative approach called forecast entropy28 is chosen because it can not only measure the predictability for IBFP but also provide the important information on the delay time and the number of parameters needed for reconstructing the original dynamical system. The latter will be very valuable for building a chaotic predictive model for IBFP. The rest of this paper is organized as follows: Section 2 introduces briefly the ironmaking blast furnace process. The detailed algorithm for measuring the predictability of a time series is provided in Section 3. Section 4 gives the experimental data considered in this work. The results and discussion are presented in Section 5. Finally, Section 6 concludes this paper. 2. Ironmaking Blast Furnace Process A typical BF consists of a tall, vertical, chimneylike shaft furnace. It is narrow at the top, increasing in diameter downward but becoming narrow again suddenly almost at the bottom. The zone of a BF from the top to the bottom can be roughly divided into five sections called throat, stack, belly, bosh, and hearth, respectively. A representative schematic diagram of a BF is shown in Figure 1, in which the main chemical reactions are listed. In general, a BF can serve continuously for 6-10 years, producing ∼3 000-10 000 tons of hot metal each day.

Very often, investigating the predictability of a time series is implemented by measuring its spatial complexity. A commonly used approach for measuring the complexity of a distribution is to calculate the Boltzmann-Gibbs-Shannon (BGS) entropy, which stems from information theory and can be expressed as the expected value of a surprise function over all states.29 Here, the surprise function is defined as the minus logarithmic value of the probability of finding a point in a specific neighborhood. In the practical calculations, denote by pi the probability of finding a point in a neighborhood i, then the surprise function is defined to be -ln pi. Thus, the BGS entropy can be written as

S)-

∑i pi ln pi

(1)

In general, the larger S is, the more complex the distribution is, i.e., the more difficult it is to predict the next point. Despite its popularity in the complexity measurement, the BGS entropy often fails in some very complex cases because it just considers the information in regular space. Therefore, in this paper, an alternative technique called forecast entropy is used for measuring the predictability of a time series instead of BGS entropy. The main idea of forecast entropy is to calculate the summation of weighted BGS entropy of a distribution on different scales. 28 Compared with the BGS entropy, the advantage of forecast entropy lies in its considering the information not only in regular space but also in tangent space. In what follows, a detailed algorithm for calculating forecast entropy is given.28 Ndat Suppose a one-dimensional finite time series {xi}i)1 ) {x(t0 N + (i - 1)∆t)}i)1dat sampled from a dynamical system, where t0 is the initial sampling time, ∆t is the sampling interval, and Ndat is the size of time series satisfying Ndat ) 2m for the convenience of calculations. Here, m is a positive integer. Denote the maximum and minimum values of this time series as xmax and xmin, respectively, then the forecast entropy is defined as m

F)

∑ aiSi i)1

(2)

Ind. Eng. Chem. Res., Vol. 47, No. 9, 2008 3039

where Si is the BGS entropy on scale i and ai is the associated weighting coefficient, which is dependent only on the size of the time series, Ndat. The partition of scales is often executed by separating the group into subgroups, and then into smaller subgroups according to the numbers of points. For instance, for one, the time series is partitioned into two subgroups at the center of [xmin, xmax ]. Denote nL and nR as the numbers of points on the left and right subgroups, then the probabilities of finding the next point in the left and right subgroups can be expressed to be pL ) nL/(nL + nR) and pR ) nR/(nL + nR), respectively. Hence, according to eq 1, it is easy to obtain the BGS entropy on scale 1 as

S1 ) -pL ln pL - pR ln pR

(3)

For scale 2, the left and right subgroups are again partitioned into smaller subgroups at their centers, and the numbers of points in the new subgroups are denoted as nLL, nLR, nRL, and nRR, respectively. Then the values of probabilities in each regime are pLL ) nLL/(nLL + nLR), pLR ) nLR/(nLL + nLR), pRL ) nRL/ (nRL + nRR), and pRR ) nRR/(nRL + nRR). Finally, obtaining the BGS entropy on scale 2 is

S2 ) -pLL ln pLL - pLR ln pLR - pRL ln pRL - pRR ln pRR (4) The above procedure is iterated to obtain more BGS entropies at larger scales until the scale is m. Besides the BGS entropies, there are yet weighting coefficients ai, i ) 1, ‚‚‚, m, needed to be determined in order to calculate forecast entropy. In ref 28, the determination of these parameters is implemented by defining the forecast entropy of an ideal random system as 1. More details about the ideal random system and how to determine these parameters will be described in Appendix A. Here, direct substitution of the results of aj into eq 2 can get the forecast entropy to be

F)

m

1 m

2 ln 2

(2S1 +

∑ Si) i)2

(5)

When the forecast entropy is used to evaluate the predictability of a distribution, there are two extreme cases for the values of forecast entropy, i.e., zero, representing any completely predictable system such as periodical systems, and 1, representing a totally unpredictable system such as an ideal random system.28 Usually the forecast entropy of a practical sequence lies between the two extreme values. The larger the forecast entropy is, the more difficult it is to predict the future behavior. The above procedure for calculating the forecast entropy is only suitable for a one-dimensional system. In the case of multidimensional systems, the forecast entropy, for simplicity, may be roughly defined as the average of the forecast entropies of the system along each coordinate.28 A slightly complex and more accurate algorithm for measuring the forecast entropy of multidimensional systems that utilizes the information in tangent space will be described in Section 5. 4. Experiment Data For the current research, the data sets of time series of silicon content in BF hot metal collected from a medium-sized BF with a volume of 2000 m3 and two pint-sized BFs with volumes of 750 and 380 m3 are taken as the samples, labeled (a), (b), and (c), respectively, as shown in Figure 2. Obviously, the values of silicon content in BF (a) are universally lower than those in BFs (b) and (c). This implies a higher furnace temperature and

Figure 2. Time series of observed silicon content in BF hot metal measured from BF (a), (b), and (c).

larger energy consumption in BFs (b) and (c). Also, the stronger fluctuation of silicon content emerges in BFs (b) and (c). It is also a possible result that the pint-sized blast furnace can exhibit more complex behavior than the medium-sized one, which may cause larger control difficulties for the former. The main motivation of selecting these three kinds of BFs is, on the one hand, to gain deeper insight into in-furnace phenomena with typical size and, on the other hand, to capture the possible difference in dynamical behavior of BFs with different sizes. The amounts of silicon in hot metal depicted in Figure 2 are measured by a technique often used in a laboratory. In each tapping, a sample of liquid hot metal is ladled and then sent to the laboratory for compositions measurement through spectroscopic analysis. The tapping frequency is ∼1 h for BF (a) and ∼2 h for BFs (b) and (c), and the size of the time series is 2 048 ) 211. 5. Results and Discussion 5.1. Calculation of Forecast Entropy for Original Sequences. To calculate the forecast entropy of a blast furnace system only from the time series of silicon content in BF hot metal, the dynamics of this system must be reconstructed from the observations at first. Ndat stand for the actual measured value of silicon Let {yisignal}i)1 content in BF hot metal exhibited in Figure 2. Under the condition of choosing suitable delay time τ and embedding dimension d, the patterns signal signal T Yisignal ) (ysignal , yi-τ , ‚‚‚, yi-(d-1)τ ) ∈ Rd, i

i ) (d - 1)τ + 1, ‚‚‚, Ndat - 1, Ndat (6) are thought to be capable of tracing out the orbits of the reconstructed system.30 Furthermore, under generic conditions, i.e., τ, d , Ndat, such patterns have one-to-one correspondence to the states of the original physical system, and the reconstructed dynamics is topologically equivalent to the underlying dynamics provided that d g 2dA + 1, where dA is the box-counting dimension of the chaotic attractor.26,31 The calculation of forecast entropy of this d-dimensional reconstructed system, denoted by Fd in this context, is performed according to28

Fd )

2

d-1

∑ (d - j)Fj

d(d - 2) j ) 1

(7)

3040

Ind. Eng. Chem. Res., Vol. 47, No. 9, 2008

Figure 3. Forecast entropy of original silicon series versus the delay time, τ, of BF (a), (b), and (c) in the case of (I) n ) 2d+2 and (II) n ) 16 when the embedding dimension, d, changes from 2 to 9 increasingly.

where Fj represents the forecast entropy of the distribution on signal ) plane. Detailed descriptions about the calcuthe (yisignal, yi-jτ j lation of F may be found in Appendix B. As a sound method, forecast entropy should be consistent when used in practice, i.e., it should be able to handle the randomness of the data and improve its accuracy as the amount of data increases. Therefore, the consistency of forecast entropy should be validated before it is applied to actual measured silicon series. Appendix C exhibits the detailed process of validating the consistency of forecast entropy. Utilizing the algorithm described in Appendix B, the forecast entropy of the measured silicon sequences can be achieved readily under different reconstructing parameters. Results are given in Figure 3. Here, two cases are considered in the procedure of calculating the forecast entropy. One is the number of neighbors, n (see step 1 in Appendix B), changing with the embedding dimension according to the relationship of n ) 2d+2, the other is the fixed n, i.e., n ) 16, denoted as cases I and II, respectively. The former can give the optimal embedding dimension and delay time for prediction purposes, while the latter can give the best candidate of dynamical dimension and delay time for reconstructing the dynamics of the system producing the sequences.28 A look at Figure 3 might suggest the following information. In case I, (i) except at τ ) 1, Fd of BF (a) increases basically as d changes from 2 to 6, then decreases after d g 7, while Fd of BFs (b) and (c) interlace with one another when d changes from 2 to 4, then decrease after d g 5; (ii) at a given d, Fd of the three BFs will remain almost unchanged as τ increases after a quick increase from τ ) 1 to τ ) 2; (iii) because of a limit to Ndat ) 211, only embedding dimensions of up to d ) 8 are included in this case of n ) 2d+2, and for these considered d and τ, the minimum Fd of the three BFs all appear at τ ) 1 and d ) 3. While in case II, (i) Fd of the three BFs decrease monotonously as d changes from 2 to 6, then remain almost unchanged afterd g 7; (ii) similar to that in case I, Fd also increase quickly in the starting part of τ, then remain almost unchanged as τ increases; (iii) the minimum forecast entropies of the thee BFs appear at τ ) 1 and d ) 7. As an alternative method, forecast entropy can provide important information on the delay time and the number of parameters needed for reconstructing system dynamics. This information seems to be obvious for these selected BFs from the above results. However, the above results should be

interpreted with great caution. The main reasons are that, on the one hand, the omnipresence of noise such as process noise and measurement noise may add the randomness of actual silicon series, and on the other hand, the size of the selected data is only 2 048, being rather short compared with the usually used size of data for nonlinear time series analysis, which may be not long enough to make the results approach the true values as close as they possibly can. The validity of the above results may be reduced due to the presence of random noise and the short size of the data sets. Therefore, it is necessary to verify the validity of the above results even if the algorithm of forecast entropy exhibits good performance in handling the data contaminated by noise.28 To this task, attention is turned below. 5.2. Calculation of Forecast Entropy for Filtered Sequences. 5.2.A. Preprocessing of Experimental Data. As every realistic measurement, especially like the signal observed from complex IBFP, is inevitably contaminated by noise, such as process noise and measurement noise, the preprocessing must be undertaken before the measured signal is used for analysis. To reduce the noise limitation on detecting the inherent characteristics of the data, an extremely simple nonlinear method of noise reduction is adopted here.32 The method states the } to be following: assume a finite-scale time series {ysignal i composed of a clean signal yiclean and a noise signal ηi such that

yisignal ) yiclean + ηi

(8)

then the clean signal yiclean can be obtained through replacing each yisignal by the average value of this coordinate over points in a suitably chosen neighborhood. The selecti on of neighborhood is based on the widely used method of delay coordinates, also called embedding theorem.26,30 The detailed procedure is as follows:32 First, construct multivariate state vectors involving the past and future coordinates as signal , ‚‚‚, yisignal ) Ui ) (yi-k +l

(9)

where k and l are fixed positive integers; then, find all the neighboring points of Ui in a small ball of radius  through the l∞ norm, i.e., signal signal |Uj - Ui|∞ t max{|yj-k - yi-k |, ‚‚‚, |yjsignal signal signal |, ‚‚‚, |yj+l - yi+l |} <  (10) ysignal i

Finally, calculate the mean value of the (k + 1)th entry of the neighboring points according to

yclean ) i

1

∑

|Ω(Ui)|Uj ∈ Ω(Ui)

ysignal j

(11)

where Ω(Ui) represents the neighbor set of vector Ui and |Ω(Ui)| means the number of elements within the set Ω(Ui). The calculated mean value is taken as the requested clean signal at time i. There are three factors that may affect the amount of noise reduction, i.e., k, l, and . Seen from eq 10, the selection of the l∞ norm instead of the usual Euclidean norm can greatly reduce the effect of the values of k and l on the acquirement of the neighboring points. Thus, the effect of their values is ignored here, and for the sake of simplicity, let them be the same value of 1. The effect of  may be found by calculating the dependence of the absolute noise level σ2 ) 〈η2〉 on the value of . Shown in Figure 4 are the results.

Ind. Eng. Chem. Res., Vol. 47, No. 9, 2008 3041

Figure 4. Dependence of the absolute noise level on the size of the chosen neighborhood for BF (a), (b), and (c). Figure 6. Forecast entropy of filtered silicon series versus the delay time, τ, of BF (a), (b), and (c) in the case of (I) n ) 2d+2 and (II) n ) 16 when the embedding dimension, d, changes from 2 to 9 increasingly. Table 1. Analysis of Results of the Forecast Entropy for the Original Series (Typed with Normal Font) and the Filtered Series (Typed with Italic Font)

From Figure 4, it can be observed that the absolute noise amplitudes of three groups of measurements will graduate into increasing slowly after a quick increase for a period of time with the change of . The optimal value of  should be the starting point when the value of σ2 increases slowly.32 The arrows point out the optimal values of  used for noise reduction, ∼0.35 for BF (a), 0.45 for BF (b), and 0.60 for BF (c). It is also noted that, in the optimal places for noise reduction, the order of magnitude of noise of the measurements is the largest for BF (c), the second largest for BF (b), and the smallest for BF (a). This indicates that a stronger disturbance of noise is encountered in the measurements of the pint-sized BF, which will make larger trouble for understanding the dynamical behavior of the pint-sized BF. The presence of a stronger disturbance may be the source indicating that the inner phenomena taking place inside the pint-sized BF is more complex than that inside the medium-sized BF. These intuitive results can further be confirmed according to the calculation of forecast entropy in what follows. Substituting the optimal  obtained from Figure 4 into eqs 10 and 11 can easily yield the approximate clean signal of the measurements from three BFs. Figure 5 gives the same data points as shown in Figure 2 after noise reduction. Clearly, the filtered signals become much steadier than the original ones. Moreover, during some intervals, the signal is nearly the same. Remarkably, this method for reducing noise is not a well-chosen one. The randomness existing in the original signal may be filtered through this method. Therefore, the filtered signal may be not the true value of silicon content in BF hot metal. However, it is always a puzzle to obtain

case II, for reconstruction purpose

number of BF

d

τ

FEa

d

τ

FE

BF (a) BF (b) BF (c)

3 3 3

1 1 1

0.2211(0.2121) 0.2911(0.2914) 0.2550(0.2944)

7 7 7

1 1 1

0.0097(0.0095) 0.0120(0.0098) 0.0113(0.0103)

a

Figure 5. Time series of filtered silicon content in BF hot metal measured from BF (a), (b), and (c).

case I, for prediction purpose

FE represents the forecast entropy.

the effective signal from the observed signal contaminated by noise with unknown statistical features. A discussion of optimal noise-reduction method for the measurements is beyond the scope of this work. As a reference of original signal, it is still of practical significance to analyze this filtered signal. Naturally, the filtered signal is approximately taken as an example in which randomness of data is reduced for the subsequent research. 5.2.B. Calculation of Forecast Entropy. The same procedure for calculating forecast entropy is made on the filtered signal. Shown in Figure 6 are the results. Compared with the value of forecast entropy in Figure 3, the forecast entropy in Figure 6 is a little larger at the same d and τ. This is reasonable because the presence of noise will add the randomness of the original silicon series, which thus leads to the appearance of larger forecast entropy. A closer look at Figure 6 suggests that, in case I, (i) Fd of BF (a) increases as d increases from 2 to 6 except at τ ) 1 where F2 > F3, while Fd of BFs (b) and (c) interlace with each other when d ) 2 and d ) 3; then Fd of the three BFs increase monotonously after d > 3; (ii) at every given d, the same phenomena can be observed as that in Figure 3, i.e., Fd of the three BFs will remain almost unchanged as τ increases after a quick increase in the starting part of τ; (iii) the minimum of Fd of the three BFs all appear at τ ) 1 and d ) 3; while in case II, the same information is found as that in Figure 3. By comparing Figure 3 with Figure 6, it can be observed that there is a difference in the values of forecast entropy between the original series and the filtered series at every given d and τ. However, the properties of the results are the same. Namely, the minimum of forecast entropy of two groups of series appear at d ) 3 and τ ) 1 in case I and at d ) 7 and τ ) 1 in case II. Thus, the results obtained from Figure 3 should be valid. These results may serve as guidelines for determining

3042

Ind. Eng. Chem. Res., Vol. 47, No. 9, 2008

the information on delay time and the number of needed parameters for these selected BFs, both for prediction purpose and for reconstruction purpose. A detailed analysis of these results is reported in Table 1, where FE represents the forecast entropy. From Table 1, it can be obtained that a 3-dimensional embedding space with τ ) 1 is enough to do a prediction of silicon content for BFs (a), (b), and (c), and the delay time and the number of parameters needed for reconstructing their dynamics are τ ) 1 and d ) 7, respectively. The practical significance of the result for prediction purpose is that a predictor of the form

Yi+1 ) f(Yi)

(12)

where Yi ) (yi, yi-1, yi-2)T ∈ R3 can have the optimal performance if one uses the historical data of silicon sequences to do a prediction; the practical significance for reconstruction purpose is that the minimum number of macroscopic variables needed for modeling ironmaking BF dynamics is 7. As a result, the state-space model of the form

[

xi+1 ) g(xi, ui, wi) yi ) h(xi, Vi)

]

(13)

in which g and h are projects; x ∈ R7, u ∈ Rp (p is a positive integer), and y ∈ R1 are the state, input, and output variables; and w and V are the process noise and measurement noise, respectively, may be constructed to describe the IBFP. Here, y may be the value of silicon content in hot metal, u may be the blast volume, blast temperature, pulverized coal injection rate, oxygen enrichment, ore-to-coke rate, etc., and w and V are generally assumed to be the normally distributed white noise. A further look at Table 1 might suggest that the number of parameters and the delay time needed for building a predictive model or for reconstructing the BF dynamics is independent of the size of the BF. This result is easily elucidated because, in theory, the mechanism of IBFP remains unchanged; thus, the same reactions and transfer phenomena will take place during the IBFP whatever the size of the BF. In addition, Table 1 also shows that the forecast entropies of every BF, a little more than 0.2, all depart from 1 greatly. This implies that the evolution of silicon sequences does not follow the stochastic rule, and the prediction of silicon sequences is much easier than the ideal random sequence. Therefore, the predictability of silicon sequences can, in turn, act as an indicator of its determinism. Although the selection of embedding dimension and delay time is independent of the size of the BF for these selected cases, there is still a slight difference in the forecast entropies of the three BFs. As one might expect, larger forecast entropy appears in the pint-sized BF, which implies the larger complexity inside it. The possible reasons for this may be that (i) there is a stronger noise disturbance in the pint-sized BF; (ii) the raw materials used in every BF, such as iron ore, coke, pulverized coal, etc., contain different components, or different operating conditions, even different melters, can cause different complexity; and (iii) the pint-sized BF can exhibit more complex behavior than a medium-sized BF just because of its small size. More work needs doing to validate these possible reasons. 6. Conclusions From the above discussion, the following conclusions can be reached: (i) A deterministic mechanism that governs the dynamics of the evolution of silicon series is apparently indicated through measuring the predictability, and the prediction

of silicon content is much easier than the ideal random series. (ii) There is no difference in selecting the delay time and the number of parameters needed both for prediction purpose and for reconstruction purpose for BFs with different sizes, being approximately τ ) 1 and d ) 3 for the former and τ ) 1 and d ) 7 for the latter. (iii) Larger difficulty will be encountered in predicting the silicon content of the pint-sized BF because of the presence of the larger forecast entropy. (iv) The stronger noise disturbance is found in the silicon series of the pint-sized BF, which may be the source of producing the larger complexity. It should be noted that, in spite of knowing the number of degrees of freedom and the delay time used to model dynamics of the evolution of silicon sequences, there is a lot of work needed to be done in finding candidate parameters for modeling and the proper projects appearing in eqs 12 and 13. The corresponding work is being done now. Also, it should be noted that the selection of optimal  for noise reduction is through the direct observation on dependence of absolute noise level on the size of the chosen neighborhood, which is rule-of-thumb and inevitably subjective. This may cause filtering of the signals to be too much or insufficient. Obtaining the true signals from the real signals contaminated by noise with unknown statistical features is always a challenging topic, and much work needs to be done on this in the future. In conclusion, through the calculation of the forecast entropies for selected BFs, a deterministic mechanism is apparently indicated in the evolution of silicon sequences. This may serve as a guideline for finding a novel tool, such as chaos, to predict and control the complex IBFP. Furthermore, the delay time and the number of parameters needed to construct the predictive model or the dynamical model of silicon content in BF hot metal is obtained and found to be independent of the size of the selected BFs. It is hoped that the obtained knowledge can give some help to investigate the complex IBFP both in theory and in practice. Appendix A The main purpose of defining an ideal random system (IRS) is to provide a standard of unpredictability by setting its forecast entropy as 1, the upper limit of forecast entropy. As an IRS, it is supposed to have a uniform distribution in its possible range Ndat is generated by of values.28 For example, if the series {xi}i)1 the IRS, xi can be any value in [xmin, xmax] in equal probability. Thus, if the area [xmin, xmax] is partitioned into Ndat equal bins, there is just a single point in every bin, like the pattern

where the digit 1 means the number of points in the bin. Hence, according to eq 1, the BGS entropy of the IRS is ln 2 at scale 1, 2 ln 2 at scale 2, ‚‚‚ , and 2m-1 ln 2 at scale m. Substituting these results into eq 2 yields the forecast entropy of IRS to be m

FIRS ) ln 2

ai2i-1 ) 1 ∑ i)1

(A-1)

The determination of the parameters ai(1 e i e m) is performed by considering the forecast entropy of another distribution having the pattern

Ind. Eng. Chem. Res., Vol. 47, No. 9, 2008 3043

Figure B1. Schematic projected diagram from d-dimensional space to two-dimensional.

which is half of FIRS. However, from the point of view of BGS entropy, the latter has the same BGS entropy as that of IRS at all other scales except scale m; Sm ) 0 for the latter. Thus, another expression for parameters ai is

0.5FIRS ) ln 2

m-1 ai2i-1 ∑i)1

(A-2)

Combining eqs A-1 and A-2 can yield

am )

1 m-1

2

m-1 ai2i-1 ∑i)1

(A-3)

Similar relations between ai and a1, a2, ‚‚‚, ai-1 (2 e i e m) can be obtained through considering the forecast entropy of the distributions

and

system onto the (yl, yl-τ) plane is the same as that on the (yi, yi-τ) plane, where l ) i - τ, i - 2τ, ‚‚‚, i - (d - 2)τ, so that the forecast entropies of the distributions on these planes are considered to be the same. Therefore, there are just (d - 1) planes needed projecting onto. The projected schematic diagram is exhibited in Figure B1. • Step 3. Calculate the angle of every neighboring point projected on the (yi, yi-jτ) plane. For example, Yk1 ) (yk1, yk1-τ, ‚‚‚, yk1-(d-1)τ), the angle is arctan{yk1+1 - yk1}/{yk1-jτ+1 - yk1-jτ} if yk1+1 - yk1 and yk1-jτ+1 - yk1-jτ are positive. If the sign of yk1+1 - yk1 or yk1-jτ+1 - yk1-jτ is minus, the angle must be adjusted to make it lie in the area [0, 2π]. • Step 4. Divide the area [0, 2π] into n equidistant bins with every distance of interval 2π/n, and count the numbers of points in every bin according to the angles calculated in Step 3. • Step 5. Calculate the local forecast entropy, Fkj, of the distribution of the angles using eq 5. • Step 6. Update k ) k + 1 and repeat steps 1-5 until k ) Ndat. • Step 7. Take the average of these local forecast entropies as Fj, i.e.,

F )

1

j

. Thus, the ultimate relations among these parameters are

am ) am-1 ) ‚‚‚ ) a3 ) a2 ) a1/2

1 2

m-1

ln 2

∑

Ndat - (d - 1)τ k)(d-1)τ+1

Fkj

(B-1)

(A-4)

Note that the norm adopted to search the neighbors in step 1 is the l∞ norm instead of the usual Euclidean norm to reduce the computation time.

(A-5)

Appendix C

Substituting these results into eq A-1 can yield a1 to be

a1 )

Ndat

Finally, eq 5 can be obtained from eqs 2, A-4, and A-5. Appendix B The calculation of Fj can be outlined as the following steps: 28

• Step 1. Find n nearest-neighboring points of Yk(k ) (d 1)τ + 1) in regular space, denoted as Yk1, Yk2, ‚‚‚, Ykn, respectively. • Step 2. Project these neighboring points onto some planes constructed by any two entries of the patterns appearing in eq 6, such as (yi, yi-τ), (yi, yi-2τ) plane, etc. In principle, there are in total d(d - 1)/2 planes needed to be projected onto. However, the phase diagram projected by the d-dimensional reconstructed

The consistency test of forecast entropy is implemented by considering the chaotic Lorenz system33

[

x˘ ) σ(y - x) y˘ ) x(r - z) - y z˘ ) xy - bz

]

(C-1)

which is a three-dimensional vector filed. Let the parameter values be σ ) 16.0, r ) 45.92, and b ) 4.0, for which the well-known “butterfly” attractor exists, the initial values be x0 ) -3.16, y0 ) -5.13, and z0 ) 13.31, and the step size be 0.01; then we can obtain a time series (such as x series) using the fourth-order Runge-Kutta method. To check the consistency of forecast entropy, we first calculate the forecast entropy of the obtained clean series under

3044

Ind. Eng. Chem. Res., Vol. 47, No. 9, 2008

yisignal ) measured silicon content in hot metal at the time of t0 + (i - 1)∆t yiclean ) actual silicon content in hot metal at the time of t0 + (i - 1)∆t d ) embedding dimension τ ) delay time Fd ) forecast entropy of d-dimensional reconstructed system n ) number of neighbors Literature Cited

Figure C1. Dependence of forecast entropy of Lorenz series contaminated by noise on the amount of data. The dashed part represents the forecast entropy of clean Lorenz series with the amount of data to be 1024.

the condition of reconstructed parameters d ) 3, τ ) 1, the number of neighbors n ) 16, and the amount of data Ndat ) 1 024 ) 210, partly motivated by Yao et al.28 and partly for the sake of convenience. The forecast entropy result is ∼0.0417 and shown in Figure C1 as the dashed part. This result is assumed to be the “true” forecast entropy value of the Lorenz series. Then, under the same condition as above, we calculate the forecast entropy of the Lorenz series contaminated by noise, which is a pseudorandom series with uniform distribution in [0, 1], and the signal noise rate is ∼30. Afterward, according to the idea of Monte Carlo, the amount of contaminated Lorenz series is continuously increased to Ndat ) 219 with the increment of last data size every time, and the same procedure of calculating forecast entropy is made for these series. Finally, all forecast entropy values of the contaminated Lorenz series at different data sizes are obtained. Figure C1 shows the variety of forecast entropies as the amount of data increases. Seen from Figure C1, the forecast entropy of contaminated Lorenz series will converge to 0.0417 as the amount of data increases. This indicates that the forecast entropy value of the contaminated Lorenz series will approach the “true value” as the data size increases. Therefore, the algorithm of forecast entropy can be thought to be approximately consistent. Acknowledgment This work has been financially supported by Zhejiang Provincial Natural Science Foundation of China under Grant No. Y107110, Research Fund for the Doctoral Program of Higher Education of China (for new teachers) under Grant No. 20070335161, Joint Funds of NSFCsGuangdong under Grant No. U0735003, and Natural Science Foundation of China under Grant No. 60604029. The authors would like to thank Dr. Weiguang Yao for his helpful discussion about the calculation of forecast entropy and the anonymous reviewers for many insightful comments. Nomenclature S ) Bolzmann-Gibbs-Shannon entropy Ndat ) size of time series t0 ) initial sampling time ∆t ) sampling interval m ) satisfying the equation of Ndat ) 2m F ) forecast entropy Si ) Bolzmann-Gibbs-Shannon entropy on scale i Ri ) weighting coefficient on scale i

(1) Tsuchiya, N.; Tokuda, M.; Ohtani, M. Kinetics of silicon transfer from a gas phase to molten iron. Tetsu to Hagane. 1972, 58, 1927. (2) Tsuchiya, N.; Tokuda, M.; Ohtani, M. The transfer of silicon from gas-phase to molten iron in blast-furnace. Metall. Trans. B 1976, 7, 315. (3) Ozturk, B.; Fruehan, R. J. Kinetics of the reaction of SiO(g) with carbon saturated iron. Metall. Trans. B 1985, 16, 121. (4) Ozturk, B.; Fruehan, R. J. The reaction of SiO (g) with liquid slags. Metall. Trans. B 1986, 17, 397. (5) Batra, N. K. Adiabatic flame temperature and silicon transfer in a blast furnace. Trans. Indian Inst. Met. 1997, 50, 169. (6) Sugawara, K.; Morimoto, K.; Sugawara, T.; Dranoff, J. S. Dynamic behavior of iron forms in rapid reduction of carbon-coated iron ore. AIChE J. 1999, 45, 574. (7) Liu, D. Y.; Wang, G. X.; Litster, J. D. Unsaturated liquid percolation flow through nonwetted packed beds. AIChE J. 2002, 48, 953. (8) Matsui, Y.; Mori, S.; Noma, F. Kinetics of silicon transfer from pulverized coal injected into blast furnace under intensive coal injection. ISIJ Int. 2003, 43, 997. (9) Gustavsson, J.; Andersson, A. M. T.; Jo¨nsson, P. G. A thermodynamic study of silicon containing gas around a blast furnace raceway. ISIJ Int. 2005, 45, 662. (10) Singh, H.; Sridhar, N. V.; Deo, B. Artificial neural nets for prediction of silicon content of blast furnace hot metal. Steel Res. 1996, 67, 521. (11) Radhakrishnan, V. R.; Mohamed, A. R. Neural networks for the identification and control of blast furnace hot metal quality. J. Process Control 2000, 10, 509. (12) Chen, J. A predictive system for blast furnaces by integrating a neural network with qualitative analysis. Eng. Appl. Artif. Intell. 2001, 14, 77. (13) Zheng, D. L.; Liang, R. X.; Zhou, Y.; Wang, Y. A chaos genetic algorithm for optimizing an artificial neural network of predicting silicon content in hot metal. J. UniV. Sci. Technol. Beijing 2003, 10, 68. (14) Helle, M.; Pettersson, F.; Chakraborti, N.; Saxe´n, H. Modeling noisy blast furnace data using genetic algorithms and neural networks. Steel Res. Int. 2006, 77, 75. (15) Pettersson, F.; Chakraborti, N.; Saxe´n, H. A genetic algorithms based multi-objective neural net applied to noisy blast furnace data. Appl. Soft Comput. 2007, 7, 387. (16) Waller, M.; Saxe´n, H. Application of nonlinear time series analysis to the prediction of silicon content of pig iron. ISIJ Int. 2002, 42, 316. (17) Bhattacharya, T. Prediction of silicon content in blast furnace hot metal using partial least squares (PLS). ISIJ Int. 2005, 45, 1943. (18) Hao, X. J.; Shen, F. M.; Du, G.; Shen, Y. S.; Xie, Z. A blast furnace prediction model combining neural network with partial least square regression. Steel Res. Int. 2005, 76, 694. (19) Martı´n, R. D.; Obeso, F.; Mocho´n, J.; Barea, R.; Jime´nez, J. Hot metal temperature prediction in blast furnace using advanced model based on fuzzy logic tools. Ironmaking Steelmaking 2007, 34, 241. (20) Waller, M.; Saxe´n, H. Time-varying event-internal trends in predictive modeling-methods with applications to ladlewise analyses of hot metal silicon content. Ind. Eng. Chem. Res. 2003, 42, 85. (21) Wales, D. J. Calculating the Rate of Loss of Information from Chaotic Time Series by Forecasting. Nature 1991, 350, 485. (22) Gao, C. H.; Zhou, Z. M.; Qian, J. X. Chaotic identification and prediction of silicon content in hot metal. J. Iron Steel Res. Int. 2005, 12, 3. (23) Gao, C. H.; Qian, J. X. Time-dependent fractal characteristics on time series of silicon content in hot metal of blast furnace. ISIJ Int. 2005, 45, 1269. (24) Zhou, Z. M. Measurement of time-dependent fractal dimension for time series of silicon content in pig iron. Physica A 2007, 376, 133. (25) Nightingale, R. J.; Dippenaar, R. J.; Lu, W. K. Developments in blast furnace process control at Port Kembla based on process fundamentals. Metall. Mater. Trans. B 2000, 31, 993.

Ind. Eng. Chem. Res., Vol. 47, No. 9, 2008 3045 (26) Taken, F. Detecting strange attractors in fluid turbulence. Lect. Notes Math. 1981, 898, 336. (27) Shiner, J. S.; Davison, M.; Landsberg, P. T. Simple measure for complexity. Phys. ReV. E 1999, 59, 1459. (28) Yao, W. G.; Essex, C.; Yu, P.; Davison, M. Measure of predictability. Phys. ReV. E 2004, 69, 066121. (29) Shannon, C. E. A mathematical theory of communication. Bell. Syst. Tech. J. 1948, 27, 623. (30) Packard, N. H.; Crutchfield, J. P.; Farmer, J. D.; Shaw, R. S. Geometry from a time series. Phys. ReV. Lett. 1980, 45, 712.

(31) Sauer T.; Yorke J. A.; Casdagli M. Embedology. J. Stat. Phys. 1991, 65, 579. (32) Schreiber, T. Extremely simple nonlinear noise-reduction method. Phys. ReV. E. 1993, 47, 2401. (33) Lorenz, E. N. Deterministic nonperiodic flows. J. Atmos. Sci. 1963, 20, 130.

ReceiVed for reView June 27, 2007 ReVised manuscript receiVed February 4, 2008 Accepted February 4, 2008 IE070879S