arXiv CSopen access

Barrier-enforced multi-objective optimization for direct point and sharp interval forecasting

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

Barrier-enforced multi-objective optimization for direct point and sharp interval forecasting ∗ Worachit Amnuaypongsa1 , Yotsapat Suparanonrat1 , Pana Wanitchollakit2 , and Jitkomut Songsiri†1 1 2

Department of Electrical Engineering, Faculty of Engineering, Chulalongkorn University, Bangkok, Thailand Department of Computer Engineering, Faculty of Engineering, Chulalongkorn University, Bangkok, Thailand e-mail: [email protected], [email protected], [email protected], [email protected]

arXiv:2604.18492v1 [cs.LG] 20 Apr 2026

April 21, 2026

Abstract This paper proposes a multi-step probabilistic forecasting framework using a single neural-network based model to generate simultaneous point and interval forecasts. Our approach ensures non-crossing prediction intervals (PIs) through a model structure design that strictly satisfy a target coverage probability (PICP) while maximizing sharpness. Unlike existing methods that rely on manual weight tuning for scalarized loss functions, we treat point and PI forecasting as a multi-objective optimization problem, utilizing multi-gradient descent to adaptively select optimal weights. Key innovations include a new PI loss function based on an extended log-barrier with an adaptive hyperparameter to guarantee the coverage, a hybrid architecture featuring a shared temporal model with horizon-specific submodels, and a training strategy. The proposed loss is scale-independent and universally applicable; combined with our training algorithm, the framework eliminates trial-and-error hyperparameter tuning for balancing multiple objectives. Validated by an intra-day solar irradiance forecasting application, results demonstrate that our proposed loss consistently outperforms those in current literature by achieving target coverage with the narrowest PI widths. Furthermore, when compared against LSTM encoder-decoder and Transformer architectures–including those augmented with Chronos foundation models–our method remains highly competitive and can be seamlessly adapted to any deep learning structure.

1 Introduction In time series forecasting applications for a variable y, predictions are generally categorized into point forecasts and interval forecasts. A point forecast, denoted as ŷ, provides a single estimate for a future time step, typically representing the most likely value or the mean of the predictive distribution. In contrast, an interval forecast–commonly as a form of prediction interval (PI) quantifies the uncertainty of the estimate. It is defined by a range [ˆl, û] that is expected to contain the true future value of y with a pre-specified confidence level. In renewable energy (RE) applications, the random nature of RE poses significant challenges to grid stability and energy management system eﬀiciency. Point forecasts are used across multiple timescales in unit commitment and economic dispatch to maintain the balance between supply and demand while ensuring grid stability and economic eﬀiciency. Shifting to risk-based decision making, PIs provide the necessary bounds for determining operating reserve, and robust optimizations in grid management [EMMS12, ZWS22, AWS25]. The cost of point and probabilistic forecast errors is primarily driven by the magnitude and sign of the deviation alongside the price spread between markets, as well as the PI quality. These inaccuracies lead to significant grid operational expenses, ineﬀicient dispatch, and the economic burden of maintaining reserves to cover uncertainty [GZK+ 24]. Given the importance of PI estimation, an extensive literature on quantile regression (QR) [Koe17] and QR neural network (QRNN) variants using pinball loss has emerged; however, producing high-quality intervals remains challenging due to the inherent trade-off between reliability (how well the PI covers the actual samples) and sharpness (how tight the interval is). These two metrics are typically measured by the Prediction interval coverage probability (PICP) and the PI width. These objectives are naturally conflicting: a high PICP–indicating superior coverage–often results in a wider interval. In practice, excessive width leads to over-conservatism in uncertainty-aware system planning; in power systems, this incurs higher operational costs as more resources must be reserved for worst-case scenarios. To address this, several studies propose specialized loss functions to balance these goals by penalizing PI width. Commonly, PI loss functions merge PICP and width into a scalar value, such as the multiplicative CWC [KNCA11] or additive DIC [ZWW+ 15]. Recent work has shifted toward gradient-based formulations like the likelihood-based QD loss [PBZN18], which uses smooth PICP ∗ This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible. † Corresponding author.

approximations. This framework has been refined in recent variants [SLG24]. Reducing PI over-conservatism is a primary focus of [AWS24], which introduced a large-width penalty using the sum-of-K largest function. This was evolved into the Sum-k loss [AWS25] used with NN-based models to achieve even tighter intervals. In contrast to the QRNN approach, these methods explicitly incorporate PI width as an optimization objective rather than a secondary outcome. Nevertheless, they require hyperparameters to balance PICP and PI width. Since achieving the target PICP is the priority, performance depends on carefully tuning these parameters. Risk-aware decision-making requires both point forecasts and uncertainty quantification through PIs. Consequently, research has shifted toward integrated NN-based models capable of generating both outputs simultaneously. Literature in this direction proposed a loss function that merges point and PI objectives into a single value. Based on quantile loss, [KJL20] provided a framework for NNs to produce a point prediction and valid PI by leveraging a quantile-based approach and split conformal prediction. While the paper addresses valid coverage, it does not fully optimize for maintaining the narrowest possible intervals. Several approaches introduced complex loss functions that require careful hyperparameter tuning to balance competing objectives. For example, extended from QD loss, the QD+ loss [PBZN18] added MSE as a point loss and a penalty for coverage violations. IPIV [SKR22] constructed point forecasts via a convex combination of PI bounds, using a weighted sum of PI and regression objectives, which balances PICP and width via a weight parameter. A midpoint approach [LSH+ 22] defined the point forecast as the PI midpoint and replaced the Gaussian log-likelihood variance term with a width penalty, though it lacks explicit PICP guarantees and requires extra tuning. Further complexity is introduced in the EMVE loss [SLZA25] which combined modified pinball loss, ordering quantile penalties, and coverage/width terms to address extreme tail estimation, significantly increasing complexity. Finally, the EMVE loss [SLR+ 25] integrates width, Gaussian likelihood, and quantile-violation penalties; however, because these terms use disparate units and rely on a Gaussian assumption, the framework may not generalize well to non-Gaussian data. Other methodologies focused on specific neural network architectures or training strategies. For instance, PI3NN [LZLZ22] employed three separate networks for ˆl, ŷ, and û–each trained independently with MSE–which increases model complexity and computational overhead. Alternatively, DualAQD [MS23] utilized a two-stage process where a point forecasting model is pre-trained before estimating PIs using an adaptively weighted sum of width and coverage penalties. However, because these loss terms are not properly scaled, controlling their relative importance during training remains a challenge. In a general context, forecasting quality depends on model architecture, the loss function, training mechanisms, and input selection. For multi-step point and probabilistic forecasting, the primary objective is producing a valid PI that covers the point prediction and achieves the target PICP, as insuﬀicient coverage–such as in RE reserve preparation–can lead to a higher deficit penalty like the value of lost load than other penalties [ZWS21]. Once the PICP target is reliably met, the secondary goals are to minimize PI width to reduce over-conservatism and maintain acceptable deterministic forecasts. Balancing these goals across multi-step horizons is particularly challenging, as interval width naturally expands with increasing lead times. While existing literature addresses these goals to some extent, the main challenge remains the simultaneous optimization of multiple objectives. Using a scalarized loss requires a proper balance of weights, which becomes diﬀicult to tune when objective terms overlap in function or are not properly scaled. This paper addresses these critical points by proposing the following contributions: 1. Structural enforcement of valid PIs: Unlike existing literature that relies on complex penalty terms to prevent bound crossing or to force ŷ within the PI, we structurally design the network output layer to satisfy these constraints. This architectural approach ensures valid PIs by construction and reduces the need for complex penalties. 2. Barrier-enforced composite loss: We propose a novel loss function combining regression loss for point estimation, PI width penalties, and a new PICP loss. The latter utilizes an extended log-barrier function to strictly penalize coverage violations. All terms are designed to be scale-independent for universal applicability. 3. Multi-objective training algorithm: By treating point and PI estimation as dual objectives, we apply a training algorithm that identifies a common gradient descent direction to mutual improvement in both objectives. This approach eliminates the weight tuning typically required by scalarized loss functions in existing literature. 4. Seamless architecture integration: The proposed smooth loss and training mechanisms are easily integrated into advanced deep learning frameworks, demonstrated by augmenting foundation models into tailor-made probabilistic forecasting architectures. These improvements are integrated into our proposed framework, SolarPointPI, consisting of a horizon-specific NN architecture, a novel loss function, and a targeted training strategy. While the framework is universal, we demonstrate its eﬀicacy through solar irradiance forecasting, with implementation details tailored specifically for this application. Existing literature on point and interval forecasting of solar irradiance often lacks a direct mechanism for constructing PIs as primary model outputs or incorporating PI widths as a design objective. Most frameworks generate PIs through a secondary stage: for instance, [FKF18] parameterized Beta distributions to estimate CDFs, while others rely on postprocessing techniques like KDE on XGBoost forecasts [LMC+ 22], Gaussian processes [HZW+ 20], or data clustering with k-means [SSP16]. Even modern hybrid architectures that merge statistical methods with NNs, such as those in [GLS23] and [SDS23], typically focus on quantile production rather than width control. This two-stage approach–seen in [WYP+ 17], 2

Future regressor

Auto-lag regressor

Exogeneous-lag regressor Prediction interval construction -step prediction sub-models

Figure 1: A proposed point and PI forecasting model architecture. where deep CNN point forecasts are subsequently converted into quantiles–results in sharpness being an outcome of the statistical distribution rather than a primary design constraint. While advancements in solar energy forecasting could alternatively rely on complex input selection, or high-capacity deep learning models [PTSN+ 23] requiring vast datasets, our framework focuses on practical improvements of PI quality that bypass the prohibitive costs of extensive data requirements.

2 Methodology This section presents the methodology for generating point and PI forecasts in a multi-step-ahead forecasting framework. It consists of three parts: a model architecture, a loss function, and a training algorithm. The model architecture explains how the point forecast and its corresponding upper and lower bounds are structured to ensure that the point estimate remains within the PI. A loss function is proposed to firmly achieve PIs with a desired coverage probability and with a special modification for solar irradiance data that have a diurnal cycle. The training algorithm details how the multi-objective optimization problem is addressed using the multi-gradient descent algorithm (MGDA).

