Spatio-Temporal Kriging Based Economic Dispatch Problem Including Wind Uncertainty.

: The incorporation of wind generation introduces challenges to the operation of the power system due to its uncertain characteristics. Therefore, the development of methods to accurately model the uncertainty is necessary. In this paper, the spatio-temporal Kriging and analog approaches are used to forecast wind power generation and used as input to solve an economic dispatch problem, considering the uncertainties of wind generation. Spatio-temporal Kriging takes into account the spatial and temporal information given by the database to enhance wind forecasts. We evaluate the performance of using the spatio-temporal Kriging, and comparisons are carried out versus other approaches in the framework of the economic power dispatch problem, for which simulations are developed on the modiﬁed IEEE 3-bus and IEEE 24-bus test systems. The results show that the use of Kriging-based spatio-temporal models in the context of economic power dispatch can provide an opportunity for lower operating costs in the presence of uncertainty when compared to other approaches. forecasting (FCST), power planning studies (PLNG), and unit commitment or economic dispatch. approaches. Scenario-based approaches are found to have a lower cost than the deterministic method under uncertainty.


Introduction
Renewable generation has received interest over the last decades due to its environmentally sustainable and cost-effective operation relative to traditional energy sources. The integration of this intermittent energy into the power grid presents a range of challenges in the operation and management of these systems because energy can not dispatch in the conventional sense [1][2][3][4]. The scheduling of wind energy in power systems is a decision-making challenge due to the uncertainty inherent in the availability of wind. Addressing the uncertainty of renewable generations is crucial to increase economic benefits and enforce the criteria of reliability for decision-making [1].
Wind depends on several factors, such as temperature, humidity, direction, and other aspects, which can lead to unpredictable behavior. Besides, it is non-stationary and typically has strong diurnal and seasonal trends. It is also temporally and spatially auto-correlated and demonstrates heteroscedasticity when the variance of the disturbances is not consistent in the observations. It is, therefore, a physical phenomenon in which modeling is challenging due to its characteristics [2]. Consequently, there is an urgent need to develop more accurate, scalable, and realistic approaches to modeling short-term renewable generation uncertainties in order to accommodate higher penetration of renewable energy sources into the grid [1].
In the literature, there is a broad range of methods for forecasting wind speed and converting these values to wind power forecasts using power curves [1,2,5]. A widely used approach to represent renewable energy uncertainties is the creation of a set of time-series called scenarios for forecasting the future [2]. Scenario approach has played an essential role in several stochastic and robust optimization problems, such as unit commitment [6,7], economic dispatch [8][9][10], energy trading policy [6], and storage capacity sizing [11].
Although commonly used, generating scenarios that are capable of accurately informing system operators of power generation uncertainties remains a challenging problem [1]. The complexity of modeling unknown stochastic processes related to renewable energy is one of the most critical challenges in the generation of scenarios. Researchers have performed extensive work on the use of probabilistic and statistical models. In [12] Gaussian copula functions are used to generate statistical scenarios that take into consideration the complete temporal dependency structure of wind speed. [13] examines tail dependency when choosing the copula function for wind speed modeling. [4] uses Auto-Regressive Moving Average (ARMA) time series and Monte-Carlo simulation to generate a set of wind speed and load forecast error scenarios. A comprehensive modeling approach for multi-site wind power generation scenarios is presented in [2]. The model uses time series and principal components techniques, along with data preprocessing techniques to capture wind characteristics from multiple sites.
Nevertheless, most of the methods listed above required a significant amount of site-specific expertise for renewable resource modeling. The intermittent and time-varying nature of renewables, the complex spatial and temporal interactions, make it difficult for most of these approaches to be applied and challenging to implement in practice [1,14]. Consequently, the scenarios generated may not reflect the inherent patterns and realistic time-series of actual historical observations of renewable energy resources. Notwithstanding their computational benefits, the limited capacity of mentioned models in capturing intricate patterns of correlations in real-world processes has encouraged the scientist community in the last decade to develop groups of models that capture spatial-temporal interactions.
In spatio-temporal modeling, both spatial and temporal variations of wind speed are taken into account in one model, and the prediction is made accordingly. Table 1 highlights a summary of the works found in the literature on the use of spatio-temporal models in power systems, in particular for forecasting (FCST), power planning studies (PLNG), and unit commitment or economic dispatch. Notice that there is no single method to model the spatio-temporal relationship of the variable to be forecasted. Also, Table 1 specifies the spatial-temporal variables (ST-VAR), whether a point forecasting is made (P-FCST) or scenarios (Scen) are generated, time frame (TF), and finally the model used. The proposed model is included in Table 1 to observe the differences with respect to the cited literature.
In [17], for three separate weather variables, namely irradiance, ambient temperature, and wind speed, the forecasting is performed using a Compressive Spatio-Temporal Forecasting model. Then, it estimates the output of the PV power plant applying a power conversion model. A spatio-temporal Markov chain model for wind power forecasting is proposed in [18]. It extends the traditional discrete-time Markov chain and incorporates off-site reference information to improve the forecast performance of regional wind farms. In [19], the spatial-temporal correlation between member wind farms and aggregate wind power is modeled on a joint distribution model based on copula theory. Inverse sampling was applied to the joint distribution conditional (CDF) to generate aggregated probabilistic forecasts. In [20], the authors propose a probabilistic framework for predicting wind speed from measured spatial-temporal data. The frame is based on the decomposition of spatial-temporal covariance and simulation using this decomposition. This framework allows temporal extrapolation and spatial interpolation. In [21], the authors proposed a probabilistic spatial-temporal model. The model is based on the quantile regression adapted to integrate the LASSO (Least Absolute Shrinkage and Selection Operator) variable selection process. This makes it particularly suitable for situations in which data from a high number of PV installations is available since it is capable of dealing directly with the resulting high dimensionality and over-fitting issues that characterize the use of such large amounts of data. A deep convolution neural network topology for wind speed forecasting with spatio-temporal correlation (PDCNN) is presented in [22]. The model is a unified framework, integrating convolutional neural networks (CNNs) and a multi-layer perceptron (MLP). PDCNN generates the forecasted wind speed by using the spatial-temporal correlations learned. In [23], the authors propose a covariance modeling approach that considers spatial-temporal anisotropy in the spatial-temporal kriging model. Spatial and temporal covariance functions are individually designed to generate synthetic wind power profiles to simulate the effects of fluctuations in power. The formulation of the Wiener cyclo-stationary filter is designed to predict wind based on its past values and the spatial-temporal correlation between the data measured in [24]. An iterative, stochastic gradient predictor that uses a multichannel LMS algorithm has also been suggested. The Wiener cyclo-stationary filter and the LMS provide speed and direction predictions. Two non-linear methods for the production of short-term spatial-temporal wind speed forecasts are presented in [25]. A kernel least mean square algorithm and a kernel least-square recursive algorithm are introduced and used to generate wind speed forecasts at different locations. In [26], a graph-based convolution network is proposed for spatial-temporal wind speed forecasting. The proposed architecture is further improved by using the Rough Set Theory to capture deep wind data interval features by incorporating upper and lower bound parameter approximations in the model. In [27], the authors propose a spatial-temporal model that produces PV power forecasts for a power plant 4 of 25 using measurements from other nearby plants. The method makes it possible to take into account local weather conditions in the spatial-temporal model. The model also incorporates the Least Absolute Shrinkage and Selection Operator to improve the selection of variables.
Spatio-temporal correlations of irradiance data are analyzed for forecasting irradiance in [28,29]. Spatial and spatio-temporal methods are used for model irradiation at an arbitrary point based on a given time series irradiation at specific locations. In [30], the authors use an empiric variogram to analyze the spatial correlation of different locations on a wind farm. They use universal kriging and a Bayesian dynamic model for hierarchical wind farm modeling. Besides, the Gibbs sampling approach is used to analyze the model and to develop wind speed predictions. [32] proposes two Bayesian spatio-temporal models to obtain full probabilistic wind power forecasts. The models are based on a stochastic partial differential equation (SPDE) approach to spatio-temporal modeling, which enables rapid inference through integrated nested Laplace approximations (INLA) as well as dimensional reduction. [33] proposes a multi-channel linear forecast model ARMA or MARMA for short-term wind speed forecasting using neighboring wind speed measurements around the target location. A comprehensive modeling methodology for multi-site wind power generation scenarios for forecasting power generation is presented in [2]. The model uses time series and Principal Component techniques along with pre-processing data techniques to directly capture the main wind characteristics from multiple locations, such as distinct diurnal and seasonal trends, non-Gaussianity, and spatial and temporal correlations.
The works mentioned above have only been used to predict or study the production of renewable energy. The implementation of the approaches and methods used for the operational planning of electrical systems has not been published.
The papers listed below have been identified in the literature on the use of spatio-temporal models for operational planning, especially in ED. In [15,16] a Trigonometric Direction Diurnal Model is used, which integrates the spatio-temporal wind information into the forecasting model. It also formulates an ED problem incorporating the available spatio-temporal wind power forecast data for the short term. Authors in [31] propose a finite-state Markov chain model for wind farm generation based on spatial and temporal characteristics. The Point-forecast of wind farm generation is derived from the Markov chains and integrated into the power system for ED. A scenario generation method based on the Kriging model, along with the importance sampling method, is proposed in [36] to describe the wind uncertainty in the ED of the Energy Internet. The Kriging model is used to to estimate the objective function value of the wind energy problem corresponding to each scenario. In [34], the authors use dynamic uncertainty sets as a modeling technique for the dynamic relationship between uncertainties at the decision-making level. Uncertainty sets explicitly define the temporal and spatial model correlations in wind speed. They also propose a robust two-stage adaptive optimization model for multi-period ED. In [38], the authors propose a rotating regime-switching space-time diurnal wind speed forecasting model (RRSTD) that allows the forecast regime to vary with the dominant wind direction and with the seasons. Moreover, they formulate an ED model that takes into account the information on the space-time wind forecast modeled by the RRSTD. In [39], the authors use the Kriging approach in combination with the Vector Auto-Regressive models for wind forecasting and apply it to dynamic line rating studies. A modified fuzzy adaptive PSO assisted by a Kriging model is proposed in [40] to solve the optimization problem in active distribution networks with the incorporation of multiple controllable resources. The Kriging model is used to calculate the power flow of the active distribution networks approximately during the evolutionary search process.
Based on the above, few studies have shown that spatio-temporal models are used to forecast and create power wind scenarios for a particular location. Besides, there are limited references that consider spatio-temporal Kriging models in the scope of scenario generation, unit commitment, and economic power dispatch (Table 1). In this work, we focus mainly on the Kriging and analogous models to generate wind power scenarios. Kriging or Gaussian Process Regression is a common kernel-based regression model capable of modeling complex functions [41]. It is often referred to as the Best Linear Unbiased Estimator (BLUE) and, as the name implies, is a linear estimator [42]. It determines the values at points of interest as a weighted sum of values at other points, and its variations are primarily based on assumptions about the underlying mean function of the data distribution [42]. The interpolation error is minimized by studying and modelling the spatio-temporal distribution of the points already measured. This spatio-temporal distribution or spatio-temporal variability is represented in the form of an experimental variogram [41]. The variogram is the basis for the application of the Kriging. Generally, Kriging models are just spatial models but can also extend to spatio-temporal problems by considering time as an additional dimension [23]. Compared to other regression methods, the Kriging approach has the following advantages: 1) it provides not only the estimate of the value of a function, but also the mean squared error of the estimation, the so-called Kriging variance [41]; 2) it does not assume the data follow any probability distribution and can fully use both temporal and spatial information; 3) it can use a limited set of sampled data points to estimate the value of a variable over a continuous spatial field.
In this paper, a wind forecasting approach is proposed. It uses spatial and temporal correlations existing in the database of a wind farm to predict the power generation of specific locations using the spatio-temporal Kriging method. The method discusses in this paper differs from the work developed in [36], as it uses a spatio-temporal Kriging model and analogous models to generate wind power scenarios. [36] uses only a spatial Kriging model. Besides, in [36] the model is used to estimate the value of the ED objective function corresponding to each scenario generated by the sampling method. Concerning the work described in [2], the principal components technique and time series are used to create wind scenarios for a security assessment application. Both [36] and [2] works require probabilistic modeling in their initial stages. Most of the works developed and shown in Table 1 are about point forecasts of renewable generation. In contrast, this paper combines the spatio-temporal Kriging method and analogs for generating scenarios. Therefore, The main contributions of this paper include: a) The development of a wind speed forecasting algorithm for renewable energy systems based on the spatio-temporal Kriging method and analog models; b) Comparison of the proposed method with commonly used methods for generating scenarios and c) Demonstration of the feasibility and practicality of the proposed model for integrating into the economic dispatch framework.
The rest of the paper is structured as follows: Section II presents an overview of the scenario-generation methods for wind forecasting, followed by the presentation of the proposed spatio-temporal wind forecast model, and the ED formulation. Section III presents the solution methodology. Numerical tests of the integration of the spatial-temporal wind forecast with economic dispatch are provided in Section IV. Conclusions and future research are discussed in Section V.

