Probabilistic Voltage Harmonic Compliance-Margin Forecasting with Adaptive Conformal Prediction at Low-Voltage Distribution Points of Common Coupling

1. Introduction

The proliferation of power-electronics-interfaced loads at the low-voltage tier of distribution networks has made harmonic distortion a routine operational concern. Grid-connected photovoltaic inverters, electric-vehicle (EV) chargers, variable-speed drives, and single-phase switching loads inject non-sinusoidal currents that propagate as voltage distortion at every connected node, producing elevated transformer losses, accelerated insulation ageing, protective-relay mis-operation, and cable derating near the thermal limit [1,2]. Taghvaie et al. [3] review the increasingly complex harmonic characteristics and estimation challenges in modern power networks, while Li et al. [4] and Kazemtarghi et al. [5] document power-quality challenges associated with distributed photovoltaic and EV charging integration.

IEEE 519-2022 [6] specifies voltage harmonic-distortion limits at the PCC, including an 8% total-voltage-distortion reference limit for systems at or below 1 kV. IEC 61000-4-30 [7] and IEC 61000-4-7 [8] define the measurement and aggregation procedures used to obtain power-quality indices. In this study, these aggregated quantities form the available telemetry modality. Methods that require waveform-level inputs [10,11,12,13] cannot be directly applied when only aggregated index telemetry is available. The present work therefore evaluates forecasting and uncertainty quantification using the retained index measurements rather than continuous waveforms.

Compliance-oriented monitoring is most useful when it provides an estimate of the remaining distance to a reference boundary before a potential exceedance. The available datasets contain no IEEE 519 exceedance events. Consequently, the present evaluation focuses on probabilistic UTHD and compliance-margin forecasting under the observed operating conditions rather than on validated exceedance-event detection. The absence of observed exceedances is therefore treated as a characteristic of the available data, not as evidence that exceedance detection has been validated.

The study addresses four methodological gaps identified in the reviewed literature. First, harmonic forecasting studies use different targets, input modalities, and uncertainty constructions [13,14,16,18]. Second, the present study evaluates multiple forecast horizons on a common operational telemetry dataset using a fixed baseline suite and statistical comparison. Third, the proposed uncertainty layer uses split conformal calibration rather than a parametric residual model, while recognising that the formal finite-sample guarantee depends on the usual exchangeability assumption [33,34]. Fourth, the study evaluates frozen cross-site transfer and quantifies feature-distribution shift rather than retraining a separate point predictor at each site.

An earlier conference pilot [17] established the feasibility of telemetry-first harmonic-state detection on a single 11-hour record at 4-second cadence; this journal paper extends that pilot to IEC 61000-4-30 10-min aggregation, distribution-free probabilistic intervals, three-site external validation, and full computational-footprint reporting.

The contributions are fourfold. (C1) A telemetry-native framework is formulated to forecast voltage UTHD and compliance margin from IEC 61000-4-30 index exports without waveform access. (C2) A multi-output HGBR produces 10-, 30-, 60-, and 120-min ahead per-phase UTHD forecasts and is compared with a fixed baseline suite using Diebold–Mariano tests and block-bootstrap confidence intervals. (C3) An adaptive split-conformal layer provides empirically calibrated prediction intervals, with the exchangeability condition and the observed coverage–width trade-off reported explicitly. (C4) Frozen-model cross-site transfer is evaluated on two independently instrumented external sites, with the point predictor kept fixed and the conformal layer recalibrated on a small labelled target subset; feature-distribution shift is quantified alongside predictive degradation.

The remainder of this paper is structured as follows. Section 2 surveys the related literature. Section 3 formalises the monitoring problem. Section 4 describes the proposed framework. Section 5 presents the experimental results. Section 6 provides discussion including limitations. Section 7 draws conclusions.

2. Related Work

2.1. Harmonic Standards and Measurement Framework

IEEE 519-2022 [6] sets compliance limits at the point of common coupling (PCC): UTHD to $8\%$ and individual voltage harmonic distortion to $5\%$ per order for systems at or below 1 kV. Current-side limits are expressed as total demand distortion (TDD) indexed by the short-circuit ratio $I_{SC}/I_L$ and additionally require $I_{SC}$ and $I_L$, which are absent from standard telemetry exports, making voltage-centric targets the practical choice under index-only data retention.