2.1 Model architecture The proposed model is designed to jointly produce the upper bound û, lower bound ˆl, and point forecast ŷ for a forecasting horizon H, ensuring that the point forecast always lies within the PI. Let x denote auto-lag regressors, xex exogenous lag regressors, z future regressors, and y the target vector (for multi-step outputs). Figure 1 illustrates that the model consists of a common model Mc (θc ) and submodels Mk (θk ) with their model parameters, θc and θk for k = 1, . . . , H. The common model, which can be any neural network model, is shared across all forecasting steps and processes historical information from the lagged inputs, x and xex . Each submodel, on the other hand, specializes in a specific forecast horizon, taking as input the hidden representation from the common model together with the future regressor z corresponding to its prediction step. For a k-step forecast, the k th submodel releases three quantities: ŷ, ∆û, and ∆ˆl, and we can construct the PI according to û = ŷ + ∆û and ˆl = ŷ − ∆ˆl. To ensure that the point forecast ŷ remains within the PI, we enforce ∆û and ∆ˆl to remain non-negative using the Softplus activation function, which is strictly positive and differentiable, to avoid the zero-gradient issue and ensure proper gradient flow for the PICP objective. A linear activation is retained for the point forecast ŷ. The proposed model in Figure 1 is flexible and can be served for any forecasting application, allowing users to substitute any model architecture for Mc (θc ) and Mk (θk ). In our experiments, we find that using an LSTM in Mc (θc ) to capture temporal dynamics and a feed-forward neural network in Mk (θk ) is suﬀicient to produce desirable results in our setting. We introduce the SolarPointPI framework, which integrates the proposed model with a new PI and point loss function, and a multi-objective training algorithm - all detailed in subsequent sections, to jointly produce PI and point forecasts.

2.2 Problem formulation In this study, we define two learning objectives: producing high-quality point forecasts and constructing reliable PIs, optimized via a point forecast loss (regression loss) and PI loss, respectively. We denote the model parameters as θ = (θc , θ1 , . . . , θH ). Conceptually, the point forecast ŷ and PI bounds (ˆl, û) are functions of θ, indexed by time t (as sample index), and k = 1, 2, . . . , H for a k-step prediction. Denoted as ŷ(t + k|t; θ), these outputs are calculated at time t for the future time t + k. Consequently, loss functions are typically aggregated across both samples and prediction steps.

Standard log-barrier

r = 0.5 r=1 r=2

4 2

0 3

0 z

=1 =2 =5 = 10 = 20

80 60

)

r(z)

