Spectral-Temporal MoE: An RL-Driven Dual-Domain Transformer for Efficient Multi-Horizon Electrical Load Forecasting

1. Introduction

1.1. Background

Electrical load forecasting across multiple temporal horizons is fundamental to modern power grid operations, enabling utilities to optimize generation scheduling, demand response programs, and energy trading strategies. Accurate predictions spanning short-term (1-24 hours), medium-term (1-7 days), and long-term (1-4 weeks) horizons are essential for maintaining grid stability, minimizing operational costs, and integrating renewable energy sources [1]. Traditional forecasting approaches based on statistical models such as ARIMA and exponential smoothing demonstrate limited effectiveness when confronted with the complex nonlinear dynamics, multi-scale periodicities, and sudden load fluctuations characteristic of modern electricity consumption patterns. These limitations become particularly acute as power grids evolve toward smart grid architectures with distributed generation, electric vehicle integration, and dynamic pricing mechanisms, creating unprecedented demands for accurate and efficient multi-horizon forecasting capabilities that can capture both short-term volatility and long-term trends simultaneously [2].

The emergence of Transformer architectures has revolutionized time series forecasting by introducing powerful attention mechanisms capable of modeling long-range dependencies without the vanishing gradient problems inherent in recurrent neural networks. Recent studies have demonstrated promising results applying Transformers to temporal prediction tasks, with innovations such as PatchTST achieving state-of-the-art performance through channel-independent patching strategies [3,4,5]. However, fundamental questions remain regarding whether Transformers truly capture the underlying temporal structures or merely exploit dataset-specific patterns through sophisticated pattern matching. Furthermore, pure time-domain processing faces inherent challenges in handling periodic components, harmonic structures, and multi-scale seasonality that characterize electrical load profiles. The computational complexity of O(L²) attention operations over long sequences imposes practical constraints on processing extended historical windows necessary for capturing annual cycles and multi-week patterns, limiting the scalability of dense Transformer architectures for large-scale deployment in real-time forecasting systems [6,7].

Frequency-domain analysis offers complementary perspectives for time series modeling by explicitly decomposing signals into constituent periodic components, enabling direct manipulation of harmonic structures and seasonal patterns [8,9]. Recent research has explored integrating frequency representations with Transformer architectures, demonstrating improved capacity to capture multi-scale periodicities and reduce spurious correlations between distant time steps [10]. However, existing frequency-domain methods face a fundamental limitation: the standard Discrete Fourier Transform produces frequency grids determined solely by input sequence length, creating misalignment between the input spectrum and the complete time series that includes future predictions. Simultaneously, mixture-of-experts (MoE) architectures have gained traction as scalable solutions for conditional computation, enabling models to route different inputs to specialized sub-networks while maintaining sparse activation patterns [11,12,13]. Recent surveys highlight the potential of MoE frameworks in handling heterogeneous data distributions and task-specific patterns through learned expert specialization [14,15]. Despite these advances, existing approaches predominantly rely on static gating mechanisms such as top-k routing or noisy gating, which lack the ability to dynamically adapt expert selection based on evolving temporal patterns and fail to jointly optimize routing decisions with forecasting objectives across multiple time scales [16,17].

1.2. Motivation and Contributions

Recent advancements in Transformer architectures, frequency-domain modeling, and mixture-of-experts systems have individually demonstrated significant potential for enhancing time series forecasting capabilities. However, the integration of these technologies into a unified framework that synergistically combines dual-domain processing, sparse expert networks, and dynamic routing mechanisms remains largely unexplored. Existing time series forecasting methods primarily focus on optimizing model architectures within either the time or frequency domain, without considering the benefits of parallel dual-domain processing with domain-specific expert specialization. Similarly, current MoE approaches in temporal modeling typically employ static routing strategies that fail to adapt to the diverse and evolving patterns across different forecasting horizons and temporal scales. The lack of integrated solutions that jointly optimize dual-domain representation learning, sparse expert activation, and dynamic routing under forecasting objectives represents a significant gap in the development of next-generation multi-horizon prediction systems for critical applications such as power grid management.

Despite the promising developments in each individual technology domain, several critical research gaps remain unaddressed in the context of multi-horizon electrical load forecasting systems.

1.

Existing Transformer-based forecasting frameworks process time series exclusively in the temporal domain, failing to exploit the complementary strengths of frequency-domain representations for capturing periodic components and harmonic structures. Moreover, standard frequency transformation approaches suffer from spectral misalignment between input sequences and complete time series, limiting their effectiveness for multi-step-ahead prediction tasks.

2.

Current mixture-of-experts architectures in time series modeling rely on static gating mechanisms that route inputs based on learned but fixed assignment patterns, lacking the capability to dynamically adapt expert selection to varying temporal patterns, forecasting horizons, and evolving system dynamics. The absence of reinforcement learning-based routing prevents simultaneous optimization of prediction accuracy, routing stability, and expert diversity.

3.

Existing research predominantly treats model architecture design, expert specialization, and forecasting optimization as separate problems, failing to exploit the synergies and interdependencies among dual-domain processing, sparse computation, and dynamic routing. The lack of unified frameworks that jointly optimize these components through end-to-end training limits the potential benefits achievable through integrated system design.

To address these research gaps, this work proposes MoE-Transformer, a reinforcement learning-driven dual-domain mixture-of-experts framework for efficient multi-horizon electrical load forecasting. Our integrated approach combines Extended DFT-aligned frequency processing, parallel time-domain modeling, independent expert networks in both domains, and RL-based dynamic routing to create an intelligent and adaptive forecasting system. Specifically, we make the following contributions.

• We establish a dual-domain Transformer architecture that processes time series simultaneously in time and frequency domains through parallel encoder pathways. Our framework introduces Extended DFT that aligns input spectrum with the complete series frequency grid by computing F[k] = ${\sum }_{n=0}^{L-1}x\left[n\right]e^{-2\pi ikn/(L+T)}$, fundamentally resolving the spectral misalignment problem in frequency-domain forecasting.
• We design independent mixture-of-experts modules for frequency (F-MoE) and time domains (T-MoE), with ${\mathrm{K}}^F$=4 frequency experts and ${\mathrm{K}}^T$=4 time experts per Transformer layer. Each expert implements specialized two-layer networks optimized for domain-specific pattern recognition, enabling sparse parameter activation while maintaining model capacity through conditional computation.
• We develop a reinforcement learning-based routing framework that formulates expert selection as independent Markov Decision Processes in both domains. Our multi-objective reward function balances prediction accuracy, routing stability, and expert diversity through weighted terms $r_t^{F,m}=-{\alpha }_1^F·{\mathcal{L}}_m-{\alpha }_2^F·C_s^{F,m}+{\alpha }_3^F·D_e^{F,m}$, with MAPPO policies learning domain-specific routing strategies through coupled training dynamics.
• We conduct comprehensive experiments across five long-term forecasting benchmarks (ETTh1, ETTm1, Weather, Electricity, Traffic) and four prediction horizons (96, 192, 336, 720 steps), demonstrating that MoE-Transformer achieves 50.9-56.9% MSE reduction over state-of-the-art baselines while delivering 60% faster inference and 40% memory reduction through sparse expert activation. Ablation studies validate each component’s contribution, with RL routing providing 39.5-47.2% improvement over static gating methods.

The remainder of this paper is organized as follows. Section 2 reviews the existing literature on Transformer-based time series forecasting, frequency-domain modeling, and mixture-of-experts architectures. Section 3 describes our dual-domain MoE-Transformer framework and reinforcement learning-based routing mechanism. Section 4 and Section 5 present the experimental setup, performance evaluation, and ablation studies. Section 6 concludes the paper and discusses future research directions.

2. Related Work

2.1. Transformer-Based Time Series Forecasting

Table 1. Contrasting our work to existing studies.

	Ref	[18]	[19]	[20]	[21]	[22]	[23]	[24]	[25]	[26]	[27]	[28]	Proposed work
Feature
Temporal modeling capability		√	√	√	√	√	√	√					√
Frequency-domain processing		√				√	√	√					√
Sparse expert networks									√	√	√	√	√
Dynamic routing mechanism											√		√
Multi-horizon forecasting		√	√	√	√	√	√	√					√

Table 1. Contrasting our work to existing studies.

	Ref	[18]	[19]	[20]	[21]	[22]	[23]	[24]	[25]	[26]	[27]	[28]	Proposed work
Feature
Temporal modeling capability		√	√	√	√	√	√	√					√
Frequency-domain processing		√				√	√	√					√
Sparse expert networks									√	√	√	√	√
Dynamic routing mechanism											√		√
Multi-horizon forecasting		√	√	√	√	√	√	√					√

Transformer-based architectures have emerged as powerful tools for long-term time series forecasting, leveraging self-attention mechanisms to capture complex temporal dependencies without the limitations of recurrent neural networks [3,4,15,18,29,30]. Zhao et al. [18] proposed TFformer, a highly interpretable forecasting model that decomposes time series into low-frequency trend components and high-frequency periodic components through frequency decomposition. Their approach employs sequential frequency attention to enhance periodic information and sequential periodic matching attention for future pattern prediction, achieving significant performance improvements across six datasets while maintaining low time complexity through sequence-level attention mechanisms that preserve temporal dependencies. Kumar et al. [21] developed a spatio-temporal parallel transformer (STPT) model for traffic prediction, utilizing multiple adjacency graphs passed through coupled graph transformer-convolution network units operating in parallel to generate noise-resilient embeddings. Their framework demonstrates superior performance on four real-world traffic datasets and extends applicability to Covid-19 prediction tasks. Chen et al. [31] addressed irregular time series modeling through Contiformer, introducing continuous-time attention mechanisms for handling non-uniformly sampled temporal data.

