Reconstruction-Aided Short-Term Power Forecasting for Data-Limited Distributed PV Stations

1. Introduction

As distributed photovoltaic (PV) penetration increases, stochastic, variable, and geographically dispersed generation makes monitoring, forecasting, and dispatch more difficult in low-voltage distribution networks [1,2]. Compared with centralized plants, distributed installations are numerous, individually small, and widely separated. Differences in module orientation, tilt, local shading, and micro-meteorological conditions therefore produce distinct power profiles. The cost of communication terminals and data acquisition prevents many sites from providing power measurements at fine time intervals or accurate local weather data; available records are often limited to installed capacity, hourly energy, and coordinates [3,4]. These missing series must be reconstructed before reliable forecasts can be made for such sites.

Three groups of methods are relevant to reconstructing power at data-limited PV stations. Neighbor-based estimation [5] includes capacity scaling [6], spatiotemporal imputation [7], and multi-site correlation models [8]. These methods can use readily available measurements, but station-specific cloud effects and differences between frequency bands may be smoothed out. Transfer-based methods [9], including feature transfer [10], autoregressive transfer [11], and spatiotemporally weighted multi-site models [12], alleviate limited training data at new stations. They do not directly resolve the recovery of a 5-min profile when only hourly energy is available as a constraint. Decomposition-based methods separate trends from short-period variations [13]. Wu et al. [14] combined similar-day clustering with dual decomposition. Other recent models use variational mode decomposition (VMD) with a Transformer [15], quadratic mode decomposition with a hybrid bidirectional gated recurrent unit (HBiGRU) [16], or wavelet-threshold denoising with VMD [17]. These studies mainly forecast power from measured series rather than reconstructing fine-resolution power from hourly energy.

Site-specific global horizontal irradiance (GHI) is also unavailable at many data-limited stations. Current approaches use numerical weather prediction and post-processing [18], satellite retrieval and forecasting [19], or data-driven reconstruction [20]. High-resolution irradiance forecasts [21] and weather-foundation-model and satellite combinations [22] provide broad, continuous coverage, but their spatial resolution may not capture local cloud shading. Satellite methods extrapolate cloud fields [23], merge imagery with ground observations [24], or couple cloud images with clear-sky models [25]; their accuracy depends on image resolution, retrieval error, and data availability. Data-driven reconstruction instead learns nonlinear spatial mappings from available irradiance, weather, or PV records. Spatial interpolation can estimate an unobserved location from nearby measurements [26], and Weyll et al. [27] evaluated ExtraTrees for GHI forecasting. For the present task, the target-cell reconstruction must preserve cloud-related fluctuations across several time scales after clear-sky normalization.

After power and micro-meteorological data have been reconstructed, measurements from nearby PV stations can supply spatial information to the forecasting model. Existing regional methods use station clustering [28], deep temporal models [29], or graph-based models [30]. Subregional forecasting [31] and BIRCH coupled with an L-Transformer [32] group stations with similar generation behavior. Other work combines XGBoost decomposition with a long short-term memory network (LSTM) [33] or uses a hybrid bidirectional gated recurrent unit (BiGRU) model [34] to learn temporal dependence. Graph-based models introduce inter-station relationships directly; Ref. [35], for example, combines graph attention and a Transformer for regional PV data-quality enhancement and ultra-short-term forecasting. These methods usually assume that the target site has a reasonably complete power or weather history. When both fine-resolution power and local irradiance are absent, those inputs must be reconstructed before spatial dependence can be used reliably.

This study addresses two linked tasks. The first is to infer 5-min target-station power from hourly energy and neighboring-station measurements. The second is to recover local GHI without discarding cloud-related variation at different time scales. The reconstructed signals are then used with spatial information for short-term forecasting.

(1)

We therefore built a three-stage framework comprising VMD-based power reconstruction, spatial GHI reconstruction, and LightGBM forecasting with spatial graph features. A data-limited station is defined here as one with hourly energy but no high-resolution power measurements. During GHI reconstruction, the target grid cell is withheld as an input, although its historical labels are available for model training. The two signals are reconstructed at 5-min resolution before the next time step is forecast. The specific contributions are:

(2)

First, the cross-station reconstruction method decomposes the target station's hourly energy and the neighboring stations' energy and power sequences into frequency-specific VMD modes. Low-frequency modes describe the intraday trend, intermediate modes retain weather disturbances, and high-frequency modes capture rapid variation. A multidimensional similarity score selects a reference mode for each target energy mode. ExtraTrees transfers the reference-station energy-to-power relationship to the target mode, after which an hourly-energy conservation constraint yields the 5-min power profile.

(3)

Second, GHI at an unobserved target cell is completed from the clear-sky indices and multiscale features of surrounding cells. Clear-sky normalization suppresses the deterministic effects of solar elevation and seasonal periodicity, and VMD then separates the remaining irradiance disturbances by scale. ExtraTrees learns the nonlinear mapping from observed reference cells to the target cell. The resulting location-matched GHI series supplies the missing micro-meteorological input for forecasting.

2. Data Reconstruction for Data-Limited Distributed PV Stations

2.1. Cross-Station Power Reconstruction Using VMD and ExtraTrees

The proposed cross-station method reconstructs high-resolution target-station power from hourly energy records and measurements at neighboring reference stations. It learns nonlinear relationships between reference-station energy modes and power modes, then transfers those relationships to the target station.

