The experimental dataset was obtained from the Vaulx-en-Velin IDMP station in France, located at 45°46′43′′ N and 4°55′21′′ E. It comprised diffuse horizontal illuminance, diffuse horizontal irradiance, zenith luminance, and atmospheric variables measured at 1 min resolution over the 2019-2024 period. For model development, only the 2024 records were used to train the ANFIS structure. The year 2024 was selected as the development year without applying any criterion based on model performance or data quality. Once the model had been developed using only the 2024 campaign, five additional complete annual measurement campaigns, corresponding to 2019-2023, were selected exclusively for external validation to assess its robustness to interannual variability. These datasets were kept completely separate from model training and parameter selection, ensuring an independent assessment of model performance across different years. Because each validation dataset covers a complete calendar year, the evaluation incorporates seasonal and meteorological variability, allowing model generalization to be assessed across the range of solar and sky conditions occurring throughout the year. Details of the measurement system, including the instrumentation, can be found in the documentation of the ENTPE CIE-IDMP station [14].

Prior to the analysis, the downloaded records were processed according to the CIE [15] quality-control recommendations. The dataset was further refined by discarding invalid samples and removing observations acquired when the solar altitude was lower than 5°.

The explanatory variables included in the correlation analysis of Table 1 were derived from the irradiance measurements and solar-position parameters according to Eqs. (1)-(5). First, the extraterrestrial normal irradiance was calculated from the solar constant corrected for the annual variation in the Sun-Earth distance:

I0=IscE0 (1)

where I_sc=1361 W/m² is the solar constant [16], and E0 is the eccentricity correction factor, obtained as shown in Eq. (2):

E0=1+0.033cos(2πdn365) (2)

with d_n representing the day number. The clearness index Kt was then computed from Eq. (3) as the ratio between the measured global horizontal irradiance and the corresponding extraterrestrial horizontal irradiance:

Kt=GHII0sin(α) (3)

where α is the solar altitude angle. In addition, two sky-condition descriptors proposed by Perez et al. [17] were considered. The sky clearness index was calculated from the diffuse and direct normal irradiance components and the solar zenith angle Z, expressed in radians, as shown in Eq. (4).

ε=(DHI+DNIDHI+1.041Z3)(1+1.041Z3) (4)

The sky brightness index, Δ, was obtained from Eq. (5), where m is the relative optical air mass.

Δ=mDHII0 (5)

As shown in Table 1, the correlation pattern for diffuse horizontal illuminance (LDFH) highlights the dominant role of the diffuse radiative component and daylight-related sky conditions. Among the analyzed predictors, DHI exhibits a strong positive association with LDFH, confirming that diffuse horizontal irradiance is the main radiometric predictor of diffuse horizontal illuminance. The sine of the solar altitude, sin(α), was also retained because it accounts for the geometric dependence of horizontal daylight quantities on the apparent position of the sun. In addition, zenith luminance was included as a final predictor because it provides direct photometric information on the brightness state of the sky vault, which is especially relevant for diffuse illuminance estimation. Although other variables, such as GHI, BHI, Kt, ε and Δ, among others, were analyzed to characterize the radiative and sky-condition context, they were not incorporated into the final LDFH model. Several of these variables are strongly linked to the same underlying irradiance balance or can be derived from the measured irradiance components, which could introduce redundant information into the ANFIS structure. Therefore, the final input set was restricted to DHI, sin(α), and zenith luminance. This selection preserves the main physical and photometric mechanisms controlling LDFH while keeping the model structure simple and reducing redundancy among predictors. The influence of the selected predictors on model performance is further assessed using reduced-input ANFIS configurations, as discussed in Section 4.

The predictive accuracy of the models was assessed using nRMSE (%), nMBE (%), and nMAE (%). These normalized indicators were adopted as the main error metrics to compare model behavior across the different datasets. Eqs. (6)-(8) define these statistics, where y_j is the measured value, ŷ_j is the corresponding model estimate, ȳ is the mean measured LDFH and n is the number of observations.

nRMSE(%)=100∑j=1n(ŷj−yj)2ny¯ (6)

nMBE(%)=100∑j=1n(ŷj−yj)∑j=1nyj (7)

nMAE(%)=1001n∑j=1n|ŷj−yj|y¯ (8)