r = min 100, Deviation of PICP (

100

r=1 r=2 r=5 r = 10 r = 25 r = 50

Extended log-barrier

40 20 3

0 z

(a) Illustration of the extended log-barrier function ψr (z).

0.25 0.00 0.25 0.50 0.75 1.00 Deviation of PICP

(b) Illustration of adaptive r formulation.

Figure 2: Penalty function components showing the extended log-barrier and the adaptive parameter r. Point estimation loss The point estimation loss can be employed using any standard regression loss ℓ and is normalized by data range RQ , the number of samples N , and the number of prediction steps. LPoint (θ) =

XX 1 1 X LPoint,k-step (θ) = ℓ(y(t + k), ŷ(t + k|t; θ)). H HN RQ t=1 H

k=1

(1)

The RQ normalization factor is is defined as the interquantile range RQ = qy (0.95) − qy (0.05), rendering the loss function unitless and independent of variable magnitude. Instead of using min-max range, RQ is less sensitive to outliers. This work applies the MAE loss, ℓ(y, ŷ) = |y − ŷ| as it offers greater robustness to outliers compared to the standard MSE loss. PI estimation loss The primary objective of PI estimation is to guarantee the target PICP across multi-step horizons, followed by minimizing interval width. In RE applications, ensuring this coverage is particularly challenging when managing reliability to avoid reserve deficit penalties [ZWS21]. To address this, let the PICP deviation from a target probability p be z = p − PICP, which is desired to be non-positive. Constrained optimization often enforces the requirement z ≤ 0 using the log-barrier function ψr (z) = − 1r log(−z), which requires z to be strictly negative. However, as a PI loss, the log-barrier is undefined for z ≥ 0, rendering the loss ill-defined if the model is initialized in an infeasible region. To resolve this issue while preserving a strong penalty on the PICP deviation, we propose to use the extended log-barrier penalty that strongly penalizes z. ( − 1r log(−z), if z ≤ − r12 , ψr (z) = (2) rz − 1r log r12 + 1r , otherwise. The extended log-barrier function with parameter r > 0 was used in [KDY+ 22] in the context of learning deep neural networks with constraints. The logarithm term in the standard log-barrier is replaced with a first-order Taylor expansion beyond a threshold −1/r2 , yielding a function that is strictly convex, twice differentiable, and defined for all z ∈ R, thereby guaranteeing continuous gradients throughout training even when the coverage constraint is violated. The hyperparameter r acts as a barrier softening factor; the larger, the higher penalty on the infeasible region as shown in Figure 2a. To dynamically tighten the penalty boundary and increase the penalty on the model when the target coverage pk is far from the desired PICP, we calculate r adaptively and inversely proportional to the PICP deviation: ρ r = min 100, , (3) |p − PICP| where a multiplier factor ρ is empirically set to 10 to ensure a stable and reliable optimization process (Figure 2b). This adaptive mechanism allows r to remain small when the model is far from the target, effectively softening the extended log-barrier function ψr (z) (Figure 2ba). This prevents aggressive gradients from forcing parameters too abruptly during early training. As the PICP approaches the target, r increases toward its cap of 100, hardening the penalty boundary to strictly enforce the coverage constraint. Consequently, the coverage penalty ψr (pk − PICPk (θ)) is applied across each k-step prediction to ensure reliable interval estimation. 4

For PI width penalty, denote an interval width as w(θ) = û(θ) − ˆl(θ) (indexed by samples and k-step predictions) and is a function of θ. We apply the concept of the Sum-k loss function [AWS25] which penalizes the sum of the k-largest interval widths. This approach effectively reduces extremely large intervals, mitigating over-conservatism when using PI bounds in robust optimization. Defining this width penalty begins with specifying a hyperparameter K (positive integer) as the number of large width samples, such as K = ⌊(0.3N )⌋ meaning that the top 30% of width samples will be penalized heavily. Then, w[i] denotes the i-th largest width according to w[1] ≥ w[2] ≥ · · · ≥ w[N ] . The Sum-k loss penalizes the large-width samples and the remaining smaller-width samples in different degrees controlled by a hyperparameter λ ≥ 0. The width penalty in the Sum-k loss is described by " # K N X 1 1 X λ W(θ) = w[i] (θ) + w[i] (θ) . (4) RQ K i=1 N −K i=K+1

Setting λ < 1, such as 0.8 in our experiment, means that we penalize more heavily on the average large width as compared to the average of smaller widths. The normalization term RQ scales the loss relative to the target variable range, making W unitless and ensuring comparability across different datasets. While the sum-of-k-largest function is non-differentiable, frameworks like PyTorch handle it effectively using a subgradient approach. Proposed PI loss. Combining the PI coverage and PI width penalty described above, the overall PI loss across all H forecasting horizons is defined as the average of the k-step prediction loss LPI (θ) =

1 X 1 X LPI,k-step (θ) = ψr (pk − PICPk (θ)) + Wk (θ). H H H

k=1

(5)

Each term in LPI is unitless and scale-independent. It requires only one user-defined hyperparameter: λ in (4). Smooth approximation of calculating PICP. Conventionally, the standard counting function, 1{y ∈ [ˆl, û]} is used to compute PICP. To preserve differentiability for backpropagation, it is replaced with a smooth tanh approximation: 1 X 1tanh {yi ∈ [ˆli , ûi ]}, N i=1 where 1tanh {y ∈ [ˆl, û]} = (1/2) max{0, tanh s(y − ˆl) + tanh s(û − y) }, and s > 0 is a smoothing factor. N

PICP =

(6)

PI coverage loss for solar energy. In solar irradiance forecasting, irradiance is zero at night. Enforcing uniform coverage across all hours complicates the trade-off between interval width and coverage, as nighttime zeros are trivially covered even by wide, unoptimized PIs. To address this, we apply the extended log-barrier function separately to daytime and nighttime samples, assigning distinct target coverage levels (p): 0.90 for day and 0.15 for night. Samples are split into day and night sets by directly thresholding irradiance values (e.g., 0.001), bypassing the need for date-time metadata in the training process. We adjust the PICP by applying a smooth mask function, max(0, tanh(·)), that returns a mask value in [0, 1] to isolate daytime and nighttime samples: PICPday = PICPnight =

N 1 X 1tanh {yi ∈ [ˆl, û]} · max 0, tanh(s(yi − threshold)) , Nday i=1

N X

Nnight i=1

1tanh {yi ∈ [ˆl, û]} · max 0, tanh(s(threshold − yi )) ,

(7)

(8)

PN PN where Nday = i=1 max{0, tanh(s(yi − threshold))} and Nnight = i=1 max{0, tanh(s(threshold − yi ))}. Accordingly, the extended log-barrier penalty for each k-step in (5) is split into separate day and night terms: ψr(0.90 − PICPday ) + ψr(0.15 − PICPnight ). We conclude this section by defining the goal of model learning as the simultaneous minimization of two distinct objectives: the point and PI estimation losses as defined in (1) and (5), respectively. Traditional approaches merge multiple goals into a single value as minimizing γ1 LPoint + γ2 LPI , while selecting a value of weights γ = (γ1 , γ2 ) yields a Pareto optimal solution. This leaves the user to find a meaningful balance between two objectives that often results in a trial and error process. Unlike this direction, we treat this task strictly as a multi-objective optimization (MOO) problem. Formally, we seek to optimize the vector-valued loss function: minimize L(θ) = (LPoint (θ), LPI (θ)) . θ

(9)

By adopting the MOO framework instead of scalarization, we avoid the biases inherent in fixed weighting. This approach enables the model to more effectively explore the trade-offs between point estimation accuracy and PI quality. The following section introduces a framework designed to achieve this balance. 5

2.3 Training algorithm We treat the proposed problem as an MOO task and first outline its fundamental optimality. We then utilize the multiple gradient descent algorithm (MGDA) [Dés12, SK18] to identify a common descent direction that simultaneously improves all objectives. This is achieved by dynamically calculating weights at each iteration based on the local geometry of the loss gradients. Finally, we summarize the training algorithm resulting from applying these concepts to our specific problem (9). Multi-objective optimization We first introduce the fundamental concepts of MOO. Consider m tasks, each associated with an objective function Li : Rn → R for i ∈ [m] where [m] denotes the index set {1, 2, . . . , m}. Unlike single-objective optimization, solutions in MOO cannot be ranked using a single scalar criterion; instead, they are evaluated through the concept of dominance. The concept of optimality can be defined as follows. Definition 1. [Mie99] A solution θ(a) ∈ Rn is said to be dominated by another solution θ(b) ∈ Rn iff Li (θ(a) ≤ Li (θ(b) ) for all i ∈ [m] and there exists at least one i ∈ [m] such that Li (θ(a) < Li (θ(b) ). A solution θ⋆ is said to be Pareto optimal if it is not dominated by another solution. A Pareto set is the set of all Pareto optimal solutions. A Pareto front is the set of all objective function values of the Pareto optimal solutions. While gradient-based methods in deep learning typically target stationary points PM for single-objective problems, MOO requires a Pareto stationary definition. Denote ∆m the simplex set: {a ∈ Rm | i=1 ai = 1, ai ≥ 0, i ∈ [m] }. Definition 2. [Dés12] The objective functions Li for all i ∈P [m] are said to be Pareto-stationary at θ⋆ iff there exists a m convex combination of the gradients, that is equal to zero: k=1 γk ∇θ Lk (θ⋆ ) = 0 where γ ∈ ∆m . Pareto stationary is a necessary condition for Pareto optimality. If all Li are convex with γi > 0, i ∈ [m], it is the Karush-Kuhn-Tucker (KKT) suﬀicient and necessary condition for Pareto optimality [Dés12]. While many studies in MOO or multi-task learning aim to identify the entire Pareto front, our objective in (9) is more targeted. We seek a Pareto optimal solution specifically within the high-reliability region where the PICP is high. Consequently, we focus on a neighborhood of the target LPI value rather than the complete front. Research directions in MOO in finding a single Pareto optimal solution generally fall into two categories: scalarization via fixed weights or dynamic weight learning [CLZ+ 25]. Among gradient-balancing methods, we select MGDA because it addresses both magnitude and directional conflict between competing objectives. Furthermore, its inherent scale-invariance eliminates the trial-and-error process of hyperparameter tuning. Multiple gradient descent algorithm MGDA extends the steepest gradient descent to MOO and is proved to converge to a Pareto stationary point. The main result is given as follows. Pm Theorem 1. [Dés12] Consider the convex hull of the objective gradients given by C = {u ∈ Rn | u = k=1 γk ∇Lk (θ), γ ∈ ∆m (the simplex) }. Let g be the minimum-norm element of C at θ0 . Then we have two possible outcomes: i) g = 0 and hence, Li (θ) for i ∈ [m] are Pareto-stationary at θ0 , and ii) when g ̸= 0 then −g is a descent direction common to all the objectives, i.e., ⟨∇Li , −g⟩ ≤ 0, i ∈ [m]. Additionally, if g is interior to C, it can be shown that ⟨∇Li , g⟩ = ∥g∥2 , ∀i ∈ [m]. From Theorem 1, if we choose our search direction as −g, specifically from the minimum-norm element in C, this does not only provide any common descent direction, the improvement for any objective i: Li (θ) − Li (θ − ηg) ≈ −η∇LTi g is at least −η∥g∥2 where η is the step size. It is an optimal descent direction in the sense of fairness. The algorithm ensures that the task with the least favorable gradient (the one most in conflict with others) still receives possible improvement allowed by the geometry. MGDA requires solving the minimum-norm problem to obtain the optimal weight γ: ⋆

γ = argmin γ

m X k=1

γk ∇Lk (θ)

subject to

γ ∈ ∆m (the simplex).

(10)

Pm The problem is to find the shortest vector within the convex hull of objective gradients. If g = k=1 γk⋆ ∇Lk (θ) ̸= 0 then −g is a descent direction that improves all objectives. For m objectives where m ≥ 2, the problem (10) is a quadratic program where [SK18] use the Frank-Wolfe algorithm to solve it. For two objectives, the problem (10) reduces to a one-dimensional problem: γ1⋆ = argmin ∥a∇L1 + (1 − a)∇L2 ∥22 subject to 0 ≤ a ≤ 1, a

γ2⋆ = 1 − γ1⋆ .

(11)

The objective of (11) is a2 ∥∇L1 − ∇L2 ∥2 + 2a(∇L1 − ∇L2 )T ∇L2 + ∥∇L2 ∥22 . The zero-gradient condition with respect T 2 −∇L1 ) ∇L2 ⋆ to a is a0 = (∇L ∥∇L1 −∇L2 ∥2 . Therefore, with constraint 0 ≤ a ≤ 1, the optimal a is obtained by clipping a0 onto [0, 1]. 6

and

Figure 3: Geometric interpretation of the minimum-norm problem [SK18]. The optimal weights γ1∗ and γ2∗ identify the shortest vector in the convex combination (black dashed line) formed by ∇L1 and ∇L2 . The faint grey arrows represent the shortest-norm vector that is a linear combination of two gradients when there is no constraint on γ1 + γ2 = 1. • Case 1: a0 ≤ 0 ⇔ ⟨∇L1 , ∇L2 ⟩ ≥ ∥∇L2 ∥2 ⇔ ∥∇L1 ∥ cos θ ≥ ∥∇L2 ∥. This gives a⋆ = 0, γ1⋆ = 0, γ2⋆ = 1. • Case 2: 0 < a0 < 1 ⇔ ⟨∇L1 , ∇L2 ⟩ ≤ ∥∇L1 ∥2 and ⟨∇L1 , ∇L2 ⟩ ≤ ∥L2 ∥2 . Hence, a⋆ = a0 , γ1⋆ = a0 , γ2⋆ = 1 − γ ⋆ . • Case 3: a0 ≥ 1 ⇔ ⟨∇L1 , ∇L2 ⟩ ≥ ∥∇L1 ∥2 ⇔ ∥∇L2 ∥ cos θ ≥ ∥∇L1 ∥. This gives a⋆ = 1, γ1⋆ = 1, γ2⋆ = 0. In conclusion, the analytical solution to (11) can be compactly expressed as (∇L2 − ∇L1 )T ∇L2 ⋆ γ1 = max min , 1 , 0 , γ2⋆ = 1 − γ1⋆ . ∥∇L2 − ∇L1 ∥22

(12)

This solution is illustrated geometrically in Figure 3. If the gradient of task i is already contained in the direction of task j (examined by ∥∇Lj ∥ cos θ ≥ ∥∇Li ∥), then by following the descent direction of task i, we are guaranteed to decrease loss j as well. The optimal weight is then assigned entirely to task i which provides the most eﬀicient common descent at that specific geometry. Conversely, when gradients are balanced or opposing, weights are distributed between both tasks to establish a balanced, common descent direction. To optimize model parameters, we integrate the common descent direction −g from MGDA into the Adam optimizer. While Adam typically tracks moments of a stochastic gradient, we utilize it to track the MGDA-derived consensus direction. This approach is justified as −g satisfies the descent condition for all objectives simultaneously. Under the framework of [LRJ23], convergence is maintained for such descent directions even under relaxed smoothness assumptions. Training procedure By combining the proposed loss, the adaptive hyperparameter rule, and MGDA, we establish the comprehensive training process shown in Algorithm 1. Denote L1 = LPoint and L2 = LPI . At the start of each epoch, the model computes PICP(train) using the full training set. This step allows the extended log-barrier parameter r to update based on current coverage levels before any weights change. Inside each mini-batch, the model calculates two different forms of loss to guide the learning process. The first form is the weighted loss, Lweighted , which is used for gradient-based updating θ with a learning rate η. This loss is a combination of the point estimation loss L1 and the PI loss L2 evaluated on the training set, balanced by the MGDA weights γ1 and γ2 . These weights represent the relative importance of each task gradient calculated from (12). By using these γ values, the model finds an update direction for model parameters with the Adam optimizer that helps both point and PI forecasts without one task harming the other. The learning rate is managed by a cosine annealing schedule that includes a linear warm-up stage to ensure stable convergence. The second form is the total validation loss, Lval , which is the sum of L1 and L2 on the validation set. This total loss acts as the monitor for the early stopping rule. While the weighted loss changes in every batch to balance gradients, the total validation loss provides a stable metric to check for overfitting. Training stops if Lval does not improve for a specific number of epochs, known as patience.

Training and validation loss

3 × 10 1

101

2 × 10 1 Loss

100

10 1 10

20 Epoch

Gradient norm

102

20 Epoch <

103 102 Inner product

101 100

Inner product

104

Point PI

Point, Point 2 PI 2

PI >

20 Epoch

Adaptive weight

0.8

101

0.6

100

0.4

10 1

0.2

10 3 0

4 × 10 2

1.0

10 2

10 1

10 1 6 × 10 2

100 0

Gradient norm

Weighted total loss

Total (training) Total (validation)

Loss

101

Loss

Total loss

Point (training) Point (validation) PI (training) PI (validation)

20 Epoch

0.0

1( Point) 2( PI)

20 Epoch

Figure 4: MGDA training process: (top) point and PI losses, total loss, and weighted loss; (bottom) gradient norms, inner product between task gradients, and adaptive weights γ1 , γ2 . Algorithm 1 Point and prediction interval learning via MGDA Require: Problem parameters: target probability p, Sum-k loss parameter λ < 1 (for solar-context, target probabilities: pday = 0.90, pnight = 0.15) Algorithm parameters: lr, min_epoch, max_epoch, patience 1: Initialize model parameters θ ▷ initial guess for weights and biases 2: L∗ ▷ initialize best validation loss to infinity val ← ∞ 3: patience counter ← 0 4: for each epoch e = 1, 2, . . . do 5: for each k-step do (train) 6: Compute PICPk on training set ▷ Equation (6) 7: Update log-barrier parameter rk based on coverage deviations using p ▷ Equation (3) 8: for each mini-batch B do ▷ MGDA inner loop 9: Infer model outputs: (ˆli , ûi ) and ŷi for all i ∈ B 10: Compute point and PI losses: L1 and L2 ▷ Equation (1) and Equation (5) 11: for each task j ∈ {1, 2} do 12: gj ← ∇ θ Lj ▷ compute task gradient 13: Solve MGDA for (γ1∗ , γ2∗ ) ▷ using (12) 14: θ ← Adam(θ, ∇θ Lweighted = γ1∗ g1 + γ2∗ g2 ) ▷ update parameters via Adam optimizer (val) (val) 15: Compute validation loss Lval ← L1 + L2 ∗ ∗ ∗ 16: if Lval < Lval , then Lval ← Lval ; θ ← θ; patience counter ← 0; else patience counter ← patience counter + 1 17: if e ≥ min_epoch and (patience counter ≥ patience or e ≥ max_epoch); then break 18: return model with parameters θ ∗ Example of training and validation loss Figure 4 illustrates the training process of the MGDA during initial epochs. The training and validation loss curves show that the PI loss (LPI ) spikes sharply in the first few epochs. This reflects the large log-barrier penalty when the initial coverage is far below the target. Meanwhile, the point forecast loss (LPoint ) stays low and stable throughout the training. Both losses converge within five epochs, and their validation curves track closely with the training curves, indicating that the model does not overfit. The total loss (LTotal ) also drops steadily after the initial steps. The weighted total loss reveals occasional sharp spikes at later epochs. These spikes coincide with the adaptive tightening of the barrier parameter (r) as 8

the coverage probability approaches its target, and they resolve quickly as the optimizer adjusts. The behavior of the adaptive weights γ is directly governed by the geometric relationship between the loss gradients. Initial convergence (epochs 0–7): the gradient of LPoint is substantially smaller than and contained within the direction of LPI . Geometrically, the projection of ∇LPI onto ∇LPoint exceeds ∥∇LPoint ∥2 . Consequently, MGDA assigns the entire weight to the point loss (γ1 = 1), prioritizing the refinement of the more restrictive task. Gradient balancing (epochs 25–35): In the later stages of training, the two gradients become more comparable in orientation. Because neither task clearly contains the other, the optimizer distributes weights between γ1 and γ2 . This distribution represents a balanced trade-off, ensuring a common descent direction that simultaneously satisfies both point accuracy and PI objectives. MGDA consistently prioritizes the point loss (γ1 ≈ 1) as a normalizing force. Under standard scalarization, the PI loss would dominate the updates because its gradient is roughly 100 times larger, effectively obscuring the point loss signal. By assigning the weight to the smaller, more restrictive point loss gradient, MGDA prevents this imbalance. Because the gradients are largely aligned, the PI loss still improves as a free rider, while the model ensures that the sensitive point estimation remains accurately minimized. Notes on multi-objective framework. The challenge of point and PI forecasting can be framed as a multi-objective optimization problem in several ways. While PICP and PI width are fundamentally conflicting, the relationship between point estimation loss and PI loss is more synergistic. Given that these losses are aggregated over a k-step horizon, we explored the following frameworks: • Explicit PICP-width splitting: We initially treated PICP and PI width as two independent objectives. However, the model struggled to satisfy both simultaneously, often converging to a trivial solution with near-zero PICP and a small PI width. • Composite-loss MGDA: We tried constructing two distinct composite losses, each representing a different weighted combination of PICP and PI width. Loss 1 was designed to penalize PI width more aggressively, while Loss 2 prioritized coverage. We applied the MGDA to dynamically determine the optimal weights for these two composites. Despite this, the resulting PICP remained below target levels. • Multi-objective k-step losses: One could define H distinct losses, where each loss k combines point and PI estimation for the k-th prediction step. However, these losses are highly correlated rather than contradictory, offering little tradeoff for a multi-objective framework to exploit. Furthermore, this approach still requires pre-defined weights to balance the point and PI objectives within each step.

3 Data description 3.1 Dataset Measurements of solar irradiance. The historical global horizontal irradiance data utilized in this study were collected from 104 measurement stations distributed across Thailand. These stations are operated by the Department of Alternative Energy Development and Eﬀiciency (DeDe), Ministry of Energy, Thailand, as well as private solar farms. All measurements were sampled at a 15-minute resolution. The temporal coverage varies by source: the DeDe dataset spans from January to December 2023, whereas the data provided by the private solar farms cover the period from January 2022 through June 2023. The geographical distribution of these 104 sites is illustrated in Figure 5a. Missing data were addressed via linear interpolation for short gaps lasting less than 6 hours; longer periods of missing observations were retained as NaN values and subsequently excluded during the batch data preparation process. Re-analysis weather data. To supplement the ground-based measurements, we incorporated historical reanalysis solar irradiance data (Icams ) from the Copernicus Atmosphere Monitoring Service (CAMS). It is crucial to note that this dataset consists of retrospective, reanalyzed historical records rather than real-time monitoring or operational forecasts. CAMS derives the surface solar irradiance by integrating satellite-based cloud parameter observations with detailed aerosol and atmospheric composition data using the Heliosat-4 method [QOB+ 17]. This reanalysis data provides a spatially and temporally consistent baseline that accounts for historical atmospheric attenuation. Satellite-derived cloud data. The cloud cover data were extracted from Himawari-8 satellite RGB imagery [CUE23], possessing a spatial resolution of 2 km × 2 km. These images were captured daily between 06:00 and 19:50 at 10-minute X−LB where X represents the intervals. The cloud index (CI) is computed by normalizing the raw pixel color values: CI = UB−LB raw pixel value, and LB = 0 and UB = 255 denote the lower and upper bounds of the 8-bit color spectrum, respectively. In this study, we specifically utilized the cloud mask (CICM ) and the red-channel cloud index (CIR ) due to their strong statistical correlation with ground-level solar irradiance attenuation. An example of the Himawari-8 cloud imagery is depicted in Figure 5b. To handle missing data, gaps shorter than one hour were interpolated, followed by rolling mean smoothing and temporal resampling to 15-minute resolution. 9

(a) Site location map.

(b) Himawari-8 imagery example.

Figure 5: Overview of geographical data sources: (a) Geographical distribution of the 104 solar measurement stations across Thailand and (b) an example of the Himawari-8 satellite RGB cloud imagery utilized in this study. Solar irradiances 800

W/sqm

600 400 200 0 2023-01-01 Actual irradiance (I)

2023-01-02

Time

2023-01-03

Clear sky irradiance (Iclr)

2023-01-04 CAMs irradiance (Icams)

Figure 6: A sample time-series of irradiance data. Clear-sky irradiance. The clear-sky irradiance (Iclr ) represents the theoretical maximum solar radiation reaching the surface under cloud-free conditions. It is estimated using an improved clear-sky model [STT+ 25, Ene24] derived from the original Ineichen–Perez formulation [IP02], with the Linke turbidity factor (TL ) recalibrated for Thailand to improve local accuracy. The model takes site-specific latitude, longitude, and altitude as inputs to estimate the ideal daily irradiance profile. From this baseline, the clear-sky index is defined as k = I/Iclr , which normalizes observed irradiance against the theoretical maximum. This index captures the attenuation of solar radiation due to cloud cover and aerosols, effectively isolating stochastic atmospheric variability from deterministic diurnal and seasonal cycles. Figure 6 presents a sample time-series profile illustrating the measured irradiance, the clear-sky irradiance, and the CAMS irradiance.

3.2 Data arrangement The experiments considered in this paper involve various model architectures with different data arrangement requirements. Consequently, we employ two distinct data splitting strategies: a stratified shuffling technique and a spatial cross-validation technique. The former is used to evaluate objective functions, while the latter is used to benchmark model architectures. Regressors. We incorporate three types of data inputs for the model: auto-regressive features consisting of historical irradiance measurements, exogenous lag regressors obtained from the cloud index, and future regressors including CAMS irradiance, clear-sky irradiance, and an encoded hour index. To preserve periodic continuity between hour 23 and hour 0, the hour index is transformed using a sine-based cyclical encoding: sin(πh/24), where h denotes the hour index.

Data for the experiment of comparing loss functions. The experiment in Section 4.2 compares the effectiveness of loss functions and utilizes the architecture shown in Figure 1 with a 4-hour historical window and a 4-hour forecast horizon. To ensure a uniform data distribution, days are categorized by sky condition (clear, partly cloudy, or cloudy) using the clear-sky index. We then apply daily-block shuffling, where entire 24-hour segments with the same location are randomly permuted. This approach preserves the intraday temporal correlations of each daily time series while randomizing their chronological order. The dataset of 4,120,032 samples is partitioned into train:validation:test using an 80:10:10 ratio. A stratified split is employed to ensure a balanced representation of sky conditions across all subsets. Data for the experiment of comparing model architectures. In the experiment comparing model architectures (Section 4.3), specific deep learning implementations require data to be fed in chronological order. Consequently, rather than shuffling temporal blocks, we employ a spatial cross-validation strategy. The dataset of 3,922,496 samples is partitioned entirely by the geographic location of solar measurement stations using an 80:10:10 ratio for train:validation:test. Each subset is curated to ensure spatial diversity, covering all regions of Thailand (North, East, South, and Central). By maintaining the full temporal sequence for each station, every split inherently captures a complete year of seasonal patterns while ensuring the model generalizes across different geographic locations.

4 Experimental results This section evaluates the performance of the proposed framework against established baselines and current literature across two key dimensions: the eﬀicacy of the loss function and the integration of deep learning architectures. Furthermore, we detail the benchmarking framework, the standard metrics employed for evaluation, and the relevant computational considerations. Python implementation is available at https://github.com/energyCUEE/PIPointForecast.

4.1 Evaluation metrics Let y and ŷ denote the ground-truth irradiance and its point estimate, respectively, with [ˆl, û] representing the prediction interval. Regressors are denoted by x. The evaluation metrics are aggregated over the test set with N samples. Typical metrics for evaluating point forecast accuracy are MAE, RMSE, and MBE, described by v u N N N u1 X 1 X 1 X |yi − ŷi |, RMSE = t (yi − ŷi )2 , MBE = (yi − ŷi ). (13) MAE = N i=1 N i=1 N i=1 We assess the trade-off between coverage reliability and interval sharpness using several PI metrics. To facilitate comparison across datasets from different studies, these metrics are normalized to ensure scale-independence. To avoid the sensitivities of a min-max range, which can be skewed by outliers, we normalize by the inter-quantile range: RQ = qy (0.95) − qy (0.05). • Prediction interval coverage probability (PICP): This metric quantifies interval reliability by measuring the proportion of observed values that successfully fall within the estimated PI. The PICP is defined as: 1 X ˆ PICP = 1{li ≤ yi ≤ ûi }, N i=1 N

(14)

where 1{E} is a counting function equal to 1 if E is true, and 0 otherwise. The PICP is expected to approach the nominal confidence level p specified in the training process. • Prediction interval normalized average width (PINAW): This metric evaluates the sharpness of PI where lower PINAW values signify narrower, sharper intervals, and are highly preferable provided the target PICP is achieved. PINAW =

N 1 X (ûi − ˆli ), N RQ i=1

(15)

Without data-scale normalization, the mean prediction interval width (MPIW) is also used in literature. • Prediction interval normalized average large width (PINALW): While average width is a standard metric, it may not fully capture the critical characteristics of PIs, particularly since large width samples can disproportionately drive the conservatism in robust design. To address this, PINALW [AWS25] computes the average of the K-largest PI widths, specifically targeting samples that exceed the τ -quantile of the width distribution. PINALW(τ ) =

K 1 X |w|[i] , KRQ i=1

(16)

where |w|[i] denotes the ith largest absolute PI width, i.e., |w|[1] ≥ |w|[2] ≥ · · · ≥ |w|[N ] , and K = ⌊(1 − τ )N ⌋ is the total number of samples exceeding the τ -quantile threshold. For example, in our evaluation, τ is set to 0.5, giving the average of the largest 50% of all widths to represent the large width samples. • Winkler score [Win72]: The Winkler score offers a comprehensive evaluation by jointly capturing sharpness and reliability. It penalizes wide intervals and imposes an additional, more severe penalty when an observation falls outside the constructed bounds. For a PI with nominal probability p, the normalized Winkler score is defined as: N i 1 X 2 hˆ ˆ ˆ Winkler = |ûi − li | + (li − yi )1{yi < li } + (yi − ûi )1{yi > ûi } . N R i=1 1−p

(17)

A lower Winkler score indicates superior model performance, implying that the constructed lower and upper bounds align tightly with the theoretical quantiles of (1 − p)/2 and (1 + p)/2. While other probabilistic metrics exist, such as CRPS [YvdMM20], they are inapplicable to direct PI forecasts. Additionally, although CWC [KNCA11] combines PICP and width into a single score, it requires a hyperparameter to balance these terms; we prefer to evaluate them individually to maintain clarity. Evaluation setting. We evaluate solar irradiance forecasting performance exclusively during daytime hours (06:00–18:00) to ensure the metrics reflect actual modeling challenges rather than trivial nighttime zeros. Since nighttime data √ accounts for half of the dataset, its inclusion would artificially lower the MAE by a factor of 1/2 and the RMSE by 1/ 2, leading 2 to an overly optimistic assessment of model accuracy [Son25]. The MAE, RMSE, and MBE are reported in W/m for comparison consistency with other solar energy studies, while the PI metrics are scale-independent and unitless.

4.2 Loss comparison This section compares the eﬀicacy of the proposed loss function against several established methods from prior research. To ensure a fair comparison, the underlying model architecture is kept consistent with our proposed design shown in Figure 1 across all experiments. Specifically, for the QR, QD+, and EMQ methods, the models are configured to directly output three distinct values: the point forecast ŷ, the upper bound û, and the lower bound ˆl. Conversely, the IPIV method strictly utilizes its proposed architectural design, which constructs the point forecast via a convex combination of the bounds. While all loss functions, ŷ, ˆl, û, PICP and widths are functions of the model parameters θ and regressors x, we omit these arguments hereafter to simplify the notation. The following benchmark methods have been adapted for our evaluation: 1. Quantile regression (QR): Quantile regression directly estimates target quantiles by minimizing the pinball function defined as ρτ (r) = max{τ r, (τ − 1)r} for a given quantile τ . To produce a PI that covers the target probability of PN 1 − α, we use the pinball loss consisting of two quantile terms: LQR = N1 i=1 ρα/2 yi − ˆl) + ρ1−α/2 yi − û). 2. Quality-driven loss (QD+) [PBZN18]: Building upon the standard quality-driven loss, QD+ incorporates point forecasting and boundary constraints with the loss function consisting of four terms: the first two represent the original QD loss that is served for width penalty and PICP, while the latter two terms signify the MSE loss and penalty on the forcing the ŷ to lie within the PI. 2 LQD+ = (1 − λ1 )(1 − λ2 )MPIWcapt. + λ1 (1 − λ2 ) max 0, (1 − α) − PICP + λ2 LMSE + ϵLPenalty . (18) where the MPIWcapt. is the mean PI width penalty that measures only the samples that y lies within the PI, and PN LPenalty = N1 i=1 max(0, ˆli − ŷi ) + max(0, ŷi − ûi ) penalizes violations of ŷ ∈ [ˆl, û]. The hyperparameter λ1 ∈ (0, 1) balances the PI width and PICP, while λ2 ∈ (0, 1) determines the weight of the point forecast error (LMSE ) relative to the PI metrics. The hyperparameter ϵ scales the penalty for constraint violations. A notable drawback of this approach is the necessity to tune λ1 , λ2 , and ϵ set by the user’s preference. Furthermore, the soft penalization does not strictly guarantee that the point forecast ŷ will remain bounded within the interval bounds. 3. Integrated prediction interval and value predictions (IPIV) [SKR22]: The method produces three outputs: ŷ, ˆl, v ∈ (0, 1) and constructs the point forecast via a convex combination: ŷ = v û + (1 − v)ˆl which helps guarantee that ŷ always lies within the PI, while positioning ŷ asymmetrically. The loss is a weighted sum of the PI and point forecast objectives: LIPIV = βLPI + (1 − β)LPoint with a weight β (default 0.5) where LPI = MPIWcapt. +