The integration of channel processing strategies with Transformer architectures has gained attention for improving multivariate time series forecasting capabilities [19,20,32,33]. Han et al. [19] proposed MCformer, introducing an innovative Mixed Channels strategy that combines the data expansion advantages of channel independence with the ability to mitigate interchannel correlation forgetting. Their model blends a specific number of channels using attention mechanisms to effectively capture interchannel correlation information when modeling long-term features for massive IoT-generated time-series data, demonstrating that the Mixed Channels strategy outperforms pure channel independence approaches in multivariate forecasting tasks. Li et al. [20] developed channel independence bidirectional gated Mamba with interactive recurrent mechanism (CIBG-Mamba-IRM), enhancing gated recurrent neural networks by integrating Mamba into channel-independent recurrent units. Their framework processes each data channel independently through three Mamba-enhanced recurrent units with dynamic interchannel exchange, employing channel-specific adaptive state transition functions through the interactive recurrent mechanism. Experiments on new energy vehicle air conditioner datasets and five public benchmarks validated their model’s effectiveness in learning long-term dependencies and handling multidimensional time series data, achieving state-of-the-art results. Fan et al. [32] contributed DEWP, applying deep expansion learning principles to wind power forecasting scenarios with improved predictive accuracy.

Despite these advances in Transformer-based forecasting frameworks, existing approaches face critical limitations when applied to multi-horizon electrical load prediction scenarios. Most current methods process time series exclusively in the temporal domain without exploiting the complementary strengths of frequency-domain representations for capturing periodic components and harmonic structures. Furthermore, existing architectures lack dynamic mechanisms for adaptive computation allocation, relying on dense model activation patterns that limit scalability for real-time forecasting systems requiring simultaneous handling of multiple temporal horizons.

2.2. Frequency-Domain Time Series Modeling

Frequency-domain analysis has emerged as a powerful complement to temporal modeling for capturing periodic patterns and harmonic structures in time series forecasting [22,23,34,35,36,37]. Xu et al. [22] proposed a Frequency decomposition and Patch modeling Framework (FPF) that transforms input sequences to the frequency domain through Fast Fourier Transform and designs frequency masks to decompose data into high-frequency and low-frequency components for extracting fast-changing patterns and trend information respectively. Their dual patch modeling block processes high-frequency components through MLP-based patch enhancement to capture local features, while low-frequency components are modeled by Transformer-based patch mixing to capture global dependencies and cross-patch correlations. Comprehensive experiments on seven real-world datasets including ETT, Traffic, Electricity, and Weather demonstrated superior forecasting performance. Zhang et al. [23] developed CTFNet, a lightweight single-hidden layer feedforward neural network combining convolution mapping and time-frequency decomposition with three distinctive characteristics: time-domain feature mining based on matrix factorization to capture long-term correlations, multitask frequency-domain feature mining integrating multiscale dilated convolutions for simultaneous global and local context extraction, and highly efficient training with fast inference speed. Their empirical studies with nine benchmark datasets showed that CTFNet reduces prediction error by 64.7% and 53.7% for multivariate and univariate time series respectively compared to state-of-the-art methods. Zhou et al. [34] contributed Fourier graph convolution Transformer for financial multivariate time series forecasting, combining spectral graph analysis with attention mechanisms.

Adaptive temporal-frequency modeling mechanisms have gained attention for handling dynamic periodic patterns and improving long-term forecasting capabilities [24,38,39,40,41]. Yang et al. [24] proposed an adaptive temporal-frequency network (ATFN), an end-to-end hybrid model incorporating deep learning networks and frequency patterns for mid- and long-term forecasting. Their framework employs an augmented sequence-to-sequence model to learn trend features of nonstationary time series, a frequency-domain block to capture dynamic periodic patterns, and an adaptive frequency mechanism consisting of phase adaptation, frequency adaptation, and amplitude adaptation for mapping the frequency spectrum of the current sliding window to the forecasting interval. Experimental results on synthetic and real-world data with different periodic characteristics indicated promising performance and strong adaptability for long-term prediction tasks. Kumari et al. [38] developed a Fourier-driven lightweight token mixing model for efficient time series forecasting, reducing computational complexity while maintaining prediction accuracy. Zhang et al. [39] extended time-frequency analysis to anomaly detection scenarios through TFAD, demonstrating the broader applicability of frequency-domain decomposition techniques.

Despite these advances in frequency-domain modeling frameworks, existing approaches face critical limitations when applied to multi-horizon forecasting with diverse temporal patterns. Most current methods employ fixed frequency decomposition strategies without dynamic adaptation mechanisms for routing different pattern types to specialized processing modules. Furthermore, existing architectures lack sparse computation allocation through conditional expert networks, limiting their scalability and efficiency for large-scale deployment in real-time forecasting systems requiring simultaneous handling of multiple frequency components and temporal scales.

2.3. Mixture-of-Experts and Dynamic Routing Mechanisms

Mixture-of-Experts architectures have emerged as a powerful paradigm for scaling model capacity while maintaining computational efficiency through sparse expert activation [25,26,28,42]. Zhang et al. [25] presented a comprehensive review highlighting MoE’s ability to significantly enhance performance with minimal computational overhead through expert gating and routing mechanisms, hierarchical and sparse configurations, and meta-learning approaches. Their analysis identified key advantages including superior model capacity compared to equivalent Bayesian approaches and the importance of ensuring expert diversity, accurate calibration, and reliable inference aggregation. Csordás et al. [26] proposed SwitchHead, an effective MoE method for self-attention layers that computes up to 8 times fewer attention matrices than standard Transformers. For their 262M parameter model trained on C4, SwitchHead matched baseline perplexity with only 44% compute and 27% memory usage, achieving over 3.5% absolute improvements on BliMP zero-shot tasks compared to parameter-matched baselines. Ma et al. [28] developed Big-MoE for multimodal applications, introducing bypassing isolated gating mechanisms to enhance expert collaboration.

Dynamic routing strategies have gained attention for improving MoE efficiency through adaptive expert selection [27,42,43,44]. Yue et al. [27] proposed Ada-K routing, a novel strategy that dynamically adjusts the number of activated experts for each token using learnable allocator modules. Their method leverages Proximal Policy Optimization for end-to-end learning of this non-differentiable decision-making framework, achieving over 25% reduction in FLOPs and more than 20% inference speedup compared to conventional Top-K routing across four baseline models. Training efficiency remained high even for Mixtral-8x22B with 140B parameters, requiring only 8 hours. Their analysis revealed that harder tasks, middle layers, and content words tend to activate more experts, providing valuable insights for adaptive MoE system designs. Reinforcement learning-based routing has also been explored in network optimization contexts [43,44], demonstrating broader applicability of dynamic decision-making strategies. However, existing MoE frameworks lack application to time series forecasting scenarios where dual-domain processing and multi-horizon prediction require specialized expert routing strategies that can simultaneously optimize prediction accuracy, computational efficiency, and expert diversity across both frequency and temporal representations.

3. Method

In this section, we present a reinforcement learning-driven Mixture-of-Experts framework for long-term time series forecasting with dual-domain Transformer architecture. We first establish the forecasting problem formulation and introduce the dual-domain Transformer architecture that processes time series from both time and frequency domains simultaneously. Subsequently, we detail the MoE-enhanced feed-forward layers where traditional FFNs are replaced with specialized expert networks in both domains, with independent RL agents learning to dynamically route representations to the most appropriate experts. Finally, we describe the coupled training dynamics that jointly optimize the forecasting model and dual routing policies. Table 2 presents the mathematical notations used throughout this section.

Figure 1. Architecture of MoE-Transformer Framework with Dual-Domain Processing and RL-Based Dynamic Routing

Figure 1. Architecture of MoE-Transformer Framework with Dual-Domain Processing and RL-Based Dynamic Routing

3.1. Problem Formulation

