Online Prediction of Hot Metal Temperature Using Transient Model and Moving Horizon Estimation

Abstract

Precise control of hot metal temperature (HMT) is crucial for achieving stable operation of a blast furnace, but it is difficult due to the sluggish dynamics caused by the huge heat capacity. To cope with such difficulty, this work aims at developing a method that can predict future HMT by adopting moving horizon estimation (MHE) based on a one-dimensional transient model. MHE is useful to successively adjust model parameters so that the undesirable influence of past disturbances on the prediction is minimized. The real application result demonstrated that the root mean square error (RMSE) of HMT of eight-hour-ahead prediction was only 11.6°C. The high-performance prediction enables operators to realize the efficient operation of the blast furnace.

1. Introduction

In the steel industry, a blast furnace smelts and reduces iron ore while consuming huge energy. The heat indices of the furnace, such as the hot metal temperature (HMT), are controlled by manipulating the coke ratio at the furnace top, the blast moisture, and the pulverized coal (PC) flow rate. Hot metal and its viscous byproduct, i.e. the slag, are drained from tap holes. HMT needs to be controlled because the slag drainage becomes extremely difficult when the heat level becomes too low; on the contrary, the problem of enormous fuel consumption and CO₂ emission occurs when the heat level becomes excessively high.

The huge heat capacity of blast furnace leads to a long response time of approximately 12 hours. Such slow dynamics makes blast furnace operation difficult. Operators frequently take excessive actions, which cause variation of the HMT. To reduce the HMT fluctuation, it is useful to provide the operators with the information on the future prediction of the furnace state. Hence, the one-dimensional (1D) transient model was developed to predict the future HMT in the previous work.¹⁾

Various blast furnace models have been developed: for example, 1D steady-state model,²⁾ 1D transient model,³⁾ 2D steady-state model,^4,5,6,7) and 3D transient model.⁸⁾ When these models are used for process control, estimation errors are caused by changes in the iron content of iron ore and the carbon content of pulverized coal. Stochastic phenomena inside the furnace such as the fluctuation of burden descent and the gas flow path make the prediction of furnace state even more difficult. These disturbances also deteriorate the prediction accuracy of HMT. To decrease the estimation errors, it is effective to update the model parameters and boundary conditions based on the measurement of process variables such as the composition of the furnace top gas.

There have been numerous state estimation methods, which can be classified into two types: the sequential type such as Kalman filter⁹⁾ and particle filter,¹⁰⁾ and the non-sequential type such as moving horizon estimation (MHE)^11,12) and the adjoint method.¹³⁾ The former modifies the model parameters to coincide the calculated process variables with the present actual data, whereas the latter changes the model parameters in the past to reproduce the actual time-series data.

In the case of the blast furnace with the large time constant, the current estimation errors stem from the accumulation of disturbances in the past. Hence, it is better to adjust model parameters and boundary conditions so that the trajectory of the estimated variables matches with actual data. For this reason, we adopted the MHE approach, which updates the model parameters in the past and features the simple algorithm.

Appropriate selection of model parameters to be adjusted by MHE is an issue. In this study, the fundamental cause of the estimation errors was extracted using principal component analysis (PCA) and the Rist diagram.^14,15,16) PCA is a well-known multivariate analysis method, and the Rist diagram graphically organizes the mass and heat balance of blast furnace and helps users to easily understand its behavior. The analysis based on PCA and Rist diagram helps to select model parameters to be adjusted for effectively reducing the estimation errors.

Figure 1 shows the outline of the future HMT prediction using the 1D transient model and MHE. The solid line indicates the online calculation, which is conducted successively. The dashed line shows the offline calculation, which is carried out only once. The rest of this paper is organized as follows. Section 2 describes the 1D transient model and the key process variables that are used for evaluating the estimation errors caused by the disturbances. Section 3 compares the key process variables estimated by the transient model with the actual ones measured in the real plant. The estimation errors are analyzed by the PCA and Rist diagram, and the appropriate model parameters to be adjusted by MHE are selected. Section 4 presents the MHE algorithm and confirms that the estimation errors are successfully reduced by MHE. This section also evaluates the prediction accuracy of eight-hour-ahead HMT by using 1D transient model and MHE, and it is shown that adjusting the model parameters in the past is necessary to predict the future HMT accurately. Concluding remarks are given in section 5.

Fig. 1.

Outline of HMT prediction using 1D transient model and MHE. (Online version in color.)

2. Transient Model and Key Process Variables

In this section, the 1D transient model¹⁾ is briefly described, and the simulation results are presented to show the internal states and the time constant of the process. Then, the key process variables are explained.

2.1. Transient Model and Simulation

The 1D transient model involves four sub-models: a gas flow model, a reaction model, an energy balance model, and a solid flow model. 12 reactions such as the CO gas reduction and the coke gasification are taken into account. Manipulated variables such as the flow rate of pulverized coal and the coke ratio at furnace top were set as the boundary conditions. The variables in the 1D transient model are summarized in Table 1.

Table 1. Variables in 1D transient model.

Symbol	Notes	Dimension
Cp,c	Specific heat of coke	[J/kg·K]
Cp,g	Specific heat of gas	[J/kg·K]
Cp,fe	Specific heat of iron	[J/kmolFe·K]
Cp,sl	Specific heat of slag (gangue)	[J/kg·K]
E	Heat exchange coefficient	[W/m3·K]
h	Overall heat transfer coefficient	[W/m3·K]
R	Reaction rate	[kmol/m3bed·s]
ug	Mass velocity of gas	[kg/m2·s]
ufe	Molar flux of liquid iron	[kmolFe/m2·s]
usl	Mass velocity of slag	[kg/m2·s]
q	Heat-loss through furnace wall	[W/m3]
T	Temperature	[K]
Xc	Volume ratio of coke	[m3coke/m3bed]
Xo	Volume ratio of ore	[m3ore/m3bed]
ΔHR	Reaction heat	[J/kmol]
η	Distribution ratio of reaction heat	[−]
ρc	Apparent coke density	[kgcoke/m3coke]
ρfe	Iron density in iron ore	[kmolFe/m3ore]
ρsl	Gangue density in iron ore	[kg/m3ore]