√

N 2 1 X ℓ vi ûi + (1 − vi )ˆli , yi N λ max 0, (1 − α) − PICP , LPoint = N i=1

(19)

The hyperparameter λ governs the balance between PICP and interval width. The point forecast loss utilizes a regression loss ℓ. To better capture epistemic model uncertainty, IPIV further incorporates an ensemble strategy across multiple networks and a z-score calibration step. 12

4. Enhanced multi-quantile loss (EMQ) [SLZA25]: The method addresses three main challenges: poor estimation of extreme tail quantiles, quantile crossing, and suboptimal balance between PI width and the target 1 − α coverage: LEMQ = Ladj. pinball + λLorder + Lcovarage + MPIWscaled .

(20)

To improve tail estimation, Ladj. pinball aggregates pinball losses across a 3-neighbor set around α/2 and 1 − α/2 quantiles, plus the median, reducing sensitivity to any single extreme quantile estimate. The Lorder term penalizes any pair of adjacent quantiles violations via an exp(max(0, ·)) function. The term Lcovarage term penalizes samples yi falling outside the PI, while MPIWscaled weights each width sample proportionally to the coverage deficiency (smaller weight) when the nominal level 1 − α is not met. The performance results are summarized in Table 1 and visually depicted in Figure 7. The evaluation is restricted to daytime data to avoid an overoptimistic evaluation as covering trivial zero-irradiance at nighttime is effortless. Prediction interval coverage probability (PICP): Table 1 clearly illustrates a major strength of the SolarPointPI loss that consistently maintains the PICP above the 0.9 nominal level across every time step. Conversely, Figure 7 demonstrates how other loss functions suffer from severe coverage degradation. QD+ and IPIV systematically fail to reach the nominal level (yielding values as low as 0.768). Furthermore, while QR initially meets the requirement, its trajectory drops below the nominal line as the forecasting step extends. Meanwhile, EMQ maintains coverage but tends to over-cover (e.g., 0.936 at the 15-minute horizon), which is compromised to unnecessarily wide intervals. We observed that EMQ and IPIV achieve the target PICP on the test set if nighttime data were included in a separate evaluation setting. Conversely, QD+ and QR fail to reach the target coverage even with nighttime data; for QD+, this is likely due to its high sensitivity to hyperparameters, which were set according to the original study. Interval sharpness (PINAW and PINALW): Among the valid models that pass PICP criterion, the SolarPointPI loss exhibits the lowest PINAW and PINALW score across all prediction steps; see Figure 7. Notably, maintaining a significant lower PINALW at longer horizons is non-trivial, as interval widths naturally expand with increasing lead times. This performance is driven by the eﬀiciency of the Sum-k width penalty in the proposed PI loss. In contrast, benchmark methods like EMQ and QR (where valid) suffer from significantly inflated PINAW and PINALW scores (reaching PINALW of 0.528 at the 4-hour horizon). Although QD+ and IPIV report low PINAW in Table 1, this is solely because they fail the 2 PICP requirement. When unnormalizing PINAW to the solar irradiance unit W/m , the average widths range in 196-296 2 2 W/m across from 15-minute to 4-hour horizons, and the average large widths range in 291-390 W/m . Winkler score: While the Winkler score provides a convenient single-index summary of both interval sharpness and quantile matching, it can be deceptive if viewed in isolation. As shown in Table 1, our SolarPointPI maintains a competitive Winkler score, even though EMQ and QR (where valid) achieve better values. However, they suffer from heavily inflated absolute widths to achieve this. Point forecast accuracy (MAE, RMSE, and MBE): The SolarPointPI achieves the lowest MAE across all prediction steps and the second lowest RMSE to QD+ (which is not valid by PICP criterion), while MBE can be relatively high. Conversely, IPIV and QD+ struggle with high residual errors, highlighting the diﬀiculty of balancing point and interval objectives. Meanwhile, valid models like EMQ exhibit visibly increasing point errors as the prediction step extends. In conclusion, the SolarPointPI loss not only provides rigorous uncertainty quantification that maintains target coverage while reducing PI width, but also ensures the point forecast remains acceptably accurate. These outcomes stem from three components: the log-barrier penalty for reliability, the Sum-k width penalty for sharpness, and the MGDA algorithm for balancing point and PI losses. Forecasting result time series: Figure 8 presents the time series plots for solar irradiance forecasts at 15-minute and 4-hour horizons under cloudy sky conditions. The results reveal a clear trade-off: while QR and EMQ maintain coverage across all horizons, they do so through excessively wide PIs. Conversely, QD+, IPIV, and SolarPointPI generate tighter PIs that closely track y; however, QD+ and IPIV fail to meet PICP requirements. During high volatility on the third day, SolarPointPI successfully captures the data with narrow intervals at the 15-minute horizon, though coverage slightly degrades at the 4-hour mark. Overall, SolarPointPI provides an optimal balance between reliability and interval sharpness.

