A Framework for Real-Time Optimal Power Flow under Wind Energy Penetration

Developing a suitable framework for real-time optimal power flow (RT-OPF) is of utmost importance for ensuring both optimality and feasibility in the operation of energy distribution networks (DNs) under intermittent wind energy penetration. The most challenging issue thereby is that a large-scale complex optimization problem has to be solved in real-time. Online simultaneous optimization of the wind power curtailments of wind stations and the discrete reference values of the slack bus voltage which leads to a mixed-integer nonlinear programming (MINLP) problem, in addition to considering variable reverse power flow, make the optimization problem even much more complicated. To address these difficulties, a two-phase solution approach to RT-OPF is proposed in this paper. In the prediction phase, a number of MINLP OPF problems corresponding to the most probable scenarios of the wind energy penetration in the prediction horizon, by taking its forecasted value and stochastic distribution into account, are solved in parallel. The solution provides a lookup table for optional control strategies for the current prediction horizon which is further divided into a certain number of short time intervals. In the realization phase, one of the control strategies is selected from the lookup table based on the actual wind power and realized to the grid in the current time interval, which will proceed from one interval to the next, till the end of the current prediction horizon. Then, the prediction phase for the next prediction horizon will be activated. A 41-bus medium-voltage DN is taken as a case study to demonstrate the proposed RT-OPF approach.


Introduction
The dramatic increase of renewable energy penetration represents a significant challenge in the operation of energy distribution networks (DNs).In particular, wind power generation is intermittent, i.e., the DN operator has to correspondingly update the operation strategy.Therefore, it is highly desired to carry out this task by an online optimization framework.However, the optimization problem to be solved is usually high-dimensional and complicated when a large network with a detailed nonlinear model as well as mixed-integer variables is considered.Thus the computation time can be much higher than required for reacting to the fast changes of the wind power generation.Even by using advanced optimization algorithms combined with modern computation facilities, the computation time can still be too high to achieve this target.Furthermore, the feasibility of the solution should be ensured within a specified sampling time.Therefore, a computation framework addressing these difficulties needs to be developed for the implementation of real-time optimal power flow (RT-OPF) under wind energy penetration.
Energies 2017, 10, 535 2 of 28 OPF has been widely used for operation planning of power networks with only conventional generators [1][2][3].OPF with renewable energy generation (REG) was taken into account in [4] and an active-reactive OPF in active DNs introduced in [5].However, these OPF problems are solved offline for planning or scheduling network operations, and therefore, intermittent wind power generation, which is randomly deviated from the pre-assumed offline situation, cannot be handled.
The concept of RT-OPF was first presented in [6], based on a linear model and a quadratic optimization method.Without considering REG, a real-time implementation of optimal reactive power flow was reported in [7] where the sampling time was between 15 min up to four hours.Using neural networks trained on several scenarios of the uncertain demand and REG, the authors in [8] presented a RT-OPF for a 23-bus radial DN with two wind generators, achieving a sampling time less than three minutes.A real-time energy management system was proposed for mitigating pulse loads in hybrid microgrids by considering REG [9] and uncertainties in plug-in electric vehicle charging parks [10].Reported in [11] a RT-OPF was implemented for a laboratory microgrid in which a switching scheme for solving the OPF problem was proposed.If the AC optimization problem does not converge, a DC optimization, based on a linearized model, is switched to facilitate the convergence.
A RT-OPF method integrating storage devices and wind generation was presented by using a linear model predictive control (MPC) [12] and tested on a 14-bus transmission network with a sampling time of 5 min.A MPC-based control approach was reported in [13] to dynamically balance, in real-time, demand and supply in a DN with integrated photovoltaic generators and energy storage systems.The objective in [13] was to alleviate the effects of fluctuations in demands and REG.The MPC formulations in [13] were modified in [14] to enable day-ahead power reference tracking functionality.In addition, reference [14] provides the theoretical guarantees on the stability of the control strategy in a simplified setting.The simulation scenarios consider the systems under real working conditions by taking into account the forecasting errors of photovoltaic plants.
Curtailment of REG becomes necessary because of transmission congestion, lack of transmission access, excess renewable generation during low load period, and voltage or interconnection issues [15].A real-time curtailment reduction scheme was implemented on real substation hardware to manage an active DN [16] by a method of constraint satisfaction and OPF.A disadvantage of a closed-loop scheme is that the system will react only when a constraint violation takes place [16].The authors in [17] proposed a risk-based RT-OPF scheme by considering a future time horizon in which the REG is described as stochastic variables.The control decision was aimed at minimizing the curtailment while ensuring network constraints within a pre-specified risk level.
Recently, a gradient algorithm for RT-OPF was proposed in [18] which is in principle based on the idea of the quasi-sequential approach proposed in [19] and improved in [20].The problem is decomposed into a simulation stage where the model equations are solved to provide values of the state variables and an optimization stage where the controls are solved by a NLP solver.Barrier terms are introduced in the objective function to penalize the inequality violations.However, it is well known that a barrier method needs a feasible guess point which may be not easy to obtain.The main advantage of the approach is that the intermediate iterates of the algorithm satisfy the power flow equations and consequently can be utilized in real-time to track evolving network conditions.Although this ensures a low computation time, it should be noted that an intermediate solution is not optimal.
Based on a linear model, a data-driven hourly real-time power dispatch was proposed by using a probabilistic optimization approach under uncertain renewable penetration [21].The dispatchable range is at first determined, based on which the operating base points for the conventional generators are computed for the next hour.Then the operation strategy will be corrected when the observations of the real renewable energy are available.
The literature on real-time power flow shows a variety of approaches to control slack bus voltage.An automatic adjustment of transformer and variable phase-shifting transformers using Newton-Raphson power flow in transmission networks was reported in [22].Incorporating the value Energies 2017, 10, 535 3 of 28 of slack bus voltage into the calculations, reference [23] presented a hybrid real-complex current injection based load flow formulation.It is noted that the discreteness of the slack bus voltage was not considered in [22,23].The obtained values of slack bus voltage are then rounded to their nearest discrete values.In contrast, a sensitivity-based current injection power flow was proposed in [24] for the simulation of local voltage controllers in a discrete manner.By calculating load flow, the studies in [22][23][24] give the state of the network based on their input scenarios.However, the determination of optimal reference values of the slack bus voltage in DNs plays an important role in planning the operation of power systems.Therefore, in this paper, it is tried to determine the optimal discrete reference values of the slack bus voltage, in addition to other decision variables.
To reduce the computation effort, a RT-OPF approach considering variable REG and demand between two consecutive scheduling intervals (e.g., 10 min) was used in [25].The intervals are divided into several subintervals (e.g., 1 min) for which the load and REG forecasts are available.OPF is carried out only once at the start of each scheduling interval aiming at optimizing the participation factors of conventional thermal generators while holding the technical constraints for all the subintervals.This was done by incorporating the generation costs of all the subintervals.The resulting controls can be considered as 'best-fit' participation factors [26,27] which are then used for modifying the controls for each subinterval according to the forecasted data for solar, wind generations and loads.Since the curtailment of REG and/or reverse power flow to an upstream network were not considered in [25], the approach is applicable for the cases that the total REG is less than the total demand plus losses in the grid, or the cases with energy storage systems like in [28,29].
In summary, many previous studies have been made on RT-OPF in DNs.However, these available methods do not simultaneously consider online optimization of the wind power curtailments of wind stations (WSs), the discrete reference values of the slack bus voltage (leading to a mixed-integer nonlinear programming (MINLP) problem) and the variable reverse power flow to an upstream network.To bridge this gap, this work extends our previous studies [30,31] and develops a new techno-economic RT-OPF framework to ensure, in real-time, the feasibility of control strategies to be realized to the grid.The contributions of this paper can be summarized as follows:

•
A novel RT-OPF framework is developed to address the conflict between the fast changing wind power and the slow optimization computation and consequently to realize an online optimization of energy systems in a very short sampling time;

•
Discrete reference values of the slack bus voltage, wind power curtailment of WSs, and reverse power flow are considered simultaneously, leading to a MINLP problem; • A scenario generation method is integrated in the RT-OPF framework to represent uncertain wind power for the prediction horizon, which leads to a set of uncoupled MINLP problems solved by parallel computing.
The scope of our prediction-realization approach is demonstrated using a 41-bus medium-voltage DN with two WSs located at different buses and positions and thus with different wind speed.The remainder of the paper is organized as follows: Section 2 describes and formulates the OPF problem.The scenario generation method is described in Section 3. The solution framework is presented in Section 4. Results and discussions of a real case study are given in Section 5.The paper is concluded in Section 6.

Problem Formulation
The ultimate goal of the RT-OPF framework proposed in this study is to compute optimal operation strategies for DNs which will autonomously be updated according to spontaneous changes of energy penetrated from WSs.Thus, the updating time interval (sampling time) T s should be kept as short as possible.However, due to its high complexity, the computation time T OPF needed to solve the optimization problem can be much higher than the sampling time.To address this conflict, we employ the forecasted data of wind energy which are available in advance of a future time horizon T p .
Energies 2017, 10, 535 4 of 28 In this paper, this forecasted time horizon is called a prediction horizon.Since the prediction horizon T p is usually higher than the sampling time T s , we can divide the prediction horizon into M sampling times, i.e.,: In this study, we assume that the total computation time T OPF is smaller than the prediction horizon T p .Under this assumption, a prediction-realization approach for RT-OPF will be developed in this paper.In the prediction phase, the optimization problem is solved in advance for a number of probable wind energy scenarios, based on the forecasted data in the prediction horizon and its probability density function (PDF), leading to a lookup table for optional optimal operation strategies.In the realization phase, the actual wind energy data are successively available from one sampling time to the next.In each sampling time, the actual data will be compared with the predefined wind energy scenarios and an optimal operation strategy corresponding to the nearest higher scenario will be selected from the lookup table and realized in the network.In this way, an online update of the operation strategy according to the spontaneously changing wind energy generation is carried out.
In this section, the optimization problem to be solved for each prediction horizon, during the prediction phase, is defined.To explain the complex problem in a clearer way, we define at first a general optimization problem for OPF.A detailed and concrete problem definition of the RT-OPF is given in Section 4. For a prediction horizon, i.e., t ∈ [kT p , (k + 1)T p ], the OPF problem of a DN with wind energy penetration is defined as: where x is the vector of state variables, u and l are the vectors of continuous and integer decision variables, respectively.Relating to the OPF under consideration, the state vector x comprises voltages of the PQ buses, active and reactive power at slack bus, and power flows in the feeders, the continuous control vector u consists of curtailment factors of each WS, and the discrete control vector l denotes the reference values of slack bus voltage.The vector ξ represents random variables of wind energy of each WS which will be generated in the prediction horizon.In this paper, these random variables are regarded as being stochastically distributed with a known PDF ρ(ξ).Therefore, the optimization problem expressed in Equation ( 2) is a MINLP problem under uncertainty.
In fact, power demand also needs to be considered as an uncertain parameter, but in comparison to the wind power, its demand value is more predictable in the application of online optimization.Many approaches have been developed for very short-term load forecasting aiming at prediction ranges of a few minutes to an hour [32][33][34][35].In [32], it is shown that the mean absolute percentage error (MAPE) can be less than 0.2% for a 120-s prediction horizon.Furthermore, MAPE values of 3.23% and 2.44% were obtained in [35] for the prediction of 30-min ahead individual household load and aggregation load, respectively.Considering the accuracy of the forecasts in the abovementioned studies, the variation of demand in short time slots (e.g., 120 s) is not considered in this paper.Thus the forecasted demand for each prediction horizon is used in our RT-OPF framework, but it may change from horizon to horizon.

