Data Reconciliation-Based Hierarchical Fusion of Machine Learning Models

Abstract

In the context of hierarchical system modeling, ensuring constraints between different hierarchy levels are met, so, for instance, ensuring the aggregation constraints are satisfied, is essential. However, modelling and forecasting each element of the hierarchy independently introduce errors. To mitigate this balance error, it is recommended to employ an optimal data reconciliation technique with an emphasis on measurement and modeling errors. In this study, three different machine learning methods for development were investigated. The first method involves no data reconciliation, relying solely on machine learning models built independently at each hierarchical level. The second approach incorporates measurement errors by adjusting the measured data to satisfy each constraint, and the machine learning model is developed based on this dataset. The third method is based on directly fine-tuning the machine learning predictions based on the prediction errors of each model. The three methods were compared using three case studies with different complexities, namely mineral composition estimation with 9 elements, forecasting of retail sales with 14 elements, and waste deposition forecasting with more than 3000 elements. From the results of this study, the conclusion can be drawn that the third method performs the best, and reliable machine learning models can be developed.

1. Introduction

Digitalization offers significant opportunities through various machine learning (ML) algorithms that can replace conventional quality assurance methods. When developing ML models for processes, only certain parameters can typically be estimated with adequate precision. These ML model estimations carry varying degrees of error. Our modeling approach integrates engineering insights and additional data to ensure that the estimates meet specific constraints or conditions. The balance equations between hierarchical levels must be satisfied in the context of modeling hierarchical systems, and these constraints should be considered when training models.

Frequently, the systems to be modeled are naturally organized in hierarchical structures, mainly in forecasting problems; for instance, the demand for a product can be recorded on different hierarchy levels, like on a store, regional, or country level [1]. The measurements and observed values at each level will add up to the higher levels, which is called “coherence” [2]. In practical solutions, model development and prediction occur independently at each hierarchy level, resulting in incoherence in the results, so the predictions do not aggregate well.

This incoherence or balance error between the hierarchy levels was handled myopically with “bottom-up”, “middle-out”, and “top-down” methods. With these methods, it is not necessary to develop a model for each hierarchy level, but only one model is developed, and then we aggregate and/or disaggregate the predictions to other hierarchy levels [3]. For instance, the “bottom-up” technique works by developing models at the most granular level of the hierarchy and then summing lower-level predictions [4].

The fact that these myopic methods do not consider some useful information on other hierarchy levels is partly solved by combination approaches, which use statistically weighted information from all hierarchical levels [5]. However, none of the mentioned approaches result in optimal reconciliation among the predictions [6].

Optimal reconciliation is based on generating independent models for all elements of the hierarchy, where the prediction results are incoherent. However, the next step is to perform optimal reconciliation to adjust the independent predictions, which will lead to the prediction being consistent with the hierarchical structure [7]. In this case, reconciliation is formulated as a regression model which projects the elementary forecasts and predictions onto a subspace where the predictions adhere to the aggregation constraints [8]. Weighted squared-error optimal reconciliation was proposed to perform optimal reconciliation forecasts, where the base predictions are adjusted minimally due to the least-squares function [9,10]. Generalized least squares and weighted least squares are also a solution that can be used to obtain the optimal reconciled forecasts and predictions, but it is really difficult to estimate the covariance matrices used in the solution. The easiest way to perform an ordinary least-squares estimate is to reconcile the base predictions [11]. These hierarchical forecasting methods worked well in different cases, like forecasting the electrical loads of a building [12], in supply chain forecasting [13] or in tourism forecasting [14].

The presented optimal reconciliation approaches only deal with the constraint due to the hierarchical system, but other system- and modeler-defined specifications and constraints can be imagined, too; for instance, the predictions must add up to a user-specified constraint, like the percentages must add up to 100 percent [15]. Data reconciliation (DR) is a general method used to implement and define both hierarchical constraints and any other necessary constraint [16]. DR enables the correction of measurements, reducing the associated uncertainty while satisfying different constraints [17].

The advantages of DR are also presented in various case studies, for instance, in analysis of measurement outcomes within the field of analytical chemistry, where the application of DR enhances predictions through real-time mass and element balancing. This strategy has been shown to decrease the standard error in predictions, eliminating the need for additional offline analyses [18]. The efficacy of the DR method in processing analytical chemistry data is clearly beneficial. In practical problems, maintaining variable consistency is crucial, and increased process automation requires rigorous monitoring. Errors are unavoidable in the measurement, processing, and transmission of signals. These errors can degrade the performance of the monitoring and control system and sometimes lead to process failures. Thus, minimizing these error effects is critical [19]. DR improves process measurement estimates by mitigating the impact of random errors [20]. The outcomes are independent and consistent, indicating that the value model at a higher hierarchical level produces results that differ from the aggregation of the lower-level models’ results. In this case study, the ML models are structured hierarchically and subject to specific constraints.