IEC 61000-4-30 [7] defines Class A measurement and aggregation from a 200-ms basic measurement window through 3-second and 10-minute reporting intervals. IEC 61000-4-7 [8] specifies the underlying discrete Fourier analysis and the prevailing-phasor convention, which was documented in field measurements by Blanco et al. [13] and validated across a population of low-voltage sites by Meyer et al. [19]. IEC 61000-3-6 [20] defines the planning levels and emission-allocation principles by which network operators apportion harmonic budget among connected installations, translating a measured exceedance into a regulatory obligation with financial and operational consequences.

2.2. Power-Quality Monitoring: From Waveform to Telemetry

Power-quality analysis has historically been anchored to continuous high-rate waveform records, from which time-frequency features are derived through Fourier analysis, synchrosqueezed wavelet decompositions, and the Hilbert-Huang transform [1,2,21,22]. Comprehensive surveys confirm that the strongest discriminative methods depend most on this waveform modality and that the gap between laboratory and field performance is largely a data-availability issue [1,2]. Field-deployed analysers of the Metrel MI 550 Poly class [9] can export aggregated power-quality indices rather than continuous records. The present study therefore focuses on the forecasting problem that arises when such aggregated telemetry is the retained data modality.

The aggregated-index modality has also received direct analytical attention. Rodríguez-Pajarón et al. [24] forecast voltage harmonic distortion from smart-meter data in residential distribution networks, demonstrating the feasibility of learning harmonic behaviour from low-rate measurements. That work is complementary to the present framework, which uses exported power-quality indices as the sole input modality and forecasts multi-horizon voltage UTHD.

2.3. Machine Learning for Power-Quality Classification

Tree-based ensembles have also been applied to power-quality event detection and harmonic-source attribution. Mishra et al. [25] combine downsampling with an optimised random forest for islanding and non-islanding power-quality-event classification, while Mella et al. [26] use an ensemble classifier to determine harmonic contribution responsibility at microgrid PCCs. Sharma et al. [27] address reduced-sensor harmonic resonance detection and compensation, and Jaramillo et al. [28] study supervised identification of distributed-energy-resource electrical profiles in low-voltage networks.

Deep architectures remain prominent in power-quality disturbance classification. Vishwanath et al. [10] combine a synchrosqueezed wavelet transform with convolutional neural networks, Khetarpal et al. [11] use a deep convolutional autoencoder, and Chiam et al. [12] combine LSTM, wavelet, and attention mechanisms. Wang et al. [29] also investigate ensemble deep learning for automated power-quality-disturbance classification. These studies illustrate the maturity of waveform-oriented learning, while the present work targets index-only telemetry for forecasting.

2.4. Harmonic Forecasting: Deterministic and Probabilistic

Al Hadi et al. [14] investigate ANN/ANFIS-based harmonic forecasting in hybrid wind–solar systems, while a related conference study addresses adaptive filtering for harmonic mitigation [15]. Garrido-Zafra et al. [13] forecast current-harmonic magnitude and phase-angle time series from 10-min measurements at the University of Córdoba data-centre PCC using standards-based harmonic quantities and a BiLSTM, with empirical uncertainty characterisation. A related 2023 conference study also uses LSTM-based current-harmonic time-series forecasting [18]. The present work instead forecasts voltage UTHD, uses exported PQ indices as inputs, applies adaptive conformal calibration, and evaluates both longer horizons and cross-site transfer.

On the probabilistic side, Bracale et al. [16] combine clustering with probabilistic harmonic forecasting using field measurements. Masood et al. [30] apply quantile LSTM and clustering to probabilistic household-load forecasting. Xie and Sun [31] formulate a piecewise probabilistic harmonic power-flow method for unbalanced residential networks, while Rodríguez-Pajarón et al. [23] address harmonic estimation in partly monitored residential networks with unknown parameters and topology. Zhao and Milanović [32] further consider harmonic-distortion prediction in sparsely monitored networks with renewable generation. These model-based and sparsely monitored approaches provide important context, but they require network information or measurement relationships beyond the index-only forecasting setting considered here.

2.5. Conformal Prediction in Power and Energy Systems

Conformal prediction is a model-agnostic framework for constructing prediction intervals with finite-sample marginal coverage under the exchangeability assumption on the calibration and test data. Shafer and Vovk [33] establish the theoretical foundations and the split-conformal construction without distributional assumptions on the residuals. Romano et al. [34] extend it to conformalized quantile regression, yielding asymmetric heteroscedastic intervals that adapt to local difficulty, which is well-suited to compliance settings where the upper exceedance risk dominates. Gibbs and Candès [35] address non-exchangeable sequences with adaptive conformal inference, which adjusts the miscoverage level online based on observed coverage errors.