Mathematical Model for Wind forecasting and Economic Dispatch
This section discusses the methods used to create wind power scenarios and also formulates the ED problem in which the scenarios generated are subsequently used. First, a review of commonly used methods such as time series (ARMA) and Monte Carlo simulations are made along with the Kantorovich distance-based scenario reduction technique. Then the spatio-temporal models based on Kriging method are presented. It begins by explaining their spatial form and then analyze the spatio-temporal model. The idea of a technique based on analogous for creating scenarios is also explained. Finally, an ED is formulated in which the objective function is to minimize both generation and reserve costs, subject to operational constraints of the system.

Overview of scenario-generation based methods
The first step is to characterize the wind speed from the database available. For this purpose, two of the most widely used methods are used to characterize wind speeds, such as Weibull distributions and ARMA models [43][44][45]. 6

of 25
The Weibull distribution is a two-parameter function commonly used to fit the frequency distribution of wind speed [44]. The Weibull probability density function can be expressed as where f (v) is the probability of observing wind speed v. Wind speed can be represented as a stochastic process whose characteristic parameters shall be estimated. The Weibull shape and scale parameters are denoted by k and c, respectively. k is dimensionless, and it indicates how peak the side under consideration is, while c has units of wind speed (m/s) and it shows how windy the site is. The cumulative distribution function is given by: The Auto-Regressive Moving Average (ARMA) statistical model is used to model random time series based on historical data, pattern identification, and parameter estimation [45,46]. This is a hybrid of auto-regressive (AR) and moving average model. The typical ARMA (p,q) model is expressed as: where Z t is the random time series to be modeled, p AR parameters φ 1 , φ 2 , ..., φ p , and q moving average parameters Θ 1 ,Θ 2 , ...,Θ q . The input ε t is normally distributed white noise with zero mean and standard deviation σ.
Once the wind speed has been characterized, scenarios are generated using some of the most commonly used methods for this purpose, such as Time-series-based ARMA models and Monte-Carlo simulations.
In Time-Series-based ARMA model, before starting the algorithm, it is necessary to fit the wind speed data to a known probability density function (PDF) and estimate the parameters of the ARMA model. Then, the number of desired scenarios N S and the period T for generating the scenarios are established. The process starts with s=1 and t=1. For each time period t, random samples for ε t are created, and the ARMA model is evaluated by obtaining random scenarios. A distribution transformation process is carried out to get wind speed scenarios [45]. The distribution transformation is described as: where F Z (Z) is the cumulative distribution function (CDF) of the randomly generated scenarios, and φ(V) is the CDF of actual historical wind speed data. Finally, the values of wind speed are transformed into power scenarios P(V) through the power curve associated with the turbine model installed in the wind farm. It is mathematically expressed as: where A is the swept area of the wind turbine rotor (m 2 ), ρ is the air density at the wind site (kg/m 3 ), C p (V) is the overall efficiency of the wind turbine, V is the wind speed (m/s), V cutin is the cut-in speed of the wind turbine (m/s), V cutout s the cut-out speed of the wind turbine (m/s), V rated is the rated speed of the wind turbine (m/s), and P max is the maximum power output of wind generator (MW). Monte Carlo simulation is a technique used to study how a model responds to randomly generated inputs. Wind speed data is used to generate scenarios of wind power. First, random sequences of the variable Y ∼ N (0, σ) are generated and then it is obtaining the CDF function F t . For horizon t, s realizations of the variable V = Φ −1 [F t ] are obtained by applying the inverse of the Weibull cumulative distribution function (CDF) Φ [45]. Finally, the scenarios are transformed to wind power by (5). The process is repeated until period T and number of scenarios N S are reached.
The uncertainty of renewable generation can be represented by a wide range of scenarios that can cause the problem of dispatch to become computationally intractable. Therefore, a mathematical technique is needed to reduce the set of scenarios that have been generated. The resulting set has statistical characteristics similar to those of the original set. A practical scenario-reduction technique based on the Kantorovich distance is used to obtain an appropriate number of scenarios [5,45]. It can be express as: where ω and ω are the scenarios, Q and Q are the probability distributions in initial scenarios set Ω and reduced scenarios set Ω . Two different algorithms can be used to solve the problem: backward reduction or the forward selection [5,45]. An interactive fast-forward selection algorithm is used to reduce the set of scenarios. It begins with an empty set, and each iteration, one by one, calculates a scenario ω k that minimizes Kantorovich distance. ω k is extracted from the initial set and included in the new set. The algorithm continues until a specified number of scenarios have been achieved. The probabilities of each non-selected scenario are then allocated to the nearest selected ω k . Finally, a reduced set with associated probabilities is created. The interactive fast-forward selection algorithm is detailed in [5,45].