Eqations (1), (2), (3) show the energy balance model that considers the reaction heat and heat exchange between the solid and the gas.   

$$\begin{array}{c}\nabla \cdot (C_{p\text{,g}}{T}_{\text{g}}{\mathbf{u}}_{\text{g}})=\sum _j{\eta }_{{\text{g}}_j}\mathrm{\Delta }{H}_{{\text{R}}_j}{R}_j+E_{\text{g,fe}}(T_{\text{fe}}-T_{\text{g}})\\ +E_{\text{g,c}}(T_{\text{c}}-T_{\text{g}})+q\end{array}$$(1)

  

$$\begin{array}{l}{\displaystyle \frac{\partial ({\rho }_{\text{fe}}{C}_{p\text{,fe}}{X}_{\text{o}}{T}_{\text{fe}})}{\partial t}}+\nabla \cdot (C_{p\text{,fe}}{T}_{\text{fe}}{\mathbf{u}}_{\text{fe}})+{\displaystyle \frac{\partial ({\rho }_{\text{sl}}{C}_{p,\text{sl}}{X}_{\text{o}}{T}_{\text{fe}})}{\partial t}}\\ +\nabla \cdot (C_{p\text{,sl}}{T}_{\text{fe}}{\mathbf{u}}_{\text{sl}})=\\ \sum _j{\eta }_{{\text{fe}}_j}\mathrm{\Delta }{H}_{{\text{R}}_j}{R}_j+E_{\text{g,fe}}(T_{\text{g}}-T_{\text{fe}})+E_{\text{fe,c}}(T_{\text{c}}-T_{\text{fe}})\end{array}$$(2)

  

$$\begin{array}{c}{\displaystyle \frac{\partial ({\rho }_{\text{c}}{C}_{p\text{,c}}{X}_{\text{c}}{T}_{\text{c}})}{\partial t}}=\sum _j{\eta }_{{\text{c}}_j}\mathrm{\Delta }{H}_{{\text{R}}_j}{R}_j+E_{\text{g,c}}(T_{\text{g}}-T_{\text{c}})\\ +E_{\text{fe,c}}(T_{\text{fe}}-T_{\text{c}})\end{array}$$(3)

where subscripts g, fe, sl, and c denote gas, iron, slag (gangue), and coke. The subscript j denotes the indices of the reactions. q is the heat flux through the furnace wall, i.e. heat-loss, described by   

$$q=-h(T_{\text{g}}-T_{\text{w}}),$$(4)

where subscript w denotes the cooling water. This heat sink is incorporated in Eq. (1) as a source term.

The heat exchange coefficient E is the product of the heat transfer coefficient and the heat transfer area. The Ranz’s equation¹⁷⁾ was used to calculate the heat transfer coefficient between the gas and the solid/liquid, regarding the liquid iron/slag as the spherical droplet whose diameter is 0.015 m. The heat exchange area between the gas and liquid was estimated by Mada’s equation.¹⁸⁾ However, the equations overestimated HMT; therefore factor 0.1 was multiplied to the heat transfer coefficient between the gas and liquid so that the calculated HMT coincides with actual HMT. The heat exchange between the coke and the iron ore was neglected, while that between the coke and the liquid iron/slag was estimated based on the Taguchi’s work.¹⁹⁾

The convergence calculation by the implicit method was conducted at each time step with the static calculation cells; that is, the gas flow, the reaction, and the heat transfer are calculated while the solid flow is suspended. Hence the convection terms regarding the solid flow do not appear in Eqs. (2) and (3), while the convection heat transfer by the liquid flow is included in Eq. (2). The time step of the calculation was set to 5 minutes, considering that the period of burden loading is between 10 and 20 minutes. The heat transfer involving the gas flow in Eq. (1) was regarded steady-state because the residence time of gas in the furnace is within 3 seconds, which is much shorter than the time step. The transient terms in Eqs. (2) and (3) are necessary since the residence time of solid is more than 6 hours, which is longer than the time step. After the convergence, the solid moves downwards by updating the calculation cells at each time step so that each cell corresponds to one set of ore layer and coke layer.

When the process variables are calculated repeatedly by the transient model with the constant manipulated variables, they converge to the periodic steady state (hereafter referred to as the steady state). Figure 2 shows a steady-state simulation result of the temperature and the gas utilization ratio. The zero point of the vertical axis corresponds to the tuyere level. Input conditions are listed in Table 2.

Fig. 2.

Calculated steady-state profiles of (a) temperature and (b) gas utilization ratio in height direction above the tuyere level. (Online version in color.)

Table 2. Input condition of 1D transient model.

Manipulated variable	Value
Blast volume	6770 Nm3/min
Coke ratio	367 kg/t-pig
PC flow rate	960 kg/min
Blast moisture	26 g/Nm3
Blast temperature	1133°C
Enrichment oxygen	284 Nm3/min
Top gas pressure	338 kPa

Figure 2(a) shows the calculated temperature profile of the furnace. The iron ore at low temperature is loaded from the furnace top and heated through the heat exchange with the gas. The simulation result indicates that a thermal reserve zone^20,21) is located between 6 and 13 m height. In addition, another thermal reserve zone can be observed between 15 and 18 m height, where the rising rate of the burden temperature decreases by the endothermic reduction of Magnetite.