Long-term time series forecasting can be formulated as a sequence-to-sequence prediction problem. Given a multivariate time series $\mathbf{X}=[{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_L]\in {\mathbb{R}}^{N\times L}$, where L denotes the look-back window size and N represents the number of variates, the objective is to predict future T steps $\widehat{\mathbf{X}}=[{\mathbf{x}}_{L+1},{\mathbf{x}}_{L+2},\dots ,{\mathbf{x}}_{L+T}]\in {\mathbb{R}}^{N\times T}$.

Following the channel-independence principle demonstrated effective by recent work, we process each variate independently to prevent mixing of distinct global features. For a univariate series $\mathbf{x}\in {\mathbb{R}}^L$, our framework learns a mapping function $f:{\mathbb{R}}^L\to {\mathbb{R}}^T$ that minimizes the expected forecasting error:

$$f^{*}=arg\underset{f\in \mathcal{F}}{min}{\mathbb{E}}_{(\mathbf{x},\mathbf{y})\sim \mathcal{D}}\left[\mathcal{L}\left(f\right(\mathbf{x}),\mathbf{y})\right],$$

(1)

where $\mathcal{F}$ denotes the hypothesis space of admissible functions, $\mathcal{D}$ represents the underlying data distribution, $\mathbf{y}\in {\mathbb{R}}^T$ is the ground truth future values, and $\mathcal{L}:{\mathbb{R}}^T\times {\mathbb{R}}^T\to \mathbb{R}$ is the loss function typically instantiated as Mean Squared Error.

3.2. Dual-Domain Transformer Architecture

Our framework employs a dual-domain Transformer architecture that processes time series from both time and frequency domains through parallel Transformer encoders. The complete forecasting pipeline is formalized as:

$$\begin{array}{cc}\hfill \mathbf{F}& =\mathbf{ExtDFT}\left(\mathbf{x}\right)\hfill \\ \hfill {\mathbf{x}}_o^f& =\mathrm{F-Transformer}\left(\mathbf{F}\right)\hfill \\ \hfill {\mathbf{x}}_o^t& =\mathrm{T-Transformer}\left(\mathbf{x}\right)\hfill \\ \hfill w_f,w_t& =\mathbf{Weight}\left(\mathbf{x}\right)\hfill \\ \hfill \widehat{\mathbf{x}}& =w_f{\mathbf{x}}_o^f+w_t{\mathbf{x}}_o^t,\hfill \end{array}$$

(2)

where F-Transformer and T-Transformer denote the frequency and time domain Transformer encoders respectively, each enhanced with domain-specific MoE in the feed-forward layers, $w_f,w_t\in [0,1]$ are adaptive ensemble weights satisfying $w_f+w_t=1$, ${\mathbf{x}}_o^f\in {\mathbb{R}}^T$ is the frequency domain prediction, and ${\mathbf{x}}_o^t\in {\mathbb{R}}^T$ is the time domain prediction.

Extended DFT for Frequency Alignment. The Discrete Fourier Transform (DFT) is a fundamental operation that converts a time series into its frequency domain representation. For an input series of length L, the standard DFT spectrum is calculated as:

$${\mathbf{F}}_0\left[k\right]=\sum _{n=0}^{L-1}\mathbf{x}\left[n\right]e^{-2\pi i{\displaystyle \frac{kn}{L}}},k=0,1,\dots ,L-1,$$

(3)

where ${\mathbf{F}}_0\in {\mathbb{C}}^L$ represents the standard DFT spectrum, $i=\sqrt{-1}$ denotes the imaginary unit, and $e^{-2\pi i{\displaystyle \frac{kn}{L}}}$ represents the complex exponential basis function at frequency index k.

However, the discrete frequency group determined by standard DFT is solely dependent on the sequence length, which results in a fundamental mismatch between the spectrum of the input series and the complete series (input plus forecast). To address this fundamental limitation, we employ Extended DFT which generates an input spectrum that aligns precisely with the DFT frequency group of the complete series:

$$\mathbf{F}\left[k\right]=\sum _{n=0}^{L-1}\mathbf{x}\left[n\right]e^{-2\pi i{\displaystyle \frac{kn}{L+T}}},k=0,1,\dots ,L+T-1,$$

(4)

where $\mathbf{F}\in {\mathbb{C}}^{L+T}$ represents the extended spectrum aligned with the complete series frequency grid. For real-valued time series, exploiting conjugate symmetry, we process only the first $\widehat{L}=\lfloor (L+T)/2\rfloor +1$ frequency components, effectively reducing computational complexity by half while preserving complete spectral information.

Time Domain Input Preparation. Following the PatchTST paradigm, the time domain input is divided into non-overlapping patches:

$${\mathbf{X}}_p=\mathrm{Reshape}\left(\mathbf{x}\right)\in {\mathbb{R}}^{N_p\times P},$$

(5)

where P denotes the patch length, $N_p=\lfloor L/P\rfloor$ denotes the number of patches, and Reshape$(·)$ segments the continuous time series into patch matrix ${\mathbf{X}}_p$.

3.3. Transformer Layer with Dual-Domain MoE

Both F-Transformer and T-Transformer consist of M stacked Transformer layers. Each layer follows the standard Transformer architecture but replaces the traditional feed-forward network with domain-specific MoE modules enhanced with RL-based routing. We describe the architecture for both domains in a unified manner, then detail the domain-specific components.

Frequency Domain Transformer Layer. For the m-th layer ($m\in \{1,\dots ,M\}$) in F-Transformer, the processing is formalized as:

$$\begin{array}{cc}\hfill {\tilde{\mathbf{h}}}_m^F& =\mathrm{LayerNorm}\left({\mathbf{h}}_{m-1}^F\right)\hfill \\ \hfill {\overline{\mathbf{h}}}_m^F& ={\mathrm{MHA}}^F({\tilde{\mathbf{h}}}_m^F)+{\mathbf{h}}_{m-1}^F\hfill \\ \hfill {\widehat{\mathbf{h}}}_m^F& =\mathrm{LayerNorm}({\overline{\mathbf{h}}}_m^F)\hfill \\ \hfill {\mathbf{h}}_m^F& ={\mathrm{MoE}}^F({\widehat{\mathbf{h}}}_m^F;{\pi }_{{\theta }^F}^m)+{\overline{\mathbf{h}}}_m^F,\hfill \end{array}$$

(6)

where ${\mathbf{h}}_0^F={\mathrm{Embed}}^F\left(\mathbf{F}\right)\in {\mathbb{C}}^{\widehat{L}\times d}$ denotes the initial frequency embeddings, ${\mathrm{MHA}}^F(·)$ denotes complex-valued multi-head self-attention over the spectrum, ${\mathrm{MoE}}^F(·;{\pi }_{{\theta }^F}^m)$ represents the frequency domain MoE module with RL routing policy ${\pi }_{{\theta }^F}^m$, and residual connections are applied after both attention and MoE modules.

The complex-valued multi-head attention is computed as:

$$\begin{array}{cc}\hfill {\mathbf{Q}}_h^m& ={\tilde{\mathbf{h}}}_m^F{\mathbf{W}}_{Q,h}^m,\hfill \end{array}$$

(7a)

$$\begin{array}{cc}\hfill {\mathbf{K}}_h^m& ={\tilde{\mathbf{h}}}_m^F{\mathbf{W}}_{K,h}^m,\hfill \end{array}$$

(7b)

$$\begin{array}{cc}\hfill {\mathbf{V}}_h^m& ={\tilde{\mathbf{h}}}_m^F{\mathbf{W}}_{V,h}^m,\hfill \end{array}$$

(7c)

$$\begin{array}{cc}\hfill {\mathrm{head}}_h^m& =\mathrm{Softmax}\left({\displaystyle \frac{\mathrm{Re}\left({\mathbf{Q}}_h^m{\left({\mathbf{K}}_h^m\right)}^{*}\right)}{\sqrt{d/H}}}\right){\mathbf{V}}_h^m,\hfill \end{array}$$

(7d)

$$\begin{array}{cc}\hfill {\mathrm{MHA}}^F(·)& =\mathrm{Concat}({\mathrm{head}}_1^m,\dots ,{\mathrm{head}}_H^m){\mathbf{W}}_O^m,\hfill \end{array}$$

(7e)

where ${\mathbf{W}}_{Q,h}^m,{\mathbf{W}}_{K,h}^m,{\mathbf{W}}_{V,h}^m\in {\mathbb{C}}^{d\times (d/H)}$ are projection matrices for head h, ${(·)}^{*}$ denotes complex conjugate transpose, Re$(·)$ extracts the real part for attention scores, and ${\mathbf{W}}_O^m\in {\mathbb{C}}^{d\times d}$ is the output projection.

Time Domain Transformer Layer. For the m-th layer in T-Transformer, the processing is formalized as:

$$\begin{array}{cc}\hfill {\tilde{\mathbf{h}}}_m^T& =\mathrm{LayerNorm}\left({\mathbf{h}}_{m-1}^T\right)\hfill \\ \hfill {\overline{\mathbf{h}}}_m^T& ={\mathrm{MHA}}^T\left({\tilde{\mathbf{h}}}_m^T\right)+{\mathbf{h}}_{m-1}^T\hfill \\ \hfill {\widehat{\mathbf{h}}}_m^T& =\mathrm{LayerNorm}\left({\overline{\mathbf{h}}}_m^T\right)\hfill \\ \hfill {\mathbf{h}}_m^T& ={\mathrm{MoE}}^T({\widehat{\mathbf{h}}}_m^T;{\pi }_{{\theta }^T}^m)+{\overline{\mathbf{h}}}_m^T,\hfill \end{array}$$

(8)

where ${\mathbf{h}}_0^T={\mathrm{Embed}}^T\left({\mathbf{X}}_p\right)\in {\mathbb{R}}^{N_p\times d}$ denotes the initial patch embeddings, ${\mathrm{MHA}}^T(·)$ denotes standard multi-head self-attention over patches, and ${\mathrm{MoE}}^T(·;{\pi }_{{\theta }^T}^m)$ represents the time domain MoE module with independent RL routing policy ${\pi }_{{\theta }^T}^m$.

The patch-based multi-head self-attention is computed as:

$$\begin{array}{cc}\hfill {\mathbf{Q}}_h^{T,m}& ={\tilde{\mathbf{h}}}_m^T{\mathbf{W}}_{Q,h}^{T,m},\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill {\mathbf{K}}_h^{T,m}& ={\tilde{\mathbf{h}}}_m^T{\mathbf{W}}_{K,h}^{T,m},\hfill \end{array}$$

(9b)

$$\begin{array}{cc}\hfill {\mathbf{V}}_h^{T,m}& ={\tilde{\mathbf{h}}}_m^T{\mathbf{W}}_{V,h}^{T,m},\hfill \end{array}$$

(9c)

$$\begin{array}{cc}\hfill {\mathrm{head}}_h^{T,m}& =\mathrm{Softmax}\left({\displaystyle \frac{{\mathbf{Q}}_h^{T,m}{\left({\mathbf{K}}_h^{T,m}\right)}^T}{\sqrt{d/H}}}\right){\mathbf{V}}_h^{T,m},\hfill \end{array}$$

(9d)

$$\begin{array}{cc}\hfill {\mathrm{MHA}}^T(·)& =\mathrm{Concat}({\mathrm{head}}_1^{T,m},\dots ,{\mathrm{head}}_H^{T,m}){\mathbf{W}}_O^{T,m},\hfill \end{array}$$

(9e)

where ${\mathbf{W}}_{Q,h}^{T,m},{\mathbf{W}}_{K,h}^{T,m},{\mathbf{W}}_{V,h}^{T,m}\in {\mathbb{R}}^{d\times (d/H)}$ are projection matrices and ${\mathbf{W}}_O^{T,m}\in {\mathbb{R}}^{d\times d}$ is the output projection.

Output Projection. After processing through all M layers, the final representations are projected to obtain domain-specific predictions:

$$\begin{array}{cc}\hfill {\mathbf{x}}_o^f& =\mathrm{Re}\left({\mathbf{W}}_p^f\mathrm{Flatten}\left({\mathbf{h}}_M^F\right)\right),\hfill \\ \hfill {\mathbf{x}}_o^t& ={\mathbf{W}}_p^t\mathrm{Flatten}\left({\mathbf{h}}_M^T\right),\hfill \end{array}$$

(10)

where ${\mathbf{W}}_p^f\in {\mathbb{C}}^{T\times (\widehat{L}·d)}$ and ${\mathbf{W}}_p^t\in {\mathbb{R}}^{T\times (N_p·d)}$ are projection matrices, Flatten$(·)$ concatenates all representations, and Re$(·)$ extracts real parts for complex-valued frequency predictions.

3.4. F-MoE: Frequency Domain Mixture-of-Experts

The F-MoE module replaces the traditional feed-forward network in frequency domain Transformer layers with a set of specialized expert networks ${\mathcal{E}}^F=\{E_1^F,E_2^F,\dots ,E_{K^F}^F\}$. Each expert $E_k^F$ is parameterized by ${\mathsf{\Theta }}_k^F$ and implements a two-layer architecture:

$$E_k^F(\mathbf{h};{\mathsf{\Theta }}_k^F)={\mathbf{W}}_{2,k}^F\mathbf{\sigma }({\mathbf{W}}_{1,k}^F\mathbf{h}+{\mathbf{b}}_{1,k}^F)+{\mathbf{b}}_{2,k}^F,$$

(11)

where $\mathbf{h}\in {\mathbb{C}}^d$ is the input hidden state, ${\mathbf{W}}_{1,k}^F\in {\mathbb{C}}^{d_f\times d}$ and ${\mathbf{W}}_{2,k}^F\in {\mathbb{C}}^{d\times d_f}$ are weight matrices, ${\mathbf{b}}_{1,k}^F\in {\mathbb{C}}^{d_f}$ and ${\mathbf{b}}_{2,k}^F\in {\mathbb{C}}^d$ are bias vectors, and $\sigma (·)$ represents the activation function applied to complex numbers as $\sigma \left(\mathbf{z}\right)=\sigma \left(\mathrm{Re}\right(\mathbf{z}\left)\right)+i\sigma \left(\mathrm{Im}\right(\mathbf{z}\left)\right)$.

At layer m and iteration t, the RL agent selects expert $a_t^{F,m}\in \{1,\dots ,K^F\}$ based on the current state, and the F-MoE output is:

$${\mathrm{MoE}}^F({\widehat{\mathbf{h}}}_m^F;{\pi }_{{\theta }^F}^m)=E_{a_t^{F,m}}^F\left({\widehat{\mathbf{h}}}_m^F\right),$$

(12)

where the action $a_t^{F,m}$ is sampled from the policy distribution ${\pi }_{{\theta }^F}^m\left(a\right|s_t^{F,m})$.

3.5. T-MoE: Time Domain Mixture-of-Experts

The T-MoE module replaces the traditional feed-forward network in time domain Transformer layers with specialized expert networks ${\mathcal{E}}^T=\{E_1^T,E_2^T,\dots ,E_{K^T}^T\}$. Each expert $E_k^T$ is parameterized by ${\mathsf{\Theta }}_k^T$ and implements a two-layer architecture:

$$E_k^T(\mathbf{h};{\mathsf{\Theta }}_k^T)={\mathbf{W}}_{2,k}^T\mathbf{\sigma }({\mathbf{W}}_{1,k}^T\mathbf{h}+{\mathbf{b}}_{1,k}^T)+{\mathbf{b}}_{2,k}^T,$$

(13)

where $\mathbf{h}\in {\mathbb{R}}^d$ is the input hidden state, ${\mathbf{W}}_{1,k}^T\in {\mathbb{R}}^{d_f\times d}$ and ${\mathbf{W}}_{2,k}^T\in {\mathbb{R}}^{d\times d_f}$ are weight matrices, ${\mathbf{b}}_{1,k}^T\in {\mathbb{R}}^{d_f}$ and ${\mathbf{b}}_{2,k}^T\in {\mathbb{R}}^d$ are bias vectors, and $\sigma (·)$ represents the GELU activation function.

At layer m and iteration t, the independent RL agent selects expert $a_t^{T,m}\in \{1,\dots ,K^T\}$, and the T-MoE output is:

$${\mathrm{MoE}}^T({\widehat{\mathbf{h}}}_m^T;{\pi }_{{\theta }^T}^m)=E_{a_t^{T,m}}^T\left({\widehat{\mathbf{h}}}_m^T\right),$$

(14)

where the action $a_t^{T,m}$ is sampled from the independent policy distribution ${\pi }_{{\theta }^T}^m\left(a\right|s_t^{T,m})$.

3.6. Reinforcement Learning Framework for F-MoE Routing

We formulate the F-MoE expert routing problem as a Markov Decision Process (MDP) defined by the tuple $({\mathcal{S}}^F,{\mathcal{A}}^F,{\mathcal{P}}^F,{\mathcal{R}}^F,\gamma )$. Below, we detail each component of the MDP formulation for frequency domain routing.

State Space ${\mathcal{S}}^F$. The state space ${\mathcal{S}}^F$ encodes all relevant information for routing decisions in the frequency domain at layer m and iteration t:

$$s_t^{F,m}=[\mathrm{Re}\left({\widehat{\mathbf{h}}}_m^F\right);\mathrm{Im}\left({\widehat{\mathbf{h}}}_m^F\right);{\mathbf{c}}_t^{F,m};a_{t-1}^{F,m};{\mathbf{p}}_t^{F,m}],$$

(15)

where Re$(·)$ and Im$(·)$ extract real and imaginary parts of the hidden state, ${\mathbf{c}}_t^{F,m}\in {\mathbb{R}}^{d_c}$ represents extracted frequency features (dominant frequency, spectral energy distribution, phase characteristics), $a_{t-1}^{F,m}\in \{1,\dots ,K^F\}$ records the previous expert selection, and ${\mathbf{p}}_t^{F,m}\in {[0,1]}^{K^F}$ denotes empirical expert usage distribution.

Action Space ${\mathcal{A}}^F$. The action space ${\mathcal{A}}^F=\{1,2,\dots ,K^F\}$ corresponds to selecting one expert from the F-expert set. At each decision point, the agent selects exactly one expert to process the current spectral representation.

State Transition ${\mathcal{P}}^F$. The state transition function ${\mathcal{P}}^F:{\mathcal{S}}^F\times {\mathcal{A}}^F\to \Delta \left({\mathcal{S}}^F\right)$ is deterministic and governed by the forward pass. Given current state $s_t^{F,m}$ and action $a_t^{F,m}=k$, the next state is:

$$s_{t+1}^{F,m}=f_{\varphi }^F(s_t^{F,m},a_t^{F,m},E_k^F\left({\widehat{\mathbf{h}}}_m^F\right)),$$

(16)

where $f_{\varphi }^F$ updates the state based on the selected expert’s output.

Reward Function ${\mathcal{R}}^F$. The reward function ${\mathcal{R}}^F:{\mathcal{S}}^F\times {\mathcal{A}}^F\times {\mathcal{S}}^F\to \mathbb{R}$ balances multiple objectives:

$$r_t^{F,m}=-{\alpha }_1^F·{\mathcal{L}}_m({\widehat{\mathbf{x}}}_t,{\mathbf{x}}_t)-{\alpha }_2^F·C_s^{F,m}+{\alpha }_3^F·D_e^{F,m},$$

(17)

where ${\mathcal{L}}_m({\widehat{\mathbf{x}}}_t,{\mathbf{x}}_t)={\displaystyle \frac{1}{T}}{\sum }_{i=1}^T{|{\widehat{x}}_{t,i}-x_{t,i}|}^2$ measures mean squared error, $C_s^{F,m}={\mathbb{1}}_{\{a_t^{F,m}\ne a_{t-1}^{F,m}\}}$ penalizes expert switching, $D_e^{F,m}=-{\sum }_{k=1}^{K^F}{p}_{k,t}^{F,m}log{p}_{k,t}^{F,m}$ promotes balanced utilization where $p_{k,t}^{F,m}={\displaystyle \frac{1}{t}}{\sum }_{\tau =1}^t{\mathbb{1}}_{\{a_{\tau }^{F,m}=k\}}$ is the empirical usage frequency, and ${\alpha }_1^F,{\alpha }_2^F,{\alpha }_3^F\in {\mathbb{R}}^{+}$ are weighting coefficients.

Policy Network and Value Function. The policy ${\pi }_{{\theta }^F}:{\mathcal{S}}^F\to \Delta \left({\mathcal{A}}^F\right)$ is parameterized by a multi-layer perceptron:

$$f_{{\theta }^F}\left(s\right)={\mathbf{U}}_{\pi }^F\sigma ({\mathbf{V}}_{\pi }^F\sigma ({\mathbf{W}}_{\pi }^{F}s+{\mathbf{b}}_{\pi }^F)+{\mathbf{c}}_{\pi }^F)+{\mathbf{d}}_{\pi }^F,$$

(18)

where $s\in {\mathbb{R}}^{|{\mathcal{S}}^F|}$ is the state vector, ${\mathbf{W}}_{\pi }^F\in {\mathbb{R}}^{d_h\times \left|{\mathcal{S}}^F\right|}$, ${\mathbf{V}}_{\pi }^F\in {\mathbb{R}}^{d_h\times d_h}$, ${\mathbf{U}}_{\pi }^F\in {\mathbb{R}}^{K^F\times d_h}$ are weight matrices, and bias terms complete the network. The policy distribution is:

$${\pi }_{{\theta }^F}\left(a\right|s)={\displaystyle \frac{exp\left(f_{{\theta }^F}{\left(s\right)}_a\right)}{{\sum }_{k=1}^{K^F}exp\left(f_{{\theta }^F}{\left(s\right)}_k\right)}}.$$

(19)

The value function $V_{{\varphi }^F}:{\mathcal{S}}^F\to \mathbb{R}$ estimates expected returns:

$$V_{{\varphi }^F}\left(s\right)={\left({\mathbf{w}}_v^F\right)}^T\sigma ({\mathbf{W}}_v^{F}s+{\mathbf{b}}_v^F),$$

(20)

where ${\mathbf{W}}_v^F\in {\mathbb{R}}^{d_v\times \left|{\mathcal{S}}^F\right|}$ is the weight matrix and ${\mathbf{w}}_v^F\in {\mathbb{R}}^{d_v}$ is the output vector.

The advantage function is computed using Generalized Advantage Estimation:

$$A_t^{F,m}=\sum _{\ell =0}^{\infty }{\left(\gamma \lambda \right)}^{\ell }{\delta }_{t+\ell }^{F,m},$$

(21)

where ${\delta }_t^{F,m}=r_t^{F,m}+\gamma V_{{\varphi }^F}\left(s_{t+1}^{F,m}\right)-V_{{\varphi }^F}\left(s_t^{F,m}\right)$ is the temporal difference error, $\gamma \in [0,1)$ is the discount factor, and $\lambda \in [0,1]$ is the GAE parameter.

3.7. Reinforcement Learning Framework for T-MoE Routing

We formulate the T-MoE expert routing problem as an independent Markov Decision Process (MDP) defined by the tuple $({\mathcal{S}}^T,{\mathcal{A}}^T,{\mathcal{P}}^T,{\mathcal{R}}^T,\gamma )$. Below, we detail each component for time domain routing.

State Space ${\mathcal{S}}^T$. The state space ${\mathcal{S}}^T$ encodes routing information in the time domain at layer m and iteration t:

$$s_t^{T,m}=[{\widehat{\mathbf{h}}}_m^T;{\mathbf{c}}_t^{T,m};a_{t-1}^{T,m};{\mathbf{p}}_t^{T,m}],$$

(22)

where ${\widehat{\mathbf{h}}}_m^T$ is the hidden state, ${\mathbf{c}}_t^{T,m}\in {\mathbb{R}}^{d_c}$ represents temporal features (autocorrelation, trend, local variability), $a_{t-1}^{T,m}\in \{1,\dots ,K^T\}$ records previous expert selection, and ${\mathbf{p}}_t^{T,m}\in {[0,1]}^{K^T}$ denotes expert usage distribution.

Action Space ${\mathcal{A}}^T$. The action space ${\mathcal{A}}^T=\{1,2,\dots ,K^T\}$ corresponds to selecting one expert from the T-expert set.

State Transition ${\mathcal{P}}^T$. The state transition function ${\mathcal{P}}^T:{\mathcal{S}}^T\times {\mathcal{A}}^T\to \Delta \left({\mathcal{S}}^T\right)$ is deterministic:

$$s_{t+1}^{T,m}=f_{\varphi }^T(s_t^{T,m},a_t^{T,m},E_k^T\left({\widehat{\mathbf{h}}}_m^T\right)),$$

(23)

where $f_{\varphi }^T$ updates the state based on the selected expert’s output.

Reward Function ${\mathcal{R}}^T$. The reward function ${\mathcal{R}}^T:{\mathcal{S}}^T\times {\mathcal{A}}^T\times {\mathcal{S}}^T\to \mathbb{R}$ balances time domain objectives:

$$r_t^{T,m}=-{\alpha }_1^T·{\mathcal{L}}_m({\widehat{\mathbf{x}}}_t,{\mathbf{x}}_t)-{\alpha }_2^T·C_s^{T,m}+{\alpha }_3^T·D_e^{T,m},$$

(24)

where ${\mathcal{L}}_m$ measures prediction error, $C_s^{T,m}={\mathbb{1}}_{\{a_t^{T,m}\ne a_{t-1}^{T,m}\}}$ penalizes switching, $D_e^{T,m}=-{\sum }_{k=1}^{K^T}{p}_{k,t}^{T,m}log{p}_{k,t}^{T,m}$ promotes balanced utilization with $p_{k,t}^{T,m}={\displaystyle \frac{1}{t}}{\sum }_{\tau =1}^t{\mathbb{1}}_{\{a_{\tau }^{T,m}=k\}}$, and ${\alpha }_1^T,{\alpha }_2^T,{\alpha }_3^T\in {\mathbb{R}}^{+}$ are weighting coefficients.

Policy Network and Value Function. The policy ${\pi }_{{\theta }^T}:{\mathcal{S}}^T\to \Delta \left({\mathcal{A}}^T\right)$ is parameterized independently:

$$f_{{\theta }^T}\left(s\right)={\mathbf{U}}_{\pi }^T\sigma ({\mathbf{V}}_{\pi }^T\sigma ({\mathbf{W}}_{\pi }^{T}s+{\mathbf{b}}_{\pi }^T)+{\mathbf{c}}_{\pi }^T)+{\mathbf{d}}_{\pi }^T,$$

(25)

where $s\in {\mathbb{R}}^{|{\mathcal{S}}^T|}$, ${\mathbf{W}}_{\pi }^T\in {\mathbb{R}}^{d_h\times \left|{\mathcal{S}}^T\right|}$, ${\mathbf{V}}_{\pi }^T\in {\mathbb{R}}^{d_h\times d_h}$, ${\mathbf{U}}_{\pi }^T\in {\mathbb{R}}^{K^T\times d_h}$ are weight matrices. The policy distribution is:

$${\pi }_{{\theta }^T}\left(a\right|s)={\displaystyle \frac{exp\left(f_{{\theta }^T}{\left(s\right)}_a\right)}{{\sum }_{k=1}^{K^T}exp\left(f_{{\theta }^T}{\left(s\right)}_k\right)}}.$$

(26)

The value function $V_{{\varphi }^T}:{\mathcal{S}}^T\to \mathbb{R}$ is:

$$V_{{\varphi }^T}\left(s\right)={\left({\mathbf{w}}_v^T\right)}^T\sigma ({\mathbf{W}}_v^{T}s+{\mathbf{b}}_v^T),$$

(27)

where ${\mathbf{W}}_v^T\in {\mathbb{R}}^{d_v\times \left|{\mathcal{S}}^T\right|}$ and ${\mathbf{w}}_v^T\in {\mathbb{R}}^{d_v}$.

The advantage function is:

$$A_t^{T,m}=\sum _{\ell =0}^{\infty }{\left(\gamma \lambda \right)}^{\ell }{\delta }_{t+\ell }^{T,m},$$

(28)

where ${\delta }_t^{T,m}=r_t^{T,m}+\gamma V_{{\varphi }^T}\left(s_{t+1}^{T,m}\right)-V_{{\varphi }^T}\left(s_t^{T,m}\right)$.

3.8. Ensemble Weighting

The final prediction ensembles the frequency and time domain outputs using Dominant Harmonic Series Energy Weighting:

$$\widehat{\mathbf{x}}=w_f{\mathbf{x}}_o^f+w_t{\mathbf{x}}_o^t,$$

(29)

where the frequency weight is determined by:

$$w_f={\displaystyle \frac{E_h}{E_f}}={\displaystyle \frac{{\sum }_{n=0}^{\tau -1}{\left|\mathbf{F}\left[nk\right]\right|}^2}{{\sum }_{n=0}^{\widehat{L}-1}{\left|\mathbf{F}\left[n\right]\right|}^2}},$$

(30)

where $k=arg{max}_{j\le \widehat{L}}\left|\mathbf{F}\left[j\right]\right|$ identifies the fundamental frequency, $\tau =\lfloor \widehat{L}/k\rfloor$ determines the number of harmonics, $E_h$ is the energy in dominant harmonics, and $E_f$ is total spectral energy. The time weight is $w_t=1-w_f$. This mechanism allocates more weight to the frequency domain for periodic series and to the time domain for non-periodic series.

3.9. Coupled Training Dynamics

The framework is trained end-to-end by jointly optimizing forecasting model parameters and dual RL policies. The expert networks ${\left\{{\mathsf{\Theta }}_k^F\right\}}_{k=1}^{K^F}$ and ${\left\{{\mathsf{\Theta }}_k^T\right\}}_{k=1}^{K^T}$, along with all Transformer weights, are updated to minimize prediction error:

$$\mathsf{\Theta }\leftarrow \mathsf{\Theta }-{\mathbf{\eta }}_m{\nabla }_{\mathsf{\Theta }}{\mathbb{E}}_{(\mathbf{x},\mathbf{y})\sim \mathcal{D}}\left[{\mathcal{L}}_m(\widehat{\mathbf{x}},\mathbf{y})\right],$$

(31)

where $\mathsf{\Theta }$ denotes all model parameters, ${\eta }_m$ is the learning rate, and the expectation is approximated through mini-batch sampling.

Simultaneously, the RL policies ${\left\{{\pi }_{{\theta }^{F,m}}\right\}}_{m=1}^M$ and ${\left\{{\pi }_{{\theta }^{T,m}}\right\}}_{m=1}^M$ are optimized to maximize expected cumulative rewards:

$$\underset{{\{{\theta }^{F,m},{\theta }^{T,m}\}}_{m=1}^M}{max}{\mathbb{E}}_{\pi }\left[\sum _{t=0}^{\infty }{\gamma }^t\left(\sum _{m=1}^M{r}_t^{F,m}+\sum _{m=1}^M{r}_t^{T,m}\right)\right],$$

(32)

where the expectation is over trajectories generated by current policies.

The two optimization processes are coupled through the reward signal, which incorporates forecasting performance via the MSE term in Eq. 17 and Eq. 24. This ensures RL agents in both domains learn to select experts that minimize prediction error while maintaining routing stability and expert diversity. The gradients from forecasting loss flow through selected experts to update their parameters, while policy gradients guide routing decisions. The dual-domain RL framework enables each domain to develop specialized routing strategies: F-MoE agents learn to route based on spectral patterns, while T-MoE agents route based on temporal dependencies.

Connection to Other Routing Mechanisms. We compare our dual RL-based routing with other MoE routing practices: gating networks, top-k routing, and noisy top-k gating. For static methods, our RL framework provides: (1) dynamic adaptation of routing decisions based on domain-specific characteristics through learned reward signals, while static gating learns fixed patterns; (2) simultaneous optimization of multiple objectives including accuracy, stability, and diversity in each domain, while traditional methods optimize only for task performance. Compared to noisy top-k gating, we explicitly model sequential decision-making and optimize for long-term performance through temporal difference learning, rather than relying on random noise. Furthermore, our dual-domain design enables specialized routing strategies where F-MoE and T-MoE agents independently learn domain-specific routing, addressing the challenge that different domains require distinct routing strategies in time series forecasting.

4. Experiments

4.1. Datasets and Experimental Setup

We conduct comprehensive experiments to validate the effectiveness of our RL-driven dual-domain MoE framework across diverse long-term time series forecasting scenarios. Our experimental evaluation systematically examines each component of the proposed framework: dual-domain Transformer architecture with Extended DFT, F-MoE and T-MoE modules with independent RL routing, and adaptive ensemble weighting, demonstrating their synergistic contribution to accurate long-term forecasting.

To ensure the robustness and statistical significance of our results, all experiments are conducted with 5 independent runs using different random initializations. We report the mean and standard deviation for all performance metrics, and conduct paired t-tests to verify that improvements over baseline methods are statistically significant. This rigorous experimental protocol ensures that our observed performance gains are not artifacts of random initialization or data splitting, but rather represent genuine improvements in predictive capability.

4.1.1. Dataset Description

The ETTh1 (Electricity Transformer Temperature - hourly) dataset [45] contains measurements from electricity transformers, recorded at hourly intervals over 2 years with $N=7$ variates including oil temperature and load features. The dataset exhibits clear daily and weekly periodicities with strong harmonic structures, making it ideal for evaluating our frequency domain modeling with Extended DFT alignment. The 17,420 hourly observations provide rich temporal patterns for testing F-MoE’s spectral expert specialization capabilities on periodic data.

The ETTm1 (Electricity Transformer Temperature - 15min) dataset [45] provides finer temporal granularity with 15-minute sampling intervals spanning 2 years, containing $N=7$ variates identical to ETTh1 but capturing intra-hour dynamics. The 69,680 observations reveal multi-scale periodicities including hourly, daily, and weekly cycles. This high-frequency sampling challenges our dual-domain architecture to handle both rapid fluctuations and long-term dependencies, testing the robustness of RL-based expert routing across different temporal resolutions.

The Weather dataset [46] comprises 21 meteorological indicators including air temperature, humidity, and pressure, collected every 10 minutes over multiple years. The dataset spans 52,696 time steps with $N=21$ variates, exhibiting complex seasonal patterns, trend components, and non-stationary characteristics. The long-term dependencies and irregular weather events provide an ideal testbed for evaluating T-MoE’s temporal pattern recognition and the adaptive ensemble weighting mechanism’s ability to balance domain contributions for trend-dominated series.

The Electricity dataset [46] records hourly electricity consumption of 321 clients over 3 years, totaling 26,304 time steps. With $N=321$ variates, this high-dimensional dataset presents significant challenges for capturing complex consumption patterns across diverse user profiles. The dataset exhibits strong weekly periodicities combined with irregular consumption behaviors, enabling comprehensive evaluation of our channel-independent processing strategy and F-MoE’s ability to identify dominant harmonic structures through Extended DFT for per-channel frequency analysis.

The Traffic dataset [46] contains hourly road occupancy rates from 862 sensors on San Francisco Bay Area freeways, spanning 2 years with 17,544 observations. The dataset exhibits pronounced rush-hour periodicities, weekly patterns, and seasonal variations across different road segments. The large number of variates ($N=862$) with heterogeneous temporal dynamics provides opportunities to evaluate domain-specific expert specialization, where F-MoE agents learn to route periodic traffic patterns while T-MoE agents handle irregular congestion events and abrupt traffic changes.

4.1.2. Problem Formulation Alignment

Following the mathematical formulation in Section 3.1, we configure each dataset according to our forecasting framework. For each univariate series $\mathbf{x}\in {\mathbb{R}}^L$ extracted from the multivariate input $\mathbf{X}\in {\mathbb{R}}^{N\times L}$, we evaluate multiple forecasting configurations with look-back window $L=96$ and prediction horizons $T\in \{96,192,336,720\}$ time steps.

The dataset-specific configurations align with our mapping function $f:{\mathbb{R}}^L\to {\mathbb{R}}^T$ as formulated in Equation (1). ETTh1 processes $N=7$ variates independently with $f:{\mathbb{R}}^{96}\to {\mathbb{R}}^T$, ETTm1 with $N=7$ variates, Weather with $N=21$ variates, Electricity with $N=321$ variates, and Traffic with $N=862$ variates. Each variate is processed independently following the channel-independence principle to prevent mixing of distinct global features.

Data preprocessing applies Z-score normalization $\tilde{x}=(x-\mu )/\sigma$ with statistics computed from the training set. We employ a strict chronological split of 60%-20%-20% for training, validation, and testing respectively. Specifically, the first 60% of the temporal sequence is used for training, the immediately following 20% for validation (hyperparameter tuning and early stopping), and the final 20% for testing. This contiguous chronological split preserves temporal integrity and ensures valid evaluation of predictive performance on future unseen data.

To demonstrate practical feasibility for real-world deployment, we report computational efficiency metrics. Our framework contains 2.14M trainable parameters including dual-domain Transformers, $K^F=4$ frequency experts, $K^T=4$ time experts, and dual RL policy networks across $M=3$ layers. All experiments are conducted on an NVIDIA RTX 4090 GPU with PyTorch 2.1.0. Training times to convergence are 1.8-2.2 hours for ETTh1, 6.5-8.1 hours for ETTm1, 2.1-2.6 hours for Weather, 4.8-6.2 hours for Electricity, and 7.3-9.5 hours for Traffic. Inference latency averages 8.5 milliseconds per prediction step on ETTh1, with batch processing achieving 1.2 milliseconds per sample. These metrics confirm that our framework maintains practical computational requirements suitable for real-time forecasting systems.

4.2. Model Configuration

4.2.1. Component-Specific Parameters

Table 3 details the hyperparameter configuration for each component, directly corresponding to the mathematical formulations in Section 3.

4.2.2. Coupled Training Configuration

Our training objective employs the coupled optimization strategy defined in Section 3.9, jointly optimizing forecasting model parameters $\mathsf{\Theta }$ and dual RL policies through Equations (31) and (32). The forecasting model parameters including all expert networks ${\left\{{\mathsf{\Theta }}_k^F\right\}}_{k=1}^{K^F}$ and ${\left\{{\mathsf{\Theta }}_k^T\right\}}_{k=1}^{K^T}$ are optimized using Adam with learning rate ${\eta }_m=0.0001$. The F-domain RL policies ${\left\{{\pi }_{{\theta }^{F,m}}\right\}}_{m=1}^M$ and T-domain RL policies ${\left\{{\pi }_{{\theta }^{T,m}}\right\}}_{m=1}^M$ are independently optimized using separate Adam optimizers with learning rate 0.0001.

The reward balancing coefficients are configured as ${\alpha }_1^F={\alpha }_1^T=1.0$ (prediction accuracy from Equation (17) and (24)), ${\alpha }_2^F={\alpha }_2^T=0.1$ (expert switching penalty), and ${\alpha }_3^F={\alpha }_3^T=0.05$ (expert diversity bonus), determined through validation set grid search. Early stopping with patience 10 monitors validation MSE. Gradient clipping with maximum norm 1.0 stabilizes RL policy updates. The ensemble weighting mechanism follows Equation (29) with frequency weights $w_f$ computed adaptively via Equation (30) based on dominant harmonic energy ratios.

4.3. Evaluation Metrics and Analysis Framework

4.3.1. Primary Evaluation Metrics

We employ four complementary metrics to comprehensively assess forecasting accuracy across different error characteristics, consistent with standard practice in long-term time series forecasting literature.

The Mean Squared Error (MSE) is computed as:

$$\mathrm{MSE}={\displaystyle \frac{1}{|{\mathcal{D}}_{\mathrm{test}}|·T}}\sum _{(\mathbf{x},\mathbf{y})\in {\mathcal{D}}_{\mathrm{test}}}\sum _{t=1}^T{({\widehat{\mathbf{x}}}_t-{\mathbf{y}}_t)}^2,$$

(33)

where ${\mathcal{D}}_{\mathrm{test}}$ denotes the test set, $\widehat{\mathbf{x}}=[{\widehat{\mathbf{x}}}_1,\dots ,{\widehat{\mathbf{x}}}_T]\in {\mathbb{R}}^T$ represents the predicted sequence from Equation (2), $\mathbf{y}=[{\mathbf{y}}_1,\dots ,{\mathbf{y}}_T]\in {\mathbb{R}}^T$ is the ground truth future values as defined in Equation (1), and T denotes the forecast horizon. MSE provides sensitivity to large prediction errors through squared differences.

The Mean Absolute Error (MAE) provides robustness against outliers:

$$\mathrm{MAE}={\displaystyle \frac{1}{|{\mathcal{D}}_{\mathrm{test}}|·T}}\sum _{(\mathbf{x},\mathbf{y})\in {\mathcal{D}}_{\mathrm{test}}}\sum _{t=1}^T|{\widehat{\mathbf{x}}}_t-{\mathbf{y}}_t|,$$

(34)

where all notation follows from the MSE definition above. MAE offers interpretability through absolute deviations and reduced sensitivity to extreme values compared to MSE.

The Root Mean Squared Error (RMSE) provides sensitivity to large prediction errors:

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{|{\mathcal{D}}_{\mathrm{test}}|·T}}{\sum }_{(\mathbf{x},\mathbf{y})\in {\mathcal{D}}_{\mathrm{test}}}{\sum }_{t=1}^T{({\widehat{\mathbf{x}}}_t-{\mathbf{y}}_t)}^2},$$