Zhang and Zhou [36] investigate conformal uncertainty quantification for industrial time series under distribution shift, providing a relevant example of adapting uncertainty estimation to changing data distributions. Mootoo et al. [37] develop distribution-free uncertainty quantification for wireless network forecasting. These studies motivate the use of adaptive uncertainty procedures, while the present implementation retains the exchangeability caveat of split conformal prediction.

In the literature reviewed for this study, we did not identify a work that combines adaptive split-conformal calibration with harmonic compliance-margin forecasting at low-voltage PCCs. Existing probabilistic harmonic studies use model-based or empirical uncertainty constructions [13,16]; the present study instead evaluates a split-conformal uncertainty layer and reports its empirical coverage and width.

2.6. Cross-Site Transfer in Power-System Forecasting

Dong et al. [38], Zhao et al. [39], and Wang et al. [40] address domain transfer for load and photovoltaic forecasting through adversarial domain adaptation, causal mechanism transfer, and federated maximum-mean-discrepancy clustering, respectively. In each case the forecast target is electrical demand or generation rather than harmonic distortion, and no published cross-site transfer study targets harmonic compliance forecasting. The present work quantifies rather than corrects covariate shift, providing the empirical deployment-oriented evidence for this problem.

3. System Model and Problem Formulation

Consider a low-voltage PCC instrumented by a Class A power-quality analyser that exports an aggregated index vector at a cadence of $\Delta t=10$ min. At discrete time t the analyser provides

$$\mathbf{x}\left(t\right)\in {\mathbb{R}}^d,$$

(1)

where the d components are the IEC 61000-4-7 exported quantities: per-order voltage and current harmonic magnitudes, total distortion scalars, fundamental magnitudes, and, where available, prevailing-phasor angles. No continuous waveform is accessible, and no quantity outside the exported set may be used.

The forecasting task is to predict, at step horizon h, the per-phase voltage total harmonic distortion

$$\mathbf{y}(t+h)=\left[U_{\mathrm{THD}}^A,U_{\mathrm{THD}}^B,U_{\mathrm{THD}}^C\right](t+h)\in {\mathbb{R}}^m,$$

(2)

with $m=3$ phase targets and $h\in \{1,3,6,12\}$ steps, corresponding to 10, 30, 60, and 120 min. UTHD on phase p is

$$U_{\mathrm{THD}}^p\left(t\right)={\displaystyle \frac{\sqrt{{\sum }_{k=2}^K{\left(U_k^p\left(t\right)\right)}^2}}{U_1^p\left(t\right)}}\times 100\%,$$

(3)

where $U_k^p$ is the RMS magnitude of the kth voltage harmonic, $U_1^p$ is the fundamental magnitude on phase p, and K is the highest harmonic order used when individual harmonic magnitudes are available. Thus, $K=11$ at the primary Córdoba site and $K=50$ at Site A. Site B exports native total-distortion indices without individual harmonic magnitudes, so its reported UTHD is treated as an instrument-exported target and is not recomputed from the harmonic sum.

Compliance is expressed as a three-state variable with exceedance threshold ${\tau }_E=8\%$ and a proportional watch threshold ${\tau }_W$:

$$s^p\left(t\right)=\left\{\begin{array}{cc}\mathrm{Exceedance},\hfill & U_{\mathrm{THD}}^p\left(t\right)\ge {\tau }_E,\hfill \\ \mathrm{Watch},\hfill & {\tau }_W\le U_{\mathrm{THD}}^p\left(t\right)<{\tau }_E,\hfill \\ \mathrm{Normal},\hfill & U_{\mathrm{THD}}^p\left(t\right)<{\tau }_W.\hfill \end{array}\right.$$

(4)

Note that ${\tau }_W$ is an operator-configurable early-warning boundary, not itself a compliance threshold. Its derivation and calibration are described in Section 4.3.

For risk-aware operation the framework must also return, for miscoverage level $\alpha$, a prediction interval ${\widehat{C}}_{\alpha }\left(\mathbf{x}\left(t\right)\right)$ satisfying

$$Pr\left(y(t+h)\in {\widehat{C}}_{\alpha }\left(\mathbf{x}\left(t\right)\right)\right)\ge 1-\alpha .$$

(5)

Three operational constraints govern the design. (C1) No waveform access: all inputs are IEC 61000-4-30 exported indices. (C2) Split-conformal marginal coverage is interpreted under the exchangeability assumption, while performance under the observed data conditions is assessed empirically. (C3) Deployable inference latency within the resource envelope of a utility edge node.

4. Proposed Framework

The framework is organised as a three-tier pipeline comprising data preparation, modelling and evaluation, and generalisation/deployment diagnostics. Figure 1 illustrates the architecture.

4.1. Data Architecture

Three independently instrumented datasets are used. The primary dataset is the University of Córdoba data-centre PCC, providing $11,088$ IEC 61000-4-30 Class A compliant samples over 77 continuous days at a 10-min cadence, with per-order magnitudes and prevailing-phasor angles to the eleventh order on three phases (A, B, C) and the neutral conductor. The dataset was sourced independently from a public repository and is also the basis of the bidirectional LSTM study by Garrido-Zafra et al. [13], enabling indirect comparison. Two anonymised external datasets, designated Site A and Site B, provide cross-site evaluation targets. Site A uses line-to-line wiring and exports harmonic magnitudes to the fiftieth order, making triplen harmonics unobservable. Site B exports only total distortion and native demand distortion without per-order magnitudes. The full Site A and Site B contain 1,001 and 965 samples, respectively, at a 10-min sampling interval; these correspond to approximately 6.95 and 6.70 days of observations when no time gaps are assumed. The held-out 150- and 144-sample test partitions therefore correspond to approximately 25 h and 24 h, respectively. The three sites therefore share only a minimal common feature subset of 18 channels, which defines the input space for the cross-site evaluation.

4.2. Feature Engineering

Three feature families are constructed from the exported indices. Sliding-window temporal descriptors capture the recent trajectory of each target channel through windowed means, standard deviations, and lagged first differences computed over windows of 10 min and 30 min. Spectrum-shape descriptors summarise the distribution of harmonic energy through an odd-order ratio, a triplen emphasis ratio (energy at orders $k\in \{3,9\}$ relative to total), and a spectral centroid. These are motivated by the known sensitivity of data-centre UTHD to switched-mode supply loading at those orders. Prevailing-phasor encodings represent harmonic phase relationships as complex features, available only at the primary site. The magnitude-only Pearson correlation between voltage and current individual harmonic distortion is weak across mid-order harmonics ($\left|\rho \right|<0.16$ for orders $k\in \{3,5,7,9\}$, rising to $+0.41$ at $k=11$), which confirms that raw cross-channel magnitudes carry limited predictive value for voltage UTHD and motivates the engineered feature design over naive channel concatenation. The enriched feature set contains 333 features; the common-only subset shared across all three sites contains 162 features.

4.3. IEEE 519 Three-State Detection Logic

The watch threshold scales proportionally with the exceedance limit:

$${\tau }_W=\kappa {\tau }_E,\kappa \in (0,1),$$

(6)

where $\kappa$ is an operator-configurable parameter. The locked operating value $\kappa =0.80$ yields ${\tau }_W=6.4\%$ for the primary site, placing the watch boundary at 80% of the 8% exceedance limit. This threshold is an early-warning boundary, not a compliance threshold; final compliance determination remains the responsibility of the utility under IEC 61000-4-30 procedures. Downward transitions from Watch to Normal apply hysteresis to suppress state chatter, while upward transitions act immediately to preserve alarm sensitivity. The sensitivity of the state machine to $\kappa$ is characterised in Section 5.6.

4.4. Multi-Output HGBR Forecaster

The forecaster is a histogram-based gradient-boosting regressor (HGBR) deployed in a multi-output configuration with one independent head per (target, horizon) pair. Each head predicts

$$\widehat{y}(t+h)={\displaystyle \sum _{b=1}^B}\eta f_b\left(\mathbf{x}\left(t\right)\right),$$

(7)

where $f_b$ is the bth regression tree fit to the negative gradient of the squared-error loss on histogram-binned features, $\eta$ is the learning rate, and B is the boosting iteration count. HGBR is selected over deep recurrent alternatives because (a) it natively handles missing values through histogram binning, which is common in telemetry exports; (b) its per-sample inference satisfies constraint C3; and (c) its accuracy on tabular telemetry matches or exceeds the recurrent baselines as demonstrated in Section 5.3.

Hyperparameters are selected by blocked walk-forward cross-validation that respects temporal order and avoids data leakage:

$$(B^{★},{\eta }^{★},D^{★},{\lambda }_2^{★})=argmin{\displaystyle \frac{1}{V}}{\displaystyle \sum _{v=1}^V}{\mathrm{MAE}}_v,$$

(8)

where V is the number of validation blocks, D is the maximum tree depth, and ${\lambda }_2$ is the L2 regularisation coefficient. The search grid spans $B\in \{500,1000\}$, $\eta \in \{0.05,0.1,0.2\}$, $D\in \{\mathrm{None},3,5\}$, and ${\lambda }_2\in \{0.0,0.1,1.0\}$, totalling 108 configurations.

4.5. Adaptive Split Conformal Prediction

The data are partitioned chronologically into training, calibration, and test sets. After fitting on the training set, absolute residuals are computed on the calibration set:

$$r_i=|y_i-{\widehat{y}}_i|,i\in {\mathcal{I}}_{\mathrm{cal}}.$$

(9)

Theorem 1 (Marginal coverage). If the calibration and test residuals are exchangeable, the split-conformal interval

$${\widehat{C}}_{\alpha }=\left[\widehat{y}-Q_{1-\alpha }\left(r\right),\widehat{y}+Q_{1-\alpha }\left(r\right)\right]$$

(10)

attains $Pr(y\in {\widehat{C}}_{\alpha })\ge 1-\alpha$, where $Q_{1-\alpha }\left(r\right)$ is the $(1-\alpha )(1+1/n_{\mathrm{cal}})$ empirical quantile of the calibration residuals [33,34].

To adapt interval width to recent changes in the residual distribution, a local quantile over the M most recent residuals is mixed with the global quantile through parameter $\lambda \in [0,1]$:

$$Q_{1-\alpha }^{\mathrm{mix}}=(1-\lambda )Q_{1-\alpha }^{\mathrm{global}}+\lambda Q_{1-\alpha }^{\mathrm{local}}\left(M\right).$$

(11)

This local–global construction is treated as an adaptive calibration strategy, not as a distribution-free guarantee under arbitrary distribution shift. Its performance under the study conditions is assessed empirically through coverage and interval width, motivated by work on adaptive conformal inference and conformal uncertainty under distribution shift [35,36].

4.6. Asymmetric Interval Construction

Under IEEE 519 the regulatory cost of an upper exceedance exceeds that of a lower miss. The interval is built from two one-sided conformal scores:

$${\widehat{C}}_{\alpha }=\left[\widehat{y}-Q_{1-\alpha }\left(r^{-}\right),\widehat{y}+Q_{1-\alpha }\left(r^{+}\right)\right],$$

(12)

where $r^{+}$ and $r^{-}$ are the signed calibration residuals above and below the forecast. This yields an asymmetric band whose upper edge is calibrated to the empirically heavier upper tail that governs exceedance alarms.

4.7. Cross-Site Transfer Protocol

The point model and feature scaler are frozen after fitting on the source domain and applied to each target domain without retraining. Only the conformal layer is recalibrated on a small labelled target calibration split (the first 15% of the target site chronologically), which adjusts interval width to the target residual scale while leaving the point predictor unchanged. Thus, the transfer experiment is not completely label-free: target labels are used only for conformal recalibration. Covariate shift is quantified by three statistics: the mean absolute standardised mean $|\overline{z}|$, the mean absolute deviation of the standardised standard deviation from unity $|\sigma -1|$, and the fraction of target samples beyond five source standard deviations.

4.8. Evaluation Metrics

Point accuracy is reported with MAE, root mean square error (RMSE), and symmetric mean absolute percentage error (sMAPE):

$$\begin{array}{cc}\hfill \mathrm{MAE}& ={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}|y_i-{\widehat{y}}_i|,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \mathrm{RMSE}& =\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{(y_i-{\widehat{y}}_i)}^2},\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \mathrm{sMAPE}& ={\displaystyle \frac{100\%}{n}}{\displaystyle \sum _{i=1}^n}{\displaystyle \frac{|y_i-{\widehat{y}}_i|}{\left(\right|y_i|+|{\widehat{y}}_i\left|\right)/2}}.\hfill \end{array}$$

