Hot Metal Temperature Forecasting at Steel Plant Using Multivariate Adaptive Regression Splines

Abstract

Steelmaking has been experiencing continuous challenges and advances concerning process methods and control models. Integrated steelmaking begins with the hot metal, a crude liquid iron that is produced in the blast furnace (BF). The hot metal is then pre-treated and transferred to the basic lined oxygen furnace (BOF) for refining, experiencing a non-easily predictable temperature drop along the BF–BOF route. Hot metal temperature forecasting at the BOF is critical for the environment, productivity, and cost. An improved multivariate adaptive regression splines (MARS) model is proposed for hot metal temperature forecasting. Selected process variables and past temperature measurements are used as predictors. A moving window approach for the training dataset is chosen to avoid the need for periodic re-tuning of the model. There is no precedent for the application of MARS techniques to BOF steelmaking and a comparable temperature forecasting model of the BF–BOF interface has not been published yet. The model was trained, tested, and validated using a plant process dataset with 12,195 registers, covering one production year. The mean absolute error of predictions is 11.2 °C, which significantly improves those of previous modelling attempts. Moreover, model training and prediction are fast enough for a reliable on-line process control.

1. Introduction

Through the last decades, the steel industry has been facing increasing challenges for improved productivity and quality with optimum cost and minimum impact on the environment. These requirements have been met by advances in production methods and process models [1]. Temperature is one of the most critical variables to be controlled along the process route and a required input for almost any control model in steelmaking operations.

Steel is mostly produced at integrated facilities, starting from iron ore that is reduced in a blast furnace (BF) into a carbon-rich molten iron. This liquid iron, called hot metal, is then transferred to a basic-lined oxygen furnace (BOF) and transformed into liquid steel by blowing oxygen and making use of scrap and other additives. Each batch of liquid steel is often called ‘a heat’. A general overview of BF and BOF processes is provided by Ghosh [2], and a detailed description of the BOF process can be found in Miller [3].

Hot metal is usually transported from the BF to the BOF in a lined vessel, called torpedo wagon or, succinctly, torpedo. During transport and other torpedo operations, the hot metal undergoes a non-easily predictable temperature drop. However, the estimation of the final hot metal temperature is required to calculate the relative quantities of hot metal, scrap, and other raw materials to be loaded into the BOF [4]. Given that these materials account for a significant part of steel cost and carbon footprint, an accurate forecasting of hot metal temperature becomes critical for an optimal BOF operation. Consequently, temperature control from BF to BOF has been receiving much attention. However, existing studies are usually focused on individual stages of the process, such as:

• Prediction and control of hot metal temperature at BF tapping using mathematical models [5,6], measurements [7], or a combination of both [8,9,10].
• Heat transfer modelling in torpedo car and its effects on the refractory lining [11,12], or on the molten iron temperature and metallurgy [13,14].
• Optimum cycling of torpedo fleet, with thermal losses affecting the cost structure [15,16].

General models for predicting the temperature evolution along the complete hot metal route such as those existing for steel ladles are by far much less common [17]. Although some mechanistic models were proposed [18,19,20], this approach seems not reliable in real plant conditions, since many relevant phenomena are difficult to measure or even characterize.

In a previous work, the authors studied the feasibility of infrared thermometry and time series forecasting for accurate prediction of the hot metal temperature at the steel plant [21]. A combination of both methods proved to be reliable and long-term stable, as well as self-adaptive to changing production scenarios. However, this research revealed the necessity to improve the regressive part of the model. The authors also explored the application of artificial neural network (ANN) to this problem but, although resulting accuracy was satisfactory, long-term stability of the model could not be guaranteed [22].

The present work aims to obtain an advanced hot metal temperature forecast at the BOF. A model based on multivariate regression can take full advantage of the information provided by the exogenous predictors already available in the process databases.

Multivariate adaptive regression splines (MARS) stands out among other multivariate regression techniques. Since its introduction by Friedman [23,24] in the 90s it has been successfully applied to a variety of fields including life, finance, industry, business, energy, and environment. However, for the moment, the application of MARS to steelmaking processes remains limited to steel solidification by continuous casting and solid state transformations by hot and cold rolling processes [25,26,27]. Despite the predictive possibilities of MARS technique, no application to BOF steelmaking has been published yet.