Spatio-temporal based method
Wind speed, and thus wind power generation, depends on time and space. Using the spatial correlation between individual locations can significantly reduce errors in point forecasts and has the advantage of developing models that can produce forecasts at locations not included in the observation samples.

Spatial Kriging
Kriging is an interpolation method in which the measured values of the surrounding region are used in a linear combination to generate a prediction at an undetermined location. It is the optimal linear predictor of Z(s 0 ) at s 0 locations D s . Z(s 0 ) is determined at the unsampled location s 0 by the Kriging prediction where λ i and n are the weight assigned to the observation point Z(s i ) and the number of neighbors considered for the estimation of Z(s 0 ), respectively. The values of λ and k are calculated by minimizing the mean squared prediction error and in this equation, the first term can be written as where c ≡ (Cov(s 0 , s 1 ), ..., Cov(s 0 , s n )) T and Σ is a nxn matrix whose (i, j)th element is Cov(s i , s j ). Differentiating (10) with respect to λ and setting the result equal to zero gives the optimal coefficients λ * = Σ −1 c, and the optimal constant term is, therefore, k * = µ(s 0 ) − c T Σ −1 µ. The optimal linear predictor is given by Finally, the minimized mean squared prediction error, also called the Kriging variance, is: Note that the Kriging predictor requires the values of the Σ and c covariances to be known. In practice, both values can be calculated by means of a variogram or a covariance function. Usually, a parametric model is chosen for the covariance function Cov(s i , s j ) = C(s i , s j ; θ), where θ is a vector of unknown parameters to be determined from the data. Depending on how the expectation E(Z(s)) is modeled, different versions of predictors can be developed: Simple Kriging, Ordinary Kriging, and Universal Kriging.