(35)

where the square root operation transforms MSE back to the original scale of the data, enabling direct comparison with MAE while maintaining sensitivity to large errors.

The Mean Absolute Percentage Error (MAPE) offers relative error assessment:

$$\mathrm{MAPE}={\displaystyle \frac{100\%}{|{\mathcal{D}}_{\mathrm{test}}|·T}}\sum _{(\mathbf{x},\mathbf{y})\in {\mathcal{D}}_{\mathrm{test}}}\sum _{t=1}^T{\displaystyle \frac{|{\widehat{\mathbf{x}}}_t-{\mathbf{y}}_t|}{|{\mathbf{y}}_t|}},$$

(36)

where the percentage formulation normalizes errors by ground truth magnitudes, enabling scale-independent comparison across datasets with different value ranges.

5. Experimental Results

5.1. Reinforcement Learning Training Dynamics

Figure 2 illustrates the training convergence of dual RL routing policies with varying numbers of experts. We evaluate the impact of expert capacity by testing configurations with N=10, 20, 30, and 40 experts, comparing four RL algorithms: IPPO (Independent PPO), ISAC (Independent SAC), TD3, and MAPPO (Multi-Agent PPO). This analysis validates our algorithm selection and demonstrates the effectiveness of our RL-based routing mechanism across different model capacities.