Scenario Generation
It is necessary to describe the uncertain vector ξ in Equation (2) in advance of each prediction horizon.To do this, a set of wind power scenarios for the prediction horizon representing the stochastic behaviors of the wind power need to be generated, for which we need the PDF.The wind power scenarios are generated within the range [0, 1] pu where 1 pu corresponds to the rated power value Energies 2017, 10, 535 5 of 28 P w.R (n w ).N s scenarios are generated for each WS.We define N s − 1 intervals for the wind power P w (n w , n s ), where Pr is the probability operator.In this way, an equal probability between two adjacent scenarios is ensured.It is noted that (3) can be applied to any type of continuous bounded distribution as far as the area under its PDF curve equals one.Beta distribution is suggested to be highly suitable to represent the forecast errors of wind power ( [36][37][38]).Although Beta distribution cannot model the fat tail of the forecast errors perfectly [38], due to its variable kurtosis [38], it is still more suitable than the Gaussian distribution and gives reasonably accurate results [38,39].Beta PDF has been used in many recent studies [40][41][42][43] and therefore is chosen in this paper to represent wind power forecast errors.The PDF of the Beta distribution is defined as [36]: where α(n w ), β(n w ) are the first and second shape parameters of the Beta distribution.The corresponding probability distribution function is expressed as: As a result, the probability in the interval between 0 and scenario n s is n s −1 N s −1 and then the scenarios can be generated by: The parameters α(n w ), β(n w ) can be determined by [38]: where P w.M (n w ) and σ w (n w ) are the values of the mean and standard deviation of the wind power generation, respectively.For the RT-OPF, we take the forecasted values available before each prediction horizon, as P w.M (n w ).It is noted that the value of σ w (n w ) cannot be forecasted, but it can be evaluated from historical data.Figure 1 illustrates the generated scenarios for three forecasted wind power values when N s = 7 and σ w = 0.1.It can be seen that this scenario generation method leads to more scenarios near the mean value.The scenarios generated are the boundaries of intervals.In this way, the scenarios cover the whole range [0, 1] (i.e., from zero to the rated power).Consequently, our RT-OPF can deal with any actual wind power generation (also see Section 4.2).It is noted that the scenarios here are not those from the Latin hypercube sampling method where the scenarios are randomly generated inside each interval [44].In addition, it is noted that this scenario generation method is a significant improvement to that used in [30,31] where the scenarios were generated based on a constant width of intervals which cannot cover the whole range when a limited number of scenarios is chosen.For a power system with N w WSs, the total number of scenario combinations N c is: Thus, we need to define the power scenarios for individual WSs ( , ), 1, , , and the combinations of the scenarios for all WSs, ( ), 1, , , each combination is a vector), respectively.According to (9), when the number of WSs increases, the number of scenario combinations increases exponentially.In this case, parallel computing seems not reasonable to address the computational problem which is expectedly solved by the next generation of hardware and software, considering the rapid advancement of the computer technology.All scenarios finally generated for the power system is given in Table 1.The rules based on which the scenario combinations are sorted from the first row to row c N is further cleared by graphical examples in Figure 2. It is noted that the scenarios are listed from the highest to the least wind power values, due to the reason described in Section 4.2.To the best of the authors' knowledge, the integration of such a scenario generation method into a RT-OPF framework has not yet been considered.Thus, we need to define the power scenarios for individual WSs P w (n w , n s ), n s = 1, • • • , N s , and the combinations of the scenarios for all WSs, P w (n c ), , each combination is a vector), respectively.According to (9), when the number of WSs increases, the number of scenario combinations increases exponentially.In this case, parallel computing seems not reasonable to address the computational problem which is expectedly solved by the next generation of hardware and software, considering the rapid advancement of the computer technology.All scenarios finally generated for the power system is given in Table 1.The rules based on which the scenario combinations are sorted from the first row to row N c is further cleared by graphical examples in Figure 2. It is noted that the scenarios are listed from the highest to the least wind power values, due to the reason described in Section 4.2.To the best of the authors' knowledge, the integration of such a scenario generation method into a RT-OPF framework has not yet been considered.Thus, we need to define the power scenarios for individual WSs ( , ), 1, , , and the combinations of the scenarios for all WSs, ( ), 1, , , each combination is a vector), respectively.According to (9), when the number of WSs increases, the number of scenario combinations increases exponentially.In this case, parallel computing seems not reasonable to address the computational problem which is expectedly solved by the next generation of hardware and software, considering the rapid advancement of the computer technology.All scenarios finally generated for the power system is given in Table 1.The rules based on which the scenario combinations are sorted from the first row to row c N is further cleared by graphical examples in Figure 2. It is noted that the scenarios are listed from the highest to the least wind power values, due to the reason described in Section 4.2.To the best of the authors' knowledge, the integration of such a scenario generation method into a RT-OPF framework has not yet been considered.