Temperature forecasting with MARS has been mainly focused on natural systems [28,29]. Only very recently, Krzemień [30] has applied this technique to predict the temperature evolution in underground coal gasification. For this problem, the syngas chemical composition, flowrate, and temperature were used to predict the temperature of syngas 1 h ahead. The obtained absolute error was less than 15 °C in 95% of the cases for a temperature range of about 200 °C.

Some properties of MARS make it an interesting choice for this research [31]. Firstly, it is more flexible than linear regression, allowing to model nonlinearities and interactions between variables. Additionally, input data are automatically partitioned, limiting the effect of the undetected outliers that can be expected from a large industrial dataset. Moreover, it automatically includes important variables in the model and excludes unimportant ones; this helps model stability even in changing process conditions. MARS is also suitable for handling fairly large training datasets at low computing cost and then provides very fast predictions. This characteristic is very valuable for a process control model that is going to be continuously trained and repeatedly executed. Finally, MARS provides interpretable models in which the effect of each predictor can be clearly interpreted. This is true not only for additive models but also when interactions between variables are present. This strength should not be underrated, since model interpretability strongly affects process improvement and industrial users’ satisfaction.

One drawback of the method is the non-differentiability of the model at discrete points. This limitation can be overcome by locally smoothing the resulting piecewise model in a post-processing phase or by using local higher order splines [23,32]. These variants can be interesting, for example, for global sensitivity purposes [33]. However, usually local non-differentiability does not affect the prediction performance of the model.

The feasibility of MARS technique for hot metal temperature forecasting is investigated in the present work. The effect of model hyperparameters on predictions accuracy is assessed. Moreover, a novel approach for continuous model training is adopted to ensure long-term adaptive operation of the model [34,35]. The resulting improved temperature forecast, is used as an input for the BOF charge model and results in environmental, productivity, and cost improvements [36,37,38].

2. Materials and Methods

2.1. Explanatory Variables

The hot metal process from the BF to the BOF is summarized in Figure 1 and

Table 1. Model variables.

Description	Symbol	Min	Max	Unit
Initial temperature	X1	1400	1540	°C
Total holding time	X2	2	20	h
Pre-treatment duration	X3	0	40	min
Empty torpedo duration	X4	1	16	h
Empty ladle duration	X5	0	8	h
Final temperature	Y	1200	1420	°C

. Five explanatory variables, X_i, were selected for predicting the final hot metal temperature at BOF, Y. The initial temperature of the hot metal, X₁, is measured with disposable thermocouples in the iron runner of the blast furnace during tapping [39,40]. Three measurements are usually taken by cast: Firstly, just after drilling the tap hole; secondly, when the slag arises; and finally, approximately at the end of the cast [5]. Hot metal temperature for each torpedo is calculated by time interpolation between consecutive thermocouple measurements.

The total holding time, X₂, that is to say, the time between temperature measurements in the BF and in the BOF, was taken as the effective transport time. It comprises torpedo car operations, pouring of hot metal into the ladle, and ladle transport. Since the hot metal pouring may extend over a significant period of time, there is not a clear time limit between torpedo and ladle. Considering that the heat losses in torpedoes and ladles were found to be similar [21], both holding times can be grouped together without excessive simplification.

The pre-treatment duration, X₃, accounts for the thermal losses during this phase. It is assumed that the mass flowrates of desulfurizing agents and inert gases are constant and therefore, the phase duration is the main aspect to be considered. Possible differences between different desulfurizing agents are neglected assuming that hot metal stirring causes the main effect on temperature.

The time during which the torpedoes and ladle stay empty after the previous cycle, X₄ and X₅, respectively, where chosen as a convenient way of describing their initial thermal condition. Other aspects, such as the actual lining thickness or the amounts of slag and iron solidified inside the torpedo, cannot be accurately measured and are not considered in this model. Moreover, lining pre-heating is not considered because burners efficiency was not fully described in the available data. Considering that the number of pre-heated torpedoes and ladles was less than 5%, these cases were left outside the scope of the model.