Spatio-temporal Kriging
A Kriging predictor, as described above, is obtained to generate spatio-temporal predictions. For this case, an additional dimension, the temporal dimension t, is included. This section discusses how spatial Kriging prediction is generalized for use in spatio-temporal contexts.
The problem setting is as follows: Consider a spatio-temporal set of available data Z(s i , t i ) = z 1 (s 1 , t 2 ), ..., z n (s n , t n ) where z i is the value of variable z at location s and time t. The n points are the turbine sites in a wind farm and the values at (s i , t i ) are the wind speed measurements at these locations [47]. The objective is to make predictions at (s 0 , t 0 ), where no measurements are taken. The estimatorẐ(s, t) is define by: where λ i are weights for Z(s i , t i )), i = 1, ..., n. The kriging predictorẐ(s, t) is a linear combination of the observed points. The weights are based on the covariances among points in the sample and the covariances between sample points and the point to be predicted. It incorporates the covariance structure among the Z(s, t) into the weights for predictingẐ(s, t). The objective of the estimator is to minimize the error variance under the constraint of unbiasedness: As mentioned above, the kriging estimator is the optimal estimator, since that it minimizes the error variance. This quantity is also known as the mean squared prediction error (MSPE).
In order to compute estimates, we will need the covariances among all points and between each of the observed points and the point to be predicted. The usual way to obtain these is through a covariance function.
where (h s + h, h t + t) D s and (h s , h t ) T.
As mentioned above, the kriging weights can be calculated by means of a semivariogram γ(s, t) or a covariance function C(s, t). The semivariogram is more useful than the covariance function given that it can be calculated experimentally. A semivariogram is a visual representation of the covariance Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 26 October 2020 doi:10.20944/preprints202010.0513.v1 between each pair of points in the data sampled. For each pair of points in the data sampled, the semivariance is plotted against the distance or lag between them. The experimental semivariogram γ(s, t) can be written as: where the set N(h s , h t ) consists of the points that are within spatial distance h s and time lag h t of each other. The relationship between covariance function and semivariogram is: When the experimental semivariogram has been developed, it is fitted to a theoretical model of covariance C(s, t). The theoretical or model semivariogram is the distributional model that best fits the data. Spatio-temporal covariance function C(s, t) has a variety of forms that can be applied to the data [48,49], such as linear, spherical, exponential, Gaussian, or logistic models. The choice of a semivariogram model is user-defined, although statistical tools can help identify best-fitting models using several approaches, including the least square, maximum likelihood, and Bayesian methods. Before spatio-temporal covariance model has been identified, we can use it to derive semivariances at all locations and find the kriging weights.
Once the covariance matrices among points in the sample and between sample points and the prediction point have been calculated, the weights of the kriging estimator can be found. The idea is to find the estimatorẐ(s, t) = λ Z such that it satisfies the problem (14)(15). This is satisfied if ∑ λ i = 1 and the mean is stationary. To solve this optimization problem, Lagrange multipliers are used: where the last term in this sum guarantees that the unbiasedness constraint ∑ λ i = 1 is met through the Lagrange multiplier λ. Then partial derivatives of L are developed with respect to the weights λ i , set these equal to zero, and resolved. The weights can be found from the following system of linear equations (kriging equations): where , and λ is a Lagrange multiplier that appears due to the unbiasedness constraint ∑ λ i = 1. When the weights are calculated, the predicted value of Z(s, t) is estimated at an unsampled spatio-temporal location (so, to) by (13). The value of the predicted point (so, to) is equal to the sum of the value of each sampled point times that unique weight of the point (s, t). The covariance matrix used to calculate the weights differs slightly depending on the type of kriging being performed.
The analog method consists of examining the wind speed database to look for similar information to predict future occurrences. Each data in the historical archive is considered to be a possible scenario [50,51]. The method is to compare the current day forecast with the database in order to mine similar days or analogs as scenarios.