Prediction Phase
The task of the prediction phase is to solve the OPF problems corresponding to the N c scenario combinations of wind power (see Table 1) for each prediction horizon.The active and reactive power demand at bus i, denoted as (P d (i), Q d (i)) as well as the active and reactive energy prices (Price P , Price Q ) are assumed to be constant in the short prediction horizon.But they may change from horizon to horizon.In addition, the power system model/structure is considered to be as in [5,45] and fixed in the prediction phase.The OPF problem defined in Equation ( 2) is formulated here in more detail for each scenario combination n c , ( where: The objective function in Equation ( 10) aims to maximize the total revenue from the wind power f 1 (n c ), and meanwhile to minimize the total cost of the active energy losses in the grid f 2 (n c ), the cost of the active energy at slack bus f 3 (n c ), and the cost of the reactive energy at slack bus f 4 (n c ).
Here, P loss (n c ) is the grid total active power losses [5] for scenario combination n c .P w (i, n c ) is the active power of WS located at bus i for scenario combination n c .P S (n c ) and Q S (n c ) are the active and reactive power injected at slack bus, respectively (i.e., the imported active and reactive energy from an upstream high-voltage (HV) network).The vector of discrete decision variables l, in Equation ( 2), consists of slack bus voltage V S (n c ), representing the controller reference of tap positions of the on-load tap changer.The vector of continuous decision variables u includes the curtailment factors of wind power for each WS, β w (i, n c ).Here, (0 ≤ β w (i, n c ) ≤ 1), where β w = 1 when no curtailment and β w < 1 otherwise [5].
The objective function of Equation ( 10) is subject to the active and reactive power flow equations at the buses: where f P (n c ) and f Q (n c ) are the network active and reactive power functions [5] for scenario combination n c , respectively.P S (n c ) and Q S (n c ) are the active and reactive power terms included only for slack bus, respectively.In addition, the following inequality constraints should be held: Bounds of active and reactive power at slack bus: Voltage bounds of buses: Feeder sections limits: and the limits of the curtailment factors: In Equations ( 18) and ( 19), the parameters γ Ps,rev and γ Qs,rev define the percentage of the allowable reverse active and reactive power to an upstream HV network.For instance, if γ Ps,rev = 1, active power exported to the HV network is fully allowed and if γ Ps,rev = 0, no active power export is allowed.In addition, based on Equations ( 21) and ( 22), for scenario combination n c , an optimal value of slack bus voltage V S (n c ) is obtained by selecting the best ∆V S (n c ) which is a discrete variable.Therefore, the formulation of Equations ( 10)-( 24) leads to a high-dimensional, MINLP problem for each scenario combination n c .The solution of this problem is obtained by using a MINLP solver.Parallel computing can be easily carried out because each scenario is independent of the other scenarios, i.e., multiple processors can be used each of which are responsible for a number of the MINLP OPF problems.The solutions of the MINLP OPF problems lead to a lookup table providing options of operation strategies, one of which will be selected for the grid operation in the realization phase.

Realization Phase
The lookup table provides N c solutions corresponding to the scenario combinations generated based on forecasted wind power P w.M (n w ).The actually generated wind power values of the WSs are available at each sampling interval m.For each sampling interval, one of the solutions in the lookup table will be selected and the corresponding control values realized to the network.The selection is made by comparing the actual wind power P w.A (n w , m) with the wind power scenarios of each WS P w (n w , n s ), based on the following Algorithm 1: This means, for each WS n w , if P w.A (n w , m) is not equal to P w (n w , (n s − 1)), then consider it to be the nearest higher scenario.This is because a higher wind power associates with a lower curtailment factor, leading to a feasible operation strategy.Since the scenarios generated in the prediction phase cover the whole range [0, 1], for any actual value of wind power, there is a higher or equal scenario corresponding to an optimal operation strategy.It should be noted that the decision made in this way also has a certain degree of built-in conservativeness, so that feasibility of the selected solution to be realized to the grid can be ensured.This conservativeness will be decreased if more scenarios for each WS are used.But then the number of MINLP problems to be solved will be increased correspondingly.

Implementation of the Real-Time Optimal Power Flow Framework
The RT-OPF framework of the proposed prediction-realization approach is implemented as shown in Figure 3  Steps 1-5, as shown with solid lines in Figure 3, correspond to the prediction phase, while steps 6-8, presented with dashed lines, denote the realization phase.Step 8 means that, with an optimal value of the slack bus voltage, an optimal amount of wind power is penetrated to the grid in the current sampling time m.When m = M, the computation proceeds to the next prediction horizon.correspondingly.

Implementation of the Real-Time Optimal Power Flow Framework
The RT-OPF framework of the proposed prediction-realization approach is implemented as shown in Figure 3 where the execution steps are shown by the numbers.It consists of following eight steps: (1) For the current prediction horizon, provide the forecasted active Steps 1-5, as shown with solid lines in Figure 3, correspond to the prediction phase, while steps 6-8, presented with dashed lines, denote the realization phase.Step 8 means that, with an optimal value of the slack bus voltage, an optimal amount of wind power is penetrated to the grid in the current sampling time m .When m M  , the computation proceeds to the next prediction horizon.To online realize the computation framework described above, it is necessary to ensure all 8 steps to be completed inside the time slot of p T .The implementation is illustrated in Figure 4, showing two consecutive prediction horizons.In Figure 4, the red lines indicate the execution of the eight tasks along the time.It means when the optimal strategies for the current prediction horizon are realized to the network, the solutions for the upcoming prediction horizon are prepared.
In Figure 4  To online realize the computation framework described above, it is necessary to ensure all 8 steps to be completed inside the time slot of T p .The implementation is illustrated in Figure 4, showing two consecutive prediction horizons.In Figure 4, the red lines indicate the execution of the eight tasks along the time.It means when the optimal strategies for the current prediction horizon are realized to the network, the solutions for the upcoming prediction horizon are prepared.
In Figure 4, T r denotes the length of the time reserved for data management (for our case study T r = 2 s).For steps 1, 3, 5, 6, and 8, T r means the communication (i.e., sending/receiving data) time.In steps 2, 4, and 7, T r means the time for processing data after receiving the corresponding inputs.T OPF is the time to solve the N c OPF problems, which takes the largest part of the time horizon.
T p is the prediction horizon for which the forecasted data is available and its length should be the summation of the lengths of all the tasks (including OPF computations as shown in Figure 4) in the prediction phase: Here, the greatest time slot is allocated to the computation of the OPF problems corresponding to the generated scenarios (i.e., T OPF ).It means that, at the end of the prediction phase, the OPF results must be ready for the selection algorithm in the realization phase.
Here, the greatest time slot is allocated to the computation of the OPF problems corresponding to the generated scenarios (i.e., OPF T ).It means that, at the end of the prediction phase, the OPF results must be ready for the selection algorithm in the realization phase.The wind power spectrum has been investigated for a long time.Using short term and long term records, the Van der Hoven spectrum [46] shows the diurnal and turbulence effects of wind power.It confirms that there are considerable wind power discrepancies in short term (e.g., 1 min) periods.The spectrum has been also studied to evaluate and improve wind power prediction [47,48].However, in this study, we assume that the values of wind power prediction are provided by a forecast center.The time period over which the forecasted value is updated, is defined as the prediction horizon.In our case study, the forecasted wind power value is assumed to be available every two minutes.Thus we define the length of the prediction horizon as two minutes.The forecasted wind power profiles for one day are taken from [31] for the case study.Indeed, this is the reason that the blocks 'Wind power data' and 'Demand power data' are not included in the '8 steps' of the framework in Figure 3.This means that the way how the wind and power demand power data is obtained is not considered in this study.
As mentioned in Section 3, parallel computing is used for solving the individual MINLP problems by multiple processors, each of which solves an equal number of the optimization problems.Therefore, the time needed for the solutions to be available is the maximum time taken by the processors.But we need to allocate OPF T large enough to ensure that none of the processors The sampling time T s in which the actual values of wind power are available, could be very short (even less than a second).However, since the realization is also performed at every T s , the length of T s should be realistic in order not to damage devices by too frequent control actions in very short time intervals.On the other hand, based on Equation (1), T s theoretically can be equal to T p (i.e., M = 1).However, long sampling intervals could not response the intermittency of the wind power.
The wind power spectrum has been investigated for a long time.Using short term and long term records, the Van der Hoven spectrum [46] shows the diurnal and turbulence effects of wind power.It confirms that there are considerable wind power discrepancies in short term (e.g., 1 min) periods.The spectrum has been also studied to evaluate and improve wind power prediction [47,48].However, in this study, we assume that the values of wind power prediction are provided by a forecast center.The time period over which the forecasted value is updated, is defined as the prediction horizon.In our case study, the forecasted wind power value is assumed to be available every two minutes.Thus we define the length of the prediction horizon as two minutes.The forecasted wind power profiles for one day are taken from [31] for the case study.Indeed, this is the reason that the blocks 'Wind power data' and 'Demand power data' are not included in the '8 steps' of the framework in Figure 3.This means that the way how the wind and power demand power data is obtained is not considered in this study.
As mentioned in Section 3, parallel computing is used for solving the individual MINLP problems by multiple processors, each of which solves an equal number of the optimization problems.Therefore, the time needed for the solutions to be available is the maximum time taken by the processors.But we need to allocate T OPF large enough to ensure that none of the processors exceeds this limit.The proposed RT-OPF framework is further described by a flowchart in Figure 5.The flowchart shows the prediction and realization phases.
exceeds this limit.The proposed RT-OPF framework is further described by a flowchart in Figure 5.The flowchart shows the prediction and realization phases.

Network and Input Data
The network for the case study is taken from [49,50].It is a 41-bus, 27.6 kV typical rural DN, as shown in Figure 6.The peak power demand is 16.25 MVA [5] and the substation rating is 20 MVA.Two WSs each with .

MW
w R P  and unity power factors are located at buses 2 and 16, as shown in Figure 6.The WSs are subject to different wind profiles which make the problem more complicated.The bus number 1 is considered as the slack bus with zero voltage angle [49,50].The active and reactive energy prices are adapted and fixed with 1.67 $/MW.p T and 0.4 $/Mvar.p T [31,49], respectively.We assume that the forecasted and measured values of wind power generation are available in every 120 and 20 s, respectively.Therefore, the length of the prediction horizon p T and the length of the sampling time s T are taken as 120 and 20 s, respectively.For generating wind scenarios, we take MINLP OPF problems.To implement parallel computing, we use seven processors to solve the 49 problems (i.e., seven scenario combinations are allocated to each processor).The computations are carried out on a server with 2 Intel Xeon X5690, CPU 3.47 GHz (six cores, 12 threads) and 64 GB random access memory (RAM).The problem is formulated and coded in the general algebraic modeling system (GAMS) [51] framework and the resulting problem solved by the MINLP solver BONMIN [52], using the usual flat start [5].Based on the BONMIN manual [52], the algorithms are exact when both the objective function and the constraints are convex, otherwise they are heuristics.It is also noted that there are some possible model status messages in GAMS.In solving our MINLP problems, we always got the message No. 8 (integer solution), meaning that a feasible solution has been found to problems with discrete variables, see details in [51].However, the

Network and Input Data
The network for the case study is taken from [49,50].It is a 41-bus, 27.6 kV typical rural DN, as shown in Figure 6.The peak power demand is 16.25 MVA [5] and the substation rating is 20 MVA.Two WSs each with P w.R = 10 MW and unity power factors are located at buses 2 and 16, as shown in Figure 6.The WSs are subject to different wind profiles which make the problem more complicated.The bus number 1 is considered as the slack bus with zero voltage angle [49,50].The active and reactive energy prices are adapted and fixed with 1.67 $/MW.T p and 0.4 $/Mvar.T p [31,49], respectively.We assume that the forecasted and measured values of wind power generation are available in every 120 and 20 s, respectively.Therefore, the length of the prediction horizon T p and the length of the sampling time T s are taken as 120 and 20 s, respectively.For generating wind scenarios, we take N s = 7 which means we consider seven scenarios for each of the two WSs, leading to N c = 49 MINLP OPF problems.To implement parallel computing, we use seven processors to solve the 49 problems (i.e., seven scenario combinations are allocated to each processor).The computations are carried out on a server with 2 Intel Xeon X5690, CPU 3.47 GHz (six cores, 12 threads) and 64 GB random access memory (RAM).The problem is formulated and coded in the general algebraic modeling system (GAMS) [51] framework and the resulting problem solved by the MINLP solver BONMIN [52], using the usual flat start [5].Based on the BONMIN manual [52], the algorithms are exact when both the objective function and the constraints are convex, otherwise they are heuristics.It is also noted that there are some possible model status messages in GAMS.In solving our MINLP problems, we always got the message No. 8 (integer solution), meaning that a feasible solution has been found to problems with discrete variables, see details in [51].However, the solution achieved in this way should be considered as a local solution, since the power flow equations are nonconvex.
solution achieved in this way should be considered as a local solution, since the power flow equations are nonconvex.
) to the hourly values [31].The resulting demand trajectories for one day are shown in Figure 7a.
In the reality, the forecasted and actual wind power profiles can be acquired from environmental data centers and online measurements, respectively.The forecasted wind power profiles for one day are taken from [31] and shown as the red-dashed curves in Figure 7b,c   The active P d and reactive Q d power demand are assumed to follow the hourly IEEE-RTS fall season's days [49].However, inside each hour, the demand profiles are generated for every bus at every 120 s by adding normally distributed random values (with µ d (i) = 0 and σ d (i) = 0.01) to the hourly values [31].The resulting demand trajectories for one day are shown in Figure 7a.