Furthermore, the coefficient of determination, (R²), given in Eq. (9) and the Willmott’s index of agreement [18], defined in Eq. (10), were included as complementary goodness-of-fit indicators.

R2=1−∑j=1n(yj−ŷj)2∑j=1n(yj−y¯)2 (9)

IA=1−∑j=1n(ŷj−yj)2∑j=1n(|ŷj−y¯|+|yj−y¯|)2 (10)

The estimation of LDFH was carried out using a first-order Sugeno ANFIS. Model initialization was performed in Matlab using the Fuzzy Logic Toolbox™ [19]. Specifically, the initial FIS was obtained with the genfis function by applying grid partitioning to the selected input space. The model included three predictors, namely DHI, sin(α), and zenith luminance, and each input was represented by two Gaussian membership functions which were used to ensure continuous changes in membership degree across neighboring regions of the input space. For an input value x, the membership degree was calculated according to Eq. (11).

μ(x)=exp(−(x−c)22σ2) (11)

In the adopted ANFIS structure, the three-dimensional input domain was partitioned by assigning two fuzzy Gaussian sets to each predictor. The grid-based initialization then formed the rule base by combining the membership functions defined for DHI, sin(α), and zenith luminance. This led to 2³=8 first-order Sugeno rules, each one describing a different region of the input domain. The model was trained using all valid records from 2024 for 50 epochs with the hybrid method, which combines backpropagation and least-squares estimation. With this configuration, the model can account for the main changes in LDFH associated with diffuse irradiance, solar altitude, and sky brightness, while keeping the number of model coefficients within a controlled range. The general rule structure is expressed in Eq. (12). If DHI is A_i and sin(α) is B_j and zenith luminance is C_l, then

zk=pkDHI+qksin(α)+rkzenithluminance+sk (12)

Here, the coefficients pk, qk, rk and sk define the linear output function assigned to the k-th Sugeno rule, and zk denotes the corresponding rule output. The antecedent fuzzy subsets Ai, Bj and Cl are associated with DHI, sin(α) and zenith luminance, respectively. Once the membership degree of each input has been evaluated, the activation level of rule k, denoted by (ωk), is computed through the product operator, as indicated in Eq. (13).

ωk=μAi(DHI)μBj(sin(α))μCl(zenithluminance) (13)

Equation (14) defines the normalized firing strength of the k-th rule (ω¯k). This value is obtained by dividing the firing strength of rule (k) by the sum of the firing strengths computed for the eight rules generated from the possible combinations of the fuzzy sets indexed by (i), (j), and (l).

ω¯k=ωk∑m=18ωm (14)

The output of the first-order Sugeno FIS [20] was then obtained by aggregating the linear outputs of all rules according to their normalized firing strengths, as shown in Eq. (15). In this expression, (ŷ) represents the estimated diffuse horizontal illuminance, LDFH.

ŷ=∑k=18ω¯kzk (15)

After the initial Sugeno rule base had been generated, model training was carried out using ANFIS. The learning process adjusted both parts of the fuzzy model: the antecedent membership functions and the linear consequent equations. These parameters were estimated from the 2024 training records according to the adaptive neuro-fuzzy learning method introduced by Jang [21]. Negative output values are not physically admissible for diffuse horizontal illuminance. Therefore, when the ANFIS model produced a negative LDFH estimate, the value was set to zero before computing the final results, as indicated in Eq. (16).

LDFH=max{0,ŷANFIS,LDFH} (16)

Figure 2(a) shows that DHI has the most direct and nearly linear relationship with LDFH, confirming its role as the main radiometric predictor. However, the relationship between LDFH and sin(α) shown by Figure 2(b) presents a broader distribution, indicating that solar geometry modulates diffuse illuminance but does not explain it by itself. In turn, zenith luminance shows a clear positive association with LDFH, as shown in Figure 2(c), although with greater dispersion than DHI, suggesting that it provides additional photometric information on sky brightness. These patterns support the combined use of DHI, sin(α), and zenith luminance in the complete model. It should also be noted that, before selecting the final input structure, reduced ANFIS models based on subsets of the input variables were evaluated. The results obtained with these reduced configurations were poorer than those of the complete model, which included DHI, sin(α), and zenith luminance. Therefore, DHI, sin(α), and zenith luminance were retained as the final inputs, since this combination outperformed the simpler configurations evaluated in the preliminary comparison, as shown in the next section.