(15)

Interval quality is assessed through empirical coverage and mean width:

$$\mathrm{Cov}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}\mathbf{1}[y_i\in {\widehat{C}}_{\alpha ,i}],\overline{W}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}(u_i-l_i),$$

(16)

and the Winkler score that penalises both width and miscoverage:

$$W_{\alpha }=(u-l)+\frac{2}{\alpha }(l-y)\mathbf{1}[y<l]+\frac{2}{\alpha }(y-u)\mathbf{1}[y>u].$$

(17)

Pairwise forecast significance is assessed with the DM statistic on the loss differential $d_t=L\left({\widehat{e}}_t^1\right)-L\left({\widehat{e}}_t^2\right)$:

$$\mathrm{DM}=\overline{d}/\sqrt{2\pi {\widehat{f}}_d\left(0\right)/n},$$

(18)

with the Harvey-Leybourne-Newbold small-sample correction. Block-bootstrap confidence intervals of block length b preserve serial dependence in every reported MAE figure.

5. Experimental Results

All experiments use a global seed of 42 and a chronological 70/15/15 train, calibration, and test split. Table 1 summarises the three-site partitioning. The primary site is 99.90% Normal with no Exceedance samples. Accordingly, the reported results evaluate forecasting and uncertainty estimation within the observed operating regime rather than exceedance-event detection.