Table 1: Performance of benchmarking loss functions evaluated on daytime data (06:00 - 18:00) in the test set. PICP < 0.9 is highlighted in red. Best metrics are bolded (restricted to models with PICP ≥ 0.9). Optimal scores are the value closest to zero for MBE, and the minimum value for all other metrics. PINAW and PINALW are in %, while MAE, RMSE, 2 and MBE are in W/m . The proposed method is SolarPointPI. Objective

Metrics PICP

PINAW

PINALW

Winkler

MAE

RMSE

MBE

0.306 0.468 0.477 0.309 0.357

48.11 49.45 52.14 48.65 48.05

88.19 87.06 89.48 88.78 87.90

0.15 1.07 3.21 5.25 1.13

0.408 0.598 0.580 0.408 0.471

66.18 67.94 70.69 66.80 65.39

110.14 108.57 111.40 111.41 110.46

-1.03 -1.39 4.09 7.99 4.45

0.440 0.633 0.616 0.440 0.506

69.79 71.85 74.37 70.63 69.68

114.42 112.43 114.93 115.42 115.33

4.83 1.73 1.39 7.54 7.75

0.462 0.654 0.632 0.461 0.524

72.01 74.79 76.38 72.61 71.66

116.84 114.37 116.18 117.19 116.47

8.87 -1.11 0.42 8.39 7.08

15-minute ahead QR QD+ IPIV EMQ SolarPointPI

0.899 0.786 0.822 0.936 0.913

22.50 13.08 15.71 24.73 21.63

39.01 19.35 22.82 41.50 32.09 1-hour ahead

QR QD+ IPIV EMQ SolarPointPI

0.892 0.778 0.816 0.920 0.911

30.53 18.73 21.74 31.92 29.52

48.99 26.15 30.10 50.48 41.03 2-hour ahead

QR QD+ IPIV EMQ SolarPointPI

0.907 0.769 0.811 0.919 0.911

33.62 19.42 22.40 34.25 31.76

51.64 26.44 30.06 52.25 42.18 4-hour ahead

QR QD+ IPIV EMQ SolarPointPI

0.907 0.768 0.812 0.916 0.906

35.27 19.95 22.92 35.48 31.95