In this study, a hierarchical modeling approach with optimal DR is presented that improves the performance of ML models to facilitate optimization with digital tools. With this approach, ML models can be developed that satisfy hierarchical and other user-specified constraints. Furthermore, with the expansion of software sensors in Industry 4.0 [21], the proposed method can be implemented in edge computing devices. Various case studies with different complexities were used to evaluate the proposed methodology. The main novelties of our work are summarized as follows:

• An approach to managing hierarchical constraints was developed to enhance the performance of ML models by accounting for the prediction errors in each model (see Section 2).
• In this study, a connection was established between the summation matrix utilized in hierarchical time series (HTS) forecasting and the incidence matrix employed in traditional DR methods (see Section 2).
• The developed methods exhibit strong performance in a range of case studies with different complexities. Our tests included a three-level scenario with 9 elements in rock composition estimation from spectral signals, a three-level scenario with 14 elements in a distribution model (retail sales M5 competition), and a four-level waste deposition scenario involving more than 3000 elements (Hungarian counties, districts, and cities) (see Section 3).

2. Integrated Correction of Machine Learning Predictions Using Data Reconciliation Techniques

This section outlines the methods used in our approach to combine ML and DR techniques to improve the accuracy of ML predictions. By merely predicting each series separately, the hierarchical or grouping structure is ignored, resulting in forecasts that are not ‘coherent’, so they do not aggregate correctly and do not satisfy the hierarchical constraints.

2.1. Formulating the Integration of Machine Learning and Data Reconciliation

The goal of this study is to correct ML predictions using optimal DR considering the modeling errors in hierarchical problems. Due to the hierarchical structure used, our method deals with multivariate modeling by considering multiple target variables ($\mathbf{y}$). The basic assumption is that all modeled variables are subject to error, so the model estimates for all variables need to or should be corrected. Consider hierarchical data, the target variables (${\mathbf{y}}_t$) in the t-th sample instance can be written, in general, as follows:

$${\mathbf{y}}_t=\mathbf{f}({\mathbf{x}}_t,\mathit{\theta })+{\mathit{ϵ}}_t$$

(1)

where ${\mathbf{y}}_t$ is the vector of the target variables, ${\mathbf{x}}_t$ represents the matrix of independent variables, $\mathbf{f}$ denotes the set of ML models, $\mathit{\theta }$ represents the matrix of model parameters, and ${\mathit{ϵ}}_t$ is the prediction error vector. The model prediction ($\widehat{\mathbf{y}}$) of the target variables can be written, in general, as follows:

$${\widehat{\mathbf{y}}}_t=\mathbf{f}({\mathbf{x}}_t,\mathit{\theta })$$

(2)

Given that the measurements must adhere to the specified constraint, it is imperative that the model predictions comply as well. To achieve this, DR is used in this study. This process adjusts the predicted variables minimally to ensure compliance with a set of model constraints. It aims to reduce the discrepancy between predicted and reconciled values while considering the variance of these variables and guarantees that the reconciled parameters meet certain equality and inequality constraints [22]. Typically, the objective function for minimization is denoted as Equation (3). The general non-linear DR problem is outlined as follows:

$$\underset{{\tilde{\mathbf{y}}}_t}{\min }\left({\widehat{\mathbf{y}}}_t-{\tilde{\mathbf{y}}}_t\right){\mathbf{V}}^{-1}\left({\widehat{\mathbf{y}}}_t-{\tilde{\mathbf{y}}}_t\right)$$

(3)

subject to

$$\mathbf{h}\left({\tilde{\mathbf{y}}}_t\right)=\mathbf{0}$$

(4)

$$\mathbf{g}\left({\tilde{\mathbf{y}}}_t\right)\le \mathbf{0}$$

(5)

where ${\mathbf{V}}^{-1}$ is the inverse of the covariance matrix of errors, t is the sample instance, ${\widehat{\mathbf{y}}}_t$ is a vector of model predictions, ${\tilde{\mathbf{y}}}_t$ is a vector of reconciled values for each target variable, $\mathbf{h}$ is a vector that describes the functional form of model equality constraints, and $\mathbf{g}$ is a vector that describes the functional form of model inequality constraints [16,22].