Figure 3 shows the comparison between the initial Gaussian membership functions and those obtained after ANFIS training for the three input variables. For DHI and zenith luminance, the initial and adjusted membership functions remain nearly unchanged after training, indicating that the initial grid-partition structure was already suitable for these two inputs. DHI shows a very direct association with LDFH, while zenith luminance provides a photometric description of the sky-brightness state. In contrast, the variable sin(α) shows a more evident adaptation of its membership functions. After training, the curves associated with this variable are shifted with respect to the initial configuration, reflecting the need to readjust the fuzzy representation of solar geometry.

Figure 4 compares the experimental and predicted LDFH values for the 2024 training dataset. The points are closely distributed around the 1:1 line, with a slope of 0.9806 and an R² of 0.98, indicating strong agreement between the ANFIS estimates and the measured values. A slightly larger dispersion is observed at high LDFH values, although the general trend is well reproduced.

After developing the Type-1 Sugeno ANFIS model, an Interval Type-2 Sugeno FIS [22], [23] was implemented to represent the uncertainty associated with LDFH estimation. For this purpose, the developed Interval Type-2 Sugeno model preserves the rule structure, the centers of the membership functions, and the linear consequents obtained from the Type-1 model. Uncertainty is introduced into the Gaussian antecedents through a footprint of uncertainty, FOU, bounded by a lower membership function, LMF, and an upper membership function, UMF [22], [23]. In this study, the FOU was incorporated by asymmetrically modifying the Gaussian widths, without changing the function centers or the rule base. For each nominal Gaussian function with width (σ), the lower and upper widths were defined as shown in Eqs. (17)-(18), where (u) controls the local amplitude of the FOU, and (β) introduces asymmetry between the lower and upper bounds.

σLMF=σ(1−u) (17)

σUMF=σ(1+βu) (18)

Prediction interval models should provide accurate central estimates while generating intervals that achieve the desired coverage probability without being excessively wide. These three properties are usually quantified through mean squared error (MSE), prediction interval coverage probability (PICP), and prediction interval normalized average width (PINAW), respectively [24], [25], [26]. Previous studies have used these metrics as part of the training loss. For example, [27] combined MSE, PINAW, and a squared penalty on the difference between empirical and target coverage in order to balance central prediction error, coverage, and interval width. The objective adopted here follows this general idea but changes the treatment of the residual error. Instead of using MSE, the residual contribution is formulated as a Student-t negative log-likelihood, whe­­­­­re the error is scaled by the uncertainty derived from the IT2-Sugeno interval [28], [29]. The remaining terms control the interval width and penalize failure to reach the target coverage only when the empirical coverage falls below the prescribed level. The Student-t negative log-likelihood is obtained by taking the negative logarithm of the Student-t predictive density, resulting in a residual-dependent logarithmic term together with the scale term and the Gamma-function normalization constants [30].

For each candidate parameter vector θ=[u1,...,uM,β,ν], where M is the total number of Gaussian membership functions in the antecedent layer, the IT2-Sugeno inference mechanism propagates the footprint of uncertainty through the rule base and returns the lower and upper bounds of the predicted response, using the iterative Karnik–Mendel type-reduction algorithm [22], [23]. The genetic algorithm (GA) searches this bounded parameter space to identify the combination of u, β and ν that minimizes the proposed Student-t-based objective function, where uj(j=1...M) is the relative uncertainty amplitude associated with the j-th Gaussian membership function β is the asymmetric upper-expansion factor of the footprint of uncertainty (FOU), and ν is the number of degrees of freedom of the Student-t distribution. The search bounds are set to um∈[0.001,0.95], β∈[1,3], and ν∈[3,50], with initial values of 0.15, 1.0, and 30, respectively. The upper limit ν=50 allows the Student-t distribution to approach near-Gaussian behavior, while preventing the GA from exploring an unnecessarily large region where additional increases in ν have little influence on the distribution. For each Gaussian antecedent membership function, the asymmetric FOU is obtained by modifying the nominal Type-1 Gaussian width σm, as defined in Eqs. (17)-(18). Hence, β=1 corresponds to the symmetric FOU, whereas 1]]>β>1 introduces a larger upper uncertainty expansion.