2.1.1. Variational Mode Decomposition

VMD represents a nonstationary signal as a prescribed set of band-limited intrinsic mode functions (IMFs) with different center frequencies. This representation permits target and reference stations to be compared mode by mode at several frequency scales. Reference information is therefore chosen independently for each mode. The optimization seeks the smallest aggregate modal bandwidth while requiring the modes to reconstruct the input signal:

$$\begin{array}{c}\underset{\left\{u_k\right\},\left\{{\omega }_k\right\}}{\min }\left\{{\displaystyle \sum _{k=1}^K{‖{\partial }_t\left[\left(\delta (t)+\frac{\mathrm{j}}{\pi t}\right)\ast u_k(t)\right]{\mathrm{e}}^{-\mathrm{j}{\omega }_{k}t}‖}_2^2}\right\}\\ \mathrm{s}.\mathrm{t}.{\displaystyle \sum _{k=1}^K}{u}_k(t)=f(t)\end{array}$$

(1)

where $u_k(t)$ denotes the k-th mode component; ${\omega }_k$ gives the center frequency of the k-th mode; K specifies the number of modes; ${\partial }_t$ denotes partial differentiation with respect to t; $\delta (t)$ denotes the impulse function; * denotes convolution; and j is the imaginary unit.

PV energy and power contain a low-frequency trend together with medium- and high-frequency disturbances. The low-frequency part follows the daily irradiance cycle and regional weather background, whereas the faster components retain changes caused by moving clouds, local shading, and other short events. A quadratic penalty and Lagrange multipliers convert Equation (1) to an unconstrained problem, which is solved by the alternating direction method of multipliers. Reference modes are selected separately for each target mode rather than by assigning one station to the complete series. Figure 1 shows the VMD decomposition of hourly target-station energy.

2.1.2. Hilbert Modal Feature Extraction

For each IMF, we use the Hilbert transform to describe its local oscillation. The transform of mode k is:

$$H[u_k(t)]={\displaystyle \frac{1}{\pi }}{\displaystyle {\int }_{-\infty }^{+\infty }\frac{u_k(\tau )}{t-\tau }\mathrm{d}\tau }$$

(2)

Combining $u_k(t)$ with its Hilbert transform gives the analytic signal:

$$z_k(t)=u_k(t)+\mathrm{j}H[u_k(t)]=a_k(t){\mathrm{e}}^{\mathrm{j}{\varphi }_k(t)}$$

(3)

where, $z_k(t)$ is the analytic representation of the k-th mode;$a_k(t)$ gives its instantaneous amplitude, whereas ${\varphi }_k(t)$ gives its instantaneous phase.

Instantaneous amplitude measures the local strength of a modal oscillation, whereas instantaneous phase describes how its state evolves with time. Together, they characterize the dynamics of different frequency bands, including the rapid evolution of high-frequency disturbances.

2.1.3. Multidimensional Modal Similarity and Reference Selection

Seasonal changes affect both irradiance and correlations between stations. We therefore calculate similarity within each month, which limits comparisons across seasons but leaves enough samples for stable estimates. The matching score uses five quantities: waveform correlation, Hilbert-envelope similarity, phase consistency, ramp similarity, and modal-energy ratio.

(1)

Modal Waveform Correlation

The Pearson correlation coefficient measures waveform similarity between the target station and candidate station j for IMF k in month m:

$$S_{\mathrm{corr},j,k}^m={\displaystyle \frac{{\displaystyle \sum _{t=1}^T\left(u_{j,k}^m(t)-{\overline{u}}_{j,k}^m\right)\left(u_{0,k}^m(t)-{\overline{u}}_{0,k}^m\right)}}{\sqrt{{\displaystyle \sum _{t=1}^T{\left(u_{j,k}^m(t)-{\overline{u}}_{j,k}^m\right)}^2}\sqrt{{\displaystyle \sum _{t=1}^T{\left(u_{0,k}^m(t)-{\overline{u}}_{0,k}^m\right)}^2}}}}}$$

(4)

where $S_{\mathrm{corr},j,k}^m$ is the modal waveform correlation.${\overline{u}}_{0,k}^m$ and ${\overline{u}}_{j,k}^m$ are the mean values of the k-th IMF component for the target station and the j-th candidate station in month m, respectively. T stands for the total number of time points adopted in the calculation within this month.

(2)

Hilbert Envelope Similarity

Applying the Hilbert transform to each IMF yields its instantaneous amplitude envelope. Let $A_{0,k}^m(t)$ and $A_{j,k}^m(t)$ be the amplitude envelopes of the k-th IMF component in month m for the target station and the j-th candidate station, respectively. Their envelope similarity is calculated from:

$$S_{\mathrm{env},j,k}^m={\displaystyle \frac{{\displaystyle \sum _{t=1}^T\left(A_{j,k}^m(t)-{\overline{A}}_{j,k}^m\right)\left(A_{0,k}^m(t)-{\overline{A}}_{0,k}^m\right)}}{\sqrt{{\displaystyle \sum _{t=1}^T{\left(A_{j,k}^m(t)-{\overline{A}}_{j,k}^m\right)}^2}}\sqrt{{\displaystyle \sum _{t=1}^T{\left(A_{0,k}^m(t)-{\overline{A}}_{0,k}^m\right)}^2}}}}$$