5.1. Dataset Characterisation and Spectral Analysis

Figure 2 presents the primary-site spectral signature. Panel (a) shows per-phase individual voltage harmonic distortion by order from 2 to 11 with interquartile-range error bars. Phase A is the quietest channel across all orders; phases B and C carry elevated content at triplen orders $k=3$ and $k=9$, at even order $k=4$, and at orders $k=7$ and $k=8$, each near or above $0.9\%$. The fifth-harmonic suppression on phase B and the seventh-order asymmetry reflect the non-linear load composition at this data-centre PCC. Panel (b) traces the neutral-conductor THD over the full 77-day record: the mean is 33.7% and the maximum is 44.9%, with dominant content at triplen orders $k=3$ and $k=9$, consistent with switched-mode power supply loading throughout the monitoring period. Panel (c) shows the phase-B UTHD histogram with a kernel density overlay, centred near 4.5% with almost all mass inside the Normal band, and the watch (${\tau }_W=6.4\%$) and exceedance (${\tau }_E=8\%$) thresholds marked.

Table 2 reports the watch-threshold sweep on phase B. The false-alarm rate (FAR) is defined as $\mathrm{FAR}=N_{\mathrm{FA}}/(N_{\mathrm{FA}}+N_{\mathrm{TN}})$, where $N_{\mathrm{FA}}$ denotes samples assigned to Watch while the observed UTHD remains below the watch boundary, and $N_{\mathrm{TN}}$ denotes samples correctly assigned to Normal. The threshold sweep in Table 2 shows that the number of Watch classifications falls sharply as the watch boundary increases. The locked ${\tau }_W=6.4\%$ produces only 11 raw Watch samples in the primary record. The FAR values reported in this section are computed from the underlying classification counts used in the analysis.

5.2. Hyperparameter Selection and Model Fitting

Blocked walk-forward cross-validation over 108 configurations selects an HGBR with $B=1000$ boosting iterations, learning rate $\eta =0.05$, no explicit depth cap ($D=\mathrm{None}$), and L2 regularisation ${\lambda }_2=0.1$, yielding a mean CV MAE of $0.3235$ pp across all targets and horizons. The top five configurations lie within $0.0003$ pp of one another, demonstrating robustness to the precise hyperparameter choice inside the searched region. A single head trains in under 15 s on a commodity eight-core workstation without GPU acceleration.

5.3. Point Forecast Accuracy

Table 3 reports per-phase point accuracy at all four horizons with block-bootstrap 95% CIs. Phase A attains a 10-min MAE of $0.056$ pp; phases B and C, operating near 4.5% UTHD with less than a 2-pp margin to the watch level, are the operationally critical channels and each reach $0.274$ pp at the 10-min horizon.

Figure 3 shows MAE versus forecast horizon for all seven methods on phase B. The two HGBR variants form the lowest error band at every horizon ($0.274$ to $0.320$ pp), while persistence climbs to $0.410$ pp at 120 min, ARIMA holds near $0.71$ pp at all horizons, and ridge degrades to $1.09$ pp at 120 min. Table 4 confirms the proposed model outperforms all five baselines at $p<{10}^{-4}$ at the 10-min horizon on phase B. The DM margins satisfy $\left|\mathrm{DM}\right|>7$ against all baselines, confirming that the ranking is robust to moderate hyperparameter variation in the approximate replications.