Stochastic Dispatch Problem Formulation
Economic dispatch is the short-term determination of the optimal output of a number of electricity generation facilities, to meet the system load, at the lowest possible cost, subject to transmission and operational constraints. In this subsection, the economic dispatch problem of a power system is formulated mathematically, taking into account the following: the use of a DC optimal power flow model, the cost functions are linear, the inter-temporal constraints, such as ramping limits, are not included in the formulation, the uncertainty is presumed to be exclusively generated by wind generators, the uncertainty associated with wind generation can be effectively modeled by a finite set of scenarios and their probability of occurrence, and conventional units are considered to be entirely dispatchable from their minimum to their maximum capacity.
In this paper, the objective function of the ED model is to minimize the overall expected system costs, which consists of the day ahead dispatch cost plus the expectation of the balancing operation costs. The aim of the stochastic formulation is to assess the generation scenarios N s that have been created and to ensure that they can be included in the ED framework. Based on [52], the objective function of the model writes as follows: where P g,t , P wcut w,s,t , P dcut b,s,t , r g,s,t up , r g,s,t down , C g (.), C up g (.), C down g (.), α, β, π, G, B, W, N s , and T are the dispatched power by the deterministic generator g at t, wind curtailment for scenario s at t, load shedding for scenario s at t, up reserve by the generator g for the scenario s at t, down reserve by the generator g for the scenario s at t, the generation cost function of power plant g, the reserve costs (up and down) of power plant g, penalty cost for wind curtailment, penalty cost for load shedding, probability of occurrence of scenario s, set of conventional generators, set of buses, set of wind generators, set of scenarios, and period of time, respectively.
In order to make the solution feasible, certain constraints have to be fulfilled. Some basic constraints that need to be set are power balance equations, bounds of reserve and energy, network constraints, and bound of generators. The details of each constraint are described as follows.
• Active power balance constraints: The balancing actions must ensure a balance between generation and demand in any possible scenario s. The constraint to be satisfied can be given as where P w,s,t , δ n,s,t , δ m,s,t , x nm , Λ n , and N are the wind power scenario s of generator w at t, the voltage angle of node n under s at t, the voltage angle of node m under s at t, the reactance of the line n − m, the set of nodes directly connected to node n, and the set of nodes, respectively. We define bus 1 as the reference node by setting δ 1,s,t to 0. • Transmission capacity constraints: For physical reasons, the amount of power transmitted through a power line n − m has a limit. This limit is justified by thermal or stability considerations. Therefore, a line must be operated so that this transport limit is not exceeded under any circumstances. This is formulated as − P max nm ≤ δ n,s,t − δ m,s,t x nm ≤ P max nm , ∀n Λ n , ∀m Λ n , ∀s, ∀t (23) where P max nm is the transmission capacity of the transmission line n − m. • Generation output constraints: Each unit is designed to work between the minimum power capacity and the maximum power capacity. The following constraint ensures that the unit is within its respective rated minimum and maximum capacity.
where P min g and P max g are minimum output of generator g and maximum output of generator g. • Wind curtailment constraint: The amount of wind power production that is curtailed under each scenario s must be lower than or equal to P w,s,t : P w,s,t wcut ≤ P w,s,t , ∀w, ∀s, ∀t • Load shedding constraint: The amount of load that is shed in each scenario s has to be lower than or equal to the demand value: where γ is the maximum allowable load shedding percentage. • Reserve constraints: The reserve capacity for the balancing energy is shown below: P g,t + r g,s,t up ≤ P g max , ∀g, ∀s, ∀t (27) P g,t − r g,s,t down ≥ P g min , ∀g, ∀s, ∀t (28) r g,s,t up ≤ R g up , ∀g, ∀s, ∀t (29) r g,s,t down ≤ R g down , ∀g, ∀s, ∀t (30) where R up g and R down g are the upward reserve capacity and downward reserve capacity, respectively. • The quantity of reserves, generation and demand must be non-negative: r g,s,t up , r g,s,t down ≥ 0, ∀g, ∀s ∀t (32) P w,s,t wcut ≥ 0, ∀w, ∀s, ∀t P b,s,t dcut ≥ 0, ∀b, ∀s, ∀t (34) where Ξ = P g,t , r g,s,t up , r g,s,t down , P w,s,t wcut , P b,s,t dcut , δ n,s,t , δ m,s,t is the set of decision variables ∀g, ∀w, ∀b, ∀n, ∀m, ∀s, ∀t.
The resulting ED problem is given by: minimize (21) (35) subject to: Power Balance constraint : (22) Transmission capacity constraint : (23) The stochastic programming model (35)- (42) can be used to determine the effect of uncertain wind power generation on the expected value of total system operating costs.