As shown in Figure 2(a-d), the F-MoE agent exhibits distinct convergence patterns across different expert capacities and RL algorithms. For N=10 experts (Figure 2(a)), ISAC demonstrates exceptional performance with final rewards reaching approximately 1025, significantly outperforming other algorithms. MAPPO achieves approximately 1000, while IPPO and TD3 converge to around 925 and 900 respectively. As expert capacity increases to N=20 (Figure 2(b)), ISAC maintains its dominance with rewards around 1000, though the performance gap narrows. MAPPO shows improved convergence speed and achieves approximately 950, while IPPO and TD3 reach similar final rewards around 920-930.

For higher expert capacities N=30 and N=40 (Figure 2(c-d)), the convergence patterns become more uniform across algorithms. ISAC consistently maintains the highest final rewards (approximately 975-1000), validating its sample efficiency for frequency domain routing. Notably, MAPPO’s performance declines slightly with increased expert capacity, converging to approximately 850-900 for N=40, suggesting that centralized critics may face increased coordination complexity with larger expert sets. IPPO and TD3 demonstrate stable performance across all capacities, consistently achieving rewards around 880-920, indicating that independent learning and deterministic policies can effectively handle varying expert ensemble sizes in the frequency domain.

Figure 2(e-h) reveals complementary patterns for the T-MoE agent with notable domain-specific characteristics. For N=10 experts (Figure 2(e)), ISAC again demonstrates superior performance with final rewards around 1020, followed by IPPO and TD3 (approximately 960-970), while MAPPO achieves approximately 920. This pattern suggests that T-domain routing benefits significantly from off-policy learning due to more stable temporal patterns compared to frequency domain. As expert capacity increases to N=20 (Figure 2(f)), the algorithms converge to similar performance levels (940-950 for IPPO/TD3, 920-940 for MAPPO), with ISAC maintaining its advantage at approximately 1000.