Let us refer to the general HTS problem in Figure 1, where the relationship and subordination at the different hierarchy levels are illustrated [12]. For simplicity reasons, the measured and unmeasured variables were not labeled in the hierarchy. Every element (node) of the hierarchy (tree) is labeled as $y_{t,j}^{k,p}$, where t = 1 … n denotes the sample instance, k = 1 … K means the level of the hierarchy, p presents the node’s parent on the upper ($k-1$) hierarchy level, and j = 1 … $q_p^k$ is the number of child nodes at the k-th hierarchy level of the p-th parent.

This hierarchical structure means that the element at the root level ($y_t^0$) is the sum of the elements in the first level of the hierarchy, formally

$$y_t^0=\sum _{j=1}^{q^{1,1}}{y}_{t,j}^{1,1}$$

(6)

where $q^{1,1}$ is the number of child nodes at the first hierarchy level.

A general element on the k-th level ($y_{t,j}^{k,p}$) is the sum of its child nodes, formally written as

$$y_{t,j}^{k,p}=\sum _{i=1}^{q_j^{k+1,j}}{y}_{t,i}^{k+1,j}$$

(7)

The number of elements in the k-th hierarchical level ($q^k$) can be formalized as the sum of the number of child nodes in the upper hierarchy level, as follows:

$$q^k=\sum _{p=1}^{q^{k-1}}{q}_p^k$$

(8)

The number of nodes in the hierarchy is obtained if the number of nodes in each hierarchy level is added together.

$$q=\sum _{k=0}^K{q}^k$$

(9)

By stacking all the tree elements in a vector (${\mathbf{y}}_t$) based on Figure 1, the following is obtained:

$${\mathbf{y}}_t={[y_t^0,\underset{{\mathbf{y}}_{t,1}^T}{\underbrace{y_{t,1}^{1,1}\dots y_{t,q_1^1}^{1,1}}},\underset{{\mathbf{y}}_{t,2}^T}{\underbrace{y_{t,1}^{2,1},y_{t,q_1^2}^{2,1}\dots y_{t,1}^{2,q_1^1},y_{1,q_{q_1^1}^2}^{2,q_1^1}}},\dots ,\underset{{\mathbf{y}}_{t,K}^T}{\underbrace{\dots }}]}^T$$

(10)

where ${\mathbf{y}}_{t,K}$ represents the elements on the bottom hierarchy level (K-th level).

Every element at each hierarchy level can be calculated using a summation matrix ($\mathbf{S}$) and the bottom-level elements (${\mathbf{y}}_{t,K}$) as follows:

$${\mathbf{y}}_t=\left[\begin{array}{c}{y}_t^0\\ {\mathbf{y}}_{t,1}\\ ⋮\\ {\mathbf{y}}_{t,K}\end{array}\right]=\mathbf{S}{\mathbf{y}}_{t,K}$$

(11)

The coherence requirements within the hierarchy can be defined with the help of a summing matrix ($\mathbf{S}$), which dictates the method in which the bottom-level series aggregate at higher hierarchy levels. The element of the summation matrix ($s_{i,j}$) refers to the tree node element i and tree leaf element j. This sets the matrix element of a given node i to 1 or 0 depending on if it is an ancestor of the leaf element j or not, respectively. The dimensions of $\mathbf{S}$ are q × $q^K$, where q represents the total number of nodes in the hierarchy and $q^K$ represents the number of nodes at the bottom level.

For the optimized reconciliation of the samples, a connection must be established between the summation matrix ($\mathbf{S}$) and the equality constraint of the DR technique at each hierarchical level (Equation (4)). Due to the hierarchical structure, this equality constraint can be formalized as $\mathbf{A}{\tilde{\mathbf{y}}}_t=0$, where $\mathbf{A}$ is the incidence matrix. To create matrix $\mathbf{A}$, with dimensions q × q, the $\left[\mathbf{0}\mathbf{S}\right]$q × (q-$q^K$) null matrix and q × $q^K$ summation matrix ($\mathbf{S}$) need to be subtracted from the q × q identity matrix ($\mathbf{I}$).

$$\mathbf{A}=\mathbf{I}-\left[\mathbf{0}\mathbf{S}\right]$$

(12)

Based on this, the following relation can be obtained to formalize the linear constraint for DR using the summation matrix:

$$\mathbf{A}{\tilde{\mathbf{y}}}_t=(\mathbf{I}-\left[\mathbf{0}\mathbf{S}\right]){\tilde{\mathbf{y}}}_t=0$$

(13)

We also consider additional linear equality constraints which need to be satisfied in addition to the hierarchical constraints, as shown in Equation (14):