5.4. Conformal Prediction Interval Performance

Figure 4 shows the phase-B UTHD forecast against the 90% and 95% adaptive conformal bands over the final nine days of the test set (28 May to 5 June). The actual trace is contained within the 90% band through most of the window, and the band widens near 5 June as recent calibration residuals increase. The ${\tau }_W=6.4\%$ watch level sits comfortably above the forecast distribution, with compliance margin even at the test-set peak near 6.3%.

Table 5 reports the interval metrics. At the nominal 90% level the adaptive interval achieves 88.7% empirical coverage on phase B at 10 min with mean width 1.17 pp, and 93.7% coverage at the nominal 95% level. These are empirical test-set results; they should not be interpreted as guaranteed coverage under distribution shift. The 1.3-pp gap is controllable: as shown in Figure 5(a), coverage rises monotonically with $\lambda$ from 0.859 at $\lambda =0$ to 0.928 at $\lambda =1$. Figure 5(b) shows the M-sensitivity frontier, with coverage closing to 89.9% at $M=100$ and recovering 90.8% at $M=200$, each at a modest width cost. An operator can select any point on this characterised frontier according to operational requirements.

5.5. Cross-Site Frozen Transfer

The embedded $3\sigma$ exceedance rate in Figure 6 denotes the fraction of post-scaling feature values lying more than three source standard deviations from the source training mean, quantifying covariate shift severity at each site.

Figure 6 shows phase-B UTHD under frozen-model cross-site transfer with target-site conformal recalibration and 90% bands. Panel (a) confirms in-domain accuracy (MAE $0.274$ pp, $3\sigma$ exceedance $0.8\%$). Panel (b) shows Site A with a structural positive bias (MAE $1.976$ pp, $3\sigma$ exceedance $74.6\%$). Panel (c) shows Site B retaining substantially lower transfer error (MAE $0.478$ pp, $3\sigma$ exceedance $73.2\%$) despite similarly large absolute shift.

Table 6 and Table 7 report the quantitative cross-site results. Site A and Site B each place more than 55% of samples beyond five source standard deviations, compared with 0.1% at the source site.

5.6. Ablation Studies

Figure 7 presents two ablations. Panel (a) shows test MAE versus horizon for the four feature subsets. All four lie within $0.002$ pp at the 10-min horizon, and the common-only subset of 162 features matches or slightly outperforms the enriched 333-feature set at most horizons, validating the common-feature model as the deployable default for sites without phasor or spectrum-shape exports. Panel (b) shows FAR versus ${\tau }_W$ for all four horizons: FAR falls from $0.90$ at ${\tau }_W=4.0\%$ to near zero at $5.0\%$ and to exactly zero at the locked $6.4\%$, near-identically across all horizons, confirming the stability of the proportional threshold design.

5.7. Computational Cost and Deployability

A single HGBR head trains in approximately 14 s and the full multi-output, multi-horizon model completes within a few minutes on a commodity eight-core workstation without GPU acceleration. Per-sample inference is in the single-digit microsecond range, and the serialised model is compact relative to typical edge-computing storage budgets. These measurements support the computational feasibility of the proposed tabular model for resource-constrained deployment.

5.8. Comparison with Prior Work

The proposed HGBR and the approximate BiLSTM replication are the two strongest forecasters. The proposed model beats the BiLSTM at $p<{10}^{-4}$ on phase B at 10 min while training and serving at far lower cost. In the approximate baseline implementation used in this study, the Beta-LP reaches a phase-B 10-min MAE of $0.321$ pp against $0.274$ pp for the proposed HGBR; the baseline does not provide the same conformal uncertainty construction. A direct numerical comparison with the published results of Garrido-Zafra et al. [13] and Bracale et al. [16] is not possible because their targets differ and their hyperparameters and data splits are not reported at a reproducible level; the BiLSTM and Beta-LP entries are labelled as approximate replications. The large DM margins ($\left|\mathrm{DM}\right|>7$ in all cases) confirm that the ranking is robust to moderate baseline hyperparameter variation.

6. Discussion

6.1. Transfer Asymmetry and Feature-Space Shift

The asymmetric transfer results are consistent with substantial feature-space mismatch between the source domain and Site A. Site A uses line-to-line wiring and therefore does not expose exactly the same neutral/triplen channels as the Córdoba source site. This structural difference is consistent with the large standardised shift statistics and systematic positive prediction bias observed for Site A, although the anonymised configuration prevents direct physical attribution.