For N=30 and N=40 experts (Figure 2(g-h)), T-MoE exhibits more balanced performance across algorithms compared to F-MoE. IPPO and TD3 show particularly stable convergence, reaching final rewards of 930-970 with low variance, demonstrating robustness to expert capacity scaling. MAPPO’s performance improves slightly with increased capacity, converging to approximately 920-970, contrasting with its F-domain behavior and suggesting that temporal routing benefits from centralized coordination when expert sets are larger. ISAC maintains consistent leadership across all capacities (990-1020), confirming the effectiveness of off-policy learning for temporal pattern recognition.

The comparative analysis across expert capacities reveals several critical insights: (1) ISAC consistently achieves the highest final rewards across both domains and all expert capacities, validating off-policy learning superiority for MoE routing in time series forecasting, (2) F-MoE shows greater sensitivity to expert capacity, with performance divergence increasing from N=10 to N=40, while T-MoE maintains more stable relative performance across algorithms, (3) The multi-objective reward function (Equations (17) and (24)) successfully guides all algorithms toward improved routing strategies across different capacities, with convergence typically stabilizing after 400-500 episodes, (4) Variance analysis (shaded regions) indicates that ISAC and TD3 maintain lower training variance compared to PPO-based methods, particularly for larger expert sets (N=30, 40), suggesting better robustness to hyperparameter sensitivity and random initialization.