$${\mathbf{A}}^{*}{\tilde{\mathbf{y}}}_t={\mathbf{b}}^{*}$$

(14)

Each linear constraint can be grouped, as shown in Equation (15):

$$\begin{array}{c}\left[\begin{array}{c}\mathbf{A}\\ {\mathbf{A}}^{*}\end{array}\right]{\tilde{\mathbf{y}}}_t=\left[\begin{array}{c}\mathbf{0}\\ {\mathbf{b}}^{*}\end{array}\right]\end{array}$$

(15)

and, in a more compact form, we can rewrite this as follows:

$$\tilde{\mathbf{A}}{\tilde{\mathbf{y}}}_t=\tilde{\mathbf{b}}$$

(16)

The analytical solution of the DR (Equation (3)) with linear constraints (Equation (16)) is the following:

$${\tilde{\mathbf{y}}}_t=(\mathbf{I}-{\mathbf{V}}^{-1}{\tilde{\mathbf{A}}}^T{\left(\tilde{\mathbf{A}}{\mathbf{V}}^{-\mathbf{1}}{\tilde{\mathbf{A}}}^T\right)}^{-1}\tilde{\mathbf{A}}){\widehat{\mathbf{y}}}_t+{\mathbf{V}}^{-\mathbf{1}}{\tilde{\mathbf{A}}}^T{\left(\tilde{\mathbf{A}}{\mathbf{V}}^{-\mathbf{1}}{\tilde{\mathbf{A}}}^T\right)}^{-1}\tilde{\mathbf{b}}$$

(17)

where $\mathbf{I}$ is the identity matrix, $\tilde{\mathbf{A}}$ is the incidence matrix, and $\tilde{\mathbf{b}}$ are the constraint values. Equation (18) can be written in a more compact form as follows:

$${\tilde{\mathbf{y}}}_t=\mathbf{P}{\widehat{\mathbf{y}}}_t+\stackrel{\approx }{\mathbf{b}}$$

(18)

where $\mathbf{P}$ is the projection matrix and $\stackrel{\approx }{\mathbf{b}}$ is the correction term needed to satisfy the linear constraints.

$$\mathbf{P}=(\mathbf{I}-{\mathbf{V}}^{-1}{\tilde{\mathbf{A}}}^T{\left(\tilde{\mathbf{A}}{\mathbf{V}}^{-\mathbf{1}}{\tilde{\mathbf{A}}}^T\right)}^{-1}\tilde{\mathbf{A}})$$

(19)

$$\stackrel{\approx }{\mathbf{b}}={\mathbf{V}}^{-\mathbf{1}}{\tilde{\mathbf{A}}}^T{\left(\tilde{\mathbf{A}}{\mathbf{V}}^{-\mathbf{1}}{\tilde{\mathbf{A}}}^T\right)}^{-1}\tilde{\mathbf{b}}$$

(20)

In the case of linear regression models, Equation (21) becomes the following:

$${\tilde{\mathbf{y}}}_t=\mathbf{P}\mathit{\theta }{\mathbf{x}}_{e,t}+\stackrel{\approx }{\mathbf{b}}$$

(21)

where $\mathit{\theta }$ represents the model parameters and ${\mathbf{x}}_{e,t}$ represents the extended input features at the t-th sample instance [16].

2.2. Methods for Integrating Machine Learning and Data Reconciliation Techniques

Based on the available information, ML models can be developed in three ways, as presented in Figure 2. The first route is to develop an ML model without any DR. The second route is that if there is information about the measurement uncertainty, then the inverse covariance matrix can be defined with the standard deviation of the measurement errors, and then DR can be performed on the measured dataset. Since the measurements are independent from each other, the covariance matrix for DR contains non-zero values only in the diagonal elements. After DR, ML models can be developed using the reconciled dataset.

The third route is when there is no information about the measurement uncertainty, but the prediction errors can be used to fine-tune the model predictions.

In this case, first, the ML models are trained on the raw dataset of the measurements, and then they are reconciled with the model predictions, where the standard deviation of the prediction errors is included in the inverse covariance matrix used to perform DR. In this case, we assumed that the modeling errors are independent of each other due to the fact that an independent ML model was developed for each node. Therefore, in this case, the developed covariance matrix for DR contains non-zero values only on the diagonals. Besides this assumption, there can be cases where the modeling errors are not independent, and hence the covariance matrix for DR should be developed with care in these cases.

3. Modeling Results for Cases of Varying Complexities

Hierarchical modeling with optimal DR was performed in three different case studies with different complexity levels. The first case study is about predicting the mineral composition of different rock samples from spectral data to support oil exploration. This problem includes three hierarchy levels with nine elements.