Figure 2(b) shows the profile of the gas utilization ratio, which is the conversion ratio of CO gas into CO₂ gas. The chemical reserve zone is located between 8 and 13 m height, where the reduction reaction becomes slow and the gas composition reaches the equilibrium between Fe and FeO_1.05. Figure 3 shows the dependence of the CO/(CO+CO₂) ratio on the iron temperature. The CO/(CO+CO₂) ratio agrees with the equilibrium line between Fe and FeO_1.05 around 1000°C.

Fig. 3.

Dependence of gas composition on iron temperature. (Online version in color.)

To calculate the transient state of the process, the step responses of HMT to the three important manipulated variables for controlling HMT, i.e. the coke ratio, the blast moisture, and the PC flow rate, were investigated. A step-wise change was given to each manipulated variable after the steady state was achieved under the condition in Table 2, while the other manipulated variables were kept constant.

Figure 4(a) shows the step response of HMT to the change of the coke ratio, where the vertical axes indicate the difference from their steady-state values before the change of coke ratio. HMT increased with 6 hours delay, which corresponds to the burden descent time from the furnace top to the tuyere level. Figure 4(b) shows the step response of HMT to the blast moisture; the fast response of HMT is confirmed. Figure 4(c) shows the step response of HMT to the PC flow rate; the temperature rose gradually compared with the case of blast moisture.

Fig. 4.

Calculated step responses of HMT to three important manipulated variables: (a) coke ratio, (b) blast moisture, and (c) PC flow rate. (Online version in color.)

2.2. Key Process Variables

The process variables representing the reactions in a blast furnace need to be calculated precisely because they strongly affect the mass balance and the heat balance in the furnace.

The CO gas utilization ratio η_CO [–], which indicates the conversion ratio of CO into CO₂, is defined by   

$${\eta }_{\text{CO}}=\frac{X_{{\text{CO}}_{\text{2}}}}{X_{\text{CO}}+X_{{\text{CO}}_{\text{2}}}}$$(5)

where X_ϕ denotes the volume ratio of component ϕ in the furnace top gas. High η_CO means that the CO gas is effectively used to reduce iron ore.

The production rate (Prod) [t/min] can be calculated by the oxygen balance.   

$$Prod=\frac{V_{\text{out}}^{\text{O}}-V_{\text{in}}^{\text{O}}}{O_{\text{ore}}}$$(6)

where $V_{\text{out}}^{\text{O}}$ [kmol/min] is the flow rate of oxygen from the furnace top, $V_{\text{in}}^{\text{O}}$ [kmol/min] is the oxygen input rate to the tuyere, and $O_{\text{ore}}$ [kmol/t] is the oxygen amount in unreduced iron ore which has to be reduced to produce unit volume of hot metal. $V_{\text{out}}^{\text{O}}$ is calculated from measurable variables as follows.   

$$V_{\text{out}}^{\text{O}}=V_{\text{top}}(X_{\text{CO}}+2X_{{\text{CO}}_{\text{2}}}+X_{{\text{H}}_{\text{2}}\text{O}})$$(7)

where V_top [kmol/min] is the flow rate of top gas.

The reducing agent rate (RAR) [kg/t], which is the sum of PC rate and coke rate per iron, is evaluated by   

$$RAR=\frac{V_{\text{PC}}+V_{\text{coke}}}{Prod},$$(8)

where V_PC [kg/min] is the PC flow rate. The consumption rate of the coke V_coke [kg/min] can be estimated on the basis of the carbon balance.   

$$V_{\text{coke}}=\frac{V_{\text{out}}^{\text{C}}-V_{\text{in}}^{\text{C}}}{C_{\text{coke}}}$$(9)

where $V_{\text{out}}^{\text{C}}$ [kg/min] is the outflow rate of carbon from the furnace, $V_{\text{in}}^{\text{C}}$ [kg/min] is the carbon input rate to the tuyere, and C_coke [–] is the carbon content of the coke. $V_{\text{out}}^{\text{C}}$ depends on the reaction rates of the coke gasification and the carburization as follows.   

$$V_{\text{out}}^{\text{C}}=V_{\text{top}}(X_{\text{CO}}+X_{{\text{CO}}_{\text{2}}})+Prod\cdot [C]$$(10)

where [C] [kg/t] is the carbon content in the hot metal.

The solution loss carbon (SLC) [kg/t] is the amount of carbon consumed by coke gasification reaction, which is an endothermic reaction involving large heat absorption. SLC is also calculated from the carbon balance.   

$$SLC=\frac{V_{\text{out}}^{\text{C}}-V_{\text{in}}^{\text{O}}\cdot \raisebox{1ex}{$12$}\!\left/ \!\raisebox{-1ex}{$16$}\right.}{Prod}-[C]$$(11)

The process variables, e.g. gas composition of the top gas, top gas flow rate, and PC flow rate, are measured every 30 minutes in the actual blast furnace. Based on the measurements or the estimates of the process variables, η_CO, SLC, Prod, and RAR are obtained through Eqs. (5), (6), (7), (8), (9), (10), (11).

HMT and the four variables mentioned above, i.e. η_CO, SLC, Prod, and RAR are defined as the key process variables in this work. HMT is closely connected with RAR and SLC, since the heat balance per iron in the lower furnace is approximated by   

$$q_{\text{blast}}+RAR\cdot q_{\text{c}}\approx SLC\cdot (q_{\text{sol}}+q_{\text{c}})+q_{\text{HM}}+q_{\text{slag}}+q_{\text{heatloss}}+\zeta ,$$(12)