Finally, the actual hot metal temperature, Y, is measured in the hot metal ladle with a disposable thermocouple just before BOF loading [4,40]. The past measured values, Y, are used for model training and to assess the performance of model predictions, $\widehat{Y}$. More details for variable selection and their characteristics are given in Díaz et al. [21].

2.2. Process Dataset

The selected dataset comprises 12,195 registers, which roughly corresponds to one production year. It is composed of the actual measured hot metal temperature together with the corresponding values of the five above mentioned predictors. The first 10,000 registers were used for model training and hyperparameter testing while the last 2195 entries were kept for the final validation of the model.

Widely different variable values and scales can cause instabilities that could adversely affect the performance of the model [23]. Therefore, variables were pre-scaled to the [0,1] interval using the minimum–maximum normalization for adjusting the different dimensions and ranges to an easily interpretable standard scale:

$$x=\frac{X-\min \left(X\right)}{\max \left(X\right)-\min \left(X\right)}$$

(1)

where X represents the original variable, x its scaled version and min(X) and max(X) are the observed limits of the historic register in the whole process database, according to Table 1.

The normalized dataset is graphically displayed in Figure 2, where the minimum–maximum interval and the average value for groups of 30 heats are shown for clarity. This grouping size corresponds approximately to one production day.

It can be observed that initial and final temperatures x₁ and y, respectively, exhibit similar distributions and temporal evolution. Rather than random changes, local tendencies can also be recognized; therefore, it seems appropriate to retain and exploit the time information contained in the data.

The similarities between y and x₁ curves and histograms reflect the expected correlation between both variables. However, several features of the y curve do not match well with corresponding features of x₁, as for example the drop in y that can be observed in the vicinity of t = 10,000. This indicates that other variables are also critical.

Total holding and empty vessels durations (x₂, x₄, and x₅) have lower daily means, showing that normal production times are usually short with occasional longer times. Holding time x₂ and empty torpedo time x₄ are very similar, being dominated by torpedo logistics. Since empty torpedo movements are less critical than full ones, x₄ has more dispersion than x₂.

Empty ladle time, x₄, exhibits a multimodal distribution showing the different process situations to which the steel plant reacts with a different number of hot metal ladles in service.

Daily average of pretreatment time x₃, is much more centered as can be expected from prescribed desulfurization requirements. The bimodal histogram reveals that cases without pre-treatment (x₃ = 0) are frequent.

2.3. Multivariate Adaptive Regression Splines with Moving Training Window

The application of multivariate adaptive regression splines (MARS) to hot metal temperature forecasting at BOF shop is investigated and its performance is assessed. The following model is proposed:

$${\widehat{y}}_t={\beta }_0^{\left(w,t\right)}+{\displaystyle \sum }_{m=1}^M{\beta }_m^{\left(w,t\right)}{B}_m^{\left(w,t\right)}\left({\overrightarrow{x}}_t\right)$$

(2)

where the temperature prediction for heat number t, ${\widehat{y}}_t$, is a function of the N selected explanatory variables for that heat t: ${\overrightarrow{x}}_t=\left(x_{t1},x_{t2},\dots x_{tN}\right)$. The model is constructed as a weighed sum of M basis functions, $B_m^{\left(w,t\right)}\left({\overrightarrow{x}}_t\right)$ using constant weighing coefficients, ${\beta }_m^{\left(w,t\right)}$. The w and t indexes in coefficients and functions indicate that, for every new prediction ${\widehat{y}}_t$, both, basis functions and coefficients, are optimized by regression analysis using the w previous observations (${\overrightarrow{x}}_{t-1}$, ${\overrightarrow{x}}_{t-2}$, … ${\overrightarrow{x}}_{t-w}$). Hence, w represents the width of the moving training window.

The determination of coefficients in Equation (2) by regression analysis is straightforward, although there are multiple possible construction methods for the basis functions [41]. A non-parametric adaptive technique based on domain segmentation together with a simple sub-set of basis functions is adopted here. It was introduced in 1982 by Smith [42] for the univariate problem and extended for a multivariate setting in 1991 by Friedman [23,24] with the name of multivariate adaptive regression splines (MARS). A practical implementation of this technique, the ARESLAB toolbox for MATLAB, has been made publicly available by Jekabsons [32].