For each candidate solution θ, after type-reduction, the IT2-Sugeno model produces the prediction interval [Li,Ui], where Liand Uidenote the lower and upper interval bounds, respectively. Following the standard defuzzification of Interval Type-2 fuzzy systems [22], [23], the central prediction is computed as shown by Eq. (19):

yC,i=Li+Ui2 (19)

In this study, the GA minimizes the objective function given by Eq. (20):

JIT2=JNLL−t+λwidthJwidth+λcoverageJcoverage (20)

where JNLL−t corresponds to the dimensionless Student-t negative log-likelihood. The Student-t NLL is obtained by taking the negative logarithm of the Student-t predictive density, which leads to a residual-dependent logarithmic term together with the scale term and the Gamma-function normalization terms [30]. The dimensionless Student-t negative log-likelihood used in this study was adapted from [30] by defining the location parameter as the IT2-Sugeno central estimate and obtaining the equivalent scale from the type-reduced interval. The corresponding Student-t scale parameter derived from the IT2-Sugeno interval width is defined in Eq. (21).

σeq,i=Ui−Li2qν (21)

The Student-t quantile associated with the target coverage is defined as shown in Eq. (22).

qν=tν,1+Ctarget2 (22)

For the IT2-Sugeno optimization, the loss function used unit weight for the Student-t negative log-likelihood term, with a target coverage of 0.9545. The optimized model was subsequently evaluated using the complete training dataset. As previously mentioned, in this study, the target coverage is set to C_target=0.9545, which gives qν=tν,0.97725. The normalization factor (s) calculated as the mean experimental LDFH of the complete 2024 training dataset is defined in Eq. (23):

s=1N∑i=1Nyi (23)

Thus, the residual term in the Student-t likelihood is evaluated relative to the uncertainty scale implied by the IT2-Sugeno prediction interval, rather than through a fixed variance or a conventional squared-error term. In addition, the width term (Jwidth) shown in Eq. (20) is included to penalize excessively wide prediction intervals, following the principle that predictive intervals should be as compact as possible while achieving the prescribed coverage level [24,31]. The empirical coverage [24], [26], [27], [29] is computed as shown in Eq. (24):

Cemp=Numberofobservationssatisfying(Li≤yi≤Ui)N (24)

In addition, only failure to reach the target coverage is penalized as shown in Eq. (25). The desired coverage level is set to Ctarget=0.9545, which corresponds to the coverage associated with (k=2) in the ideal Gaussian reference case [32].

Jcoverage=max(0,Ctarget−Cemp)2 (25)

Furthermore, the term (Jwidth) shown in Eq. (26) penalizes the normalized full interval width, thereby avoiding trivially conservative solutions in which coverage is achieved solely through excessively wide FOUs and, consequently, overly broad predictive intervals [24], [31].

Jwidth=1N∑i=1NUi−Lis (26)

where Li and Ui are the lower and upper endpoints of the type-reduced interval, as shown in Eq. (27).

yC,i=Li+Ui2 (27)

Therefore, the GA searches the bounded parameter space of θ=[u1,…,uM,β,ν] and the selected solution, θ*=[u1*,…,uM*,β*,ν*], is the candidate that minimizes Eq. (20), yielding the FOU amplitudes, the asymmetric expansion factor, and the Student-t degrees of freedom that provide prediction intervals with limited width and a probabilistically coherent residual distribution. This formulation is consistent with the probabilistic forecasting principle that prediction intervals should be as sharp as possible [24], [31].

As indicated previously, the Type-1 ANFIS models were initialized using grid partitioning, with two Gaussian membership functions per input and first-order linear consequents. This configuration resulted in two rules for M1, four for M21 and M22, and eight for M3. All models were trained for 50 epochs using the valid records from 2024 and the hybrid method, which combines backpropagation to adjust the membership-function parameters and least-squares estimation for the consequents.