Test Cases
In the reality, the forecasted and actual wind power profiles can be acquired from environmental data centers and online measurements, respectively.The forecasted wind power profiles for one day are taken from [31] and shown as the red-dashed curves in Figure 7b,c.The actual wind power P w.A (n w , m) for the two WSs are generated at each 20 s using the Beta distribution with the shape parameters α(n w ) and β(n w ) corresponding to the forecasted wind power, where σ w (n w ) = 0.1 pu = 1 MW based on Equations ( 7) and (8).The resulting curves for the two WSs are shown in Figure 7b,c.
solution achieved in this way should be considered as a local solution, since the power flow equations are nonconvex.
) to the hourly values [31].The resulting demand trajectories for one day are shown in Figure 7a.
In the reality, the forecasted and actual wind power profiles can be acquired from environmental data centers and online measurements, respectively.The forecasted wind power profiles for one day are taken from [31] and shown as the red-dashed curves in Figure 7b,c

Test Cases
In the case study, possible reverse power flow is considered for the proposed RT-OPF framework.We define the forward and reverse active as well as reactive energy at the slack bus as follows: • Forward energy flow: The forward active and reactive energy from the HV network to the MV network is to be minimized based on an energy price model.

•
Reverse energy flow: The reverse power flow could have impacts on voltage profiles [53] of the upper-level network and may result in specific operational limits being exceeded at the congested primary substations [54].However, reverse flows have been considered in many studies [45,[55][56][57][58] and in reality, they are likely to happen.Therefore, in this paper we consider the cases with and without reverse power flows.
Furthermore, in the optimization problem formulation, the slack bus voltage can be either a parameter with a fixed value or a discrete free variable.Based on these considerations, we carry out RT-OPF for four different cases defined as follows: Case (1): Both reverse active and reactive power to the upstream HV network is not allowed (i.e., γ Ps,rev = γ Qs,rev = 0), and with a fixed value of the slack bus voltage (V S (m) = V S (n c ) = 1 pu).Case (2): Both reverse active and reactive power to upstream HV network is not allowed (i.e., γ Ps,rev = γ Qs,rev = 0), and with the slack bus voltage as a discrete free variable.Case (3): Both reverse active and reactive power to upstream HV network is allowed (i.e., γ Ps,rev = γ Qs,rev = 1), and with a fixed value of the slack bus voltage (i.e., V S (m) = V S (n c ) = 1 pu).Case (4): Both reverse active and reactive power to upstream HV network is allowed (i.e., γ Ps,rev = γ Qs,rev = 1), and with the slack bus voltage as a discrete free variable.

Results and Discussions
We run the RT-OPF for one day for the aforementioned four cases using the same network in Figure 6 and the input data in Figure 7.To illustrate the prediction phase and the realization phase, we firstly take the first prediction horizon (i.e., t ∈ [0, 120 s]) as an example.The forecasted values of wind power for the two WSs in the first prediction horizon are P w.M (1) = 3.8 MW and P w.M (2) = 7.05 MW.The wind power scenario combinations and the corresponding results for the four cases in the first prediction horizon are shown in Tables 2-5.
The lookup tables consist of two sections: "scenario combinations" and "optimal results".The scenario combinations are obtained as described in Section 3 and optimal results are achieved by solving the corresponding OPF problems (i.e., Equations ( 10)-( 24)).Based on Equation ( 9), there are 49 scenarios in each lookup table from which one action will be selected in each sampling interval T s .It is noted that the lookup tables are updated in each prediction horizon T p based on the new values of forecasted wind and demand power.
In the first sampling time interval (i.e., t ∈ [0, 20 s]), the real wind power are P w.A (1, 1) = 3.74 MW and P w.A (2, 1) = 4.17 MW, respectively.Therefore, the selection algorithm (shown as Step 7 in Figure 3) selects the scenario combination 27 which corresponds to the level higher than the real value of wind power for both WSs.Then the control strategy corresponding to this scenario combination will be realized to the grid.
It can be seen in Tables 2 and 3 that, for this prediction horizon, the active power at the slack bus is zero when the wind power of WSs is high.This is because of the low total active power demand (6.65 MW) and high active wind power generation.In addition, the high wind power generation leads to high curtailment (i.e., low values of the curtailment factors) to ensure feasibility.In contrast, if reverse power flow is allowed, as in Cases (3) and (4), the surplus active wind power will not be curtailed, i.e., β w (n c ) = 1, and it is exported to the upstream HV network, showing negative active power at the slack bus.The reactive power at the slack bus is always positive for all scenarios of the four cases in this time horizon, i.e., it will be imported from the upstream network.This is because of using unity power factors of the WSs and the reactive power compensation of feeder capacitive susceptance [49] (the total reactive power demand is 2.46 Mvar in this time horizon).From another perspective, for this prediction horizon, in Cases ( 2) and ( 4) where the optimization problems are formulated as MINLP, the values of slack bus voltage tend to be more than 1 pu.This is due the fact that a higher slack bus voltage results in less power losses and consequently higher values of the objective function.From the economic point of view, the optimal strategy has some degree of conservativeness due to the (limited) number of scenarios.This can be seen from the resulting active power imported from the upstream network, shown in subplot (d) in Figures 8 and 9, which means in most time the really imported value is higher than the expectedly imported value.However, when the total wind power is much higher than the forecasted value but lower than the total demand plus the losses, the really imported active power will be lower than the expected value, leading to higher total revenue for such sampling intervals as shown in subplots (f) in Figures 8  and 9.In Cases (3) and ( 4) where γ Ps,rev = γ Qs,rev = 1, i.e., the reverse power flow to upstream HV network is allowed, there will be no wind power curtailment, i.e., β w (n c ) = β w (m) = 1.Thus, high amount of energy is exported through the slack bus leading to negative values of active power as shown in subplots (d) in Figures 10 and 11.As clearly shown in subplots (e), there is no significant difference between the forecasted and actual values of reactive power at the slack bus.This is due to unity power factors of WSs which means the reactive power imported from the upstream HV network is not sensitive to wind power fluctuations and it mostly follows the same trend as the total reactive power demand in the network.However, comparing subplots (e) in Figures 8 and 9 with the ones in Figures 10 and 11, it is clearly seen that the reactive power at the slack bus is increased when the reverse power is allowed.The reason is that increasing the power flow in the network increases both the active and reactive power losses in the grid.Therefore, most of the reactive power losses must be supplied by the upstream network.
As shown in subplots (c) of Figures 9 and 11, when the slack bus voltage is taken into account as a discrete free variable, it tends to take the values higher than 1 pu which leads to decrease in power losses.At the same time, the slack bus voltage must lead to satisfying the constraints of voltage at PQ buses (i.e., V min (i) = 0.94 pu and V max (i) = 1.06 pu).When the reverse power flow is not allowed, the slack bus voltage varies between 1.06 and 1.07 as shown in Figure 9c.Allowing reverse power flow results in increasing the range to be between 1.03 and 1.07 as shown in Figure 11c.This is because power flows from buses with a higher voltage to a lower one, which means that in the case of reverse power flow, the slack bus voltage must be lower than the voltage at bus 2.
As mentioned earlier, the optimization results from the proposed approach ensure feasible operations.To illustrate this, we compare our results with those by the operation strategies based on the forecasted wind power.In the latter approach, the forecasted wind power is used for the optimization and the results are directly applied to the grid.Since no correction of the solutions is made based on the realized wind power, we call it the forecasted approach.The results of the feasibility obtained by the new approach and by the forecasted approach are compared in subplots (g) of  In Case (1), as shown in Figure 8g, using the results of the forecasted approach leads to 2066 infeasible sampling intervals.This takes place because the actual wind power is higher than the forecasted value at these sampling intervals, and thus, the forecasted curtailment factors cannot satisfy the active power balance (Equation ( 15)).In Figure 9g, there are 2110 infeasible sampling intervals for Case (2), most of which are violating the active power balance (Equation ( 15)) and the others are violating the other inequality constraints, in particular the voltage bounds at bus 2. This happens when the realized wind power is higher than the forecasted value, and meanwhile, the forecasted active power at the slack bus is greater than zero.Thus, in the reality, the active power from the slack bus is decreased and the voltage drop in the cables between the slack bus and bus 2 is reduced.Since the voltage at the slack bus is fixed, the voltage at bus 2 has to be increased which leads to violating its upper bound.This situation also occurs in Case (4) in the sampling intervals during which the forecasted active power at the slack bus is greater than zero.
Moreover, in Case (4), during the time in which the forecasted active power at the slack bus is negative, if the total actual wind power is higher than the forecasted value, the exported power to the HV network will be increased (i.e., P S (m) is increased in the reverse direction).Again, the voltage drop in the feeders between the slack bus and bus 2 is increased which leads to violating the upper bound of the voltage at this bus.In Case (4), totally 625 infeasible sampling intervals are observed, among which none of them is due to the power balance since the surplus power can be exported to the HV network.In Case (3), as shown in Figure 10g, the control strategies by the forecasted approach are also feasible, since the slack bus voltage is fixed in the middle position (i.e., 1 pu) and the surplus wind power can be exported to the HV network.
Energies 2017, 10, 535 19 of 28 exported to the HV network.In Case (3), as shown in Figure 10g, the control strategies by the forecasted approach are also feasible, since the slack bus voltage is fixed in the middle position (i.e., 1 pu) and the surplus wind power can be exported to the HV network.The feasibility of the control strategies by the proposed prediction-realization approach for all the four cases are shown with the blue line in subplots (g) of Figures 8-11.It can be seen that there is The feasibility of the control strategies by the proposed prediction-realization approach for all the four cases are shown with the blue line in subplots (g) of  It can be seen that there is no infeasible operations at all, i.e., the proposed approach can ensure the operation feasibility in any case of the wind power generation.
Subplots (h) of  give the computation time taken by the seven processors to solve the MINLP OPF problems for 49 scenario combinations in real-time.It can be seen that the computation time of each processor is less than the reserved time (T OPF = 112 s) which ensures the applicability of the RT-OPF.
Table 6 compares the results in the scope of profitability of one-day operation in the four cases.It is shown that, in Cases ( 2) and ( 4), the total yield (i.e., the objective function value) will be increased in comparison to Cases (1) and ( 3), if the slack bus voltage is considered as a variable and optimized, since, in this way, the active power losses are decreased.In contrary, when the reverse power flow is allowed (i.e., Cases (3) and ( 4)), the active power losses increase dramatically owing to higher power flow in the network.However, the effect of allowable reverse power flow is much more significant as the total yield is increased by more than three times, due to the export of the surplus wind power to the upstream HV network.

Conclusions
RT-OPF is indispensable for network operations under intermittent wind power, but its numerical implementation poses a significant challenge.In this paper, a prediction-realization approach RT-OPF is introduced for energy systems to deal with fast changing wind power.In addition, our OPF simultaneously considers curtailment factors as continuous free variables, the reference voltage at the slack bus as discrete free variables, and possible reverse power flow to an upstream network.This leads to a large-scale MINLP OPF problem with uncertain wind power generation.To address this problem, we employ the available information of wind power, i.e., forecasted value in a long time cycle (as the prediction horizon) and measured value in a short time cycle (as the sampling time interval), and developed a two-phase solution approach.In the prediction phase, most probable wind power scenarios are generated based on the Beta distribution for the prediction horizons.The corresponding MINLP OPF problems are then solved in parallel.The results are saved as a lookup table which provides a base for selecting a decision when the actual wind power value is available from one sampling interval to the next.As a result, the proposed RT-OPF framework ensures feasible solutions for the cases with and without reverse power flow to an upstream network.The solutions can be realized in a very short sampling time.The results from a case study demonstrate the applicability of the proposed RT-OPF framework.
One important point we discussed in the paper is the insurance of the operation feasibility.A series of wind power scenarios are generated according to the stochastic distribution of wind power and the corresponding optimization problems are solved in a predictive way.In the realization phase, when the measured wind power is different from the predefined scenarios, we choose the solution of the optimization problem with the scenario exactly higher than the measured value.The selection leads to a higher curtailment of the generated wind power, i.e., it guarantees the feasibility with a certain degree of conservativeness (with a lower yield).In fact, a reservation of some conservativeness to ensure feasibility is also commonly used in industrial practice.

M
Total number of sampling intervals in each prediction horizon.

N bus
Total number of buses.

N s
Total number of wind power scenarios for each WS.

N pr
Total number of processors.

N c
Total number of wind power scenario combinations.

N w
Total number of WSs.P d (i) Active power demand at bus i.Price P Price for active energy.Price O Price for reactive energy.P w.R (n w ) Rated installed wind power of WS .Q d (i) Reactive power demand at bus i. S l.max (i, j) Upper limit of apparent power flow of line between bus i and j. S S.max Upper limit of apparent power at slack bus.

Figure 2 .
Figure 2. Graphical examples of wind power scenario combinations for (a) N w = 2, N s = 7 and (b) N w = 3, N s = 4. Here, dark to light colors denote high to low wind power scenarios, respectively.

28 ( 1 )
where the execution steps are shown by the numbers.It consists of following eight steps: Energies 2017, 10, 535 9 of For the current prediction horizon, provide the forecasted active P d (i) and reactive Q d (i) demand power and wind power P w.M (n w ).(2) Generate N c wind power scenario combinations based on the Beta distribution as described in Section 3. (3) Send the generated scenarios as inputs to formulate N c MINLP OPF problems.(4) Solve the N c MINLP OPF problems with parallel computing.(5) Send the solution results as a lookup table to the selection algorithm.(6) Provide the actual wind power of WSs, P w.A (n w , m), available at the current sampling time m (for m = 1, • • • , M), to the selection algorithm.(7) Select one of the solutions from the lookup table based on P w.A (n w , m) and the selection algorithm (see Section 4.2).(8) Send the values of the controls V S (m) and β w (m) to the grid.

N
wind power scenario combinations based on the Beta distribution as described in Section 3. (3) Send the generated scenarios as inputs to formulate c N MINLP OPF problems.(4) Solve the c N MINLP OPF problems with parallel computing.(5) Send the solution results as a lookup table to the selection algorithm.(6) Provide the actual wind power of WSs, .( , ) w A w P n m , available at the current sampling time m (for 1, , m M   ), to the selection algorithm.(7) Select one of the solutions from the lookup table based on .( , ) w A w P n m and the selection algorithm (see Section 4.2).(8) Send the values of the controls

Figure 3 .
Figure 3. Framework of the proposed real-time optimal power flow (RT-OPF) for a prediction horizon.Here, HV, MV and LV denote high-voltage, medium-voltage and low-voltage, respectively.

, rT
denotes the length of the time reserved for data management (for our case study steps 1, 3, 5, 6, and 8, r T means the communication (i.e., sending/receiving data) time.In steps 2, 4, and 7, r T means the time for processing data after receiving the corresponding inputs.OPF T is the time to solve the c N OPF problems, which takes the largest part of the time horizon.

Figure 3 .
Figure 3. Framework of the proposed real-time optimal power flow (RT-OPF) for a prediction horizon.Here, HV, MV and LV denote high-voltage, medium-voltage and low-voltage, respectively.

Figure 4 .
Figure 4. Time allocation for the computational tasks of the 8 steps in Figure 3.

Figure 4 .
Figure 4. Time allocation for the computational tasks of the 8 steps in Figure 3.

Figure 5 .
Figure 5.The flowchart of the proposed RT-OPF framework.

7 sN
 which means we consider seven scenarios for each of the two WSs,

Figure 5 .
Figure 5.The flowchart of the proposed RT-OPF framework.

Figure 6 .
Figure 6.Distribution network (DN) for the case study.
. The actual wind power .( , ) w A w P n m for the two WSs are generated at each 20 s using the Beta distribution with the shape parameters ( ) (7) and (8).The resulting curves for the two WSs are shown in Figure 7b,c.

Figure 7 .
Figure 7. Trajectories of one day: (a) total active (blue-solid) and reactive (red-dashed) power demand; (b) forecasted (red-dashed) and actual (blue-solid) wind power of the first WS; and (c) forecasted (red-dashed) and actual (blue-solid) wind power of the second WS.

Figure 6 .
Figure 6.Distribution network (DN) for the case study.

Figure 6 .
Figure 6.Distribution network (DN) for the case study.
. The actual wind power .( , ) w A w P n m for the two WSs are generated at each 20 s using the Beta distribution with the shape parameters ( ) (7) and (8).The resulting curves for the two WSs are shown in Figure 7b,c.

Figure 7 .
Figure 7. Trajectories of one day: (a) total active (blue-solid) and reactive (red-dashed) power demand; (b) forecasted (red-dashed) and actual (blue-solid) wind power of the first WS; and (c) forecasted (red-dashed) and actual (blue-solid) wind power of the second WS.

Figure 7 .
Figure 7. Trajectories of one day: (a) total active (blue-solid) and reactive (red-dashed) power demand; (b) forecasted (red-dashed) and actual (blue-solid) wind power of the first WS; and (c) forecasted (red-dashed) and actual (blue-solid) wind power of the second WS.

Figure 8 .
Figure 8. Trajectories of one day for Case 1.(a) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for first WS; (b) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for second WS; (c) Forecasted (red-dashed) and actual (blue-solid) values of voltage at the slack bus; (d) Forecasted (red-dashed) and actual (blue-solid) slack bus active power; (e) Forecasted (red-dashed) and actual (blue-solid) slack bus reactive power; (f) Forecasted (red-dashed) and actual (blue-solid) total objective function value; (g) Feasibility status of the deterministic (red-dashed) and prediction-realization (blue-solid) approaches when applying actual wind power.Here, 1 denotes feasible and 0 denotes infeasible solution; and (h) computational time of the seven processors.

Figure 8 .
Figure 8. Trajectories of one day for Case 1.(a) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for first WS; (b) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for second WS; (c) Forecasted (red-dashed) and actual (blue-solid) values of voltage at the slack bus; (d) Forecasted (red-dashed) and actual (blue-solid) slack bus active power; (e) Forecasted (red-dashed) and actual (blue-solid) slack bus reactive power; (f) Forecasted (red-dashed) and actual (blue-solid) total objective function value; (g) Feasibility status of the deterministic (red-dashed) and prediction-realization (blue-solid) approaches when applying actual wind power.Here, 1 denotes feasible and 0 denotes infeasible solution; and (h) computational time of the seven processors.

Figure 9 .
Figure 9. Trajectories of one day for Case 2. (a) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for first WS; (b) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for second WS; (c) Forecasted (red-dashed) and actual (blue-solid) values of voltage at the slack bus; (d) Forecasted (red-dashed) and actual (blue-solid) slack bus active power; (e) Forecasted (red-dashed) and actual (blue-solid) slack bus reactive power; (f) Forecasted (red-dashed) and actual (blue-solid) total objective function value; (g) Feasibility status of the deterministic (red-dashed) and prediction-realization (blue-solid) approaches when applying actual wind power.Here, 1 denotes feasible and 0 denotes infeasible solution; and (h) computational time of the seven processors.

Figure 9 .
Figure 9. Trajectories of one day for Case 2. (a) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for first WS; (b) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for second WS; (c) Forecasted (red-dashed) and actual (blue-solid) values of voltage at the slack bus; (d) Forecasted (red-dashed) and actual (blue-solid) slack bus active power; (e) Forecasted (red-dashed) and actual (blue-solid) slack bus reactive power; (f) Forecasted (red-dashed) and actual (blue-solid) total objective function value; (g) Feasibility status of the deterministic (red-dashed) and prediction-realization (blue-solid) approaches when applying actual wind power.Here, 1 denotes feasible and 0 denotes infeasible solution; and (h) computational time of the seven processors.

Figure 10 .
Figure 10.Trajectories of one day for Case 3. (a) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for first WS; (b) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for second WS; (c) Forecasted (red-dashed) and actual (blue-solid) values of voltage at the slack bus; (d) Forecasted (red-dashed) and actual (blue-solid) slack bus active power; (e) Forecasted (red-dashed) and actual (blue-solid) slack bus reactive power; (f) Forecasted (red-dashed) and actual (blue-solid) total objective function value; (g) Feasibility status of the deterministic (red-dashed) and prediction-realization (blue-solid) approaches when applying actual wind power.Here, 1 denotes feasible and 0 denotes infeasible solution; and (h) computational time of the seven processors.

Figure 10 .
Figure 10.Trajectories of one day for Case 3. (a) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for first WS; (b) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for second WS; (c) Forecasted (red-dashed) and actual (blue-solid) values of voltage at the slack bus; (d) Forecasted (red-dashed) and actual (blue-solid) slack bus active power; (e) Forecasted (red-dashed) and actual (blue-solid) slack bus reactive power; (f) Forecasted (red-dashed) and actual (blue-solid) total objective function value; (g) Feasibility status of the deterministic (red-dashed) and prediction-realization (blue-solid) approaches when applying actual wind power.Here, 1 denotes feasible and 0 denotes infeasible solution; and (h) computational time of the seven processors.

Figure 11 .
Figure 11.Trajectories of one day for Case 4. (a) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for first WS; (b) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for second WS; (c) Forecasted (red-dashed) and actual (blue-solid) values of voltage at the slack bus; (d) Forecasted (red-dashed) and actual (blue-solid) slack bus active power; (e) Forecasted (red-dashed) and actual (blue-solid) slack bus reactive power; (f) Forecasted (red-dashed) and actual (blue-solid) total objective function value; (g) Feasibility status of the deterministic (red-dashed) and prediction-realization (blue-solid) approaches when applying actual wind power.Here, 1 denotes feasible and 0 denotes infeasible solution; and (h) computational time of the seven processors.

Figure 11 .
Figure 11.Trajectories of one day for Case 4. (a) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for first WS; (b) Forecasted (red-dashed) and actual (blue-solid) curtailment factors for second WS; (c) Forecasted (red-dashed) and actual (blue-solid) values of voltage at the slack bus; (d) Forecasted (red-dashed) and actual (blue-solid) slack bus active power; (e) Forecasted (red-dashed) and actual (blue-solid) slack bus reactive power; (f) Forecasted (red-dashed) and actual (blue-solid) total objective function value; (g) Feasibility status of the deterministic (red-dashed) and prediction-realization (blue-solid) approaches when applying actual wind power.Here, 1 denotes feasible and 0 denotes infeasible solution; and (h) computational time of the seven processors.

Functionsf
Objective function.F Total value of objective function for one day.F 1 Total revenue from wind power injection for one day.F 2 Total cost of active energy losses in the grid for one day.F 3 Total cost of active energy at slack bus for one day.F 4 Total cost of reactive energy at slack bus for one day.f (n c ) Total value of objective function for scenario combination n c .f 1 (n c ) Total revenue from wind power injection for scenario combination n c .f 2 (n c ) Total cost of active energy losses in the grid for scenario combination n c .f 3 (n c ) Total cost of active energy at slack bus for scenario combination n c .f 4 (n c ) Total cost of reactive energy at slack bus for scenario combination n c .F PDF Probability distribution function.f P (n c ) Network active power function for scenario combination n c .f Q (n c ) Network reactive power function for scenario combination n c .g Equality equations.ρ Density function.Parameters L Upper bound on integer variables.

Table 1 .
The list of wind power scenario combinations for all WSs.

Table 1 .
The list of wind power scenario combinations for all WSs.

Table 1 .
The list of wind power scenario combinations for all WSs.

Algorithm 1
Comparing and selection of wind power for each WS n w = 1, • • • , N w and n s ≥ 2 If P w (n w , (n s − 1)) < P w.A (n w , m) ≤ P w (n w , n s ) then consider P w.A (n w , m) as P w (n w , n s ) end Achieve P w

Table 6 .
Comparison of results of one day for four cases.