(5)

where $S_{\mathrm{env},j,k}^m$ is the Hilbert-envelope score;${\overline{A}}_{0,k}^m$ and ${\overline{A}}_{j,k}^m$ are the mean amplitude envelopes of the k-th target component and the j-th candidate component during month m, respectively.

(3)

Phase Consistency

Let ${\phi }_{0,k}^m(t)$ and ${\phi }_{j,k}^m(t)$ denote the instantaneous phases of the k-th IMF component in month m for the target station and the j-th candidate station, respectively. Then, the phase consistency is expressed as:

$$S_{\mathrm{phase},j,k}^m={\displaystyle \frac{1}{T_m}}{\displaystyle \sum _{t=1}^{T_m}\frac{1+\cos \left[{\varphi }_{j,k}^m(t)-{\varphi }_{0,k}^m(t)\right]}{2}}$$

(6)

where $S_{\mathrm{phase},j,k}^m$ is the phase-consistency score. Larger values indicate closer phase alignment between the two stations for this mode.

(4)

Ramp Variation Similarity

To compare the direction and magnitude of changes at the two stations, the first-order difference of IMF k in month m is defined as:

$$\mathsf{\Delta }{u}_{i,k}^m(t)=u_{i,k}^m(t)-u_{i,k}^m(t-1)$$

(7)

The ramp similarity between the target and candidate stations is then calculated as:

$$S_{\mathrm{ramp},j,k}^m={\displaystyle \frac{{\displaystyle \sum _{t=1}^T\left(\mathsf{\Delta }{u}_{j,k}^m(t)-{\overline{\mathsf{\Delta }u}}_{j,k}^m\right)\left(\mathsf{\Delta }{u}_{0,k}^m(t)-{\overline{\mathsf{\Delta }u}}_{0,k}^m\right)}}{\sqrt{{\displaystyle \sum _{t=1}^T{\left(\mathsf{\Delta }{u}_{j,k}^m(t)-{\overline{\mathsf{\Delta }u}}_{j,k}^m\right)}^2}}\sqrt{{\displaystyle \sum _{t=1}^T{\left(\mathsf{\Delta }{u}_{0,k}^m(t)-{\overline{\mathsf{\Delta }u}}_{0,k}^m\right)}^2}}}}$$

(8)

where $S_{\mathrm{ramp},j,k}^m$ is the ramp-similarity score; $\mathsf{\Delta }{u}_{i,k}^m(t)$ gives the first difference of the k-th IMF component in month m;${\overline{\mathsf{\Delta }u}}_{0,k}^m$ and ${\overline{\mathsf{\Delta }u}}_{j,k}^m$ are the mean first-difference sequences for the target and j-th candidate stations during month m, respectively.

(5)

Modal Energy Ratio Similarity

The energy of IMF k in month m is defined as:

$$E_{i,k}^m={\displaystyle \sum _{t=1}^T{\left[u_{i,k}^m(t)\right]}^2}$$

(9)

The modal energy ratio similarity between the target and candidate stations is calculated as:

$$S_{\mathrm{energy},j,k}^m=1-{\displaystyle \frac{\left|E_{j,k}^m-E_{0,k}^m\right|}{E_{j,k}^m+E_{0,k}^m}}$$

(10)

where $S_{\mathrm{energy},j,k}^m$ is the modal energy-ratio score;$E_{0,k}^m$ and $E_{j,k}^m$ are the modal energies for the k-th IMF at the target and candidate stations during month m, respectively. Using these five quantities, we calculate the overall similarity between the target and j-th candidate stations for the k-th IMF during month m:

$$\begin{array}{l}{S}_{j,k}^m=F_{j,k}^m({\lambda }_1{S}_{\mathrm{corr},j,k}^m+{\lambda }_2{S}_{\mathrm{env},j,k}^m\\ +{\lambda }_3{S}_{\mathrm{phase},j,k}^m+{\lambda }_4{S}_{\mathrm{ramp},j,k}^m+{\lambda }_5{S}_{\mathrm{energy},j,k}^m)\end{array}$$

(11)

$$F_{j,k}^m=1-{\displaystyle \frac{\left|{\omega }_{j,k}^m-{\omega }_{0,k}^m\right|}{{\omega }_{j,k}^m+{\omega }_{0,k}^m}}$$

(12)

where $S_{j,k}^m$ is the overall similarity between the target and j-th candidate stations for the k-th IMF during month m;$F_{j,k}^m$ is a correction factor that downweights candidate modes whose center frequency differs markedly from the target mode;${\omega }_{0,k}^m$ and ${\omega }_{j,k}^m$ give the center frequencies of the k-th mode during month m at the target and j-th candidate reference stations; ${\lambda }_1\sim {\lambda }_5$ are the five similarity weights, respectively.

Weights for the five similarity measures are obtained with the Criteria Importance Through Intercriteria Correlation (CRITIC) procedure [36]. For measure i, its information content and assigned weight are:

$$C_i={\sigma }_i{\displaystyle \sum _{h=1}^5\left(1-r_{ih}\right)}$$

(13)