where q_blast [kJ/t] is the sensible heat of hot blast, q_c [kJ/kg] is the reaction heat when the carbon is combusted in the raceway, q_sol [kJ/kg] is the heat absorption of coke gasification reaction, q_HM [kJ/t] is the sensible heat of hot metal, q_slag [kJ/t] is the sensible heat of slag, and q_heatloss [kJ/t] is the heat loss through the furnace wall. The term ζ includes the decomposition heat of blast moisture and the pulverized coal. It also contains the reaction heat of carburization and reduction of metalloids. HMT in the transient model corresponds to the iron temperature at the tuyere level, whereas HMT in the actual furnace is measured after the hot metal is drained from tap holes. The sampling period of actual HMT varies from 15 minutes to 40 minutes, and its 30-minute average was used in this work.

3. Cause of Estimation Errors

This section provides the validation result of the transient model with actual data. The cause of the estimation errors of the key process variables is investigated with PCA and the Rist diagram. Based on the result, the model parameters to compensate the estimation errors are identified.

3.1. Estimation Results

The 1D transient model was validated by comparing the calculated process variables including HMT and the real ones. The temperature distribution, the oxidization degree of iron, and the compositions of the gas and the hot metal were calculated by running the model with the manipulated variables shown in Fig. 5, where vertical axes show the deviation from the mean values. The calculation results of the key process variables and the corresponding actual data are shown in Fig. 6, where vertical axes show the deviation from the mean values of actual data. The estimation performance was evaluated by RMSE (root mean square error). Although complicated phenomena involving the liquid flow are known to occur in the hearth region,²²⁾ the 1D transient model does not take into account the residence time of the liquid in the hearth region to make the model as simple as possible. To simplify the model for its online use, the calculated HMT is shifted by 2 hours, assuming that the residence time is 2 hours. Although estimates show the similar trend to measurements, some discrepancy between them exists during day 1 to 7 and during day 22 to 28.

Fig. 5.

Actual time-series data of manipulated variables. (Online version in color.)

Fig. 6.

Comparison between actual values and calculated values of key process variables. (Online version in color.)

3.2. Rist Diagram

In this work, the Rist diagram is used for analyzing the cause of estimation errors, therefore it is briefly explained here. Figure 7 outlines the Rist diagram, which organizes the mass and heat balance of a blast furnace. The Rist diagram shows the connections among the key process variables.

Fig. 7.

Rist diagram.

The horizontal axis (O/C) represents the degree of oxidization of carbon, which is 1 for CO gas and 2 for CO₂. The upper part of the vertical axis (O/Fe) represents the degree of oxidization of Fe, which is 1.5 for unreduced iron ore (FeO_1.5), 1.05 for partially reduced Wustite (FeO_1.05), and 0 for reduced iron (Fe). The slope of segment AE represents RAR. The point W is the fixed point where the oxidization degree of iron is 1.05, and the oxidization degree of carbon is corresponding to the CO/CO₂ ratio at the equilibrium between Fe and FeO_1.05. The reduction efficiency HR/HW is the degree of reaching the equilibrium.

The segment TV includes the heat output from the lower furnace: the sensible heat of the hot metal and the slag, and the heat loss. The segment UV must pass the point P, which is the intersection of the segments SQ and AE, under the constraint of the heat balance in Eq. (12).

To improve the estimation accuracy of HMT, the behavior of point V in the actual furnace has to be reproduced by the model. Since the segment SQ is determined by blast temperature, the point V is uniquely determined by the segment AE. Hence, it is necessary to calculate the movement of the segment AE accurately.

The segment AE has two degrees of freedom. When the two of four variables, i.e. η_CO, SLC, Prod, and RAR are fixed, the rest two are determined. These constraints are useful when analyzing the cause of the estimation errors.

3.3. Interpretation of Estimation Errors with PCA and Rist Diagram

Figure 8 shows the correlation between the estimation errors of the key process variables that are shown in Fig. 6. ΔX denotes the difference between the calculated values and the actual values of variable X; calculated minus actual. The strong correlations in Figs. 8(a) 8(b) 8(c) derive from the fact that only two degrees of freedom exist among the four variables: η_CO, SLC, RAR, and Prod.

Fig. 8.

Correlation between estimation errors of key process variables. (Online version in color.)

The errors of SLC and RAR cause that of HMT, as indicated by Eq. (12). When the calculated RAR is larger than the actual RAR, the calculated HMT becomes higher than the actual HMT because more heat source per iron is consumed in the calculation. When the calculated SLC is larger than the actual SLC, the calculated HMT gets lower than the actual HMT because more heat absorption occurs in the lower furnace in the calculation.

Since the estimation errors of key process variables have strong correlations, PCA was conducted to extract the root cause behind the correlations. The principal components, which best describe the data distribution in multi-dimensional space, can be obtained by PCA. First, the estimation error of each key process variable was normalized, and the data matrix $\tilde{\mathrm{X}}$ was prepared as follows.   

