Long-Term Electricity Demand Prediction Using Non-negative Tensor Factorization and Genetic Algorithm-Driven Temporal Modeling

Electricity demand forecasting has been studied extensively using statistical, machine learning, and hybrid approaches. Comprehensive reviews by Hong et al. [1] and Lago et al. [2] summarize advances in probabilistic modeling, neural networks, and ensemble methods. Many of these rely on external variables such as temperature, population, or macroeconomic indicators to enhance short- and long-term forecast accuracy [3, 4].

Other studies have focused on incorporating sociopolitical or regulatory dimensions, such as policy-driven demand shifts and renewable energy adoption [5, 6]. Time series forecasting using grey-Markov models or seasonal ARIMA has been applied in data-scarce environments [7, 8], showing that consumption-only models can still yield useful predictions.

Short-term load forecasting using ARIMA and neural networks has also been demonstrated on pure consumption data [9]. Fan and Hyndman’s Monash Electricity Forecasting Model (MEFM) [10] uses semiparametric additive models with external factors for national-scale long-term forecasting, but its adaptability to variable-scarce scenarios remains unclear.

Tensor decomposition methods—particularly Non-negative Tensor Factorization (NTF)—have been used for pattern extraction and clustering in electricity data [11, 12, 13], and have also proven effective in domains like image processing and hyperspectral analysis [14, 15]. Sparse and incremental extensions of NTF have been developed for recommendation systems and completion tasks [16, 17], though these do not address time series forecasting.

Xu et al. [18] provide a comprehensive survey of tensor-based forecasting. Zarnaz et al. [19] demonstrate that NTF can improve generalization in financial time series forecasting, highlighting the potential of low-rank latent modeling in structured domains.

In terms of forecasting algorithms, classical ARIMA models typically rely on manually tuned parameters based on domain expertise or statistical criteria [20]. Genetic Algorithms have been adopted to automate parameter search in hybrid models and epidemic prediction tasks [21, 22]. Their flexibility makes them well-suited for non-convex, high-dimensional search spaces like those arising in time series modeling.

While deep learning methods have become increasingly prominent in short-term forecasting [23], their application to long-term national-scale prediction remains limited due to their data and compute requirements, and lack of interpretability.

In this context, our work integrates NTF and Genetic Algorithm-based time series modeling to support interpretable, scalable, and data-efficient long-term electricity demand prediction, without requiring external covariates.

The dataset used in this study is derived from publicly available electricity consumption records published by the Federation of Electric Power Companies of Japan. It includes annual consumption data from 1963 to 2015, recorded across ten regional electric utilities and detailed industry classifications.

The industry labels consist of 36 categories organized hierarchically. These are aggregated into 10 major industrial sectors: Mining, Manufacturing, Chemical Industry, Petroleum and Coal Products Manufacturing, Rubber Product Manufacturing, Ceramic and Stone Product Manufacturing, Steel Industry, Non-ferrous Metal Manufacturing, Machinery and Equipment Manufacturing 1, and Machinery and Equipment Manufacturing 2.

The full dataset is modeled as a third-order non-negative tensor $\mathcal{X}\in\mathbb{R}^{I\times J\times K}$, where $I=10$ (electric utilities), $J=10$ (industries), and $K=53$ (years). This tensor serves as the input to the Non-negative Tensor Factorization (NTF) process described in Section 1.2.1.

To evaluate forecasting performance, the dataset is split into training, validation, and test sets along the year mode. Two test cases are constructed:

• •

  Case A2010: Training data from 1963–2003, validation data from 2004–2009, test data from 2010–2015.

• •

  Case A2000: Training data from 1963–1993, validation data from 1994–1999, test data from 2000–2005.

The training tensor $\mathcal{X}_{\text{train}}$ is decomposed using the canonical polyadic (CP) NTF method. The number of update iterations is fixed at 100, and the rank $R$ is selected by evaluating reconstruction accuracy (as defined in Section 1.2.1) across different candidate ranks.

Figure 1 shows the reconstruction accuracy for each rank. Based on a trade-off between reconstruction quality and computational cost during the Genetic Algorithm (GA) optimization, we choose $R=8$ for all experiments.