$$r_{ih}={\displaystyle \frac{{\displaystyle \sum _{n=1}^{Ns}\left(S_{n,i}-{\overline{S}}_i\right)\left(S_{n,h}-{\overline{S}}_h\right)}}{\sqrt{{\displaystyle \sum _{n=1}^{Ns}{\left(S_{n,i}-{\overline{S}}_i\right)}^2}}\sqrt{{\displaystyle \sum _{n=1}^{Ns}{\left(S_{n,h}-{\overline{S}}_h\right)}^2}}}}$$

(14)

$${\lambda }_i={\displaystyle \frac{C_i}{{\displaystyle \sum _{i=1}^5{C}_i}}}$$

(15)

where ${\sigma }_i$ is the standard deviation of the i-th metric across different candidate mode pairs;$r_{ih}$ is the correlation coefficient between the i-th metric and the h-th metric; $N_s$ is the number of candidate mode pairs participating in the correlation calculation; $S_{n,i}$ and $S_{n,h}$ represent the values of the i-th and h-th similarity metrics for the n-th candidate mode pair, respectively;${\overline{S}}_n$ and ${\overline{S}}_h$ represent the mean values of the corresponding similarity metrics, respectively.

The modal energy ratio quantifies each IMF's relative contribution to the original sequence. The energy contribution of IMF k is defined as:

$${\beta }_k^m={\displaystyle \frac{E_{j,k}^m}{{\displaystyle \sum _{k=1}^K{E}_{j,k}^m}}}$$

(16)

The M highest-ranking candidate modes form the reference set ${\mathsf{\Omega }}_k^m$ for the target mode k in month m. The similarity scores are normalized to determine the contribution of each reference mode:

$$w_{j,k}^m={\displaystyle \frac{S_{j,k}^m}{{\displaystyle \sum _{r\in {\mathsf{\Omega }}_k^m}{S}_{r,k}^m}}}$$

(17)

where $w_{j,k}^m$ is the contribution coefficient of the j-th reference mode corresponding to the k-th IMF component of the target station in month m, satisfying ${\displaystyle \sum _{j\in {\mathsf{\Omega }}_k^m}{w}_{j,k}^m=1}$.

Equations (4)–(17) combine waveform correlation with envelope similarity, phase consistency, ramp similarity, and modal energy information. This provides a multidimensional measure of similarity between target and reference modes. Figure 2 summarizes the calculation procedure.

2.1.4. Modal Power Reconstruction Using ExtraTrees

ExtraTrees is a nonlinear regression ensemble comprising randomized decision trees, as illustrated in Figure 3. Here, it relates hourly energy modal features and reference power modal features to high-resolution power components, enabling reconstruction of the target station's 5-min power profile.

The training set for IMF k is defined as:

$$D_k=\{\left({\mathit{X}}_{1,k},{\mathit{Y}}_{1,k}\right),\left({\mathit{X}}_{2,k},{\mathit{Y}}_{2,k}\right),\dots ,\left({\mathit{X}}_{N,k},{\mathit{Y}}_{N,k}\right)\}$$

(18)

where ${\mathit{X}}_{i,k}$ is the input feature vector of the i-th sample under the k-th mode; ${\mathit{Y}}_{i,k}$ is the corresponding target power modal component for this sample; N is the number of samples.

The regression output of decision tree b is expressed as:

$${\widehat{\mathit{Y}}}_k^{(b)}=T_{b,k}\left({\mathit{X}}_k\right)$$

(19)

The reconstructed k-th power mode is obtained by averaging the outputs of all trees:

$${\widehat{\mathit{Y}}}_k={\displaystyle \frac{1}{B}}{\displaystyle \sum _{b=1}^B}{\widehat{\mathit{Y}}}_k^{(b)}$$

(20)

A separate ExtraTrees model is fitted for each IMF. The reconstructed components are arranged chronologically and summed to obtain the target station's 5-min power series.

2.2. Spatial Irradiance Reconstruction at Target Grid Cells

A spatial reconstruction method is used to estimate GHI at target grid cells where high-resolution local irradiance inputs are unavailable. Like the power reconstruction stage, it combines VMD features with nonlinear ExtraTrees regression. The power stage transfers relationships between hourly energy and 5-min power across stations, whereas the irradiance stage maps contemporaneous conditions at reference grid cells to a target grid cell.

2.2.1. Description of the Spatial Reconstruction Problem

Let there be a total of $N_G$ GHI grids within the study area, where the reference grid set is ${\Omega }_r$ and the target grid set is ${\Omega }_q$, expressed respectively as:

$${\Omega }_r=\{r_1,r_2,\dots ,r_{Nr}\}$$

(21)

$${\Omega }_q=\{q_1,q_2,\dots ,q_{Nq}\}$$

(22)

where $N_r$ and $N_q$ are the number of reference grids and the number of target grids, respectively, satisfying $N_r+N_q=N_G$.

For any target grid $q\in {\Omega }_q$, the method uses the GHI states of the reference grids at the same time and their multiscale modal features to reconstruct the GHI sequence at the target grid. This process can be expressed as:

$${\widehat{G}}_{q,t}=F(X_{{\Omega }_r,t},z_{q,t})$$

(23)

where ${\widehat{G}}_{q,t}$ is reconstructed GHI for target grid q at time t; $X_{{\Omega }_r,t}$ collects irradiance predictors from the observed reference grids at time t; $z_{q,t}$ contains time, location, and clear-sky radiation descriptors for the target grid; and $F(\cdot )$ denotes the spatial reconstruction function.