Proposed Methodology
The proposed scenario creation methodology consists of a three-step process: first, the spatial covariance structure of the sampled points is determined by fitting a semivariogram; second, weights derived from this covariance structure are used to interpolate values for unsampled points; and third, the forecast values are used to generate power scenarios by the analogous method. Figure 1 summarizes the proposed methodology.

Empirical variogram calculation and parametric fitting
Initially, assuming that the database has a set of spatio-temporal data measured in T instants of time on a spaced grid, as shown in the Figure 2. The possible sets of spatial distances h s between the pairs of points considered are calculated, whereas the lag times h t are t 1 , t 2 , ..., t n . Thus, combining the distance and the time lags, spatio-temporal distances are obtained. Also, N(h s , h t ) is calculated for such distances. For this purpose, it is considered: Computing the semivariances for each spatio-temporal distance h s , h t , the empirical spatio-temporal semivariogramγ(s, t) the empirical spatio-temporal semivariogram is obtained by (17). After calculating the empirical semivariogram, it is important to model it using a parametric function which will be used in the Kriging method. The most common covariance (and semivariogram) functions are: 1. The separable covariance model assumes that the spatio-temporal covariance function is represented as the product of a spatial and temporal term: 2. The product-sum covariance model with k > 0 3. The spatio-temporal metric covariance model 4. The sum-metric covariance model: 5. The exponential covariance model: The coefficients of the parametric functions are calculated by applying the weighted least squares (WLS) method and minimizing the difference between the empirical and parametric semivariograms.

Estimating the weights and derivation of the Kriging estimator
The weights are based on the covariances among points in the sample and the covariances between sample points and the point to be predicted. With the function fitted in the previous step, these covariances are calculated. The covariance matrix among the observed points C ij is: The covariances C ij are considered to be parameters, which are estimated by the modeled semivariogram. The points in the lower triangular portion of this matrix are omitted due to the symmetry of the covariance matrix. The covariances between sample points and the prediction point C i0 are: Recall that in the matrix form, the kriging equations can expressed by (20), where after inverted the partitioned matrix, the weights can be written as: The predicted value of Z at the prediction point is:

Scenario generation
The analog approach searches back in time to find similar forecasts and uses matching observations from such analogous dates directly as scenarios. This research proposes and tests a simple analog approach based on the spatio-temporal wind forecasts created by the Kriging method. The method is based on the analysis of distance measurements between the current forecast and those available in the database [50,51].
whereω t , ω t , and T are wind speed forecast vector, historical archive of observation vectors, and period respectively. Also, i = 1, 2, ..., columns( ω s ) Once the evaluation has been carried out, data with less distance are selected as scenarios. The scenarios are transformed to wind power by (5). The process is repeated until period T and number of scenarios N S are reached.

Database description
The data used in this analysis consists of one year of spatio-temporal measurements in five randomly selected turbines on an at a wind farm, between 2010 and 2011 [47]. The data contains precise hourly wind speeds and turbine coordinates. To convert wind speed into wind power, the