Site B exhibits lower transfer error than Site A despite large feature-distribution shift. This indicates that aggregate feature alignment, rather than a single physical explanation, is the more defensible interpretation of the cross-site results. With only two external sites, the shift statistics should therefore be treated as candidate deployment diagnostics rather than validated predictors of transfer success.

6.2. Neutral Harmonic Characteristics

The primary site exhibits elevated neutral-conductor THD with dominant triplen components at orders $k=3$ and $k=9$. Such harmonic structure is consistent with nonlinear single-phase loading, but the available data do not contain independently labelled load classes. The present study therefore treats this observation as a spectral characteristic rather than as a validated non-intrusive load-identification mechanism.

6.3. Coverage Gap and Operational Controllability

The 1.3-pp coverage gap at the locked operating point ($\lambda =0.3$, $M=50$, $\alpha =0.10$) is controllable rather than fixed. Figure 5 shows a fully characterised frontier: $M=100$ closes the gap to 0.1 pp (89.9% coverage) and $M=200$ recovers 90.8%. An operator can select any point on this frontier by trading interval width against coverage tightness, with the trade-off quantified in advance of deployment.

6.4. Compliance Framing

The framework provides standards-referenced early warning relative to IEEE 519-2022 voltage distortion limits. Final compliance assessment remains subject to the full measurement and aggregation procedures of IEC 61000-4-30 and applicable utility practice; the framework estimates compliance margin with calibrated uncertainty and does not replace formal compliance determination. This framing is consistent with the regulatory intent of IEC 61000-3-6, which allocates emission budget prospectively rather than adjudicating compliance retrospectively.

6.5. Limitations

Six limitations are noted. First, none of the three datasets contains IEEE 519 exceedance events, and all three sites represent commercial or IT-intensive load profiles; generalisation to sites with persistent exceedances or to residential and industrial load profiles remains an open question. Second, the external test partitions contain only 150 and 144 samples covering approximately 25 h and 24 h, respectively, which is insufficient for seasonal performance characterisation; the transfer results are preliminary evidence of cross-site portability rather than universal generalisation. Third, the point predictor is frozen during transfer but the conformal layer uses a labelled target-site calibration subset; the experiment therefore does not represent completely label-free deployment. Fourth, ${\tau }_W=6.4\%$ is operator-configured for the Córdoba dataset; practitioners should re-calibrate $\kappa$ in Equation (6) for sites with different UTHD operating points. Fifth, Site A and Site B are anonymised, so no topology or load-composition data are available to establish the physical causes of the observed transfer asymmetry. Sixth, the BiLSTM and Beta-LP entries are approximate replications that preserve the broad published architectures but not original hyperparameter tuning; the reported DM margins provide ranking-robustness evidence within the chosen baseline settings, not a substitute for an exhaustive baseline optimisation study.

7. Conclusions

This paper evaluated a telemetry-first framework for probabilistic voltage harmonic compliance-margin forecasting at low-voltage points of common coupling. The framework forecasts per-phase UTHD from retained power-quality indices without continuous waveform access and evaluates four forecast horizons using a multi-output HGBR and a fixed baseline suite. For the primary phase-B target, the proposed model achieved a 10-min MAE of $0.274$ pp and significantly outperformed persistence in the reported Diebold–Mariano comparison.

The adaptive split-conformal layer produced empirical prediction intervals whose coverage and width could be adjusted through the local-window and mixing parameters. At the locked operating point, the 90% nominal interval achieved 88.7% empirical coverage; the sensitivity analysis increased observed coverage to 90.8% at $M=200$. These results demonstrate empirical uncertainty calibration under the study conditions, while the formal marginal coverage statement remains conditional on exchangeability.

Frozen-model transfer to two independently instrumented external sites showed that prediction error changes substantially with feature-distribution shift. Because the external records are short and the conformal layer is recalibrated using labelled target observations, these results should be interpreted as preliminary evidence of cross-site portability and as motivation for deployment-oriented shift diagnostics, not as proof of long-term or label-free generalisation.

The main limitations are the absence of observed IEEE 519 exceedance events, the short external-site records, the use of anonymised target-site metadata, and the approximate reproduction of selected baselines. Future work should evaluate the framework on datasets containing documented exceedance events, longer multi-season external records, and additional network types, and should investigate site-adaptive modelling and closed-loop operational interventions.