2.2.2. Normalized GHI Representation Using the Clear-Sky Index

Direct spatial mapping of raw GHI does not adequately represent differences caused by cloud shading and meteorological variability. We therefore normalize GHI by the clear-sky index. For grid i at time t, let the observed GHI be $G_{i,t}$ and its clear-sky envelope be $G_{i,t}^{cs}$. The normalized index is:

$$c_{i,t}={\displaystyle \frac{G_{i,t}}{\max (G_{i,t}^{cs},\epsilon )}}$$

(24)

where $\epsilon$ is the effective clear-sky radiation threshold, used to prevent abnormal values caused by an excessively small denominator during nighttime or periods of extremely low radiation.

The clear-sky index is bounded to limit the influence of extreme samples, including anomalous irradiance values and cloud enhancement:

$$c_{i,t}^{\ast }=\min \left[\max (c_{i,t},0),c_{\max }\right]$$

(25)

where $c_{\max }$ is the upper limit of the clear-sky index.

2.2.3. Multiscale Irradiance Features at Reference Grid Cells

Since the clear-sky index sequence contains meteorological perturbations across different time scales, this study uses the VMD methodology from Section 2.1 to perform multiscale decomposition on the clear-sky index sequence of the known reference grids. This extracts the multiscale states of GHI perturbations for the reference grids, serving as input features for the spatial mapping model. Let the clear-sky index sequence of the known reference grid be$r_i$ be $c_{r_i}(t)$:

$$c_{r_i}(t)={\displaystyle \sum _{k=1}^K{u}_{i,k}(t)}$$

(26)

The Hilbert transform described in Section 2.1.2 provides the instantaneous amplitude and phase of each mode. These quantities characterize local oscillations. The multiscale features for a reference grid cell at time t are expressed as:

$$U_{{\mathsf{\Omega }}_r,t}=\left[\begin{array}{c}{u}_{r_1,1}(t),\dots ,u_{r_1,K}(t),\dots ,u_{r_{N_r},1}(t),\dots ,u_{r_{N_r},K}(t)\end{array}\right]$$

(27)

$$H_{{\mathsf{\Omega }}_r,t}=\left[\begin{array}{c}{a}_{r_1,1}(t),{\varphi }_{r_1,1}(t),\dots ,a_{r_{N_r},K}(t),{\varphi }_{r_{N_r},K}(t)\end{array}\right]$$

(28)

where $U_{{\mathsf{\Omega }}_r,t}$ collects VMD modes for the observed reference grid;$H_{{\mathsf{\Omega }}_r,t}$ contains the associated Hilbert features; and $a_{r_i,k}(t)$ and ${\varphi }_{r_i,k}(t)$ give the amplitude and phase of the k-th mode for reference grid i at time t, respectively.

2.2.4. Spatial Mapping Features

The target-grid inputs describe the regional irradiance state, multiscale perturbations, temporal periodicity, and spatial statistics. The reference-grid clear-sky index vector at time t is:

$$c_{{\mathsf{\Omega }}_r,t}=\left[\begin{array}{c}{c}_{r_1,t},c_{r_2,t},\dots ,c_{r_{N_r},t}\end{array}\right]$$

(29)

This vector represents the contemporaneous regional irradiance field relative to clear-sky conditions. Its mean, standard deviation, maximum, and minimum describe spatial variability across reference grid cells:

$${\mu }_{{\mathsf{\Omega }}_r,t}={\displaystyle \frac{1}{N_r}}{\displaystyle \sum _{i=1}^{N_r}{c}_{r_i,t}}$$

(30)

$${\sigma }_{{\mathsf{\Omega }}_r,t}=\sqrt{{\displaystyle \frac{1}{N_r}}{\displaystyle \sum _{i=1}^{N_r}{\left(c_{r_i,t}-{\mu }_{{\mathsf{\Omega }}_r,t}\right)}^2}}$$

(31)

$$c_{{\mathsf{\Omega }}_r,t}^{\max }=\underset{r_i\in {\mathsf{\Omega }}_r}{\max }{c}_{r_i,t},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em} c_{{\mathsf{\Omega }}_r,t}^{\min }=\underset{r_i\in {\mathsf{\Omega }}_r}{\min }{c}_{r_i,t}$$

(32)

Temporal features retain information about diurnal and annual periodicity:

$$T_t=\left[\begin{array}{c}{T}_{d,t}^{\sin },T_{d,t}^{\cos },T_{y,t}^{\sin },T_{y,t}^{\cos }\end{array}\right]$$

(33)

$$\begin{array}{l}{T}_{d,t}^{\sin }=\sin \left(2\pi {\displaystyle \frac{h_t}{24}}\right),\hspace{0.33em}\hspace{0.33em}\hspace{0.33em} T_{d,t}^{\cos }=\cos \left(2\pi {\displaystyle \frac{h_t}{24}}\right)\\ T_{y,t}^{\sin }=\sin \left(2\pi {\displaystyle \frac{d_t}{365}}\right),\hspace{0.33em}\hspace{0.33em}\hspace{0.33em} T_{y,t}^{\cos }=\cos \left(2\pi {\displaystyle \frac{d_t}{365}}\right)\end{array}$$