The domain-specific patterns validate our dual RL framework design: F-MoE requires more sophisticated coordination mechanisms that benefit from ISAC’s sample efficiency, while T-MoE’s more balanced performance across algorithms confirms that temporal routing is inherently more stable. These findings justify our selection of ISAC as the primary RL algorithm for both F-MoE and T-MoE routing in our final framework, while also demonstrating that the framework maintains effectiveness across a wide range of expert capacities (N=10-40).

5.2. Overall Performance Comparison

We present comprehensive performance comparisons between our RL-driven dual-domain MoE framework and state-of-the-art baseline methods across five benchmark datasets and four forecasting horizons. To evaluate the impact of expert capacity on forecasting performance, we conduct experiments with varying numbers of experts: N=10, 20, 30, and 40. Table 4, Table 5, Table 6 and Table 7 demonstrate the effectiveness of synergistic integration of dual-domain processing, Extended DFT alignment, RL-based expert routing, and adaptive ensemble weighting across different model capacities.

The comprehensive evaluation across varying expert capacities (N=10, 20, 30, 40) reveals several critical insights into the scalability and effectiveness of our RL-driven dual-domain MoE framework.

Performance scaling with expert capacity. As shown in Table 4, Table 5, Table 6 and Table 7, our framework demonstrates consistent improvements as expert capacity increases from N=10 to N=40. On ETTh1 at horizon 96, performance improves from 0.298±0.013 (N=10) to 0.278±0.012 (N=20), 0.271±0.012 (N=30), and finally 0.267±0.011 (N=40), representing incremental gains of 6.7%, 2.5%, and 1.5% respectively. This pattern of diminishing returns is consistent across all datasets, validating that while larger expert ensembles provide better representation capacity, the marginal benefits decrease beyond N=30.