In the second case, a Walmart commercial dataset included in a Kaggle competition was investigated, where reconciliation was performed between time series predictions on three hierarchy levels with 14 elements. The third case also includes an HTS analysis dataset used to predict regional waste deposition in Hungary on four hierarchy levels with more than 3000 elements.

3.1. Mineral Composition of the Rock Samples

The first case study uses the results of data generated during oil exploration and production as input data. The traditional approach to determining the mineral composition of rock samples has been to take field samples/drill cores collected during field investigations to the laboratory. Analytical data are prepared using traditional laboratory measurement methods, which are costly and time consuming, may include the use of hazardous substances, and are challenging due to the exploration of large areas [23]. From the point of view of technology development, the presented methodology is suitable for obtaining an accurate picture of the analysis during oil field drilling. This is a complex task for geologists in oil exploration and production. Geological analyses classify the different storage rocks into separate categories, which is particularly important from the point of view of geological interpretations [24]. Furthermore, knowing the detailed composition of the elements is also essential before executing the various efficiency-enhancing procedures of the wells that are already in operation.

Fourier-transform infrared spectroscopy (FT-IR) measurements can also determine the composition and properties of rocks. In these cases, chemometric (ML) models are necessary to analyze spectra [24]. These models offer the possibility of installing software sensors in drill heads, which allows for a detailed analysis of the composition of rocks in the field [25].

The mineral composition of the rocks used in the modeling is represented by mass percentage values, which must total 100%. Our task involves a hierarchical regression process to predict the specific mineral composition of the rocks, ensuring their sum equals 100%. The aims of this methodology are to develop an algorithm that outperforms traditional models while meeting the initial conditions. In this study, ML models were constructed for determining the mineral composition of rock samples. The target (dependent) variables $\mathbf{y}$ of the model are derived from X-ray diffraction (XRD) measurements, while the independent variables $\mathbf{X}$ are derived from FT-IR spectra.

In the first method, a partial least-squares (PLS) regression model was developed, creating a distinct model for each element’s composition. In the second method, DR was applied considering the uncertainty associated with traditional XRD measurement errors. For the third method, DR was performed on the model predictions using the uncertainty of these predictions. The training dataset consisted of 618 samples and applied 10-fold cross-validation to fine-tune the PLS parameters. The refined models were then evaluated on 305 unknown samples.

The dependent variables are the mineral compositions obtained from different rocks. Despite these traditional laboratory measurements being performed with an XRD device, the total of the mineral compositions does not add up to 100, as expected. The DR method is also applicable to conventional laboratory measurements. To calculate the standard deviation of the XRD measurements, the measurement uncertainty associated with traditional measurements reported in the literature was taken into account. Since 2000, the Reynolds Cup has been held every two years, allowing various institutions and companies to participate in lap sample measurements [26]. From these measured data, an aggregated result is compiled in

Table 1. Results of the XRD measurements for the 2006 Reynolds cup. * The results contain the average result of the first three places of the 2006 Reynolds cup (RC) 3-2 sample on the kidney [27]. ** SD: standard deviation.

	Actual	Submitted *	SD **
Quartz	29.9	30.67	0.77
K-feldpar and Plagioclase	8.6	8.0	0.6
Calcite	4.6	4.47	0.13
Kaolinite	15.0	15.7	0.7
Total clay without Kaolinite	20.2	19.17	1.03
Dolomite, Magnesite, Hematite Aragonite, Fluorite, Apatite …	21.7	21.73	0.03

, which shows the standard deviation of the XRD measurements for the six mineral compositions [27].

The hierarchical levels related to the rock case study and the relationship between parent and child nodes are presented in Figure 3. This figure presents the hierarchical structure at the highest level (level 0), where the total of the compositions must always be equal to 100. At the subsequent level (level 1), ${\mathbf{y}}_{t,1}^{1,1}$ presents a silicate, and ${\mathbf{y}}_{t,2}^{1,1}$ presents carbonate rock types. The lowest level (level 2) includes the variables of the six mineral compositions. In this case study, siderite is represented by quartz, K-feldspar and plagioclase, calcite, and kaolinite. The total clay without kaolinite, along with dolomite, magnesite, and hematite, constitutes the carbonate group, and the t index in this case study is the number of samples, not a timestamp.

The hierarchical relationship in the mineral composition is shown in detail in the summation matrix ($\mathbf{S}$).

$$\mathbf{S}=\left[\begin{array}{cccccc}1& 1& 1& 1& 0& 0\\ 0& 0& 0& 0& 1& 1\\ 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]$$

(22)