Table 2, Table 3, and Table 4 show the performance of the three reduced ANFIS models for both the 2024 training dataset and the independent test years 2019-2023. M1, based only on DHI, shows the largest errors, with nRMSE values ranging from 11.40% to 13.03% and nMAE values from 9.17% to 10.30%. The incorporation of sin(α) in M21 consistently improves performance in all years analyzed, reducing nRMSE relative to M1 by 10.7% during training and, on average, by 13.4% over the test years. M21 also achieves the highest IA and R² values among the reduced models in all test datasets, with nRMSE values ranging from 9.84% to 11.70%. The incorporation of zenith luminance in M22 also improves the results relative to M1, although to a lesser extent, with an average nRMSE reduction of 5.6% over the independent years. M22 nevertheless shows slightly lower nMAE values than M21 during training and in 2023. Taken together, these results show that both additional variables provide useful information compared with the use of DHI alone, but that, when considered individually, sin(α) provides a more consistent improvement in predictive performance than zenith luminance.

Figure 5 shows the three-dimensional distribution of the experimental LDFH and the LDFH estimated by the ANFIS model for the 2024 training dataset. In both cases, the x-axis represents DHI, the y-axis represents sin(α), and the z-axis represents zenith luminance. The color of each point indicates the LDFH value in klux. To improve the readability of the point cloud, the plotted values were averaged within cells defined by discrete intervals of DHI, sin(α), and zenith luminance. Therefore, when several records fall within the same discretization cell defined by these three variables, the mean LDFH value is represented instead of all individual points. The comparison between both plots shows that the ANFIS model adequately reproduces the global structure observed in the experimental data. In particular, the progressive increase in LDFH with DHI is preserved, as well as its joint dependence on solar elevation and zenith luminance. The similarity between the two distributions indicates that the model captures the main behavior of diffuse horizontal illuminance during the training year.

The LDFH distribution shows that it is mainly controlled by DHI. The lowest LDFH values are concentrated at low DHI, especially below approximately 0.05-0.10 kW m⁻², where the color scale remains in the lower range. As DHI increases, LDFH increases progressively, reaching the highest values in the upper DHI region, around 0.3-0.5 kW m⁻². These maximum LDFH values appear mainly when zenith luminance is also high, indicating that large diffuse illuminance is associated with the combined effect of high diffuse irradiance and strong sky brightness near the zenith. Figure 5 shows that maximum LDFH occurs under conditions of simultaneously high DHI, high zenith luminance, and moderate-to-high solar elevation.

Table 5 shows that the Type-1 Sugeno ANFIS model for the direct estimation of LDFH achieves, for the training dataset, an nRMSE of 8.55%, an nMAE of 6.10%, an IA of 0.9951, and an R² of 0.9806, indicating a good fit to the training data. For the independent test years 2019-2023, the performance remains stable, with nRMSE values between 8.37% and 10.59%, IA values above 0.991, and R² values higher than 0.967. Although these results confirm the suitability of the Type-1 ANFIS model as a deterministic estimator, they only quantify point-prediction accuracy and do not provide information about the uncertainty associated with each estimate, which motivates the use of an Interval Type-2 Sugeno model. In practical terms, the training nRMSE of 8.55% for 2024 indicates that the root mean square error corresponds to approximately 8.55% of the mean measured LDFH for the dataset. This value should therefore be interpreted as an overall measure of prediction error relative to the mean diffuse illuminance level, rather than as an 8.55% error for each individual observation. Together with the low nMBE and the high IA and R² values, it indicates that the model adequately reproduces both the magnitude and variability of LDFH.

Compared with the reduced models, the full M3 model shows the best performance across all datasets analyzed. For the 2024 training dataset, M3 reduces nRMSE to 8.55%, compared with 11.40% for M1, 10.18% for M21, and 10.43% for M22, corresponding to reductions of 25.0%, 16.1%, and 18.1%, respectively. This improvement is maintained over the independent test years: the mean nRMSE of M3 is 9.26%, compared with 12.04% for M1, 10.43% for M21, and 11.36% for M22. Thus, relative to the reduced models, M3 decreases the mean nRMSE over the test years by approximately 23.1% compared with M1, 11.2% compared with M21, and 18.5% compared with M22. The same pattern is observed for nMAE, with mean reductions over the test years of 33.1%, 20.4%, and 23.3%, respectively. In addition, M3 consistently achieves higher IA and R² values. These results show that, although DHI is the main predictor and the individual incorporation of either sin(α) or zenith luminance improves the estimation, the combined use of both additional variables provides a further and consistent improvement, thereby justifying the increased complexity of the M3 model. However, the M21 model could also be used when zenith luminance measurements are not available.