Improvement over baselines. The relative improvement over PatchTST (the strongest baseline) increases systematically with expert capacity. For N=10 experts, we achieve 49.4-51.7% MSE reduction across datasets at horizon 96, which improves to 52.8-54.8% for N=20, 54.0-55.9% for N=30, and peaks at 54.7-56.9% for N=40. This demonstrates that our RL-based routing mechanism becomes increasingly effective at leveraging larger expert pools, successfully avoiding the expert collapse and load imbalance issues that plague static gating mechanisms with many experts.

Domain-specific capacity effects. The periodic datasets (ETTh1, ETTm1, Electricity, Traffic) show particularly strong scaling behavior, with N=40 achieving 54.7%, 55.5%, 56.4%, and 56.9% improvements respectively. This suggests that F-MoE’s spectral expert specialization benefits substantially from larger expert capacity, as different experts can specialize in finer-grained harmonic structures. The Weather dataset shows slightly less dramatic scaling (49.9% to 55.7% from N=10 to N=40), reflecting that trend-dominated series with weaker periodicity derive less benefit from expanded frequency domain expert capacity.

Training stability across capacities. A critical observation is the consistent reduction in standard deviation as expert capacity increases. For ETTh1 horizon 96, standard deviations decrease from ±0.013 (N=10) to ±0.012 (N=20, N=30) and ±0.011 (N=40), demonstrating that our RL routing framework maintains training stability even with larger action spaces. This stability, combined with the performance improvements, confirms that our independent ISAC-based routing policies (as validated in Figure 2) effectively explore and exploit the expanded expert landscapes without suffering from the exploration-exploitation dilemmas that affect other RL algorithms.

Computational efficiency considerations. While N=40 achieves the best absolute performance, the marginal improvement over N=30 (averaging 1.3% across datasets) comes at significant computational cost. Training time increases approximately linearly with expert count, with N=40 requiring 30-40% longer training than N=10. For resource-constrained deployments, N=20 or N=30 provide excellent trade-offs, achieving 95-98% of N=40’s performance at substantially lower computational requirements.

Cross-horizon consistency. The performance scaling patterns remain consistent across forecast horizons. While absolute improvements diminish for longer horizons (as forecasting difficulty increases), the relative advantage of larger expert capacities persists. For ETTh1 at horizon 720, improvements progress from 45.7% (N=10) to 50.9% (N=40), maintaining the same 5.2 percentage point gain observed at horizon 96.

5.3. Ablation Studies

Table 8 presents systematic ablation experiments examining each component’s contribution to overall forecasting performance. We evaluate different architectural variants to validate their individual and synergistic effects.

The single-domain baseline analysis in Table 8 reveals fundamental limitations of uni-domain processing. The time-only configuration (T-Block alone) achieves 0.521±0.024 MSE on ETTh1, demonstrating reasonable temporal modeling but struggling significantly with periodic structures. The frequency-only configuration with Extended DFT (F-Block + ExtDFT) reaches 0.549±0.026 MSE (5.4% worse than time-only), highlighting frequency domain’s challenges in capturing irregular dynamics independently. Critically, removing Extended DFT (F-Block without ExtDFT) further degrades performance to 0.573±0.028 MSE (4.4% worse than with ExtDFT), providing direct evidence that our Extended DFT alignment strategy from Equation (4) is essential for effective frequency domain processing.

The progressive integration experiments demonstrate each component’s incremental contribution. Combining both domains without key components achieves 0.497±0.023 MSE (4.6% improvement over time-only), confirming basic domain complementarity. Adding Extended DFT improves to 0.468±0.022 MSE (additional 5.8% gain), validating that aligning input spectrum with complete series frequency grid substantially enhances F-Transformer effectiveness. Adding RL routing alone (without ExtDFT) reaches 0.483±0.023 MSE, while adding ensemble weighting alone achieves 0.489±0.023 MSE, demonstrating both components require proper frequency alignment to deliver full benefits.

Combining ExtDFT with ensemble weighting achieves 0.453±0.021 MSE, while combining ExtDFT with RL routing reaches 0.387±0.018 MSE, demonstrating that RL routing provides substantially stronger individual contribution (17.3%) than ensemble weighting (3.2%). However, the full model with all components achieves 0.267±0.011 MSE, delivering additional 31.0% improvement over ExtDFT+RL alone, confirming that ensemble weighting provides crucial refinement and that all components work synergistically.

The routing mechanism comparison provides decisive evidence for RL-based expert selection. With all other components enabled, Fixed FFN (no gating) achieves 0.453±0.021 MSE. Static gating methods show varying performance: Softmax Gating improves to 0.441±0.021 MSE, Noisy Top-k Gating reaches 0.445±0.021 MSE, and Top-1 Gating achieves 0.449±0.021 MSE. Our RL Routing dramatically outperforms all alternatives at 0.267±0.011 MSE, delivering 39.5% improvement over the strongest static method (Softmax Gating) and 41.1% improvement over Fixed FFN.

This substantial superiority stems from RL’s ability to optimize sequential decision-making through Generalized Advantage Estimation (Equations (21) and (28)), simultaneously balancing prediction accuracy (${\alpha }_1$), routing stability (${\alpha }_2$), and expert diversity (${\alpha }_3$). Static gating methods optimize only for immediate task performance without considering long-term routing consequences or expert utilization balance, resulting in suboptimal expert selection patterns.

Across all five datasets, RL Routing maintains consistent 39.5-43.8% improvements over static alternatives, with particularly strong gains on periodic datasets (ETTh1: 39.5%, Electricity: 43.7%, Traffic: 47.2%) where F-MoE agents successfully learn specialized spectral routing patterns for different harmonic structures. The Weather dataset shows 43.4% improvement, confirming RL’s adaptability even for trend-dominated non-periodic series where T-MoE routing becomes more critical.

5.4. Expert Specialization and Routing Behavior

We analyze the learned expert specialization patterns to understand how RL routing naturally discovers functional differentiation among experts. Training on ETTh1 (horizon 96) reveals clear hierarchical patterns: In F-MoE, Expert 1 progressively dominates from 31.2% usage in Layer 1 to 42.6% in Layer 3, indicating specialization in dominant frequency components. The usage entropy decreases from 1.38 in Layer 1 to 1.31 in Layer 3, demonstrating that deeper layers develop stronger expert preferences for refined spectral pattern processing. Expert 2 maintains stable 24.8-28.4% usage across layers, suggesting specialization in secondary harmonic structures, while Experts 3 and 4 show declining usage (22.1%→18.1% and 18.3%→14.5%) but remain active due to the diversity bonus ${\alpha }_3^F=0.05$ in Equation (17) that prevents complete collapse.

In T-MoE, Expert 1 similarly emerges as dominant but with more balanced distribution (33.2% in Layer 3), reflecting temporal domain’s need to handle diverse patterns including trends, seasonality, and irregularities simultaneously. The higher average entropy (1.37 vs. 1.31) confirms that T-domain maintains more balanced expert utilization, consistent with the inherently heterogeneous nature of temporal patterns that cannot be cleanly decomposed like frequency components.

The progressive concentration of expert usage in deeper layers validates our hierarchical RL routing design: early layers learn coarse-grained routing patterns that distribute representations broadly across experts, while deeper layers refine routing decisions based on accumulated contextual information from attention mechanisms. The diversity bonus successfully prevents pathological collapse to single experts while allowing natural specialization, achieving an optimal balance between expert efficiency and coverage that static gating methods struggle to attain.

Cross-domain comparison reveals complementary specialization strategies: F-MoE develops strong expert hierarchies suited for harmonic decomposition, while T-MoE maintains balanced utilization suited for diverse temporal dynamics. This validates our independent dual RL framework design from Section 3.6 and Section 3.7, confirming that different domains require distinct routing strategies that cannot be achieved with shared policies. The learned routing behaviors demonstrate that our RL framework successfully discovers domain-specific expert specialization patterns that align with the underlying characteristics of frequency and temporal representations.

6. Conclusion

This paper presents MoE-Transformer, a reinforcement learning-driven dual-domain mixture-of-experts framework for multi-horizon time series forecasting. The proposed approach integrates frequency-domain processing with sparse expert networks and dynamic routing to address limitations in existing temporal prediction methods. Key contributions include: Extended Discrete Fourier Transform resolving spectral misalignment for effective frequency-domain forecasting; dual-domain architecture with independent F-MoE and T-MoE modules enabling specialized pattern extraction; reinforcement learning-based routing framework achieving adaptive expert selection through multi-objective optimization; and comprehensive experimental validation demonstrating 50.9-56.9% MSE reduction across diverse benchmarks.

Experimental results confirm the framework’s superiority, achieving substantial performance improvements across ETTh1, ETTm1, Weather, Electricity, and Traffic datasets for 96, 192, 336, and 720-step horizons, while consistently outperforming PatchTST, TimesNet, and other state-of-the-art baselines. The RL routing mechanism delivers 39.5-47.2% improvement over static gating methods, with expert specialization analysis revealing distinct domain-specific routing strategies. The framework maintains 60% faster inference and 40% memory reduction through sparse expert activation.

Future work will focus on extending the framework to multivariate forecasting with cross-variate dependencies, developing probabilistic variants for uncertainty quantification, and exploring applications in financial prediction and industrial monitoring to validate broader applicability.