$${\displaystyle \begin{array}{l}\tilde{\mathrm{X}}=\\ \left[\begin{array}{ccccc}\mathrm{\Delta }{\eta }_{\text{CO}}(1)/{\sigma }_{\mathrm{\Delta }{\eta }_{\text{CO}}}& \mathrm{\Delta }SLC(1)/{\sigma }_{\mathrm{\Delta }SLC}& \mathrm{\Delta }RAR(1)/{\sigma }_{\mathrm{\Delta }RAR}& \mathrm{\Delta }Prod(1)/{\sigma }_{\mathrm{\Delta }Prod}& \mathrm{\Delta }HMT(1)/{\sigma }_{\mathrm{\Delta }HMT}\\ \mathrm{\Delta }{\eta }_{\text{CO}}(2)/{\sigma }_{\mathrm{\Delta }{\eta }_{\text{CO}}}& \mathrm{\Delta }SLC(2)/{\sigma }_{\mathrm{\Delta }SLC}& \mathrm{\Delta }RAR(2)/{\sigma }_{\mathrm{\Delta }RAR}& \mathrm{\Delta }Prod(2)/{\sigma }_{\mathrm{\Delta }Prod}& \mathrm{\Delta }HMT(2)/{\sigma }_{\mathrm{\Delta }HMT}\\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \mathrm{\Delta }{\eta }_{\text{CO}}(N-1)/{\sigma }_{\mathrm{\Delta }{\eta }_{\text{CO}}}& \mathrm{\Delta }SLC(N-1)/{\sigma }_{\mathrm{\Delta }SLC}& \mathrm{\Delta }RAR(N-1)/{\sigma }_{\mathrm{\Delta }RAR}& \mathrm{\Delta }Prod(N-1)/{\sigma }_{\mathrm{\Delta }Prod}& \mathrm{\Delta }HMT(N-1)/{\sigma }_{\mathrm{\Delta }HMT}\\ \mathrm{\Delta }{\eta }_{\text{CO}}(N)/{\sigma }_{\mathrm{\Delta }{\eta }_{\text{CO}}}& \mathrm{\Delta }SLC(N)/{\sigma }_{\mathrm{\Delta }SLC}& \mathrm{\Delta }RAR(N)/{\sigma }_{\mathrm{\Delta }RAR}& \mathrm{\Delta }Prod(N)/{\sigma }_{\mathrm{\Delta }Prod}& \mathrm{\Delta }HMT(N)/{\sigma }_{\mathrm{\Delta }HMT}\end{array}\right]\end{array}}$$(13)

ΔX and σ_ΔX are the estimation error of variable X and its root mean square error, respectively, and N is the number of data points. σ_ΔX was calculated by   

$${\sigma }_{\mathrm{\Delta }X}=\sqrt{\sum _{i=1}^N\mathrm{\Delta }X{(i)}^2/N}.$$(14)

The principal components are the eigenvectors of the covariance matrix ${\tilde{\mathbf{X}}}^T\tilde{\mathbf{X}}/(N-1)$. Figure 9 shows the first two principal components; the vertical axis indicates the loadings, i.e. the elements comprising the principal components.

Fig. 9.

Loadings of 1st and the 2nd principal components derived from the estimation errors of key process variables. (Online version in color.)

The first principal component shows the positive loadings of RAR and HMT and the negative ones of η_CO and Prod. When the calculated RAR becomes higher than the actual one, η_CO and Prod in the model become lower than the actual ones. Then, it causes the increase of calculated HMT compared with the actual one. Figure 10(a) shows the interpretation of the first principal component with the Rist diagram, where the increase of the slope of AE in the calculation moves point A to the right, and moves point B, point E, and point V downwards. Although the Rist diagram has been originally used for the steady-state analysis, it is useful to obtain the qualitative interpretation of the principal components under the transient states.

Fig. 10.

Interpretation of 1st and 2nd principal components with Rist diagram. (Online version in color.)

The second principal component shows the positive loadings of η_CO and HMT and the negative one of SLC, while that of RAR is nearly zero. When the calculated reduction efficiency increases under the condition of constant RAR, η_CO increases and SLC decreases in the model, and it leads to the increase of the calculated HMT compared with the actual one. Figure 10(b) shows the interpretation of the second principal component with the Rist diagram, where the downward movement of the segment AE, while keeping its slope moves point A to the right, and moves point B, point E, and point V downwards.

3.4. Model Parameters to Reduce Estimation Errors

It was found that the estimation errors can be explained by the first two principal components. Hence, it is necessary to identify the adjustable model parameters which correspond to these principal components so as to reduce the estimation errors by MHE.

The cause of the fluctuation of RAR, which seems to correspond to the first principal component, is investigated. RAR depends on the carbon content of the coke and the pulverized coal, and the iron content of the iron ore. These values are not measured continuously and subject to analysis errors. Since the coke is the dominant reducing agent in blast furnaces, the coke ratio at the furnace top is multiplied by the correction parameter α to take the effect of these disturbances into consideration.

Figure 11(a) shows the response of the key process variables when α is changed from 1.0 to 1.01. The solid line in Fig. 12 shows the key process variables at 72 hours in Fig. 11(a). All variables were normalized by the RMSE of each variable, which was obtained by Eq. (14). Since a similar pattern with the first principal component in Fig. 9 was obtained, α is one of the effective model parameters to reduce the estimation errors.

Fig. 11.

Responses of key process variables to (a) coke ratio correction parameter, (b) reduction equilibrium parameter, and (c) heat-loss parameter. (Online version in color.)

Fig. 12.

Responses of the key process variables to coke ratio correction parameter, reduction equilibrium parameter, and heat-loss parameter at steady state. (Online version in color.)

The cause of the fluctuation of the reduction efficiency, which seems responsible for the second principal component, is considered. Figure 13 shows the correlation between RAR and η_CO, where all variables are mean centered. The variance of actual η_CO under the condition of constant RAR is larger than the calculated one. Stochastic phenomena in the actual furnace, such as the change of gas flow path and the variation in the reducibility of iron ore, may affect the reaction speed of the iron ore reduction by CO gas.

Fig. 13.

Correlation between RAR and η_CO: (a) calculated data and (b) actual data. (Online version in color.)

To take the stochastic phenomena which fluctuate the reduction efficiency into the model, a parameter which adjusts the equilibrium of iron ore reduction is introduced. It was assumed that the equilibrium of the CO gas reduction is described by the equations in Table 3, where β is a model parameter. The reason why we selected not the reaction rate coefficient but the equilibrium as the parameter is as follows. If the reaction rate coefficient is decreased to reproduce the drop of the reduction efficiency, CO/(CO+CO₂) around 1000°C in Fig. 3 deviates from the equilibrium line. This deviation causes the chemical reserve zone to disappear. Since the existence of chemical reserve zone is suggested in blast furnaces, the equilibrium is used as the parameter to reproduce the fluctuation of the reduction efficiency while maintaining the characteristics of the model that involve the chemical reserve zone.