Results of scenario generation and reduction
First, the data is adjusted by regression to the Weibull distribution function. The estimated scale and shape parameters for the Weibull distribution are 10,927 and 2,14067, respectively. The distribution parameters are used to create scenarios based on the Monte Carlo method. The set of available data is then divided into two parts. One part is used to estimate the ARMA model parameters, and the other part is used for validating the ARMA model. The data are used to estimate the set of parameters and the structure of the model. The resulting model is used for scenario generation. The model parameters for the selected structure are c = 0.6282, p = 0.9178, and q = −0.1058.
Concerning spatial-temporal modeling, a set of points (triangle) and a set of wind speeds is considered, as shown in the Figure 3. Also, the value of wind speed at the new location s 0 (circle) is required to predict. To compute the estimates, the spatial and temporal distances to be used for the  construction of the semivariogram are calculated. Next, the experimental semivariogram based on equation (17) is developed in order to establish a theoretical model for the covariance function between the spatial-temporal data. Once the semivariogram has been constructed, the data is adjusted to the theoretical model. Based on the adjustment made, it was decided that the most acceptable model was the exponential curve. Figure 4 depicts the experimental semivariogram, and the model fitted. The  corresponding exponential covariance function has the form that can be defined by: where c 0 , c 1 , k, and a are parameters obtained by the fitting and h and u are spatial and temporal data, respectively.
Once the model of the semivariogram is adjusted, the kriging weights are calculated by solving equation (51): Note that the weights sum to 1.
Now the predicted value of Z at the unobserved point can be calculated by (52). Figure 5 depicts the prediction developed (black line) at the unobserved point. The next step is to use the Monte-Carlo, ARMA, and Spatio-temporal methods to generate several scenarios for each hour of the next 24 hours. Once the wind scenarios are generated, the Kantorovich distance algorithm is used to reduce the number of realizations. The sum of probabilities of all generated and reduced scenarios is always one at any given time. Figure 6 depicts the wind speed scenarios reduced (only 10 scenarios are shown). From this set, wind power scenarios are derived. Equation (5) is used to calculate output power.
Time series metrics are used to discuss scenario quality: scenario mean, mean absolute error and root mean square error are evaluated. Table 2 gives forecast error statistics for the scenarios. We measured the overall error characteristics using two values: the mean absolute error (MAE) and the root-mean-square error (RMSE). In the Table 2, it can be seen that scenarios based on spatio-temporal model have lower error values.  Concerning computational cost, the three methods were used to generate 300 scenarios. All the numerical experiments were performed using Matlab 9.2 on a desktop with a Core-i7 processor, 2.60 gigahertz, 12 gigabytes RAM and GPU GeForce GTX 960M. In short, the ARMA-based method generates scenarios from a prediction model derived from the database, and also uses the parameters k and c of the associated Weibull distribution. Average computational time for this algorithm is 1.182961 seconds. The Monte Carlo method generates scenarios using the parameters k and c of the Weibull distribution. Average computational time for this algorithm is 0.29834 seconds. Finally, the spatio-temporal method generates scenarios through analogous in the database. Average computational time for this algorithm is 8.9813 seconds.
The results show that the proposed algorithms generate scenarios that fit appropriately to the database used. Based on the error analysis, it can be verified that the scenarios based on spatio-temporal and analog method have a better fit. Finally, the Monte Carlo method presents lower computational cost and reduced number of parameters for the generation of scenarios.

ED Results
The modified IEEE-3 bus and IEEE-24 test systems are proposed for the study of the ED problem with the use of spatio-temporal Kriging and common methods for modeling the generation of wind.
The ED model is implemented in the Julia programming environment and the simulations were performed on a laptop with a Core-i7 processor, 2.60 gigahertz, 12 gigabytes RAM, and GPU GeForce GTX 960M.
In order to analyze and compare the impact of different scenarios, a test was performed to determine the appropriate number of scenarios to use in the dispatch problem by increasing the number of scenarios and observing the behavior of the objective function. The corresponding numerical results obtained are shown in Figure 7, where the value of the cost function for each method and the corresponding number of scenarios are depicted. For a set of scenarios less than 120, the cost function has variations, while for values greater than 120, it can be verified that the changes in the value of the objective function are negligible. Based on results of the figure, it can be established that a set of 120 scenarios is adequate for uncertainty assessment due to renewable generation.  The methods are compared in two aspects: costs of dispatch and computational time.

Case 1: IEEE-3 bus test system
The modified IEEE-3 bus test system is used to demonstrate the effectiveness of proposed method. The system has three conventional power plants located in bus 1, bus 2, and bus 3 and one wind generator located in bus 2. All technical specifications for generation units and the test system are reported in [54]. The installed capacity of the wind farm is 30 MW for this case. The penalty coefficient α (underestimation cost coefficient) for wind curtailment is 80$/MW and the penalty coefficient β (overestimation cost coefficient) for load shedding is 160$/MW [10]. The maximum allowable load shedding percentage γ is set to be 5% [10].
The total cost of the methods after running the optimization problem described in section 2.3 by equations (35)- (42), for an arbitrary selected day, are listed in Table 3. The table details the cost of dispatch before and after the actual wind power. As shown, the cost of the deterministic approach is lower than that of the other methods before the real wind power is dispatched. However, with actual wind power, it is found that the cost of the Spatio-temporal method is lower than the Deterministic, Monte Carlo, and ARMA costs. It demonstrates that scenario-based methods perform better than deterministic ones in the presence of uncertainty. In general, the proposed method has better overall savings. the generators is 150 MW, 150 MW, and 80 MW, respectively. It is observed that the most economical generators P 1 and P 3 are those that participate mainly in the dispatch for all cases. In all three cases, generator P 2 , which is the most expensive, dispatches its minimum capacity. Also, generator P 3 , which is the cheapest, dispatches its maximum capacity. The tables also present the amount of wind power used to complete the dispatch. For this case study, when using the ARMA method, 4.9% of the total demand is met with wind, in the case of Spatio-temporal 5%, and for Monte Carlo 3.8%. Table 4. Demand of the system The computational time of methods are listed in Table 8. For the DO method, the computation time is about 18.975s, the Spatio-temporal is about 33.221s, for Monte Carlo is about 31.720s, and for the ARMA model is about 30.6883s. The computation time of the Spatio-temporal method is longer than that of the deterministic, ARMA and Monte Carlo based method.