The incidence matrix ($\tilde{\mathbf{A}}$) is

$$\tilde{\mathbf{A}}=\left[\begin{array}{cccccccc}1& 0& -1& -1& -1& -1& 0& 0\\ 0& 1& 0& 0& 0& 0& -1& -1\\ 1& 1& 0& 0& 0& 0& 0& 0\end{array}\right]$$

(23)

In this case study, the constraint that the sum of the silicate and carbonate content should be 100 is included in matrix $\tilde{\mathbf{b}}$.

$$\tilde{\mathbf{b}}=\left[\begin{array}{c}0\\ 0\\ 100\end{array}\right]$$

(24)

$\mathbf{V}$ is a 8 × 8 matrix that includes the standard deviation of all variables in the diagonal of the matrix. So, the standard deviation of the measurements in the second method is taken from Table 1, and in the third method, it is the standard deviation of the ML prediction errors.

In the following, the results of the development of ML models for each method are presented through the average prediction errors and through the hierarchical balance error. The results of different hierarchical levels and different methods were compared with the root-mean-squared error of the prediction ($RMSE$) value, which indicates the accuracy of each estimate. At the first level, the third method performed the best for silicate ($RMSE$ = 4.577%), and the first method was the best ($RMSE$ = 4.154%) for carbonate, but here, the balance error was not met. The most significant difference is seen between the kaolinite $RMSE$ values at the two levels, where the first method performed the best ($RMSE$ = 1.584%), but in this case, the balance errors were not met either.

The results for the rock samples are summarized in the bar plots in Figure 4.

Table 2. Prediction results of the first, second, and third methods for the rock case study.

level 0	Total
y_real	100
$\widehat{y}$_1st	103.24
$\tilde{y}$_2nd	102.88
$\tilde{y}$_3rd	100
level 1	Silicate				Carbonate
y_real	79				21
$\widehat{y}$_1st	66.55				36.69
$\tilde{y}$_2nd	66.20				36.68
$\tilde{y}$_3rd	65.22				34.78
level 2	Quartz	K-feldp. and Plagio.	Kaolin.	Clay without Kaolin.	Calcite	Dolom. Magne. Hemat. …
y_real	51	14	2	12	11	10
$\widehat{y}$_1st	37.57	12.81	1.67	12.19	19.32	14.33
$\tilde{y}$_2nd	37.56	12.73	1.60	12.72	19.36	14.02
$\tilde{y}$_3rd	37.70	13.04	2.16	12.32	19.78	15.00

presents a chosen subset from the dataset of 305 rock samples, illustrating the actual composition of the samples and the outcomes achieved by each method. At level 0, the DR error is 3.24 for method 1 and 2.88 for method 2, while the third method exhibits no DR error. At level 1, method 1 has an error of 3.24, method 2 has an error of 2.88, and method 3 shows no error. At the lowest level (level 2), method 1 has a DR error of 2.11, method 2 has an error of 2.01, and method 3 has no error. These results clearly demonstrate the advantage of the third method, as the initial conditions are also satisfied, considering the errors.

The developed method helps to better predict the rock composition of oil fields. It can also be used to balance the material flows of products produced during chemical processes. It is also suitable for monitoring and controlling chemical and separation technology processes. Incorporating the correlations between predicted values facilitates the embedding of ML models into an industrial environment. In addition, the DR technique increases the usability of the models.

3.2. Retail Sales Forecasting

The second case study discussed in the article involves an HTS analysis conducted on data from the Walmart shopping chain, part of a data analysis competition announced in 2020 as M5. The competition aimed to improve the accuracy of forecasting and empirically compare various forecasting methods [28]. The input data consist of unit sales for 3049 products in the US.

These products are sold in 10 shops across three states (CA, TX, WI).

Locations such as states and shops were selected carefully, representing different waiting habits, purchasing dynamics, durable consumer goods, and fast- and slow-moving products. The M5 dataset is an excellent choice due to its meaningful hierarchies and cross-sectional levels. The daily data span from 29.01.2011 to 19.06.2016 (1969 days), with the training set data spanning from 29.01.2011 to 24.04.2016 (1913 days), a 27-day validation period spanning from 25.04.2016 to 22.05.2016, and a 28-day test period spanning from 23.05.2016 to 19.06.2016. Explanatory variables include calendar information, sales prices, and promotional activities [29]. Figure 5 illustrates a hierarchical diagram of the M5 case study, and the related incidence matrix is shown in Equation (24).

In this case, the hierarchy of levels 0, 1, and 2 has ten elements at the lowest level (level 2: shop), three elements in a 4-3-3 split at the middle level (level 1: state) above, and one element at the top (level 0: country) in $\mathbf{S}$.