(34)

where $T_t$ is the temporal periodicity feature vector corresponding to time t;$h_t$ is the hour of the day; and $d_t$ is the day of the year.

The complete input vector for target grid cell q at time t is:

$$x_{q,t}=[c_{{\mathsf{\Omega }}_r,t},U_{{\mathsf{\Omega }}_r,t},H_{{\mathsf{\Omega }}_r,t},P_{{\mathsf{\Omega }}_r,t},T_t,D_{q,{\mathsf{\Omega }}_r}]$$

(35)

where $P_{{\mathsf{\Omega }}_r,t}$ is the spatial statistical feature of the reference grid set ${\mathsf{\Omega }}_r$ at time t; and $D_{q,{\mathsf{\Omega }}_r}$ is the spatial position-related feature between the target grid q and the reference grid ${\mathsf{\Omega }}_r$, which is composed of variables such as latitude and longitude differences, spatial distance, or relative orientation.

The feature vector in Equation (35) is passed to ExtraTrees to learn the nonlinear relationship between regional multiscale irradiance conditions and the target-grid clear-sky index. The model, described in Section 2.1.4, outputs the reconstructed clear-sky index for grid cell q at time t.

The training set for target grid cell q is defined as:

$$D_q={\{(x_{q,t},c_{q,t})\}}_{t=1}^{N_{\mathrm{tr}}}$$

(36)

where $x_{q,t}$ contains predictors constructed from the observed reference grids;$c_{q,t}$ is the historical observed clear-sky index of the target grid and is used only as supervision during training; and $N_{\mathrm{tr}}$ gives the number of training samples.

3. LightGBM Power Forecasting with Spatial Graph Features

3.1. LightGBM Model

LightGBM builds a regression model by adding decision trees sequentially. At each iteration, a new tree approximates the negative loss gradient, and the accumulated tree outputs predict target power from the input features. The approach handles heterogeneous tabular predictors and nonlinear effects with modest training cost. Its structure is shown in Figure 4.

We write the training data as:

$$D_{\mathrm{p}}={\{(x_i,y_i)\}}_{i=1}^{N_{\mathrm{p}}}$$

(37)

where $x_i$ contains the predictors for sample i;$y_i$ is its PV power label; and $N_{\mathrm{p}}$ gives the number of training cases..

Adding the outputs of J regression trees yields the forecast:

$${\widehat{y}}_i=F_J(x_i)={\displaystyle \sum _{j=1}^J{f}_j(x_i)}$$

(38)

where ${\widehat{y}}_i$ denotes forecast power; J counts the trees, and $f_j(x_i)$ is the prediction from tree j.

The training objective balances empirical loss and model complexity:

$$L={\displaystyle {\displaystyle \sum _{i=1}^N}l\left(y_i,{\widehat{y}}_i\right)}+{\displaystyle {\displaystyle \sum _{j=1}^J}C\left(f_m\right)}$$

(39)

where $l\left(y_i,{\widehat{y}}_i\right)$ is the loss function, with the squared error loss employed for the power forecasting task in this paper; and $C\left(f_m\right)$ is the complexity regularization term of the j-th tree, used to mitigate model overfitting.

3.2. Leaf-Wise Growth Strategy

A level-wise tree splits all eligible nodes at one depth before proceeding deeper, which produces a balanced tree even when some splits offer little loss reduction. LightGBM instead expands a single leaf per iteration and selects the leaf that gives the greatest gain. With an equal split budget, this leaf-wise rule can achieve a larger decrease in training loss.

Figure 5. Comparison of level-wise and leaf-wise tree growth strategies.

Figure 5. Comparison of level-wise and leaf-wise tree growth strategies.

The model receives reconstructed power, reconstructed GHI, time indicators, and aggregated neighboring-station power. Their association with the next power value is strongly nonlinear around cloud passages and cloud-edge enhancement. By choosing the largest-gain leaf, LightGBM allocates more splits to these locally complex regions and can represent sharp changes without expanding the entire tree.

3.3. Histogram-Based Decision Tree Algorithm

Most predictors in the forecast are continuous, so testing every observed value as a split point would be expensive in both time and memory. LightGBM discretizes each continuous predictor into a fixed set of bins. Split statistics are accumulated at the bin level, which reduces the search cost. Figure 6 depicts this operation.

Let the continuous feature $x_i$ be mapped to the b-th feature bin, as follows:

$$bin(x_i)=b,\hspace{0.33em}\hspace{0.33em}\hspace{0.33em} b\in \{1,2,\dots ,N_{\mathrm{bin}}\}$$

(40)

where $bi{n}_b$ represents the sample set falling into the b-th bin.

3.4. Power Forecasting with Spatial Graph Features

LightGBM maps reconstructed power, reconstructed GHI, temporal features, and spatial graph aggregates to the next-step power of the target distributed PV station. The input vector is:

$$x_t=\left[\begin{array}{c}{x}_t^P,x_t^{GHI},x_t^T,x_t^{sp}\end{array}\right]$$

(41)

Spatial graph features summarize how power at nearby reference stations varies with target-station power. Denote the reference stations in the target neighborhood by ${\mathsf{\Omega }}_S$, the per-unit power of neighbor j at time t by $p_{j,t}^{pu}$, and its graph weight by $w_j^{sp}$. Their weighted aggregate is:

$$P_t^{sp}={\displaystyle \sum _{j\in {\mathsf{\Omega }}_{nb}}{w}_j^{sp}{p}_{j,t}^{pu}}$$

(42)

where $w_j^{sp}$ is the spatial graph weight, and satisfies ${\displaystyle \sum _{j\in {\mathsf{\Omega }}_S}{w}_j^{sp}=1}$.

4. Forecasting Workflow

The framework combines cross-station power reconstruction, spatial irradiance reconstruction, and LightGBM forecasting with spatial graph features. Figure 7 summarizes the four stages.

Step 1: Power reconstruction. The hourly energy sequences from the target and candidate reference stations are separated into VMD modes. A multidimensional score orders the reference modes, and normalized weights are assigned to the top three. This small reference set avoids dependence on one station while excluding many weakly related inputs. ExtraTrees estimates each target power mode from the corresponding energy and reference-power modes. After summation, the profile is adjusted so that every hourly total matches the recorded energy, producing a 5-min target-station series.

Step 2: Spatial irradiance reconstruction. Regional GHI, coordinates, and time information are collected. Clear-sky normalization reduces the influence of solar elevation and seasonal cycles. VMD extracts multiscale perturbations from the reference-grid clear-sky index. These features, together with temporal, spatial, and statistical inputs, are passed to ExtraTrees to reconstruct irradiance at the target grid cell.

Step 3: Power forecasting. Reconstructed power, reconstructed GHI, and temporal features form the basic inputs. Spatial distance and power correlation define relationships with neighboring stations, whose weighted power values provide spatial graph features. LightGBM combines these inputs through successive regression trees to predict target-station power at the next time step.

Step 4: Result analysis. Power reconstruction, spatial GHI reconstruction, and power forecasting are evaluated separately. Performance is assessed with ${\epsilon }_{\mathrm{NMAE}}$, ${\epsilon }_{\mathrm{NRMSE}}$, and $R^2$.

5. Case Study

5.1. Computing Environment and Data Sources

The reconstruction, spatial reconstruction, and forecasting models were implemented in Python with PyCharm. The case study used measured 5-min power from distributed PV stations in a UK region in 2021, together with station coordinates and installed capacities. GHI was obtained from a 7 × 7 Copernicus Atmosphere Monitoring Service (CAMS) grid covering the same region and aligned to 5-min intervals. Data from January–September 2021 were used for training, and data from October–December were reserved for testing. Computations were performed on Windows 11 with an RTX 4060 GPU, an R9-7940HX CPU, and 16 GB of RAM.

5.2. Evaluation Metrics

We quantify error with normalized mean absolute error ${\epsilon }_{\mathrm{NMAE}}$, normalized root mean square error ${\epsilon }_{\mathrm{NRMSE}}$, and measure goodness of fit with the coefficient of determination R².

$${\epsilon }_{\mathrm{NMAE}}={\displaystyle \frac{1}{N_{\mathrm{ev}}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{ev}}}\frac{|o_i-{\widehat{o}}_i|}{o_{\max }}\times 100\%}$$

(43)

$${\epsilon }_{\mathrm{NRMSE}}={\displaystyle \frac{\sqrt{\frac{1}{N_{\mathrm{ev}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{cv}}}{(o_i-{\widehat{o}}_i)}^2}}}{o_{\max }}}\times 100\%$$

(44)

$$R^2=1-{\displaystyle \frac{{\displaystyle \sum _{i=1}^{N_{\mathrm{ev}}}{(o_i-{\widehat{o}}_i)}^2}}{{\displaystyle \sum _{i=1}^{N_{\mathrm{ev}}}{(o_i-\overline{o})}^2}}}$$

(45)

where $N_{\mathrm{ev}}$ denotes the length of the evaluation series;$o_i$ is the observed value at index i; ${\widehat{o}}_i$ is the corresponding estimate;$\overline{o}$ and $\overline{\widehat{o}}$ are their respective series means; and $o_{\max }$ is the maximum observed value.

5.3. Power Reconstruction Performance

5.3.1. VMD Decomposition and Modal Matching

To determine the number of VMD modes, sensitivity analysis is conducted within the range of $K=3\sim 10$ as shown in Table 1.

Increasing the mode count lowers the reconstruction residual and inter-modal correlation but raises the share of low-energy modes (Table 1). At $K=6$, the mean residual ratio is 0.164, low-energy modes account for 33.33%, and mean absolute modal correlation is 0.023. Relative to $K=5$, the residual ratio and absolute modal correlation fall by 16.475% and 8.995%, while the low-energy share rises by 4.160%. For $K\ge 7$, low-energy modes exceed 44.792%; further increases add low-energy components and enlarge the subsequent matching and ExtraTrees models. Considering reconstruction error, retained modal energy, and computational size, we use $K=6$. Figure 8 compares the target and reference energy modes. VMD separates the original energy series into components with distinct center frequencies, and mode-wise matching assigns reference components with similar fluctuations to each target frequency band. These matched components are then used for power reconstruction.

5.3.2. Reconstruction Results