Following tensor decomposition, we apply time series forecasting to each column of the annual factor matrix $\mathbf{C}\in\mathbb{R}^{K\times R}$. Each time series $\mathbf{c}_{r}$ is modeled using an ARIMA($p_{r},d_{r},q_{r}$) process, where hyperparameters are optimized via a Genetic Algorithm as outlined in Section 1.2.2.

The predictive annual factor matrix for the validation period is denoted as:

The reconstructed tensor for the validation period is then:

The GA selects the best set of ARIMA parameters by minimizing the mean squared error (MSE) between the sum of predicted consumption and the ground-truth validation tensor.

To evaluate the effectiveness of our NTF+GA-based method, we define two baseline methods:

• •

  w/o NTF: Apply ARIMA(3,2,3) directly to each of the 100 time series corresponding to (utility, industry) combinations, bypassing tensor decomposition.

• •

  w/o GA: Apply fixed ARIMA(3,2,3) to each latent annual factor obtained from NTF without performing GA-based hyperparameter search.

These settings allow us to isolate the impact of NTF and GA independently.

We define four evaluation patterns depending on the test case and fitness criterion used during GA optimization:

• •

  Pattern A2010: GA optimized using MSE on total consumption (Case A2010).

• •

  Pattern B2010: GA optimized using reconstruction accuracy (Frobenius norm) (Case A2010).

• •

  Pattern A2000: GA optimized using MSE (Case A2000).

• •

  Pattern B2000: GA optimized using reconstruction accuracy (Case A2000).

The ARIMA hyperparameters used in fixed settings are shown in Table 1.

We evaluate the proposed method and the two baseline methods (w/o NTF, w/o GA) across four experimental patterns described in Section 2.5. Tables 2–5 summarize the MSE and reconstruction accuracy for each method. Figures 2–5 illustrate the annual total electricity demand predicted by each method compared with the ground truth.

To assess the stability of the proposed method, we examine how reconstruction accuracy varies with different random seed values used during NTF initialization. Tables 6–9 show reconstruction accuracy on validation data for five different seeds: 10, 20, 30, 42, and 45. The corresponding prediction results, averaged over the top three reconstruction runs, are plotted in Figures 6–9.

Comparing Patterns A2010 and B2010, where the only difference lies in the GA fitness function (MSE versus reconstruction accuracy), we observe contrasting behavior. As shown in Tables 2 and 3, optimizing for reconstruction accuracy leads to higher structural fidelity (0.825 vs. 0.811), while optimizing for MSE yields better forecasting performance in terms of true error (lower MSE).

Interestingly, in Patterns A2000 and B2000 (Tables 4 and 5), both MSE and reconstruction accuracy improve when the GA is optimized for reconstruction accuracy. This suggests that the optimal fitness criterion may depend on the data range, underlying consumption trends, or the degree of seasonality and noise.

To understand the generalization ability of each model, we analyze the total number of degrees of freedom (DOF) associated with each configuration. Table 10 summarizes the DOF for each experimental pattern in the proposed method, computed as:

We find that configurations with fewer DOF (e.g., B2010 and B2000) generally result in better reconstruction accuracy and sometimes better MSE, supporting the hypothesis that model simplicity contributes to generalization.

Notably, the w/o NTF baseline exhibits the highest reconstruction accuracy in several patterns. However, this is likely due to its increased effective DOF per tensor element, as ARIMA is applied independently to each utility–industry pair. In contrast, our method enforces structural constraints through NTF, which may sacrifice reconstruction accuracy but improves generalization by capturing shared latent patterns.

As shown in Tables 6–9, the NTF-based method exhibits variability in reconstruction accuracy depending on the random seed used during initialization. This is consistent with known behavior of non-convex matrix/tensor factorization methods, where local minima can lead to different decompositions.

Nonetheless, Figures 8 and 9 show that averaging over top-performing initializations provides stable predictions. This suggests that incorporating multiple initialization trials or ensemble methods could further enhance robustness.

By reducing the dimensionality of the original data through NTF, the proposed method achieves interpretable factor matrices that reflect underlying temporal trends. The use of GA further automates the selection of appropriate time series model complexity for each factor, allowing for adaptive modeling without manual tuning.

Moreover, the framework is scalable and extensible. It can be adapted to include additional modes (e.g., geographic regions, customer classes) or applied to other domains such as water or gas demand, provided a structured tensor representation is available.