52.89 26.42 29.89 52.84 41.59

PICP

PINAW PINAW (%)

0.925 0.900 PICP drops from nominal

PICP

0.875

35 30 25 20 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Prediction step

0.850

PINALW

0.825 PINALW (%)

0.800 0.775

40 30 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Prediction step

Winkler

MAE MAE (W/sqm)

Winkler

0.6 0.5 0.4 0.3

60 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Prediction step

RMSE 110

MBE (W/sqm)

RMSE (W/sqm)

100 90 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Prediction step QR

QD+

IPIV

MBE

10.0 7.5 5.0 2.5 0.0

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Prediction step

EMQ

SolarPointPI (Proposed)

Nominal 0.90 target

Figure 7: Overall performance metrics in objective function comparison.

15-min ahead prediction QR W/sqm

1000 500 0 00:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00

QD+

W/sqm

1000 500

0 00:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00

IPIV W/sqm

1000 500 0 00:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00

W/sqm

EMQ 1000 500 0 00:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00

SolarPointPI W/sqm

1000 500 0 00:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00

Point forecast (y)

Upper/lower bounds (u, l)

Prediction interval (PI)

Ground truth (y)

4-hour ahead prediction QR W/sqm

1000 500 0 00:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00

QD+ W/sqm

1000 500 0 00:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00

IPIV W/sqm

1000 500 0 00:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00

EMQ W/sqm

1000 500 0 00:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00

SolarPointPI W/sqm

1000 500 0 00:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00 06:00 12:00 18:00 24:00

Point forecast (y)

Upper/lower bounds (u, l)

Prediction interval (PI)

Ground truth (y)

Figure 8: Solar irradiance forecasts (W/sqm) comparing loss functions at 15 minute (top) and 4-hour (bottom) horizons under cloudy-sky conditions.

4.3 Model comparison We compare our proposed architecture against state-of-the-art models. These include a zero-shot foundation model Chronos 2.0, LightGBM, LSTM and Transformer models configured in an encoder-decoder architecture. To ensure consistency across all evaluations, every model is configured to release a point forecast and PI with a confidence level of 0.9. Regarding the input window, all models except LightGBM use a 2-day lag, while LightGBM uses a 4-hour lag; both settings target a 4-hour forecasting horizon. All models except the Chronos and LightGBM are trained according to the proposed loss and Algorithm 1. 1. Chronos (Zero-shot) [ASK+ 25]: Chronos-2 is a pretrained foundation model based on an encoder-only transformer architecture, designed for zero-shot generalization across universal time series forecasting tasks. It introduces a group attention mechanism that facilitates in-context learning by sharing information across related time series, allowing the model to handle univariate, multivariate, and covariate-informed tasks without the need for task-specific training. 2. LightGBM: LGBM employs a leaf-wise growth strategy and built-in feature selection for modeling eﬀiciency. To generate multi-step forecasts, we train 16 independent models, each optimized via pinball loss for the median (the 0.5-quantile) and PI bounds (the 0.05- and 0.95-quantiles). Hyperparameter tuning (learning rate, depth, and leaf constraints) is performed on the 1-step and 9-step models with the resulting parameters shared across the 1-8 and 9-16 horizons, respectively. 3. LSTM: A recurrent neural network architecture designed for sequence-to-sequence modeling. As shown in Figure 9a, the LSTM encoder processes historical inputs (lag regressors) and passes its final hidden state to an LSTM decoder. The decoder recursively incorporates future regressors to generate step-wise hidden states. These states are passed through an output head consisting of linear, GELU (Gaussian Error Linear Unit), and dropout layers. The final layer outputs a central point estimate alongside softplus-activated width adjustments (∆ˆli , ∆ûi ) to strictly formulate the lower and upper interval bounds. 4. Transformer: An architecture utilizing attention mechanisms to capture long-range dependencies. As shown in Figure 9b, the encoder processes positionally-encoded lag regressors via multi-head self-attention. The decoder utilizes multi-head cross-attention to integrate the encoder’s continuous representations with positionally-encoded future regressors. The aggregated features are routed through the output head (same as the one used in LSTM) to yield ŷ, ∆ˆli , ∆ûi . 5. Proposed SolarPointPI: Our primary proposed LSTM architecture. As shown in Figure 9c, the core LSTM network processes the sequence of lag regressors, and the final state is normalized and activated via batch normalization (BN) and ReLU layers. Rather than using a sequential decoder, this encoded representation is concatenated step-wise with the respective future regressors. The combined vector is processed by step-specific independent submodels (consisting of stacked Linear, BN, and ReLU blocks) acting as specialized heads. Finally, these submodels output the point forecast and boundary adjustments (∆ˆli , ∆ûi ) to construct the prediction interval. 6. Chronos as feature generators: We augment the Chronos foundation model as a feature extractor. Chronos processes the lag regressors to generate auxiliary future representations as the 0.05, 0.5, and 0.95 quantile outputs. These extracted features are then concatenated with the future regressors before being fed into each of the LSTM, Transformer, and SolarPointPI models. • Chronos + LSTM: Chronos outputs are applied into the decoder of the LSTM. • Chronos + Transformer: The Chronos outputs are fed into the decoder module of the Transformer. • Chronos + SolarPointPI: The Chronos outputs are concatenated with the outputs of the shared model Mc and future regressors before being fed into the k-step submodels.

(a) LSTM.

(b) Transformer.

Figure 9: Benchmarking model architectures.

Table 2: Model and algorithm hyperparameters for model comparison experiment. Model

Hyperparameters

LGBM

Submodel 1-8: learning rate: 8.8e-3, num leaves: 454, max depth: 14 min data in leaf: 477, feature fraction: 0.63, bagging fraction: 0.58, min gain to split: 0.77 Submodel 9-16: learning rate: 1.1e-2, num leaves: 456, max depth: 14 min data in leaf: 11, feature fraction: 0.98, bagging fraction: 0.76, min gain to split: 0.86 hidden size: 64, num encoder/decoder layer: 2,2, batch size: 32678, learning rate: 1e-4, optimizer: AdamW d model: 64, num encoder/decoder layer: 3,3, num heads: 8, feed forward dim: 128, batch size: 16384 (Gradient accumulation), learning rate: 3e-4, optimizer: AdamW num lstm cell: 1, lstm hidden size: 70, submodel neuron: [100,100], batch size: 32678, learning rate: 3e-4, optimizer: AdamW

LSTM Transformer SolarPointPI

PICP

0.90 0.88

PICP drops from nominal

0.86

PINALW (%)

PICP

0.92

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Prediction step

80 70 60 50 40 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Prediction step

MBE 10

100

MBE (W/sqm)

MAE (W/sqm)

MAE

PINALW

80 60

5 0 5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Prediction step Chronos (Zero-shot) LGBM

LSTM Chronos + LSTM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Prediction step

Transformer Chronos + Transformer

SolarPointPI (Proposed) Chronos + SolarPointPI

Nominal 0.90 target

Figure 10: Performance metrics comparison for model architecture benchmarking The performance evaluation metrics are reported in Table 3 and plotted across all 16 forecasting steps in Figure 10. Overall performance. As illustrated in Figure 10, the Chronos and LGBM models frequently yield PICP values below the 0.9 target coverage while producing relatively large PINALW. In contrast, remaining models trained with the proposed loss function satisfy the coverage requirement while achieving comparable PINALW and MAE. Although Chronos (a zero-shot forecaster that has not seen this dataset) provides moderately accurate point forecasts at the 15-minute horizon, its MAE increases sharply as the horizon extends, accompanied by a high MBE. While the LSTM and its Chronos-augment achieve an MBE near zero, SolarPointPI exhibits a slight residual bias. Since all deep learning models satisfy the coverage target, we compare their PI widths in Table 3 and more closely. SolarPointPI variants achieve the lowest PINALW at the 15-minute horizon. At longer horizons, the Chronos-augmented LSTM yields the narrowest intervals, followed by the Transformer variants, which consistently rank as the second best across all lead times. PINAW results follow a similar pattern: SolarPointPI is superior for the 15-minute and 1-hour lead times, while LSTM-based models show improved results for extended horizons. Regarding the Winkler score, MAE, and RMSE, the Chronos-augmented versions of LSTM, Transformer, and SolarPointPI alternately achieve the best performance across different forecasting horizons. While the metrics remain relatively competitive between these models, Chronos+SolarPointPI tends to yield the lowest scores at shorter horizons We find that attention mechanisms do not significantly outperform LSTM, supporting claims that Transformers struggle with strict temporal ordering [ZCZX23] and inter-variate impact [CCB25]. Although variants of Transformer and recurrent structure as in [KOPJ24, CJL+ 25] exist for solar power, this experiment focuses on baseline comparisons rather than architectural optimization.

Table 3: Performance metrics of benchmarking models evaluated on daytime data (06:00 - 18:00) in the test set. PICP < 0.9 is highlighted in red. Best metrics are bolded (restricted to models with PICP ≥ 0.9). Optimal scores are the value closest to zero for MBE, and the minimum value for all other metrics. PINAW and PINALW are in %, while MAE, RMSE, 2 and MBE are in W/m . The proposed method is SolarPointPI. Model

Metrics PICP

PINAW

PINALW

Winkler

MAE

RMSE

MBE

0.376 0.316 0.393 0.364 0.374 0.368 0.351 0.352

55.59 47.32 49.00 48.14 52.54 51.63 46.74 46.56

97.37 87.51 88.14 86.43 88.02 87.04 87.18 86.32

12.39 3.13 3.30 -1.32 -7.10 0.16 4.90 3.33

0.563 0.429 0.497 0.464 0.485 0.461 0.475 0.462

85.90 66.20 66.60 64.04 68.42 66.40 65.52 63.88

141.27 112.56 111.33 107.67 111.26 108.09 112.61 108.99

9.65 6.07 -0.94 1.12 -2.23 3.28 8.02 6.69

0.646 0.466 0.523 0.485 0.507 0.481 0.509 0.485

101.25 70.15 70.04 67.02 71.49 69.32 69.32 67.43

161.52 117.20 115.27 111.37 114.83 111.64 116.98 113.11

10.07 7.19 0.27 2.10 2.19 4.12 9.66 7.95

0.712 0.479 0.537 0.493 0.518 0.485 0.525 0.499

113.26 71.97 72.37 68.45 73.62 70.60 71.42 69.03

176.93 118.78 117.01 112.82 116.67 112.24 118.66 114.99

9.35 6.63 1.99 2.41 3.41 3.83 8.71 8.22

15-minute ahead Chronos (zero-shot) LGBM LSTM Chronos+LSTM Transformer Chronos+Transformer SolarPointPI Chronos+SolarPointPI

0.877 0.898 0.918 0.916 0.923 0.921 0.913 0.913