Table 3. Equilibrium of CO gas reduction reaction.

Substance	Equilibrium
FeO1.33/Fe	${\left(\frac{P_{{\text{CO}}_{\text{2}}}}{P_{\text{CO}}}\right)}_{\text{eq}\text{.}}=exp\left(\frac{603+500(\beta -1)}{T}-0.862\right)$
FeO1.33/FeO1.05	${\left(\frac{P_{{\text{CO}}_{\text{2}}}}{P_{\text{CO}}}\right)}_{\text{eq}\text{.}}=exp\left(\frac{-3\mathrm{\hspace{0.17em} }900+500(\beta -1)}{T}+4.464\right)$
FeO1.05/Fe	${\left(\frac{P_{{\text{CO}}_{\text{2}}}}{P_{\text{CO}}}\right)}_{\text{eq}\text{.}}=exp\left(\frac{-2\mathrm{\hspace{0.17em} }104+500(\beta -1)}{T}-2.657\right)$

The gain 500 in the equations of Table 3 was determined so that the variations of the equilibrium which have been reported^23,24,25) can be explained by changing β from 0.5 to 1.5. Figure 14 shows the variation of equilibrium lines when β is 0.8, 1.0, and 1.2. The adjustment of β is valid because it is well known that the equilibrium depends on the type of calcium ferrite^24,25) contained in the sintered iron ore, which is affected by the sintering process.²⁶⁾ Figure 11(b) shows the response of the key process variables when β is changed from 1.0 to 1.1. The dashed line in Fig. 12 shows the key process variables at 72 hours in Fig. 11(b), which is similar to the second principal component in Fig. 9.

Fig. 14.

Equilibrium of CO gas reduction reaction with different reduction equilibrium parameter β. (Online version in color.)

In addition to the fluctuation of RAR and reduction efficiency, the heat-loss through the furnace wall is known to change depending on the gas flow path. The gas flow tends to concentrate near the wall region and the biased flow causes the increase in the heat-loss, especially when unburned char is accumulated in the dead-man due to the excessive pulverized coal injection.²⁷⁾ The heat-loss is multiplied by the factor γ to incorporate this disturbance into the model. Figure 11(c) shows the response of the key process variables when γ is increased from 1.0 to 1.1. The key process variables at 72 hours are shown in Fig. 12 with the dotted line.

From the viewpoint of the HMT estimation, the relationship between RAR and HMT is crucial since RAR indicates the heat source per unit mass of iron. The increase of γ involves the rise in RAR and the decrease in HMT, whereas the increase of α causes the rise in RAR and HMT. The response of RAR at steady state is nearly zero in the case of β. Hence, γ is an effective parameter to explain the correlation between the estimation error of RAR and that of HMT.

4. Moving Horizon Estimation

In this section, the calculation procedure and results of MHE are presented. In the present work, the three model parameters, i.e. the coke ratio correction parameter α, the reduction equilibrium parameter β, and the heat-loss parameter γ, are successively updated so that the calculated key process variables coincide with the actual ones. After we confirm that MHE can decrease the estimation errors of the key process variables, the eight-hour-ahead prediction of HMT is performed to justify adjusting the model parameters in the past for the accurate prediction.

4.1. Calculation Procedure of MHE

In this work, the calculation process with the transient model is expressed by the update of the state variables   

$$\mathbf{x}\left(t+1\right)=f\left(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right)\right),$$(15)

where x(t) and u(t) are the state variables and the manipulated variables at time step t, respectively. The nonlinear function f is the developed 1D transient model. One time step is 30 minutes, which was determined by taking account of the process dynamics and calculation time. There are 7 manipulated variables listed in Table 2 and 15 state variables including the temperature, the oxidation degree of iron, and the gas composition in each cell. The number of cells is dynamically updated within the range between 30 and 60; thus, the number of states is 15 times 30 (=450) through 15 times 60 (=900).

The observation operator is defined as   

$$\mathbf{y}(t)=C\left(\mathbf{x}(t)\right),$$(16)

where y(t) = [y₁(t), y₂(t), y₃(t), y₄(t), y₅(t)]^T ≡ [η_CO(t), SLC(t), RAR(t), Prod(t), HMT(t)]^T is the key process variables at time step t, and C is the function to obtain the key process variables from the state variables as described in section 2.2.

Figure 15 shows the procedure of MHE, which is comprised of 6 steps as described below.

Fig. 15.

Outline of moving horizon estimation. (Online version in color.)

Step 0: Set α=1, β=1, γ=1, and x(t−A)=x₀, where x₀ is the initial condition of the state variables.

Step 1: Calculate the state variables and the key process variables during the past A time steps with Eqs. (17) and (18).   

$$\mathbf{x}(t-k+1)=f\left(\mathbf{x}(t-k),\mathbf{u}(t-k)\right)$$(17)

  

$$\mathbf{y}(t-k+1)=C\left(\mathbf{x}(t-k+1)\right)$$(18)

where k takes the values from A to 1. Actual data of manipulated variables are used, and the model parameters are kept constant.

Step 2: Preserve x(t−A+1), which is used as the initial condition at the next iteration.

Step 3: Derive the estimation errors of the key process variables during the past A time steps.   

$$\mathbf{e}(t-k)={\mathbf{y}}_{\text{act}}(t-k)-\mathbf{y}(t-k)$$(19)

where the suffix act means the actual values.

Step 4: Calculate the correction amount of the model parameters, i.e. Δα, Δβ, and Δγ, so that the estimation errors during the past A/2 steps are expressed by the combination of the responses, which are shown in Fig. 11, as accurately as possible. In other words, the sum of squared errors is minimized.   