To provide a reference for evaluating the added value of the ANFIS approach, regression-based models were also developed using the same combinations of input variables considered in the reduced and full ANFIS models. Thus, M1 uses only DHI, M21 combines DHI and sin(α), M22 combines DHI and zenith luminance, and M3 incorporates all three input variables. All regression models were fitted using the 2024 training dataset and subsequently evaluated for the independent test years 2019-2023 using the same performance metrics employed for the ANFIS models.

The regression baselines consisted of quadratic regression models fitted by ordinary least squares with an intercept. They included linear terms, squared terms, and, where applicable, all two-factor interactions. The resulting M1, M21, M22, and M3 models are shown in Eqs. (28)-(31).

M1:LDFĤ=860.428875+127.195952*DHI−0.018958*DHI2 (28)

M21:LDFĤ=329.587465+103.278083*DHI+5373.394564*sin(α)+20.907082*DHI*sin(α)−0.017816*DHI2−2256.414168*sin2(α) (29)

M22:LDFĤ=537.833481+154.130915*DHI−0.877290*zenithluminance−0.000215*DHI*zenithluminance−0.0616765*DHI2+0.000050*zenithluminance2 (30)

M3:LDFĤ=−174.649876+147.943535*DHI−1.8108435*zenithluminance+6995.932617*sin(α)+0.0006745*DHI*zenithluminance−9.527160*DHI*sin(α)+1.412972*zenithluminance*sin(α)−0.056380*DHI2+0.000018*zenithluminance2−3615.003370*sin2(α) (31)

As can be seen in Table 6, Table 7, Table 8, and Table 9, the four quadratic regression models show improved performance as additional variables are incorporated. M1, based only on DHI, presents the largest errors, with a mean nRMSE of 12.08% and a mean nMAE of 9.95% over the independent test years. The incorporation of sin(α) in M21 produces a clear and consistent improvement, reducing these mean values to 10.60% and 8.49%, respectively, corresponding to reductions of 12.2% in nRMSE and 14.7% in nMAE relative to M1. M22, which incorporates zenith luminance, also improves upon M1, although to a lesser extent, with a mean nRMSE of 11.69% and a mean nMAE of 9.09%. Therefore, when considered individually, sin(α) provides a greater improvement than zenith luminance. The best performance is obtained with M3, which combines all three input variables, with mean nRMSE and nMAE values of 9.84% and 7.33%, respectively. Compared with M1, these values represent reductions of 18.6% and 26.3%, while relative to M21, the additional reductions are 7.2% and 13.7%. M3 also generally presents the highest IA and R² values. Taken together, the results show that the combined use of DHI, sin(α), and zenith luminance provides the best performance among the quadratic regression models, although the largest individual additional contribution relative to DHI comes from sin(α).

The comparison between the ANFIS models and the quadratic regressions shows that the largest difference is observed for M3, for which ANFIS reduces the mean nRMSE from 9.84% to 9.26%, a reduction of approximately 5.9%, and the mean nMAE from 7.33% to 6.64%, a reduction of approximately 9.4%. For the M3 training dataset, the reductions reach 6.8% in nRMSE and 10.7% in nMAE, highlighting the greater ability of ANFIS to represent the nonlinear relationships among DHI, sin(α), zenith luminance, and LDFH. The Type-1 ANFIS structure also provides the basis for the Interval Type-2 extension used in this study, which generates sample-specific predictive intervals that can be evaluated in terms of their empirical coverage and width.

Figure 6 analyzes the residual behavior of the Type-1 ANFIS model for LDFH using the 2024 training dataset. The residual histogram in Figure 6(a) shows that most errors are centered around zero, indicating no substantial systematic bias in the deterministic ANFIS predictions. However, the residual distribution departs from a Gaussian distribution, with a sharper central peak and heavier tails than the Gaussian fit. Likewise, Figure 6(b) compares the empirical quantiles of the standardized residuals with the theoretical standard Student-t quantiles and shows good agreement across the central quantiles, although deviations remain in the distribution tails. Figure 6(c) provides a complementary representation of the residuals as a function of the predicted LDFH, showing an increase in residual dispersion at higher predicted illuminance levels. From a modeling perspective, Figure 6(a) indicates that relatively large prediction errors occur more frequently than would be expected under a Gaussian assumption. This represents a limitation of a purely Gaussian description of the residual uncertainty and supports the use of the Student-t likelihood in the IT2 optimization because of its greater robustness to heavy-tailed errors. However, the remaining deviations in the tails of the Student-t Q-Q plot also indicate that this distribution does not perfectly reproduce all residual behavior.