24.00 23.18 24.75 22.31 24.22 24.00 20.81 20.99

39.15 38.15 32.46 30.47 31.88 32.08 30.55 30.37

1-hour ahead Chronos (zero-shot) LGBM LSTM Chronos+LSTM Transformer Chronos+Transformer SolarPointPI Chronos+SolarPointPI

0.846 0.893 0.911 0.909 0.909 0.911 0.908 0.908

33.08 31.69 30.58 28.15 29.38 28.96 28.78 28.00

54.31 48.20 39.23 37.34 37.69 38.32 39.32 38.49

2-hour ahead Chronos (zero-shot) LGBM LSTM Chronos+LSTM Transformer Chronos+Transformer SolarPointPI Chronos+SolarPointPI

0.851 0.890 0.908 0.910 0.909 0.910 0.908 0.907

40.40 33.88 31.22 29.22 30.70 29.87 30.37 29.40

65.82 50.13 39.45 38.30 38.91 39.19 40.14 39.77

4-hour ahead Chronos (zero-shot) LGBM LSTM Chronos+LSTM Transformer Chronos+Transformer SolarPointPI Chronos+SolarPointPI

0.884 0.889 0.902 0.910 0.909 0.914 0.908 0.908

50.36 34.97 30.97 29.76 31.49 30.96 31.41 30.36

80.46 50.93 38.58 38.89 39.49 40.28 40.71 40.46

PICP PICP

0.92 0.91 0.90

PINALW (%)

0.89

PINALW

40 30 20 10 0

MAE

MAE (W/sqm)

80 70 60 50 40

MBE

MBE (W/sqm)

10 5 0 5 15min LSTM Chronos + LSTM

30min

Transformer Chronos + Transformer

SolarPointPI (Proposed) Chronos + SolarPointPI

4h Nominal 0.90 target

Figure 11: Performance comparison of base models with and without Chronos 2.0 Effect of Chronos 2.0 feature augmentation. Figure 11 illustrates the impact of utilizing Chronos 2.0 as a feature generator across the three base architectures. All model variants maintain a valid PICP above the 0.90 target level. Notably, augmenting SolarPointPI with Chronos features reduces the PINALW, MAE, and MBE across all prediction steps, suggesting that Chronos features facilitate narrower prediction intervals and more accurate point forecasts without compromising reliability. For the LSTM model, the inclusion of Chronos leads to a significant reduction in all metrics across nearly all horizons, with the exception of the 4-hour horizon, where the PINALW and MBE experience slight degradation. Conversely, the benefit of Chronos for the Transformer architecture is primarily reflected in the MAE and MBE; unexpectedly, the PINALW for the Chronos+Transformer configuration increases relative to the base model. Ultimately, whether the improvement gained by adding Chronos is statistically or practically significant depends heavily on the specific base architecture, the choice of metrics, and the forecast horizons. Consequently, selecting a superior model may require considering additional factors, such as computational complexity and numerical implementation requirements.

15-min ahead prediction W/sqm

Chronos (Zero-shot)

4-hour ahead prediction

LGBM

Chronos (Zero-shot)

LGBM

1000 500 0 06:00 12:00 06:00 12:00 18:00 06:00 12:00 06:00 12:00 18:00 06:00 12:00 06:00 12:00 18:00 06:00 12:00 06:00 12:00 18:00

W/sqm

1000

LSTM

Chronos + LSTM

LSTM

Chronos + LSTM

500 0 06:00 12:00 06:00 12:00 18:00 06:00 12:00 06:00 12:00 18:00 06:00 12:00 06:00 12:00 18:00 06:00 12:00 06:00 12:00 18:00

W/sqm

1000

Transformer

Chronos + Transformer

Transformer

Chronos + Transformer

500 0 06:00 12:00 06:00 12:00 18:00 06:00 12:00 06:00 12:00 18:00 06:00 12:00 06:00 12:00 18:00 06:00 12:00 06:00 12:00 18:00

W/sqm

1000

SolarPointPI

Chronos + SolarPointPI

SolarPointPI

Chronos + SolarPointPI

500 0 06:00 12:00 06:00 12:00 18:00 06:00 12:00 06:00 12:00 18:00 06:00 12:00 06:00 12:00 18:00 06:00 12:00 06:00 12:00 18:00

Point forecast (y)

Upper/lower bounds (u, l)

Prediction interval (PI)

Ground truth (y)

Figure 12: Solar irradiance forecasts (W/sqm) comparing model architectures at 15-minute and 4-hour prediction horizons under cloudy-sky conditions. Example of forecasting result. Figure 12 illustrates the performance of various models at 15-minute and 4-hour horizons under cloudy sky conditions. The zero-shot Chronos model exhibits a very high uncertainty estimate, resulting in extremely wide PIs at both horizons. Similarly, the LGBM model produces relatively large intervals to capture sudden irradiance spikes. In contrast, the LSTM, Transformer, and SolarPointPI models show a tighter fit to the ground truth, with visually comparable PI widths and point forecast accuracy. All methods produce wider PIs at the 4-hour horizon, as uncertainty naturally increases as the prediction horizon extends.

4.4 Computational complexity All benchmarks were conducted on an Ubuntu 24.04 LTS system equipped with an Intel Core i7-14700 processor (20 cores, 28 threads), an NVIDIA GeForce RTX 5070 Ti GPU with 16 GB of VRAM, and 64 GB of RAM, using Python 3.13.11 within PyTorch implementation [PGM+ 19]. To evaluate model complexity, the ptflops library [Sov24] was used to calculate the Multiply-Accumulate operations (MACs) and total parameter count, with the resulting MACs multiplied by approximately two to estimate total FLOPs. Inference latency was benchmarked over 100 iterations following a 10-iteration warmup phase to stabilize hardware caches. All inference measurements were performed on the test set with a batch size of 1, reflecting the practical deployment scenario for solar forecasting, where the model predicts a single sample at a time. As shown in Figure 13 and Table 4, the computational complexity analysis demonstrates the eﬀiciency of the proposed SolarPointPI model benchmarking with other model architectures. SolarPointPI achieves the lowest FLOP count at 0.642 MFLOPs. The LSTM records a higher count of 1.708 MFLOPs, due to its encoder-decoder structure that requires two sequential passes through the network. SolarPointPI avoids this by using a single-pass architecture, and its lightweight output heads further contribute to the reduction in arithmetic operations. The Transformer has the highest FLOP count among the compact models at 4.370 MFLOPs, a result of its O(n2 ) attention complexity where n is the sequence length. This quadratic scaling becomes particularly costly at 23

0.31

0.25

0.8

Table 4: Computational complexity comparison.

0.64

0.4

0.79

0.6 1.71 0.10 0.65

Metric value

1.0

1.59

4.37

1.2

0.2 0.0

LSTM MFLOPs

Transformer Parameters (M)

Model

MFLOPs

Parameters (M)

Latency (ms)

Chronos (zero-shot) LSTM Transformer SolarPointPI

3246.600

119.478

7.323

1.708 4.370 0.642

0.104 0.254 0.306

0.649 1.590 0.789

SolarPointPI Latency (ms)

Figure 13: Computational complexity. batch size 1, where the overhead of the attention mechanism is not distributed across samples. SolarPointPI has more parameter counts but fewer FLOPs than LSTM because it pairs a single LSTM with separate, submodules for each output step. In contrast, the LSTM (encoder-decoder) reduces parameters by sharing a single output head, but its dual-LSTM encoder-decoder architecture significantly increases overall computational complexity. In terms of inference latency, SolarPointPI completes a forecast in 0.789 ms, which is comparable to the LSTM at 0.649 ms and substantially lower than the Transformer at 1.590 ms. These results confirm that the eﬀiciency gains from the single-pass design translate directly to reduced runtime. The inference time increases significantly when considering the Chronos (zero-shot) baseline, which requires 3,246.6 MFLOPs and 119.478 M parameters. In Chronos-augmented variants, inference time is dominated by this pretrained backbone; even with frozen weights, the full parameter set must be evaluated at each step. Nevertheless, total inference time remains in the millisecond range, making real-time deployment feasible. Given the resulting performance gains, the additional latency introduced by Chronos is a justifiable trade-off.

5 Conclusions This paper demonstrates the effectiveness of a unified NN framework for simultaneous point and interval forecasting, driven by a novel PI loss and a multi-objective training strategy using multiple gradient descent directions. We successfully addressed the fundamental trade-off between target coverage (reliability) and interval width (sharpness) by employing an extended log-barrier function for the former and the Sum-k loss for the latter. Experimental results in intra-day solar irradiance forecasting–a domain characterized by high fluctuation–confirm that our method consistently achieves the target PICP with superior sharpness compared to existing scalarized loss functions. Furthermore, our approach eliminates the intensive hyperparameter tuning typically required in literature. The framework’s versatility was validated by integrating it with state-of-the-art architectures, including LSTM encoder-decoders, Transformers, and Chronos-augmented models, showing competitive performance across all configurations. While this work focuses on prediction intervals, the framework is architecturally flexible and can be adapted to different input features or more complex models depending on data availability for real-time deployment. Future research will explore extending this multi-objective optimization approach to other forms of uncertainty quantification, such as cumulative distribution functions (CDFs) or multi-quantile estimation, to support a broader range of risk-based decision-making tasks.

6 Acknowledgment This research project is financially supported by the Ratchadaphiseksomphot Endowment Fund, Chulalongkorn University. During the preparation of this work the authors used Gemini 3 Flash in order to improve readability and language. After using this tool, the authors reviewed and edited the content as needed and took full responsibility for the content of the published article.

References [ASK+ 25]

Abdul Fatir Ansari, Oleksandr Shchur, Jaris Küken, Andreas Auer, Boran Han, Pedro Mercado, Syama Sundar Rangapuram, Huibin Shen, Lorenzo Stella, Xiyuan Zhang, Mononito Goswami, Shubham Kapoor, Danielle C. Maddix, Pablo Guerron, Tony Hu, Junming Yin, Nick Erickson, Prateek Mutalik Desai, Hao Wang, Huzefa Rangwala, George Karypis, Yuyang Wang, and Michael Bohlke-Schneider. Chronos-2: From univariate to universal forecasting, 2025.

[AWS24]

Worachit Amnuaypongsa, Wijarn Wangdee, and Jitkomut Songsiri. Probabilistic solar power forecasting using multi-objective quantile regression. In 2024 18th International Conference on Probabilistic Methods Applied to Power Systems (PMAPS), pages 1–6. IEEE, 2024.

[AWS25]

Worachit Amnuaypongsa, Wijarn Wangdee, and Jitkomut Songsiri. Neural network-based prediction interval estimation with large width penalization for renewable energy forecasting and system applications. Energy Conversion and Management: X, 27:101119, 2025.

[CCB25]