The model is built in two phases: Forward selection and backward deletion. In the forward pass, a pair of new basis functions is added at each step, resulting in an increasingly complex model that will not generalize well to new data. In the backward pass this overfitted model is pruned by removing the least effective term at each step.

2.3.1. Forward Phase

The forward phase starts with a model which consists of just the intercept term, that is, $B_0^{\left(w,t\right)}\left({\overrightarrow{x}}_t\right)=1$. At each additional step a new pair of basis functions is added to the model. Each new pair of basis functions consists of a basis function already existing in the model multiplied by a pair of functions of the form ${\left(x_i-c\right)}_{+}$ and ${\left(c-x_i\right)}_{+}$, where the subscript “+” means the positive part, that is:

$$\begin{array}{c}{\left(x_i-c\right)}_{+}=\{\begin{array}{c}{x}_i-c,\mathrm{if}{x}_i>c,\\ 0,\mathrm{otherwise},\end{array}\\ {\left(c-x_i\right)}_{+}=\{\begin{array}{c}c-x_i,\mathrm{if}{x}_i<c,\\ 0,\mathrm{otherwise}\end{array}\end{array}$$

(3)

where x_i is one of the input variables and c is one of the values of that input in the dataset (i.e., an observation of x_i), which is usually referred to as the knot of the pair of basis functions. In other words, a new pair of functions from the collection

$${\left\{{\left(x_i-c\right)}_{+},{\left(c-x_i\right)}_{+}\right\}}_{\begin{array}{c}c\in \left\{x_{t-1,i},x_{t-2,i},\dots ,x_{t-w,i}\right\}\\ i=1,2,\dots ,N\end{array}}$$

(4)

is multiplied by a basis function already in the model to obtain the new pair of basis functions. The algorithm compares all the combinations of existing terms, variables, and values for each variable in the data set.

One additional restriction is applied to the formation of basis functions: Each variable can appear at most once in a product (i.e., maximum self-interaction order, S = 1). This prevents the formation of higher-order powers of a single variable, which increase or decrease too sharply near the boundaries of the domain. Such powers can be approximated in a more stable way with piecewise linear functions. It is also possible to set an upper limit on the order of interaction between different variables, I. For example, a limit of two will allow pairwise products of linear functions, but not three or higher order products. An upper limit of one results in an additive model with no interactions between variables.

The coefficients β_m are estimated by minimizing the residual sum-of-squares, that is, by standard linear regression finding the pair of basis functions that gives the maximum reduction in error. This pair is added to the model and this process of adding terms continues until the number of terms in the model reach a prescribed limit.

2.3.2. Backward Phase

At the end of the forward phase, a large model of the form in Equation (2) is obtained containing M = M_F basis functions. This model typically overfits the data and will not generalize well to new data. Therefore, a backward deletion pass is applied by iteratively deleting the basis function whose removal causes the smallest increase in residual squared error until the model has only the intercept term. An estimated best model, ${\widehat{f}}_{\lambda }$ for each model size λ, where λ is the number of terms in the model, is obtained at each step. Cross-validation could be used to estimate the optimal model size, but for computational savings the generalized cross-validation (GCV) is preferably used and works well in practice [23,43,44]:

$$GCV=\frac{\frac{1}{w}{{\displaystyle \sum }}_{i=t-1}^{t-w}{\left(y_i-{\widehat{f}}_{\lambda }\left({\overrightarrow{x}}_i\right)\right)}^2}{{\left(1-\frac{C\left(\lambda \right)}{w}\right)}^2}$$

(5)

The GCV criterion is the average-squared residual of the fit to the data, increased by a penalty factor that accounts for the increased variance associated with increasing model complexity C(λ) [32]:

$$C\left(\lambda \right)=\lambda +d\left(\lambda -1\right)/2$$

(6)

where λ is the number of basis functions in the model including the intercept term, and d is the penalty. Since (λ − 1)/2 is the number of knots, C(λ) penalizes the model not only for its number of basis functions but also for its number of knots. Some mathematical and simulation results suggest that a price of three parameters for selecting a knot in a piecewise linear regression is a convenient choice (d = 3) for the general case, and two (d = 2) for an additive model [44]. However, d can always be selected as a model hyperparameter to be optimized.

The resulting GCV is an estimator of the mean squared error of the model in prediction mode, therefore the minimum value of GCV reflects a compromise between fit and model complexity, resulting in a better generalization of the model.

The forward pass adds basis functions in pairs at each step, while the backward pass deletes one basis function, the less relevant, at a time. Consequently, terms are often not seen in pairs in the final model.

In a previous modelling work, time series were applied to solve this problem. It was shown that the introduction of exogenous predictors into an auto-regressive integral moving average (ARIMA) model improved its performance [21]. Conversely, the potential benefit of introducing the L last observations of the predicted variable as input variables will be analyzed in the present work. Accordingly, for L > 0, ${\overrightarrow{x}}_t$ in Equation (2) is expanded as ${\overrightarrow{x}}_t=\left(x_{t1},x_{t2},x_{t3},x_{t4},x_{t5},y_{t-1},y_{t-2},\dots y_{t-L}\right)$.

2.3.3. Model Hyperparameters

As previously discussed, six possible model hyperparameters have been considered: S, I, M_F, d, L, and w. The first four are typical settings of the standard MARS technique, while the last two (L and w) are the original of this study. In particular, w is a consequence of the novel continuous training approach adopted here which is required to ensure long-term performance of the model. The basic configuration of these six hyperparameters, together with the ranges to be tested, and the finally adopted configuration can be found in

Table 2. Model hyperparameters with their initial values, explored ranges, and final configuration.

Description	Symbol	Base	Min	Max	Final
Maximum self-interaction order	S	1	1	1	1
Maximum interaction order	I	1	1	5	1
Maximum functions in the forward phase	MF	21	3	49	21
Penalty for model complexity	d	2	0	12	2
Width of moving training window	w	1000	10	2000	2000
$y_{t-1},y_{t-2}\dots y_{t-L}$ as additional predictors	L	0	0	6	4

.

An initial hyperparameter set of {S = 1, I = 1, M_F = 21, d = 2, L = 0, w = 1000} was taken as the base case. Self-interactions of variables where not allowed (S = 1) in order to avoid singularities near the boundaries of the domain [23,32]; this criterion was maintained along the study. The remaining parameters where varied to assess their effect on model performance. Interactions between variables where not allowed for the base case (I = 1) in order to start with a simple additive model. The maximum number of base functions in the forward phase, M_F, was taken initially as maximum (21, 2N + 1), where N is the number of input variables, as suggested by Milborrow [45]. The GCV penalty per knot was chosen as d = 2, as suggested by Hastie [44] for additive models. Finally, the initial choices for w and L, were based on previous experiences of the authors [21,22].

3. Results

3.1. Computation

The proposed moving MARS model was implemented on MATLAB using the ARESLab toolbox [32]. The different cases were solved on a regular 64-bit laptop with a 2.6 GHz processor and 4 GB RAM. Computation cost was mainly due to the training of the model and it increased with the number of observations used to train de model, w, and the number of input variables, N = 5 + L. However, the resulting training time of a single case was always below 4 s. Considering that the interval between heats is typically 15–30 min, the computation time is far to be a limiting factor for the real time application of the model.

The model was developed using a test dataset formed by 10,000 registers. The initial 2000 registers were reserved for initializing the training window with a maximum of 2000 registers. Consequently, each study case involved testing 8000 different models, from t = 2000 to t = 10,000. The overall study required the evaluation of more than 2.4 × 10⁶ different models. Results were compared in terms of the mean absolute error (MAE). This is a convenient choice since the error of hot metal temperature forecast is proportional to the excess or defect in hot metal consumption [38]. Therefore, MAE reduction has a straight translation to economic savings and environmental improvements [36,37,38].

3.2. Study of Hyperparameters

The influence of different parameter configurations on model performance is shown in Figure 3. It can be seen that w has a notable effect on model accuracy: The greater the w, the lower the MAE provided that the other hyperparameters are equal. This effect is more important for lower w; while for w > 1000 the reduction in MAE is lower than 0.0003. Additionally, the effect of any other hyperparameter tends to be lower as w increases. As a conclusion, it is interesting to set w ≥ 1000 not only to improve model accuracy but also to reduce its sensibility with regard to the other hyperparameters.

As shown in Figure 3a, if any order of interaction between input variables (I > 1) is allowed, the MAE increases. This effect is dramatic for w < 500 and much less pronounced for w > 500. This result suggests that the actual process performs in this way and that increasing I, only provides unnecessary degrees of freedom to the model, with the side effect of overfitting.

Figure 3b shows that the initial choice of M_F = 21 for the maximum number of functions in the forward phase is close to its optimum value when w ≥ 1000. In fact, the best MAE for w = 1000 and 2000 is obtained when M_F = 11 but the improvement in MAE is only around 0.0003 with regard to M_F = 21. This finding suggests that, for this particular problem, M_F = 2N + 1 is a good choice, but max(21, 2N + 1) = 21 can also be a reasonable selection for a large dataset [45]. The small effect of an additional increase of M_F is logical, since the number of basis functions in excess at the end of the forward phase will be pruned in the backward phase of the method. However, M_F still has some influence because GCV is a good proxy—but only a proxy—of the forecasting capabilities of the model. Consequently, increasing M_F above the optimum gives slightly overfitted models.

The lowest MAE is obtained when the penalty for model complexity, d, is between 4 and 6 as can be seen in Figure 3c. For w < 500 a clear minimum is reached whereas for a larger w the MAE remains around the minimum for any d ≥ 6. Again, the effect in MAE is below 0.0004 for w ≥ 1000, indicating the robustness of generalization of the model when it is trained with enough data.

The introduction of the L previous observations of the predicted variable as additional predictors causes very different effects depending on the value of w. Figure 3d shows that for w = 100 it has an adverse effect on MAE, regardless the value of L. For w between 200 and 500 the improvement is modest (0.0003 to 0.0008), reaching a minimum for 2 < L < 4. Finally, for w ≥ 1000 an evident reduction of 0.0012 in MAE is attained for 3 < L < 5. This finding is in line with the results previously obtained for a simpler model based on moving average smoothing (MAS) [21]. In that case, in which the only predictors were the L previous observations of the predicted variable, the optimum value was 4 < L < 6. It can be inferred from both models that up to four or five previous observations give some information of the temperature tendency of the process. Older observations seem to be affected by changes in the process and they are not correlated with present conditions.

Finally, Figure 4 shows the effect of the most significant hyperparameters (M_F, d, and w) when L = 4. It can be seen that the incorporation of four lagged observations not only improves MAE but also makes it less dependent on the other hyperparameters. For w = 2000, any M_F > 15 gives the lowest MAE of 0.0506. Similarly, for w = 2000, any d ≥ 2 gives also the minimum MAE. It can be concluded that any model configuration within {S = 1, I = 1, M_F > 15, d ≥ 2, L = 4, w = 2000} gives the best results.

3.3. Model Validation

Based on the analysis described above, the final configuration of the model hyperparameters, indicated in Table 2, {S = 1, I = 1, M_F = 21, d = 2, L = 4, w = 2000} was applied to the validation dataset for assessing the actual forecasting performance of the method. The model errors from t = 10,001 to t = 12,195 are shown in Figure 5.

The absolute errors are always lower than 0.25 and the daily mean error is always lower than 0.0500. The global MAE for the validation dataset is 0.05076 which is consistent with the value of 0.0506 obtained during hyperparameter optimization, confirming that GCV is a convenient estimator for the forecasting capability of the model [23,43,44].

4. Discussion

As indicated before, the application of the final MARS model {S = 1, I = 1, M_F = 21, d = 2, L = 4, w = 2000} to the validation dataset provided 2195 evaluations of the method from t = 10,001 to t = 12,195. The resulting basis functions and coefficients in Equation (2) at t = 11,755 are shown in

Table 3. Model equations for MARS {S = 1, I = 1, M_F = 21, d = 2, L = 4, w = 2000} at t = 11,755. The mean squared error (MSE) and the generalized cross validation (GCV) of the model predictions for the training data are 0.00408 and 0.00418, respectively. The columns MSE and GCV indicate the new values obtained when the basis function is removed from the model.

Basis Function		Coefficient		MSE	GCV
$B_0$	1	${\beta }_0$	0.4411	-	-
$B_1$	(x2 − 0.2672)+	${\beta }_1$	−0.5320	0.00541	0.00553
$B_2$	(0.6998 − x4)+	${\beta }_2$	0.3201	0.00526	0.00537
$B_3$	(x1 − 0.4643)+	${\beta }_3$	0.3777	0.00510	0.00521
$B_4$	(0.2672 − x2)+	${\beta }_4$	0.4428	0.00455	0.00465
$B_5$	(0.4643 − x1)+	${\beta }_5$	−0.5489	0.00430	0.00439
$B_6$	(0.6409 − x6)+	${\beta }_6$	−0.1124	0.00414	0.00423
$B_7$	(x5 − 0.0312)+	${\beta }_7$	−0.0683	0.00412	0.00421
$B_8$	(x6 − 0.6409)+	${\beta }_8$	0.1408	0.00412	0.00421
$B_9$	(x7 − 0.5818)+	${\beta }_9$	0.0862	0.00411	0.00420
$B_{10}$	(x3 − 0.4571)+	${\beta }_{10}$	−0.0526	0.00410	0.00419
$B_{11}$	(x8 − 0.6818)+	${\beta }_{11}$	0.0928	0.00409	0.00418

. The model at this point is taken as an example to discuss the features of the model. In this case, the model comprises 12 basis functions including the intercept term. Considering that 21 functions where allowed in the forward phase (M_F = 21), it is inferred that nine functions with the lowest contribution to GCV where pruned in the backward phase.

As can be deducted from Table 3, the most important basis functions are those containing x₁, x₂, and x₄, that is, the initial hot metal temperature, the total holding time, and the empty torpedo time, respectively. Empty ladle time and pre-treatment time have less importance; in fact, the actual temperature of the previous heat x₆ = y_t−1, is a better predictor. Moreover, the MARS method automatically excludes the less important variables. For this particular heat (t = 11,755), the actual temperature four heats before, x₉ = y_t−4, is not included in the model. This indicates that its contribution to model performance is negligible or even adverse. This is not the general case for every heat; for example, at t = 10,656 a basis function is also included for x₉. In general, it is positive to include x₉ in the model to improve its performance, as indicated in Figure 3c for curve w = 2000 at L = 4. It can be seen that some basis functions appear in pairs, as is the case for x₁, x₂, and x₆. Other basis functions are individual, indicating that the corresponding symmetric functions were removed in the backward phase. A basis function without its symmetric variant indicates that either the involved variable is relevant only above a particular value (this is the case of x₃, x₅, x₇, and x₈), or in the lower part of its range (as for x₄). As can be seen, the adaptive knot location of MARS method succeeds in capturing the nonlinearities of the data using segmented linear regression.

The model equations obtained for other heats have very similar shape with slight changes. The most important features are always present in the model. Other characteristics arise or go away depending on the history of the previous w heats.

The actual shape of the basis functions can be better understood combining Table 3 and the graphical representation of the model in Figure 6. The line plots located in the diagonal of the mosaic show the model output, that is, the predicted hot metal temperature, $\widehat{y}$, against each individual input variable, x_i. The rest of the input variables are kept at their mid-points (x_k = 0.5 for k ≠ i).

The line plot of ${\widehat{y}}_t$ versus x₁ indicates that the effect of the initial temperature tends to damp as it increases. This is a coherent result, since higher thermal losses are foreseen from higher initial temperatures in all the phases of the process. A similar reasoning can be applied to the empty torpedo time, x₄, considering that the thermal losses are expected to be higher for shorter times, since the lining temperature is higher. The effect of holding time, x₂, and empty ladle time, x₅, was found to be almost linear within the considered ranges, and with the anticipated slope, as can be seen in the related plots.

As predictable, the higher the temperature of the previous heats, the higher the expected temperature for the actual heat. Moreover, this effect seems to be only relevant, or more relevant, for higher temperatures. Although it can be hardly perceived in Figure 3 (but can be confirmed in Table 3), the slope of the line plot for ${\widehat{y}}_t$ versus x₆ is smaller for x₆ < 0.64. This behavior is more evident for x₇ and x₈ which are relevant only above 0.58 and 0.68, respectively.

The response surface for all the pair-wise combinations of the input variables are represented in the contour plots located above the diagonal in Figure 6. Again, the input variables outside the considered pair are kept at the mid-point of their ranges (x_tk = 0.5 for k ≠ i and k ≠ j). It can be seen that the adaptive knot location succeeds in representing the non-linear features of the multivariate process dataset.

Finally, the response surface around the actual value of the input at heat number t = 11,755, ${\overrightarrow{x}}_t$ is shown in the elements below the diagonal in Figure 6. Here, the color scale gives the predicted temperature (${\widehat{y}}_t$ = 0.056) and the circle marks the input vector.

The proposed model is dynamically trained for every new prediction to be made. This dynamic behavior is better appreciated in Figure S1, an animated version of Figure 6 showing model representation from t = 10,004 to t = 12,195. Figure S1 is included in the Supplementary Materials.

The good predictive performance of the proposed model can be better appreciated by comparison with other previously developed models as shown in Figure 7. A common reference model is the moving average smoothing (MAS) which forecasts temperature as the average of the w previous observations. The best MAE for this method is 0.0721 for w = 5; higher w values provide no benefit. However, for w > 100, time-series auto-regressive integral moving average with eXogenous predictors (ARIMAX) models can be applied providing a MAE below 0.0690 for w > 1000. On the other hand, for w > 20, any MARS configuration provides a major improvement over simple reference models but also over the more advanced ARIMAX models. The initial base configuration for MARS {S = 1, I = 1, M_F = 21, d = 2, L = 0, w = 1000} gives a MAE of 0.0518 giving a 25% error reduction with reference to ARIMAX. The best performing MARS {S = 1, I = 1, M_F = 21, d = 2, L = 4, w = 2000} provides 0.0506 of MAE, representing a 27% and a 30% of error reduction with regard to ARIMAX and MAS, respectively. This model configuration is the new benchmark for this problem.

As demonstrated in the previous section, the introduction of lagged terms as additional predictors further reduces the MAE of the model. However, it is worth noting that the higher the number of input variables, the higher the w required to obtain an optimum performance of the model, as can be seen in Figure 3d and Figure 7. The base configuration of MARS {S = 1, I = 1, M_F = 21, d = 2, L = 0} requires at least ten previous observations of the five input variables, as illustrated in Figure 7. Configurations with additional input variables would require more than one hundred previous observations of all the inputs. This limitation poses a potential problem when some registers are missing in the process database. However, a judicious implementation of the model should relieve this problem, as indicated by the shaded region in Figure 7. This region delimits the lowest MAE that can be achieved by applying the best configuration for the available data-window at execution time.

5. Conclusions

An improved version of the standard MARS technique has been proposed for hot metal temperature forecasting at the steel plant. It employs five robust process variables and the four previous observations of the predicted variable as model inputs. A moving window approach is used for model training in order to avoid the need for periodic manual tuning.

There is no precedent for the application of MARS techniques to BOF steelmaking and a comparable temperature forecasting model of the BF–BOF interface has not been published before.

The effect of model hyperparameters in this problem has been explored and understood. The most influencing factor is the training window width, being the greater the best. However, widths above 1000 provided only marginal improvement in MAE. Similarly, the introduction of the four previous temperatures as additional predictions gave a small additional reduction in MAE. However, any little refinement in hot metal temperature forecasting is considered relevant, provided the mega-scale effect of the steelmaking industry [21,36,37,38].

The optimum hyperparameter configuration gave a MAE of 0.05076 (11.2 °C) on the validation phase. This result outperforms those of previous modelling attempts, with a MAE ranging from 0.083 to 0.069 (18.3 to 15.2 °C) [21].

The required computing time was lower than 4 s in total for training plus forecasting, the former thanks to GCV approximation, and the latter due to the simple mathematical formulation of the model. This computing time ensures real time operation, allowing reiterative execution for forecast updating.

The resulting model succeeds in capturing the non-linearities of the process but it is still easy to interpret. The model adapts to changes along time, since the relevant input variables, the number of basis functions and, of course, the fitting parameters are re-calculated for every new prediction. In a future work the adaptability performance of the model will be tested and optimized by introducing synthetic changes in the process.