Case 2: IEEE-24 bus test system
The modified IEEE-24 bus test system is used with twelve conventional generators and one wind generator located in bus 3. All technical specifications for generation units and the test system are reported in [55]. The installed capacity of the wind farm is 500 MW for this case. The penalty coefficient Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 26 October 2020 doi:10.20944/preprints202010.0513.v1 α (underestimation cost coefficient) for wind curtailment is 80$/MW and the penalty coefficient β (overestimation cost coefficient) for load shedding is 160$/MW [10]. The maximum allowable load shedding percentage γ is set to be 5% [10]. The total the total cost of the methods after running the optimization problem described in section 2.3 in equations (35)- (42), for an arbitrary selected day, are listed in Table 9. As shown, the cost of the deterministic approach is lower than that of the other methods before the real wind power is dispatched. However, with actual wind power, it is found that the cost of the Spatio-Temporal method is lower than the Deterministic, Monte Carlo, and ARMA costs. The computational time of methods are listed in Table 10. For the DO method, the computation time is about 23.2292s, the Spatio-temporal is about 429.401s, for Monte Carlo is about 410.833 s, and for the ARMA model is about 448.2883s. The computation time of the Spatio-temporal method is longer than that of the deterministic and Monte Carlo method and shorter than that of ARMA based method. Finally, an estimation of the behavior of the total cost of the ED problem is carried out over 30 days. Figure 8 depicts the cost behaviour for DO (black), Monte Carlo (green), ARMA (magenta), and Spatio-temporal (red) approaches. Scenario-based approaches are found to have a lower cost than the deterministic method under uncertainty.  Figure 8 depicts that the Spatio-temporal method has the lowest costs most of the days tested. However, there are also points, as in t 10 , where the ARMA method has the lowest cost. Table 11 reports the average total cost of each case and also indicates, in percentage, the cost decrease obtained when compared to the deterministic method. As seen in the Table, with the spatio-temporal approach, the highest percentage of cost reduction (1.2301%) is achieved.

Conclusions
This research examined three wind scenario creation methods. The methods allow the operator to leverage historical data to create many scenarios. Also, this paper described a wind prediction model using a spatio-temporal model based on the Kriging method. Besides, from the spatio-temporal prediction, wind scenarios were generated using analogs. These scenarios were then converted into wind power, considering the power curve of the wind turbines. An economic dispatch model was then established so that the generated scenarios can be integrated into the formulation. Finally, the problem was implemented and evaluated on modified IEEE 3-bus and IEEE 24-bus systems. The following conclusions were established based on the studies: • The proposed approach for spatio-temporal modeling is computationally efficient and can be applied to predict values on unobserved points, from data measured around the unmeasured point. In this work, the data of the wind speed and the location of these measurements are used. • Because the weights of the kriging estimator depend on the modeled semivariogram, kriging is very sensitive to the mis-specification of the semivariogram model. • In general, interpolation accuracy by kriging is limited if the number of observations sampled is small. The data are limited in spatial scope, or the data are not broadly spatially correlated. In these cases, it is challenging to generate an experimental semivariogram.
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 26 October 2020 doi:10.20944/preprints202010.0513.v1 • The number of scenarios is an important parameter and can have a significant impact on cost function of the ED problem. A good value for N s is likely to depend on the method of scenario creation. For the presented approach, the stabilization of the objective function occurred in values close to 120 scenarios, as can be seen in Figure 7. • An economic dispatch model that considers wind generation in its formulation is presented in this research. Wind generation is created using a Kriging-based spatio-temporal model and other methods commonly used for this purpose. The proposed approach's effectiveness is demonstrated by the simulation of the modified IEEE 3-bus and the IEEE-24 bus systems. • The results show that spatio-temporal and analogs models for generating scenarios can be integrated efficiently in the ED problem framework and reduce costs.The reduction in operating costs when using spatio-temporal scenarios is 1.2301%.
Scenario creation and ED problem under uncertainty remains an open area for research. Many of the methods presented in this paper could be improved upon to produce better scenarios and to attain ED solutions with lower costs in less time. Scenario creation methods have the potential to be applied to other sources of uncertainty in the power system. Besides, we can further take advantage of spatio-temporal and analog-based approaches and provide different strategies for modeling uncertainties on both the supply and demand sides. Finally, savings due to ED under uncertainty are likely to vary from system to system and with wind energy penetration. The number of scenarios used may play a large role in the dispatch costs. Evaluations of the potential for savings should be made for other systems, especially those with flexibility requirements. Further work will focus on the application in a real power grid.