$$\begin{array}{l}\underset{\mathrm{\Delta }\alpha ,\mathrm{\Delta }\beta ,\mathrm{\Delta }\gamma }{min}\sum _{q=1}^5\sum _{n=0}^{\frac{A}{2}}(e_q\left(t-\frac{A}{2}+n\right)-R_q^{\alpha }\left(\frac{A}{2}+n\right)\mathrm{\Delta }\alpha \\ {\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }-R_q^{\beta }\left(\frac{A}{2}+n\right)\mathrm{\Delta }\beta -R_q^{\gamma }\left(\frac{A}{2}+n\right)\mathrm{\Delta }\gamma )}^2\end{array}$$(20)

where e_q(t) denotes the estimation error of the q-th key process variable, and $R_q^p(s)$ is the response of the q-th key process variable to the model parameter p at time step s; it is just like the s-th step response coefficient.

Step 5: Update the model parameters.   

$$\alpha \leftarrow \alpha +\mathrm{\Delta }\alpha$$(21)

  

$$\beta \leftarrow \beta +\mathrm{\Delta }\beta$$(22)

  

$$\gamma \leftarrow \gamma +\mathrm{\Delta }\gamma$$(23)

Step 6: Update the time step; t←t+1. Go to step 1.

In the present work, the parameter A is set to 144 time steps, i.e. 72 hours, because Fig. 11 shows that it takes about 3 days to reach the new steady-states. Further discussion on how to determine A is described in section 4.3. The x₀ in step 1 was set to the steady-state values of state variables under the typical operating condition shown in Table 2. $R_q^{\alpha }(s)$, $R_q^{\beta }(s)$, and $R_q^{\gamma }(s)$ for 144 time steps are shown in Figs. 11(a), 11(b), and 11(c), respectively. It takes more than 36 hours for the key process variables to settle the new steady states when the model parameters are changed; that is, $R_q^p(s)$ for the first 36 hours is smaller than the steady-state value. Hence, the correction amounts of the model parameters may be estimated excessively large if the initial period of the response is included. Thus, the recent A/2 time steps were set as the fitting period in step 4.

4.2. Estimation Results with MHE

Figure 16 shows the estimation result with the transient model and MHE, i.e. y(t) in step 1. Vertical axes show the deviation from the mean values of actual data. Compared with Fig. 6, the estimation errors of the key process variables decreased, except for SLC. MHE adjusts the model parameters so that the sum of the estimation errors of the five key process variables is minimized. Although the slight increase of the estimation error of SLC is confirmed, the estimation errors of η_CO, RAR and HMT reduced remarkably. In particular, the estimation error of HMT reduced by 7.8°C. This result shows that MHE allows us to use the 1D transient model persistently for the online control of HMT.

Fig. 16.

Estimation result of key process variables using 1D transient model and MHE. (Online version in color.)

The solid lines in Fig. 17 show the trend of the estimated model parameters. The dashed lines show the fixed model parameters which were employed in the calculation without MHE in section 3.1. The α, β, and γ were adjusted to 0.988, 0.985, and 1.015, respectively, so that the estimation errors are minimized. The estimated coke ratio correction parameter α in Fig. 17(a) indicates the fluctuation of the iron content of iron ore or the carbon content of the coke by more than 1%. Since Fig. 11(a) shows that the change of α by 1% causes the deviation of HMT by 15°C at steady state, the correction of α is crucial for the online estimation of HMT.

Fig. 17.

Result of parameter estimation using 1D transient model and MHE. (Online version in color.)

4.3. Prediction of Future HMT

The discussion so far was concerning the estimation accuracy at present. From the viewpoint of HMT control, the prediction accuracy of the future HMT is more important than the estimation accuracy of the current HMT.

The residence time of iron inside the furnace is roughly 8 hours, therefore the key process variables are predicted up to 10 hours ahead using the following model.   

$$\mathbf{x}(t+m+1)=f\left(\mathbf{x}(t+m),\mathbf{u}(t)\right)$$(24)

  

$$\mathbf{y}(t+m+1)=C\left(\mathbf{x}(t+m+1)\right)$$(25)

To predict the key process variables up to 10 hours (20 time steps) ahead, m takes the values from 0 to 19. It was assumed that the manipulated variables u(t) are kept constant in the future at the current values.

Figure 18 shows the prediction result of the key process variables including HMT. The manipulated variables are shown in the left column, and the key process variables are in the right column. The origin of the horizontal axis is the moment when the prediction is performed. The actual HMT in the future coincides with the predicted HMT even though the actual manipulated variables in the future were not used in the calculation. The heat capacity of the furnace is so large that the future HMT in 10 hours is largely influenced by the accumulation of manipulations in the past.

Fig. 18.

Prediction result of the key process variables using 1D transient model and MHE. (Online version in color.)

Predicting the future HMT with the fixed (constant) manipulated variables is important to achieve stable operation because it is useful to grasp the effect of past manipulations and determine appropriate manipulations in the actual operation. In other words, such future HMT prediction helps operators understand what will happen in the future if they keep the current operating condition (manipulated variables).

Then, the prediction accuracy of eight-hour-ahead HMT was evaluated with long-term data. Considering that the target control accuracy is 10°C and the accuracy of the temperature measurement is 5°C, the target prediction accuracy was set to 15°C. The predicted HMT change in 8 hours (16 time steps) at time step t is obtained by   

$$\mathrm{\Delta }HM{T}_{\text{pre}}(t)=y_5(t+16)-y_5(t).$$(26)