Three reconstruction schemes were evaluated at 5-min resolution: capacity scaling (Scheme 1), spatiotemporally weighted multi-station estimation (Scheme 2), and VMD–ExtraTrees (Scheme 3). The target station's measured 5-min series was excluded from every input. Only hourly energy, installed capacity, and coordinates were supplied to the models; the withheld measurements were used exclusively for evaluation. Table 2 and Figure 9 give the results.

Scheme 1 reproduces the daily profile in Figure 9, but its estimates depart from the measurements during local fluctuations and at several peaks. This method scales a reference profile only by installed capacity; it therefore misses differences in installation geometry and local weather. Scheme 2 uses several reference stations and reduces part of this error. Its weighted average still suppresses some fast variation around peaks and ramps. Relative to Schemes 1 and 2, Scheme 3 reduces ${\epsilon }_{\mathrm{NMAE}}$ by 70.654% and 50.735%, respectively, decreases ${\epsilon }_{\mathrm{NRMSE}}$by 68.747% and 51.042%, respectively, and increases $R^2$ by 0.275 and 0.094, respectively. With VMD–ExtraTrees, the daily shape and shorter variations are retained because matching and nonlinear ExtraTrees regression are performed separately for each frequency band.

5.4. Spatial Irradiance Reconstruction Performance

Figure 10 compares spatial interpolation, direct ExtraTrees, and VMD–ExtraTrees at the target GHI cell. All methods reproduce the daytime cycle, but their behavior differs around peaks, troughs, and rapid changes associated with cloud motion. Distance-based interpolation is smoother than the reference series and misses several abrupt transitions. Direct ExtraTrees follows more of the local structure, although errors remain at some peaks. The VMD–ExtraTrees trace is closest to the reference across both slow and fast variations. Clear-sky normalization removes much of the deterministic solar cycle, leaving VMD to characterize perturbations by scale before nonlinear spatial mapping. This separation improves the representation of regional irradiance and local cloud variability.

5.5. Power Forecasting Accuracy

Five models were compared to evaluate the effects of input features and forecasting architecture: baseline LightGBM (M1), LightGBM with irradiance features (M2), LSTM with spatial graph features (M3), BiGRU with spatial graph features (M4), and the proposed LightGBM with spatial graph features (M5). Table 3 and Figure 11 present their forecasting performance.

Table 3 shows that the ${\epsilon }_{\mathrm{NMAE}}$ , ${\epsilon }_{\mathrm{NRMSE}}$, and $R^2$ of the M1 model are 0.928%, 3.546%, and 0.915, respectively. After incorporating irradiance features into the M2 model,${\epsilon }_{\mathrm{NRMSE}}$ decreases to 3.521%, $R^2$ increases to 0.916, while ${\epsilon }_{\mathrm{NMAE}}$ is nearly unchanged from M1, indicating a modest improvement in these metrics. Overlap between the information in reconstructed power and irradiance may explain the limited gain. The M3 and M4 models use recurrent neural networks to extract temporal variation features of the power sequence, with their $R^2$ being 0.914 and 0.916, respectively, and overall performance is similar to M1 and M2. The BiGRU model M4 achieves a lower ${\epsilon }_{\mathrm{NRMSE}}$ through bidirectional temporal information extraction; however, in the 5-min forecasting scenario, lagged power, differences, and rolling statistics already describe much of the short-term variation, which may limit the additional benefit of recurrent models. Compared with the M1, M2, M3, and M4 models, the ${\epsilon }_{\mathrm{NMAE}}$ of the M5 model decreases by 59.375%, 59.419%, 63.362%, and 61.570%, respectively; its ${\epsilon }_{\mathrm{NRMSE}}$ decreases by 67.061%, 66.828%, 67.228%, and 66.846%, respectively; and its $R^2$ increases by 0.075, 0.074, 0.076, and 0.074, respectively.

In Figure 11, M1–M4 reproduce the overall PV profile but miss several peaks and fast transitions. M5 stays closer to the measured curve, including the rising and falling ramps shown in the inset. In addition to reconstructed power, GHI, and calendar variables, M5 receives simultaneous power from neighboring stations. Those measurements contain current information about nearby cloud conditions that is absent from the target-station history alone. This advantage applies when synchronous neighboring-station data are available.

6. Conclusions

We developed a forecasting framework for distributed PV stations that have hourly energy records but lack high-resolution power and local irradiance. It reconstructs power across stations, recovers GHI spatially, and supplies both signals and spatial graph features to LightGBM. The case study yielded three main findings:

(1)

VMD–ExtraTrees reconstructs 5-min target-station power from hourly energy records and neighboring-station measurements. Relative to spatiotemporally weighted multi-station estimation, the coefficient of determination increases by 0.094. The reconstructed profiles capture both diurnal variation and local fluctuations.

(2)

Clear-sky normalization and multiscale VMD features support nonlinear GHI reconstruction with ExtraTrees. Compared with spatial interpolation and direct ExtraTrees regression, the reconstructed irradiance more closely follows peaks, troughs, and short-term fluctuations, providing a local irradiance input for the target station.

(3)

The final LightGBM model with spatial graph features uses reconstructed power, reconstructed GHI, temporal periodicity, and correlations among neighboring stations. Its coefficient of determination R² reaches 0.990. This result supports combining reconstructed signals with spatial information for short-term forecasting in this case.

Future work will examine how changing weather and cloud movement affect inter-station relationships and develop dynamic spatial models for different weather conditions.