The hierarchical relationship in the sales dataset is shown in detail in the summation matrix ($\mathbf{S}$).

$$\mathbf{S}=\left[\begin{array}{cccccccccc}1& 1& 1& 1& 1& 1& 1& 1& 1& 1\\ 1& 1& 1& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 1& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 1& 1\\ 1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]$$

(25)

The incidence matrix ($\tilde{\mathbf{A}}$) is

$$\tilde{\mathbf{A}}=\left[\begin{array}{cccccccccccccc}0& 1& 0& 0& -1& -1& -1& -1& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0& -1& -1& -1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& -1& -1& -1\end{array}\right]$$

(26)

In this case study, the sum of the lower levels should be output to predict the upper level. Matrix $\tilde{\mathbf{b}}$ can be written in the following form:

$$\tilde{\mathbf{b}}=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$$

(27)

Figure 6 illustrates the performance of several forecasting models along with their hierarchical results. The metric used is the $RMSE$, displayed on the y axis. Different DR techniques are shown on the x axis, with colours indicating the hierarchical levels. In general, the middle-out method has the highest $RMSE$ (2.421), while the bottom-up method has the lowest (2.241). Despite the small overall differences in the $RMSE$, the bottom-up method is the best performer. At the country level, the top-down method has the highest $RMSE$ (5.408), whereas the bottom-up method, with a significantly lower $RMSE$ (4.822), is more efficient and suitable for country-level forecasting. At the state level, the top-down method has the highest $RMSE$ (2.911), while the bottom-up method has the lowest (2.487). In this category, the bottom-up method is the most effective. At the shop level, there are slight differences among all four methods. The top-down method has the highest $RMSE$ (1.656), while the bottom-up method has the lowest (1.569).

With the introduction and development of the analytical DR solution in this article, the overall $RMSE$ from the sales data is 2.305, which is in the mid-range compared to the other methods. At the country level, the analytical DR solution’s value is 4.848, ranking second after the bottom-up method, indicating strong performance. At the state level, the DR analytical solution’s value is 2.839, also showing good performance relative to the other methods. At the shop level, the DR analytical solution’s value is 1.633, which is significantly lower than that of the other methods, though the RMSE of the bottom-up method is even lower. In this case study, the initial dataset satisfies the hierarchical balance and does not require reconciliation. Consequently, the second method is less pertinent. Only the first and third methods are suitable in this scenario. By emphasizing one of the forecasted outcomes on 18.05.2016, the results of the first and third methods are summarized in

Table 3. Prediction results of the first and third methods in the sales case study.

level 0	Country
y_real	22.0
$\widehat{y}$_1st	14.102
$\tilde{y}$_3rd	13.938
level 1	CA				TX			WI
y_real	11.0				1.0			10.0
$\widehat{y}$_1st	6.034				1.569			4.416
$\tilde{y}$_3rd	6.278				2.124			5.536
level 2	CA_1	CA_2	CA_3	CA_4	TX_1	TX_2	TX_3	WI_1	WI_2	WI_3
y_real	2.0	4.0	1.0	4.0	0.0	1.0	0.0	1.0	3.0	6.0
$\widehat{y}$_1st	1.656	1.314	1.685	0.432	0.244	0.052	0.656	2.627	1.872	1.885
$\tilde{y}$_3rd	1,729	1.376	1.879	1.294	0.458	0.544	1.122	2.442	1.501	1.593

. The results show that the third method gives similar results to the first method, but the third method involves comparing the data. At the national level, the third method provides the best results, but bottom-up is the second worst method. At the state level, the bottom-up method gives the best result, with the third method in the middle. There are slight differences at the shop level, with the lowest error obtained by the DR techniques given by the bottom-up method and the second lowest by the third method.

When comparing the results tested on the data of 18 May 2016, the analytical solution of the HTS analysis using the DR technique (third method) showed little difference from the results of the first methodology. Although the third methodology was among the top performers at the country level, it ranked in the lower half for the state and shop cases.

3.3. Waste Management Hierarchical Time Series Prediction with Data Reconciliation