Figure 19 shows the prediction accuracy of the HMT change of 8 hours ahead under various fitting periods; the parameter A is set to 18, 36, and 144. The data period is the same as Fig. 16, which includes 1400 data points. It is confirmed that the longer fitting period leads to better prediction accuracy. A was set to 144 time steps in section 4.2. Comparing the case with MHE of A=144 time steps and the case without MHE, the prediction accuracy improved by 0.4°C in RMSE. The prediction accuracy in the case of A=18 steps is inferior to that without MHE, which indicates that too short fitting period results in the deterioration of the prediction accuracy. It is because the estimation error at present is affected by the past disturbances. Figure 11 suggests that it takes more than 36 hours for the effect of the disturbances to fully appear. Thus, it is necessary to modify the model parameters more than 36 hours before and recalculate the state variables for incorporating the effect of disturbances into the model calculation. These results demonstrate the effectiveness of the prediction method of HMT with the transient model in which the model parameters in the past are successively corrected so that the calculated trajectory of the process variables agrees with the actual data.

Fig. 19.

Prediction accuracy of HMT change based on 1D transient model and MHE under various fitting periods. (Online version in color.)

5. Conclusion

In this research, the prediction of HMT was performed by the 1D transient model and MHE for the thermal control of a blast furnace. To use the transient model for process control persistently, this work aimed at reducing the estimation errors of the model caused by disturbances in the actual furnace. The cause of the estimation errors was analyzed by PCA and the Rist diagram. As a result, it was found that the effective model parameters to reduce the errors were the coke ratio correction parameter, the equilibrium parameter of CO gas reduction, and the heat-loss parameter. These model parameters were successively updated by MHE, and the estimation error of HMT was successfully reduced by 8°C. By using the 1D transient model and MHE, the prediction accuracy of eight-hour-ahead HMT was 11.6°C in RMSE. It was found that adjusting the model parameters in the past is necessary to predict the future HMT accurately.

References

• 1)   Y.  Hashimoto,  Y.  Kitamura,  T.  Ohashi,  Y.  Sawa and  M.  Kano: Control Eng. Pract., 82 (2019), 130.
• 2)   J.  Yagi and  I.  Muchi: Trans. Iron Steel Inst. Jpn., 10 (1970), 392.
• 3)   M.  Hatano,  Y.  Misaka,  Y.  Matoba and  K.  Otsuka: Trans. Iron Steel Inst. Jpn., 22 (1982), 524.
• 4)   D.  Fu,  Y.  Chen,  Y.  Zhao,  J.  D’Alessio,  K. J.  Ferron and  C. Q.  Zhou: Appl. Therm. Eng., 66 (2014), 298.
• 5)   X. F.  Dong,  A. B.  Yu,  S. J.  Chew and  P.  Zulli: Metall. Mater. Trans. B, 41 (2010), 330.
• 6)   Y.  Shen,  B.  Guo,  S.  Chew,  P.  Austin and  A.  Yu: Metall. Mater. Trans. B, 47 (2016), 1052.
• 7)   K.  Yang,  S.  Choi,  J.  Chung and  J.  Yagi: ISIJ Int., 50 (2010), 972.
• 8)   K.  Takatani,  T.  Inada and  Y.  Ujisawa: ISIJ Int., 39 (1999), 15.
• 9)   G.  Ueno,  T.  Higuchi,  T.  Kagimoto and  N.  Hirose: SOLA, 3 (2007), 5.
• 10)   T.  Higuchi: J. Stat. Comput. Simul., 59 (1997), 1.
• 11)   T.  Otsuka: Trans. Soc. Instrum. Control Eng., 35 (1999), 1253.
• 12)   J. B.  Rawlings and  B. R.  Bakshi: Comput. Chem. Eng., 30 (2006), 1529.
• 13)   M.  Kano,  S.  Miyazaki,  K.  Ito and  K.  Hirahara: Earth Planets Space, 65 (2013), 1575.
• 14)   A.  Rist and  N.  Meysson: J. Met., 19 (1967), 50.
• 15)   A.  Spanlang,  W.  Wukovits and  B.  Weiss: Chem. Eng. Trans., 52 (2016), 973.
• 16)   Y.  Ujisawa,  K.  Nakano,  Y.  Matsukura,  K.  Sunahara,  S.  Komatsu and  T.  Yamamoto: ISIJ Int., 45 (2005), 1379.
• 17)   W. E.  Ranz and  W. R.  Marshall: Chem. Eng. Prog., 48 (1952), 141.
• 18)   J.  Mada,  H.  Shinohara and  M.  Tsubakihara: Chem. Eng. Jpn., 27 (1963), 978 (in Japanese).
• 19)   S.  Taguchi,  H.  Kubo,  N.  Tsuchiya,  K.  Ichifuji and  K.  Okabe: Tetsu-to-Hagané, 68 (1982), 2303 (in Japanese).
• 20)   K.  Okabe,  T.  Hamada and  A.  Watanabe: Tetsu-to-Hagané, 55 (1969), 764.
• 21)   H.  Itaya,  T.  Fukutake,  K.  Okabe and  T.  Nagai: Tetsu-to-Hagané, 62 (1976), 472 (in Japanese).
• 22)   L.  Shao and  H.  Saxen: ISIJ Int., 51 (2011), 228.
• 23)   D.  Noguchi,  K.  Ohno,  T.  Maeda,  K.  Nishioka and  M.  Shimizu: ISIJ Int., 53 (2013), 1350.
• 24)   T.  Usui,  H.  Kawabata,  T.  Fujimori,  I.  Fukuda and  Z.  Morita: Tetsu-to-Hagané, 78 (1992), 982 (in Japanese).
• 25)   T.  Maeda and  Y.  Ono: Tetsu-to-Hagané, 75 (1989), 416 (in Japanese).
• 26)   Y.  Iwami,  T.  Yamamoto,  T.  Higuchi,  N.  Oyama,  M.  Sato and  H.  Sodani: Tetsu-to-Hagané, 101 (2015), 177 (in Japanese).
• 27)   R.  Murai,  A.  Murao,  K.  Goto,  M.  Sato and  T.  Ariyama: Tetsu-to-Hagané, 88 (2002), 249 (in Japanese).