Figure 7 compares the initial and GA-optimized FOUs of the antecedent membership functions in the IT2-Sugeno model proposed for LDFH estimation. The initial FOUs are constructed from the antecedent membership functions of the previously trained Type-1 ANFIS model. The final FOUs are obtained by optimizing the uncertainty parameters through the genetic algorithm, using the proposed loss function as the optimization criterion.

As can be observed in Figure 7, the optimized FOUs introduce different uncertainty levels across the input variables and membership functions, with larger or smaller expansions reflecting the uncertainty assigned to each antecedent membership function during optimization. The membership function centers are kept fixed from the trained Type-1 ANFIS in order to preserve the learned fuzzy partition, rules, and nominal input-output relationship. Thus, the GA optimization adjusts the FOU parameters without retraining the antecedent structure or the consequents. Uncertainty is incorporated at the antecedent level, allowing the IT2-Sugeno model to generate prediction intervals while retaining the structural basis of the Type-1 ANFIS model.

In addition, Table 10 reports the GA-optimized FOU parameters of the antecedent membership functions used in the LDFH IT2-Sugeno model. The centers and nominal Gaussian widths of the underlying Type-1 membership functions remain unchanged after optimization, thereby preserving the fuzzy partition learned by the Type-1 ANFIS. The GA adjusts the FOU parameters (u and β), together with the Student-t degrees of freedom parameter (ν). The optimized results show that the largest uncertainty amplitudes are associated with the DHI membership functions and the second membership function of sin(α).

Figure 8 shows the spatial distribution of the uncertainty estimated by the IT2-Sugeno model for LDFH over the 2024 training dataset. The highest uncertainty values are mainly concentrated at low solar elevation and at low to moderate zenith luminance levels, whereas the uncertainty tends to decrease at higher sin(α) and higher zenith luminance. This pattern is consistent with the three-dimensional LDFH distributions shown previously for the measured and ANFIS-estimated values. As shown in Figure 5, the largest LDFH values occur mainly for high DHI, high sin(α), and high zenith luminance. However, Figure 8 shows that these high-illuminance regions are not necessarily associated with the largest predictive uncertainty. Therefore, the IT2-Sugeno interval model does not merely increase the interval width as LDFH increases. Instead, it captures the local uncertainty structure of the input-output relationship, with wider intervals in the lower-solar-elevation and lower-zenith-luminance regions.

Table 11 reports the empirical coverage and uncertainty of the IT2-Sugeno interval model for LDFH. The number of samples (N) is large in all cases, ranging from 95,438 observations in the 2024 training dataset to more than 118,000 observations in each external validation year, where N represents the number of records evaluated rather than the number of statistically independent observations. In the 2024 training dataset, the empirical coverage reaches 95.10%, which is very close to the nominal target coverage of 95.45%. This indicates that the GA-optimized IT2-Sugeno intervals provide an appropriate uncertainty level for the data used during model development.

A similar behavior is observed in the external validation years 2019 and 2020, where the empirical coverage is 95.69% and 96.24%, respectively. The years 2021 and 2022 also remain close to the target level, with coverages of 94.43% and 94.76%. The mean uncertainty remains relatively stable throughout the years analyzed, ranging from 3.59 to 4.14 klux. The median values are slightly lower than the corresponding means, which suggests that most samples are associated with moderate interval widths, while a smaller fraction of observations leads to wider intervals.

The lowest coverage is obtained in 2023 (93.00%), which also presents the largest mean uncertainty (4.14 klux), indicating that the lower coverage is not simply due to narrower prediction intervals. In addition, 2023 shows the highest nRMSE (10.59%) and the most negative nMBE (-3.48%) among the independent test years. Since the model was trained using only the 2024 dataset and subsequently applied without retraining to each independent year, these differences may reflect interannual variations in the relationship between the input variables and LDFH. Such variations could therefore contribute to the lower coverage observed in 2023.