Waste management is a critical issue in Hungary due to the increasing amount of solid waste generated yearly. Effective waste management strategies require the accurate forecasting of waste generation, which can be achieved through time series analysis. Much of the literature focuses on predicting waste quantities as accurately as possible because precise solid waste forecasts are crucial for the circular economy, aiming to maximize recycling and enhance energy efficiency [30]. An integral part of these forecasts involves predicting waste amounts at collection points, which is vital from a workload perspective [31]. Eryganov et al. employed reconciliation in HTS forecasting, which is known to improve the quality of initial forecasts [32]. Moreover, HTS and DR technologies were applied in analyzing hazardous waste, utilizing a tree structure that mirrors the organizational layout of a region (including regions, micro-regions, and their sections), with the autoregressive integrated moving average (ARIMA) algorithm being used for predictions [33]. However, waste data are often incomplete and inconsistent, which results in unreliable predictions. HTS analysis has proven to be a powerful method for forecasting when dealing with multiple time series at various levels of aggregation. Additionally, DR techniques can enhance accuracy by resolving data inconsistencies. In this case study, HTS analysis and DR techniques were investigated to forecast solid waste generation within Hungary’s hierarchical waste management system.

The solid waste data for Hungary were sourced from the Hungarian Central Statistical Office database. This comprehensive dataset includes the volume of solid waste at the national level for all 19 counties, including Budapest (the capital city), 175 districts, and 3155 settlements. Consequently, the HTS analysis incorporates levels 0, 1, 2, and 3.

The modeling process used the Holt–Winters exponential smoothing (HWES) algorithm in this research, with annual data aggregation [34]. Data on solid waste quantities were collected from 2010 to 2022. The HWES models were trained from 2010 to 2019 and forecasted until 2022, and the measured data are available for 2020, 2021, and 2022. A settlement is selected in

Table 4. Summary of the forecast results of the first and third routes of a selected settlement, as well as its district and county levels, in the case study concerning the amount of solid waste in Hungary in 2020, 2021, 2022 (tons).

	2020	2021	2022
level 0	Total
$y\_real$	3,301,482.3	3,350,245.0	3,215,457.4
$\widehat{y}$_1st	3,254,760.8	3,256,817.8	3,258,874.8
$\tilde{y}$_3rd	3,254,503.5	3,256,560.6	3,258,617.7
level 1	County
$y\_real$	119,844.3	120,933.9	116,180.1
$\widehat{y}$_1st	116,495.1	117,577.4	118,659.6
$\tilde{y}$_3rd	121,780.4	122,861.6	123,942.9
level 2	District
$y\_real$	15,558.7	16,903.2	15,751.4
$\widehat{y}$_1st	14,640.5	15,013.5	15,386.6
$\tilde{y}$_3rd	14,767.7	15,140.1	15,512.5
level 3	Settlement
$y\_real$	1341.7	1427.4	1355.0
$\widehat{y}$_1st	1037.7	1031.7	1025.7
$\tilde{y}$_3rd	1038.3	1032.3	1026.3

, and the forecast obtained by the first ML method and the third ML + DR method is displayed. The second method is not considered because in the first case, the balance error is met, so the DR of the input data of the ML models is not relevant in this case study. The reliability of the methods was examined for the years 2020, 2021, and 2022. The selected settlement is a randomly selected small town in Hungary, located in Veszprém County.

Figure 7 shows the RMSE value of the model errors, taking into account the national level, including all counties, all districts, and all settlements, in 2020, 2021 and 2022. Figure 7 clearly shows that the biggest errors are at the county level, and these errors are lower if DR is applied after ML. Taking into account the years, the highest error is seen in 2021, and the lowest is seen in 2022.

4. Conclusions

This study introduces a novel approach that combines data reconciliation techniques with machine learning. In addition to the already widespread and established methods, this research formalizes a unique analytical approach. In the case of hierarchical systems, machine learning models do not satisfy various constraints arising from hierarchical structures, and the balance equations between hierarchical levels must be satisfied when modeling hierarchical systems. In this study, three different machine learning model development methods were investigated and compared: (1) model building without data reconciliation, (2) model building with reconciled measurement data based on the measurement uncertainty, and (3) model building with direct reconciling of the machine learning model predictions based on the modeling errors of each model.

In this research, each method was examined through three case studies of different complexities. The results show that the fine-tuning of model predictions with optimal data reconciliation helps to improve the accuracy and reliability of model predictions, and all necessary constraints can be satisfied with this technique.

The presented method enables data integration from different measurement techniques, incorporating errors from different methods into the machine learning prediction process. In subsequent developments, this approach will allow us to create specific algorithms that take estimation errors into account and manage them efficiently. By remedying the imbalance, the goal was to develop more accurate, flexible, scalable, and usable models.

The limitations of the method presented in this study are that the hierarchical structure of the system, the quality of the data, and the reliability of the sources must be known precisely, and all the essential limitations of the system must be determined. If the hierarchical structure or constraints of the system change, the model structure must also be updated.