Yu Chen, Nathalia Céspedes, and Payam Barnaghi. A closer look at transformers for time series forecasting: Understanding why they work and where they struggle. In Forty-second International Conference on Machine Learning, 2025.

[CJL+ 25]

Shanlin Chen, Tao Jing, Mengying Li, Hiu Hung Lee, Ming Chun Lam, and Siqi Bu. Interpretable transformer based intra-day solar forecasting with spatiotemporal satellite and numerical weather prediction inputs. Energy and AI, page 100667, 2025.

[CLZ+ 25]

Weiyu Chen, Baijiong Lin, Xiaoyuan Zhang, Xi Lin, Han Zhao, Qingfu Zhang, and James T Kwok. Gradient-based multi-objective deep learning: Algorithms, theories, applications, and beyond. arXiv preprint arXiv:2501.10945, 2025.

[CUE23]

CUEE Satellite Receiving Station. Himawari-8 satellite RGB image archives, 2023. Accessed: 2024. Available: https://himawari.optemis.space/archives.

[Dés12]

Jean-Antoine Désidéri. Multiple-gradient descent algorithm (MGDA) for multiobjective optimization. Comptes Rendus. Mathématique, 350(5-6):313–318, 2012.

[EMMS12] Pavel V Etingov, Jian Ma, Yuri V Makarov, and Krishnappa Subbarao. Online Analysis of Wind and Solar Part I: Ramping Tool. Technical report, Pacific Northwest National Lab. (PNNL), Richland, WA (United States), 01 2012. [Ene24]

EnergyCUEE. Improved clear-sky model, 2024. Available: https://github.com/energyCUEE/improved _clearsky_model.

[FKF18]

Seyyed A Fatemi, Anthony Kuh, and Matthias Fripp. Parametric methods for probabilistic forecasting of solar irradiance. Renewable Energy, 129:666–676, 2018.

[GLS23]

Tilmann Gneiting, Sebastian Lerch, and Benedikt Schulz. Probabilistic solar forecasting: Benchmarks, postprocessing, verification. Solar Energy, 252:72–80, 2023.

[GZK+ 24]

Oktoviano Gandhi, Wenjie Zhang, Dhivya Sampath Kumar, Carlos D Rodríguez-Gallegos, Gokhan Mert Yagli, Dazhi Yang, Thomas Reindl, and Dipti Srinivasan. The value of solar forecasts and the cost of their errors: A review. Renewable and Sustainable Energy Reviews, 189:113915, 2024.

[HZW+ 20] Chao Huang, Zhenyu Zhao, Long Wang, Zijun Zhang, and Xiong Luo. Point and interval forecasting of solar irradiance with an active gaussian process. IET Renewable Power Generation, 14(6):1020–1030, 2020. [IP02]

Pierre Ineichen and Richard Perez. A new airmass independent formulation for the Linke turbidity coeﬀicient. Solar Energy, 73(3):151–157, 2002.

[KDY+ 22]

Hoel Kervadec, Jose Dolz, Jing Yuan, Christian Desrosiers, Eric Granger, and Ismail Ben Ayed. Constrained deep networks: Lagrangian optimization via log-barrier extensions. In 2022 30th European Signal Processing Conference (EUSIPCO), pages 962–966. IEEE, 2022.

[KJL20]

Danijel Kivaranovic, Kory D. Johnson, and Hannes Leeb. Adaptive, distribution-free prediction intervals for deep networks. In Silvia Chiappa and Roberto Calandra, editors, Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108, pages 4346–4356. PMLR, 26–28 Aug 2020.

[KNCA11]

Abbas Khosravi, Saeid Nahavandi, Doug Creighton, and Amir F. Atiya. Lower upper bound estimation method for construction of neural network-based prediction intervals. IEEE Transactions on Neural Networks, 22:337–346, 3 2011.

[Koe17]

Roger Koenker. Quantile regression: 40 years on. Annual review of economics, 9:155–176, 2017.

[KOPJ24]

Jimin Kim, Josue Obregon, Hoonseok Park, and Jae-Yoon Jung. Multi-step photovoltaic power forecasting using transformer and recurrent neural networks. Renewable and Sustainable Energy Reviews, 200:114479, 2024.

[LMC+ 22]

Xianglong Li, Longfei Ma, Ping Chen, Hui Xu, Qijing Xing, Jiahui Yan, Siyue Lu, Haohao Fan, Lei Yang, and Yongqiang Cheng. Probabilistic solar irradiance forecasting based on XGBoost. Energy Reports, 8:1087–1095, 2022.

[LRJ23]

Haochuan Li, Alexander Rakhlin, and Ali Jadbabaie. Convergence of Adam under relaxed assumptions. Advances in Neural Information Processing Systems, 36:52166–52196, 2023.

[LSH+ 22]

Yuandu Lai, Yucheng Shi, Yahong Han, Yunfeng Shao, Meiyu Qi, and Bingshuai Li. Exploring uncertainty in regression neural networks for construction of prediction intervals. Neurocomputing, 481:249–257, 4 2022.

[LZLZ22]

Siyan Liu, Pei Zhang, Dan Lu, and Guannan Zhang. PI3NN: Out-of-distribution-aware prediction intervals from three neural networks. In 10th International Conference on Learning Representations, ICLR 2022, 2022.

[Mie99]

Kaisa Miettinen. Nonlinear multiobjective optimization, volume 12. Springer Science & Business Media, 1999.

[MS23]

Giorgio Morales and John W. Sheppard. Dual accuracy-quality-driven neural network for prediction interval generation. IEEE Transactions on Neural Networks and Learning Systems, pages 1–11, 2023.

[PBZN18]

Tim Pearce, Alexandra Brintrup, Mohamed Zaki, and Andy Neely. High-quality prediction intervals for deep learning: A distribution-free, ensembled approach. In Proceedings of the 35th International Conference on Machine Learning, volume 80, pages 4075–4084. PMLR, 10–15 Jul 2018.

[PGM+ 19] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. arXiv, 2019. [PTSN+ 23] Quentin Paletta, Guillermo Terrén-Serrano, Yuhao Nie, Binghui Li, Jacob Bieker, Wenqi Zhang, Laurent Dubus, Soumyabrata Dev, and Cong Feng. Advances in solar forecasting: Computer vision with deep learning. Advances in Applied Energy, 11:100150, 2023. [QOB+ 17]

Z. Qu, A. Oumbe, P. Blanc, B. Espinar, G. Gesell, B. Gschwind, L. Klüser, M. Lefèvre, L. Saboret, M. Schroedter-Homscheidt, and L. Wald. Fast radiative transfer parameterisation for assessing the surface solar irradiance: The Heliosat-4 method. Meteorologische Zeitschrift, 26(1):33–57, 2017.

[SDS23]

Pardeep Singla, Manoj Duhan, and Sumit Saroha. A point and interval forecasting of solar irradiance using different decomposition based hybrid models. Earth Science Informatics, 16(3):2223–2240, 2023.

[SK18]

Ozan Sener and Vladlen Koltun. Multi-task learning as multi-objective optimization. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.

[SKR22]

Eli Simhayev, Gilad Katz, and Lior Rokach. Integrated prediction intervals and specific value predictions for regression problems using neural networks. Knowledge-Based Systems, 247:108685, 7 2022.

[SLG24]

Adnan Saeed, Chaoshun Li, and Zhenhao Gan. Short-term wind speed interval prediction using improved quality-driven loss based gated multi-scale convolutional sequence model. Energy, 300:131590, 8 2024.

[SLR+ 25]

Adnan Saeed, Chaoshun Li, Saeed Rubaiee, Mohd Danish, and Sana Anwar. Enhanced wind speed forecasting for sustainable power systems: A deep learning framework unifying deterministic predictions and uncertainty quantification. Energy, 335:137979, 10 2025.

[SLZA25]

Adnan Saeed, Chaoshun Li, Qiannan Zhu, and Belal Ahmad. Deterministic forecasts and prediction intervals for wind speed using enhanced multi-quantile loss based dilated causal convolutions. IEEE Transactions on Sustainable Energy, 16:2002–2014, 7 2025. 26

[Son25]

Jitkomut Songsiri. Solar energy forecasting in thailand: Concepts, modeling and case-study results, sep 2025. Department of Electrical Engineering, Faculty of Engineering, Chulalongkorn University, Available: http://jitkomut.eng.chula.ac.th/pdf/solarforecast2025.pdf.

[Sov24]

Vladislav Sovrasov. ptflops: a flops counting tool for neural networks in Pytorch framework, 2018-2024. Available: https://github.com/sovrasov/flops-counter.pytorch.

[SSP16]

Enrica Scolari, Fabrizio Sossan, and Mario Paolone. Irradiance prediction intervals for PV stochastic generation in microgrid applications. Solar Energy, 139:116–129, 2016.

[STT+ 25]

Suwichaya Suwanwimolkul, Natanon Tongamrak, Nuttamon Thungka, Naebboon Hoonchareon, and Jitkomut Songsiri. Deep-learning-based and near real-time solar irradiance map using Himawari-8 satellite imageries. Solar Energy, 288:113262, 2025.

[Win72]

Robert L. Winkler. A decision-theoretic approach to interval estimation. Journal of the American Statistical Association, 67(337):187–191, 1972.

[WYP+ 17] Huaizhi Wang, Haiyan Yi, Jianchun Peng, Guibin Wang, Yitao Liu, Hui Jiang, and Wenxin Liu. Deterministic and probabilistic forecasting of photovoltaic power based on deep convolutional neural network. Energy conversion and management, 153:409–422, 2017. [YvdMM20] Dazhi Yang, Dennis van der Meer, and Joakim Munkhammar. Probabilistic solar forecasting benchmarks on a standardized dataset at Folsom, California. Solar Energy, 206:628–639, 2020. [ZCZX23]

Ailing Zeng, Muxi Chen, Lei Zhang, and Qiang Xu. Are transformers effective for time series forecasting? In Proceedings of the AAAI conference on artificial intelligence, pages 11121–11128, 2023.

[ZWS21]

Changfei Zhao, Can Wan, and Yonghua Song. Operating reserve quantification using prediction intervals of wind power: An integrated probabilistic forecasting and decision methodology. IEEE Transactions on Power Systems, 36:3701–3714, 7 2021.

[ZWS22]

Changfei Zhao, Can Wan, and Yonghua Song. Cost-oriented prediction intervals: On bridging the gap between forecasting and decision. IEEE Transactions on Power Systems, 37:3048–3062, 7 2022.

[ZWW+ 15] Guoyong Zhang, Yonggang Wu, Kit Po Wong, Zhao Xu, Zhao Yang Dong, and Herbert Ho-Ching Iu. An advanced approach for construction of optimal wind power prediction intervals. IEEE Transactions on Power Systems, 30:2706–2715, 9 2015.

Record · ID 120512 · SHA-256 9ebd16c130e52adc
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.