Table 12 presents selected examples from the 2024 training dataset to illustrate the behavior of the proposed model. The corresponding 95.45% nominal predictive intervals are also reported as [ŷL, ŷR]_95.45%. The IT2-Sugeno central estimates remain close to the Type-1 ANFIS outputs, as expected, because the interval extension preserves the rule base, consequent parameters, membership-function centers, and initial Gaussian widths of the original ANFIS model. Furthermore, the predicted values are close to the corresponding experimental LDFH values. The uncertainty is introduced through the footprint of uncertainty associated with the antecedent membership functions. The selected examples show that the proposed model provides both direct LDFH estimates and sample-specific predictive intervals using DHI, zenith luminance, and sin(α) as input variables. Thus, the Type-1 Sugeno ANFIS supplies the main LDFH estimation structure, while the IT2-Sugeno extension adds local uncertainty information for each observation. From an applied perspective, these sample-specific intervals can complement the central LDFH estimate in daylighting design and control applications. For example, the lower and upper prediction bounds may be considered when the estimated daylight availability is compared with a reference illuminance level or a lighting-control setpoint. In this way, the model provides information on whether a decision based on the central estimate is associated with a relatively narrow or wide prediction interval, rather than treating all LDFH estimates as equally certain. This additional information may contribute to more robust daylight-responsive control strategies by considering the uncertainty associated with each LDFH estimate.

To assess whether uncertainty can also be quantified when zenith luminance measurements are unavailable, the same IT2-Sugeno formulation and GA optimization procedure were applied to the reduced M21 model, which uses only DHI and sin(α) as inputs. Since the methodology and optimization settings are identical to those described for M3, only the resulting FOU structure and uncertainty are discussed here.

Table 13 shows the optimized parameters of the reduced M21 model, which uses only DHI and sin(α), and may therefore be of interest for stations where zenith luminance measurements are not available.

Figure 9 shows that, as in the full M3 model, the nonuniform redistribution of uncertainty indicates that the optimization does not simply increase the uncertainty of the reduced model globally, but rather distributes it selectively across different regions of the input space according to the predictive-uncertainty patterns identified during optimization.

As shown in Table 14, over the five independent test years (2019-2023), the reduced M21 IT2-Sugeno model yielded mean nRMSE and nMAE values of 11.59% and 9.01%, respectively, compared with 9.33% and 6.71% for the full M3 model. These nRMSE and nMAE values correspond to the IT2-Sugeno central predictions. Each reported mean is the arithmetic mean of the corresponding annual values for the five independent test years. Thus, removing zenith luminance increased the mean nRMSE by 24.3% and the mean nMAE by 34.4%. Despite this reduction in point-estimation accuracy, M21 maintained a mean empirical coverage of 95.74%, compared with 94.82% for M3. The coverage obtained with M21 is very close to the nominal target of 95.45%, indicating that the proposed optimization procedure effectively adjusts the FOU to achieve a coverage level close to the prescribed target even when less input information is available. However, this is achieved by widening the prediction intervals, rather than by improving the underlying point-estimation accuracy. Accordingly, the uncertainty increased from 3.90 klux for M3 to 4.48 klux for M21, corresponding to an increase of 14.8%. These results indicate that zenith luminance improves predictive accuracy and reduces uncertainty, whereas the reduced M21 formulation remains a practical alternative for stations where zenith luminance measurements are unavailable, with the comparison against M3 quantifying the corresponding loss in point-estimation accuracy and interval precision. Similarly, Figure 10 shows that the full M3 model consistently outperforms the reduced M21 model in both predictive accuracy and uncertainty. For all independent test years, M3 yields lower nRMSE values and narrower prediction intervals.

A comparison of Table 12 and Table 15 shows the effect of removing zenith luminance on both point estimation and predictive uncertainty for the same selected data points. The full M3 model provides estimates closer to the experimental LDFH values in most the examples shown. Nevertheless, the experimental values remain within the 95.45% predictive intervals generated by M21 in all cases, showing that the reduced model still provides intervals that contain the experimental values for these selected examples despite the loss of one input variable. This is achieved, however, through generally less precise uncertainty characterization, with wider intervals in some cases, particularly for high-illuminance conditions. For the selected examples, the mean interval width increases from 7.24 klux for M3 to 11.13 klux for M21. These results illustrate how the proposed IT2 optimization compensates for the reduced input information by adapting the FOU and, when necessary, widening the predictive intervals to maintain coverage close to the prescribed target.