Probabilistic Security-Constrained Preventive Control under Forecast Uncertainties Including Volt/Var Constraints

Abstract

The continuous increase in generation from renewable energy sources, marked by correlated forecast uncertainties, requires specific methodologies to support power system operators in security management. This paper proposes a probabilistic preventive control to ensure N-1 security in presence of correlated uncertainties of renewable sources and loads. By adopting a decoupled linear formulation of the AC load flow equations, the preventive control is decomposed into two subsequent linear programming problems, the former concerning the active power and the latter the voltage/reactive power-related issues. In particular, in the active control problem, the algorithm combines Third Order Polynomial Normal Transformation, Point Estimate Method, and Cornish–Fisher expansion to model the forecast uncertainties and characterize the chance constraints in the problem. The goal is to find the optimal phase shifting transformer tap setting, conventional generation redispatching, and renewable curtailment at the minimum cost to assure the probabilistic fulfillment of N and N-1 security constraints on branch active power flows. The second stage solves another linear programming problem, which aims to minimize the adjustments to generators’ set-point voltages to avoid violations at node voltages and branch-rated limits due to reactive power flows while meeting generator reactive power constraints. Simulations performed on an IEEE test system demonstrate the effectiveness of the proposed security control method in limiting the probability of violating security limits in N and N-1 state, including voltage/reactive power constraints, in presence of correlated uncertainties.

1. Introduction

The spread of Renewable Energy Sources (RES) is profoundly changing the planning and operation of power systems. In general terms, the large-scale integration of renewables has brought a higher level of uncertainty in the electric system due to their stochastic and intermittent generation.

For this reason, the uncertainty management represents a fundamental step for the efficient and economical operation of electric systems. Therefore, specific strategies for uncertainty management are crucial to favor RES penetration and guarantee the grid’s security.

Different sources of uncertainty have to be considered according to the phase of power system management, from long term planning to operational planning and quasi-real-time control [1,2]. In an operational planning context, important sources of uncertainty regard the load demand and renewable generation forecast, as well as the outages of components, also accounting for incumbent adverse weather conditions. On the other hand, neglecting forecast uncertainties in power system management may lead to security issues [3,4].

Therefore, advanced methods are needed to address uncertainties to minimize the probability of violating system operational limits.

In this context, probabilistic methodologies based on Chance Constrained (CC) formulation [2,5,6,7,8,9], wherein the constraints are formulated in terms of violation probability of operational limits, can be very useful. Authors in [5] analyze different mathematical formulations of the chance constraints that, however, cannot consider specific properties (e.g., asymmetry) of the probabilistic distributions and their correlations. Reference [6] proposes the CC and Security Constrained Optimal Power Flow (SC-OPF), exploiting corrective actions from control devices like Phase Shifting Transformers (PSTs) and HVDC links. In [7], a CC OPF formulated as a convex conic optimization problem is presented together with an efficient cutting plane algorithm for its solution. Authors in [8] describe a real-time CC dispatch based on quadratic programming and dynamic power flow. This approach also considers constraints on conventional generator ramps on frequency deviation. In addition, the chance constraints are modeled through Cornish–Fisher expansion and cumulants’ methodology. The solution proposed in [9] incorporates generator reserve within the CC formulation. Authors in [10] propose an SC OPF based on CC that explicitly considers the probability of contingency events and potential failures in the operation of post-contingency corrective controls. It can be remarked this work is based on a full AC model solved through a heuristic approach but it does not consider uncertainties on load and RES injections. In [11], an algorithm for solving joint CC DC OPF problem is proposed, using a sample-based approach that avoids making strong assumptions on the distribution of the uncertainties and scales favorably to large problems. In [12], quadratic programming for the AC OPF is proposed for applications in smart distribution networks. This approach exploits linear constraints to represent power losses, resulting in a trade-off between accuracy, performance, and robustness. Finally, [13] presents a risk-based CC optimal power flow model by leveraging the methods of both adjustable uncertainty set and robust optimization.

The present paper proposes an innovative CC and SC control that assures the probabilistic fulfillment of N-1 security constraints in presence of forecast uncertainties due to renewable sources and loads. To this aim, it performs a set of preventive actions, i.e., the redispatch of active powers injected by conventional generators, the regulation of Phase Shift Transformers, and the renewables’ curtailment, to solve active power-related issues. However, it can also adjust the voltage references of the generating units to solve potential issues related to voltages and reactive powers.

The proposed methodology uses Point Estimate Method (PEM) and Third order Polynomial Normal Transformation (TPNT) [14,15] to obtain the raw moments (and, therefore, the cumulants) of the quantities subject to constraints (e.g., active power flows), taking into account correlated forecast uncertainties.

The control is formulated as a two-stage problem: the former stage is formulated as a linear programming problem, including a linearized model for PSTs based on distribution factors, with the final aim to minimize the redispatching cost of conventional generators and renewable curtailment. The chance constraints concern the technical limits of generators, the invariance of total generation (conventional and from RES), and the active power flow limits both in state N and N-1. The chance constraints are formulated as inequalities that involve quantiles of uncertain variables evaluated from cumulants through the Cornish–Fisher expansion. The first stage represents an upgrade of the work presented in [16], integrating a linearized model of the PSTs for active power flow management. The latter stage is additional with respect to [16], and it is formulated as a linear programming problem on a linearization of the original AC power flow equations. In this stage, the control variables are the variations of voltage reference set-points at PV generation buses, with the aim to avoid violations at PQ node voltages and at branch rated limits, limiting the relevant reactive power flows while meeting the generator reactive capability limits.

It is essential to highlight that the proposed two-stage approach has been designed consisting of two linear optimization algorithms for its future application to large electrical networks, thus, limiting the computational burden and guaranteeing globally optimal solutions. Finally, the presented methodology has been implemented in Matlab/GAMS environment, exploiting GAMS [17] as specific software for optimization problems improving the implementation and the computation efficiency with respect to the architecture used in [16] completely relying on Matlab.

In summary, this paper provides the following contributions:

• The proposal of a two-stage formulation of the control as an efficient linear programming problem, which includes both active power-related and voltage/reactive power-related constraints. The latter ones are not considered by all the works reported in the introduction, except for [10], which proposes a heuristic solution of the full AC power flow.
• The integration of a linearized model of the AC power flow equations favors the decoupling of active and reactive power control problems in the proposed security-constrained preventive control strategy.
• The design and the implementation of a linear optimization algorithm for the Volt/Var control (i.e., the second stage of the proposed CC and SC control) which is based on the above-mentioned linearized model and which receives the inputs from the active power control stage.
• The integration of a linearized model of the PSTs, based on Phase Shifter Distribution Factors (PSDFs), for the active power flow management within the first stage of the presented methodology. The adoption of PSDF-based modeling for PSTs, together with the linearization of AC power flow equations, allows one to maintain a linear programming formulation that can be solved in a very efficient way.

The paper is organized as follows: Section 2 presents the preliminary requirements in terms of power system modeling for the definition of the proposed control framework and Section 3 describes the detailed mathematical formulation. Section 4 presents the case study and the simulation scenarios and describes the results obtained from the application of the algorithm. Section 5 concludes.

2. Power System Modeling: Preliminary Requirements for an Effective Integration into the Control Framework

Before formulating the proposed control, it is important to describe the techniques available from the literature briefly and used in the present approach to effectively integrate some aspects of power system operation in the control framework. This integration has been implemented to preserve a linear formulation which is very time efficient. Specifically, the control will benefit from:

• The linearization of AC power flow equations, based on the technique in [18].
• The representation of control devices, specifically the Phase Shifting Transformers (PST), based on the application of PSDFs [19,20].
• The representation of the asymmetric distributions and the correlation of forecast uncertainties for RES generation and load absorption in the chance constraints of the control, combining the following techniques, namely the Third Order Polynomial Normal Transformation, the Point Estimate method [14,15], and the Cornish–Fisher expansion [16].

In the sequel, these techniques are recalled with the relevant references.

2.1. AC Power Flow Linearization

This subsection recalls the well-known decoupled linear formulation of the AC power flow in [18], defining the quantities and symbols used in the control framework.

To facilitate the paper comprehension, it is helpful to recap the standard equation related to the AC power flow:

$$\begin{array}{cc}\hfill P_i& =V_i\sum _{k=1}^{N_B}{V}_k\left[G_{ik}cos\left({\theta }_{ik}\right)+B_{ik}sin\left({\theta }_{ik}\right)\right]i=1,\dots ,N_B\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill Q_i& =V_i\sum _{k=1}^{N_B}{V}_k\left[G_{ik}sin\left({\theta }_{ik}\right)-B_{ik}cos\left({\theta }_{ik}\right)\right]i=1,\dots ,N_B\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill P_{ik}& =-{\tau }_{ik}{G}_{ik}{V}_i^2+V_i{V}_k\left[G_{ik}cos\left({\theta }_{ik}\right)+B_{ik}sin\left({\theta }_{ik}\right)\right]\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill Q_{ik}& ={\tau }_{ik}{B}_{ik}{V}_i^2-S{U}_{ik}{V_i^2+V}_i{V}_k\left[G_{ik}sin\left({\theta }_{ik}\right)-B_{ik}cos\left({\theta }_{ik}\right)\right]\hfill \end{array}$$

(4)

where $P_i$ is the active power at node i, $N_B$ is the number of nodes, $V_i$ is the voltage magnitude at node i, $G_{ik}$ is the real part of the element $(i,k)$ of the admittance matrix, $B_{ik}$ is the imaginary part of the element $(i,k)$ of the admittance matrix, ${\theta }_{ik}$ is the angle difference between node i and k, $Q_i$ is the reactive power injection at node i, $P_{ik}$ is the active power flow related to line $(i,k)$, $Q_{ik}$ is the reactive power flow through line $(i,k)$, ${\tau }_{ik}$ is the transformer tap ratio of line $(i,k)$, while $S{U}_{ik}$ is the half of susceptance of line $(i,k)$.

2.1.1. Active Power

In this work, the active power flows and the angles are evaluated with the following assumptions:

$$V_i=V_k=1\mathrm{p}.\mathrm{u}.$$

(5)

$$G_{ik}=0(\mathrm{line} \mathrm{resistance} \mathrm{assumed} \mathrm{zero})$$

(6)

Exploiting Equations (5) and (6), Equation (1) becomes:

$$P_i=\sum _{k=1}^{N_B}\frac{1}{X_{ik}}·{\theta }_{ik}$$

(7)

where $X_{ik}$ is the reactance of line $(i,k)$. Equation (7) can be represented also in matrix form as

$$P=C\mathsf{\Theta }$$

(8)

wherein $\mathsf{\Theta }$ is the vector of the angles at the buses with entries ${\theta }_i$, while

$$C_{ik}=-\frac{1}{X_{ik}}and{C}_{ii}=\sum _{\begin{array}{c}k=1\\ k\ne i\end{array}}^{N_B}\frac{1}{X_{ik}}$$

(9)

in which the slack busbar row and column are deleted. Inverting gives

$$\mathsf{\Theta }=C^{-1}Por{\theta }_i=\sum _{k=1}^{N_B-1}{C}_{ik}^{-1}{P}_{k}i=1,\dots ,N_B-1$$

(10)

Similarly, using Equations (5) and (6), Equation (3) can be written as

$$P_{ik}=\frac{{\theta }_i-{\theta }_k}{X_{ik}}$$

(11)

from which substituting Equation (10) gives

$$P_{ik}=\frac{1}{X_{ik}}\sum _{j=1}^{N_B-1}\left(C_{ij}^{-1}-C_{kj}^{-1}\right)P_j$$

(12)

in which if bus i is slack, $C_{ij}=0$, and therefore:

$$P_{slack}=\sum _{i=1}^{N_B}{P}_i$$

(13)

Notice that for the active power, this decoupled approach is equivalent to the DC power flow [18].

2.1.2. Reactive Power

As far as reactive power is concerned, in this formulation, the linearization is performed starting from the assumption that the voltage magnitudes are close to 1 p.u.

Under this assumption, it is possible to rewrite Equation (2) as

$$Q_i=\sum _{i=1}^{N_B}{\left(V_i+V_k\right)·A}_{ik}i=1,\dots ,N_B$$

(14)

where

$$A_{ik}=G_{ik}sin\left({\mathsf{\Theta }}_{ik}\right)-B_{ik}cos\left({\mathsf{\Theta }}_{ik}\right)$$

(15)

$$A_{ii}=-B_{ii}$$

(16)

Notice that ${\theta }_{ik}={\theta }_i-{\theta }_k$, with ${\theta }_i$ and ${\theta }_k$ calculated through Equation (10).

Equation (14) can be rearranged obtaining:

$$Q_i=\sum _{i=1}^{N_B}{A}_{ik}^{\prime }{V}_k-W_{i}i=1,\dots ,N_B$$

(17)

where

$$A_{ik}^{\prime }=A_{ik}k\ne i,A_{ii}^{\prime }=A_{ii}+W_i$$

(18)

$$W_i=\sum _{k=1}^{N_B}{A}_{ik}i=1,\dots ,N_B$$

(19)

At this point, it is possible to write Equation (17) in matrix form as

$$Q=A^{\prime }V-W$$

(20)

and also its partition in blocks:

$$\left[\begin{array}{c}{Q}_l\\ Q_g\end{array}\right]=\left[\begin{array}{cc}{M}^{\prime }& L^{\prime }\\ N^{\prime }& J^{\prime }\end{array}\right]\left[\begin{array}{c}{V}_l\\ V_g\end{array}\right]-\left[\begin{array}{c}{W}_l\\ W_g\end{array}\right]$$

(21)

where $Q_l$ is the vector of PQ nodes while $Q_g$ is the vector of PV nodes and slack bus.

From Equation (21), the following relations can be obtained:

$$V_l=({M^{\prime })}^{-1}{Q}_l-{\left(M^{\prime }\right)}^{-1}{H}^{\prime }l=1,\dots ,n_1$$

(22)

$$Q_g=D^{\prime }{Q}_l+E^{\prime }g=1,\dots ,n_2$$

(23)

where $n_1$ is the number of PQ nodes, $n_2$ is the number of PV nodes plus slack bus, and

$$H^{\prime }=W_l-L^{\prime }{V}_g$$

(24)

$$D^{\prime }=N^{\prime }{\left(M^{\prime }\right)}^{-1}$$

(25)

$$E^{\prime }=D^{\prime }{H}^{\prime }+J^{\prime }{V}_g-W_g$$

(26)

Finally, exploiting the assumption that the voltage magnitudes are close to 1 p.u within Equation (4), it is possible to write the reactive power flows as

$$Q_{ik}={\beta }_{ik}{V}_i+A_{ik}{V}_k+{\gamma }_{ik}$$

(27)

where

$${\beta }_{ik}=2{\mu }_{ik}+A_{ik}$$

(28)

$${\gamma }_{ik}=-{\mu }_{ik}-A_{ik}$$

(29)

$${\mu }_{ik}={\tau }_{ik}{B}_{ik}-S{U}_{ik}.$$

(30)

Summing up, the active power flows and the angles are obtained through a DC power flow formulation, while the voltage magnitudes at PQ nodes, the reactive power injections at PV nodes and the reactive power flows are evaluated, respectively, through Equations (22), (23), and (27).

2.2. PST Modeling

The technique exploited in this work is taken from [19,20] and is based on the phase shifter distribution factors that are similar to the Power Transfer Distribution Factors (PTDF) [21], which will be used in the proposed formulation to evaluate the contribution on the active power control.

The PST has proven to be a robust technology that can significantly contribute to regulating power flows [19,20,22].

The implemented model of the PST is approximate, but according to [19,20], the accuracy is very satisfactory as long as the PST angle is not extremely large (>30°).

For the sake of clarity, this subsection briefly recalls the well-known formulation of the PSDFs, defining the quantities and symbols used in the control framework. In particular, the notation is the same in Section 2.1. If a single PST with a phase shift $\alpha$ is located in series with a line $(s,t)$, the active power through that line can be evaluated as [19,20]:

$$P_{st}=P_{st}^{†}+\alpha ·F_{st}^{\alpha }$$

(31)

where $P_{st}^{†}$ is the value of the flow with $\alpha =0$ (i.e., without PST), while $F_{st}^{\alpha }$ is the PSDF calculated as

$$F_{st}^{\alpha }=\frac{\partial P_{st}}{\partial \alpha }=\frac{1}{X_{st}}\left[1+\frac{1}{X_{st}}\left(2C_{st}^{-1}-C_{ss}^{-1}-C_{tt}^{-1}\right)\right]$$

(32)

where $C^{-1}$ terms are computed according to Equation (9).

It is worth noting that the presence of a single PST influences the flow of the entire network. Thus, the active power flow of a line $(p,q)$ without a PST can be evaluated as

$$P_{pq}=P_{pq}^{†}+\alpha ·F_{pq}^{\alpha }$$

(33)

where again $P_{pq}^{†}$ is the value of the flow with $\alpha =0$, while $F_{pq}^{\alpha }$ is calculated as

$$F_{pq}^{\alpha }=\frac{\partial P_{pq}}{\partial \alpha }=\frac{1}{X_{pq}}\frac{1}{X_{st}}\left(C_{pt}^{-1}-C_{ps}^{-1}+C_{qs}^{-1}-C_{qt}^{-1}\right)$$

(34)

If the network is equipped with more PSTs, Equations (31) and (33) must be modified as follows [19]:

$$P_{st}=P_{st}^{†}+{\alpha }_{st}·F_{st}^{{\alpha }_{st}}+\sum _{\begin{array}{c}m,n\\ (m,n)\ne (s,t)\end{array}}{\alpha }_{mn}·F_{st}^{{\alpha }_{mn}}$$

(35)

for a line $(s,t)$ with a PST, and the presence of other PSTs in lines $(m,n)$, while

$$P_{pq}=P_{pq}^{†}+\sum _{\begin{array}{c}m,n\end{array}}{\alpha }_{mn}·F_{pq}^{{\alpha }_{mn}}$$

(36)

for a line $(p,q)$ without a PST device.

2.3. Uncertainty Modeling

The assessment of the probability distributions for the operating quantities (e.g., branch power flows) is performed starting from correlated and non-Gaussian forecast errors by applying a combination of Point Estimation Method and Third Order Polynomial Normal Transformation [14].

The former calculates the raw moments of the outputs of a generic function which is fed with m stochastic inputs with known probability distributions by evaluating the aforementioned function in specific tuples of the m variables. The latter technique expresses correlated stochastic variables with non-symmetric probability distributions as combinations of normal independent variables.

The quantiles of the probability distributions related to the operating quantities used in the chance constraints of the control are calculated based on the Cornish–Fisher expansion, which computes the quantile of a distribution starting from the statistical moments of the same distribution. In particular, the $\alpha$ quantile $q_{\alpha }$ for a stochastic variable with skewness S and kurtosis K is given in Equation (37).

$$q_{\alpha }={\zeta }_{\alpha }+\frac{({\zeta }_{\alpha }^2-1)}{6}S+\frac{({\zeta }_{\alpha }^3-3{\zeta }_{\alpha })}{24}K-\frac{(2{\zeta }_{\alpha }^3-5{\zeta }_{\alpha })}{36}{S}^2$$

(37)

where ${\zeta }_{\alpha }={\mathsf{\Phi }}^{-1}\left(\alpha \right)$ is the $\alpha$ quantile of the normal distribution with inverse CDF ${\mathsf{\Phi }}^{-1}$. More details on these techniques can be found in [16].

3. The Proposed Chance Constrained Control: Mathematical Formulation

The present section describes the probabilistic control in presence of forecast uncertainties. By adopting the decoupled linear formulation of the AC load flow equations illustrated in Section 2.1, the preventive control is decomposed into two subsequent linear programming problems, respectively, related to active power and to voltage/reactive power-related issues. Section 3.1 describes the mathematical formulation of the optimization algorithm for active power control, while Section 3.2 describes the mathematical formulation of the optimization algorithm for the Volt/Var control.

3.1. Active Power Control

The algorithm related to the active power control is represented by a security-constrained preventive redispatching of active power conventional generation, including also renewable curtailment and PST control. Furthermore, this active power control can consider the forecast uncertainties on renewable energy sources and load through the definition of a stochastic optimization algorithm based on a Chance Constrained (CC) formulation, wherein the constraints can be formulated in terms of the probability of violating operational limits.

The goal of the proposed redispatching problem is to define a power flow optimization model to minimize the redispatching costs and assure the fulfillment of probabilistic N and N-1 security constraint power flows, taking into account controlled uncertainties in renewable generation and load demand. Thus, the active control algorithm is formulated as a linear programming problem that aims to minimize the redispatching costs related to conventional generation and renewable curtailment.

The proposed optimization procedure exploits the linear approximation of AC power flow equations, described in Section 2.1.

The chance constraints concern the technical limit of the generating unit, the invariance of the global generation (composed of conventional and renewable), and active power flow limits in N and N-1 states. The chance constraints are represented by inequalities involving quantiles of uncertain quantities obtained by Cornish–Fisher expansion from cumulants.

The method presented in this section is an evolution of the one presented in [16], including the model of the PSTs in the active power flow management (see Section 2.2 for the PST representation).

The following subsections illustrate the parameters, the variables, and the formulation of the problem (including the reformulation of the chance constraints).

3.1.1. Problem Parameters

This section describes all the parameters for the active power control algorithm.

The considered transmission network is supposed to be composed of: $N_B$ buses, $N_{BR}$ branches, $N_G$ Dispatchable Generation Units (DGUs), $N_W$ Renewable Energy Sources Units (RESUs), $N_{PST}$ PSTs, and $N_D$ loads.

The data supposed to be available are:

• Network buses, branches, and topology parameters. In particular:

  −

  a PTDF matrix, defined in N state, with respect to DGUs, indicated as S (with entries $S_{b_{r}j}$, $b_r=1,\dots ,N_{BR},j=1,\dots ,N_G$) and, with respect to RESUs, indicated as $S_W$ (with entries $S_{W,b_{r}h}$, $b_r=1,\dots ,N_{BR},h=1,\dots ,N_W$), while with respect to load absorption, indicated as $S_D$ (with entries $S_{D,b_{r}d}$, $b_r=1,\dots ,N_{BR},d=1,\dots ,N_D$).

  −

  PTDF matrices computed in N-1 conditions, i.e., in presence of contingencies; specifically a set of $N_{CTG}$ is considered, and the corresponding PDTF matrices are indicated with $S^{\left(k\right)}$, $S_W^{\left(k\right)}$ and $S_D^{\left(k\right)}$, with $k=1,\dots ,N_{CTG}$. Notice that the notation used for the entries of these parameters is the same used for the previous point of this list.

  −

  The branches active power rates, indicated as $T_{lim,b_r}$ with $b_r=1,\dots ,N_{BR}$.

• PST Parameters:

  −

  PSDF matrix in N state indicated as F (with entries $F_{b_r}^{{\alpha }_y}$, $b_r=1,\dots ,N_{BR},y=1,\dots ,$$N_{PST}$).

  −

  PSDF matrices in N-1 conditions computed in presence of a line contingency, represented by $F^{\left(k\right)}$ (with entries $F_{b_r}^{{\alpha }_y\left(k\right)}$, $b_r=1,\dots ,N_{BR},y=1,\dots ,N_{PST}$, $k=1,\dots ,$$N_{CTG}$).

• Maximum and minimum dispatchable power of DGUs: given the j-th DGU, with $j=1,\dots ,N_G$, the minimum and maximum level of dispatchable power are indicated with $P_{max,j}$ and $P_{min,j}$, respectively.
• Maximum power potentially generated by each RESU: the maximum active power that the h-th RESU can generate, with $h=1,\dots ,N_W$, is indicated with $W_{max,h}$.
• Initial power dispatch of each DGU, indicated with $P_{0,j}$ ($j=1,\dots ,N_G$).
• Maximum and minimum phase shift for the PSTs: given the y-th PST, the minimum and the maximum angle are indicated with ${\alpha }_{max,y}$ and ${\alpha }_{min,y}$, respectively.
• Forecasts of the power generated by RESUs: given the h-th RESU, with $h=1,\dots ,N_W$, the forecast of the generated power indicated with $W_{0,h}$. Moreover, forecasts errors $\Delta W_{err,h}$ are modeled by a non-symmetric distribution with known moments up to the 4-th moment.
• Forecasts of the power demanded by the d-th load, with $d=1,\dots ,N_D$, indicated with $L_{0,d}$. Moreover, forecasts errors $\Delta L_{err,d}$ are modeled by a Gaussian distribution with a known standard deviation.
• Initial power flow: given the initial power dispatch and the forecasts of load demand and RES generation, an initial active power flow is defined as $T_{0,ex,l}$ for the $b_r$-th branch, with $b_r=1,\dots ,N_{BR}$.
• Initial power flow under contingency: given the k-th N-1 contingency, with $k=1:N_{CTG}$, an initial active power flow is defined as $T_{0,ex,b_r}^{\left(k\right)}$ for the $b_r$-th branch, with $b_r=1,\dots ,N_{BR}$.
• Costs of upward and downward generation of DGUs: the upward and downward costs of the j-th DGUs are indicated with $c_{+,j}$ and $c_{-,j}$ ($j=1,\dots ,N_G$), respectively.
• Costs of RES curtailment: cost of curtailment of the h-th RESU is indicated with $c_{curt,h}$ ($h=1,\dots ,N_W$).

3.1.2. Problem Variables

The optimization variables for the active power control are:

• Upward and downward power set-point variations of the j-th DGU, indicated with $\Delta P_{+,j}$ and $\Delta P_{-,j}$, respectively. These variables are supposed to be positive ($\Delta P_{+,j},\Delta P_{-,j}\ge 0$).
• Curtailment of the h-th RESU, indicated with $\Delta W_{-,h}$. Moreover, this variable is supposed to be positive ($\Delta W_{-,h}\ge 0$).
• Active power flow expected after redispatching on the $b_r$-th branch without contingencies, indicated with $T_{ex,b_r}$.
• Active power flow expected after redispatching on the $b_r$-th branch with the occurrence of the k-th contingency, indicated with $T_{ex,br}^{\left(k\right)}$.
• Phase shift of the y-th PST, indicated with ${\alpha }_y$.

In addition to the above-listed optimization variables, effectively introduced in the optimization problem, the following variables are defined to model the system uncertainties due to RES and load forecast:

• Contribution of the j-th DGU to regulation to balance forecast errors, indicated with $\Delta P_{R,j}$.
• Variation of the active power flow after redispatching on the $b_r$-th branch without contingencies, indicated with $d{T}_{b_r}$ and such that the actual active power flow $T_b{}_r$ is given by

  $$T_{b_r}=T_{ex,b_r}+d{T}_{b_r},b_r=1,\dots ,N_{BR}$$

  (38)
• Variation of the active power flow after redispatching on the $b_r$-th branch in case of the k-th contingency represented $d{T}_{b_r}^{\left(k\right)}$ and such that the actual relevant active power flow $T_{b_r}^{\left(k\right)}$ is evaluated as

  $$T_{b_r}^{\left(k\right)}=T_{ex,b_r}^{\left(k\right)}+d{T}_{b_r}^{\left(k\right)},b_r=1,\dots ,N_{BR},k=1,\dots ,N_{CTG}$$

  (39)

Finally, a set of feasibility variables is defined to relax the power flow constraints of each branch $b_r$ both with and without contingencies. They are indicated with ${\epsilon }_{+,j}$ and ${\epsilon }_{-,j}$ in the case of no contingencies, and with ${\epsilon }_{+,j}^{\left(k\right)}$ and ${\epsilon }_{-,j}^{\left(k\right)}$, in the case of contingency k-th. Each of these variables is associated with an optimization weight, indicated with $G_{+,j}$ and $G_{-,j}$ (the same for variables with and without contingencies). The higher the values of these optimization weights, the more “rigid” (or less relaxed) will be the corresponding constraints. Before introducing the active power control formulation,

Table 1. Parameters for the optimal active power control algorithm.

Description	Symbol	Size	Unit
Upward generation cost of DGUs	$c_{+,j}$	$1\times N_G$	€/MW
Downward generation cost of DGUs	$c_{-,j}$	$1\times N_G$	€/MW
Cost of RESU curtailment RESU	$c_{curt,h}$	$1\times N_W$	€/MW
Power flow relaxing weight (upward)	$G_{+,b_r}$	$1\times N_{BR}$	€/MW
Power flow relaxing weight (downward)	$G_{-,b_r}$	$1\times N_{BR}$	€/MW
PTDF matrix for DGUs	$S_{b_{r}j}$	$N_{BR}\times N_G$	p.u.
PTDF matrix for RESUs	$S_{W,b_{r}j}$	$N_{BR}\times N_G$	p.u.
PSDF matrix for PSTs	$F_{b_r}^{{\alpha }_y}$	$N_{BR}\times N_{PST}$	p.u.
Initial power flow	$T_{0,ex,b_r}$	$1\times N_{BR}$	kW
PTDF matrix for DGUs (k-th contingency)	$S_{b_{r}j}^{\left(k\right)}$	$N_{CTG}\times N_{BR}\times N_G$	p.u.
PTDF matrix for RESUs (k-th contingency)	$S_{W,b_{r}j}^{\left(k\right)}$	$N_{CTG}\times N_{BR}\times N_G$	p.u.
PSDF matrix for PSTs (k-th contingency)	$F_{b_r}^{{\alpha }_y\left(k\right)}$	$N_{CTG}\times N_{BR}\times N_{PST}$	p.u.
Initial active power flow (k-th contingency)	$T_{0,ex,b_r}^{\left(k\right)}$	$N_{CTG}\times N_{BR}$	MW
Branch power limit	$T_{lim,b_r}$	$1\times N_{BR}$	MW
Maximum phase shift for PSTs	${\alpha }_{max,y}$	$1\times N_{PST}$	rad
Minimum phase shift for PSTs	${\alpha }_{min,y}$	$1\times N_{PST}$	rad
Right quantile on power flow	${\overline{q}}_{dT,b_r}$	$1\times N_{BR}$	MW
Left quantile on power flow	$q_{dT,b_r}$	$1\times N_{BR}$	MW
Right quantile on power flow (k-th cont.)	${\overline{q}}_{dT,b_r}^{\left(k\right)}$	$N_{CTG}\times N_{BR}$	MW
Left quantile on power flow (k-th cont.)	$q_{dT,b_r}^{\left(k\right)}$	$N_{CTG}\times N_{BR}$	MW
Maximum power for DGUs	$P_{max,j}$	$1\times N_G$	MW
Minimum power for DGUs	$P_{min,j}$	$1\times N_G$	MW
Initial set-point for DGUs	$P_{0,j}$	$1\times N_G$	MW
Forecasted generation for RESUs	$W_{0,h}$	$1\times N_W$	MW
Maximum power for RESUs	$W_{max,h}$	$1\times N_W$	MW

and

Table 2. Variables for the optimal active power control algorithm.

Description	Symbol	Size	Unit
Upward power set-point variation for DGUs	$\Delta P_{+,j}$	$1\times N_G$	MW
Downward power set-point variations for DGUs	$\Delta P_{-,j}$	$1\times N_G$	MW
Curtailment for RESUs	$\Delta W_{-,h}$	$1\times N_W$	MW
Power flow relaxing variable (upward)	${\epsilon }_{+,b_r}$	$1\times N_{BR}$	MW
Power flow relaxing variable (downward)	${\epsilon }_{-,b_r}$	$1\times N_{BR}$	MW
Power flow relaxing variable (up.) (k-th cont.)	${\epsilon }_{+,b_r}^{\left(k\right)}$	$N_{CTG}\times N_{BR}$	MW
Power flow relaxing variable (down.) (k-th cont.)	${\epsilon }_{-,b_r}^{\left(k\right)}$	$N_{CTG}\times N_{BR}$	MW
Phase shift for PST	${\alpha }_y$	$1\times N_{PST}$	rad
Power flow expected after redispatching	$T_{ex,b_r}$	$1\times N_{BR}$	MW
Power flow expected after redispatching (k-th cont.)	$T_{ex,b_r}^{\left(k\right)}$	$N_{CTG}\times N_{BR}$	MW

summarize the list of the parameters and the involved variables for the first stage of the proposed procedure.

3.1.3. Problem Formulation

The active power control is formulated as a linear optimization problem whose objective is to minimize the redispatching cost of the DGUs, the cost of the RESU curtailment, and the cost of the relaxing variables:

$$\begin{array}{cc}\hfill min\sum _{j=1}^{N_G}& \left(c_{+,j}\Delta P_{+,j}+c_{-,j}\Delta P_{-,j}\right)+\sum _{h=1}^{N_W}{c}_{curt,h}\Delta W_{-,h}+\hfill \\ \hfill +\sum _{l=1}^{N_{BR}}& \left[G_{+,b_r}{\epsilon }_{+,b_r}+G_{-,b_r}{\epsilon }_{-,b_r}+\sum _{k=1}^{N_{CTG}}\left(G_{+,b_r}{\epsilon }_{+,b_r}^{\left(k\right)}+G_{-,b_r}{\epsilon }_{-,b_r}^{\left(k\right)}\right)\right]\hfill \end{array}$$

(40)

Negligible costs are assumed for PST control actions. In the next part of this section, all the constraints that describe the considered system are reported. These are subdivided into a set of equality and a set of inequality constraints.

Power Balance Equality Constraints

The generic variation of power output for the j-th DGU $\Delta P_j$ must take into account two aspects, namely the variation of the unit set-points $\Delta P_{+,j}$ and $\Delta P_{-,j}$ determined by the active power control, and the automatic load frequency response $\Delta P_{R,j}$ of the unit, which reacts to balance the forecast error of RES generation and load consumption. This means that $\Delta P_j$ can be expressed as follows:

$$\Delta P_j=\Delta P_{+,j}-\Delta P_{-,j}+\Delta P_{R,j}j=1,\dots ,N_G$$

(41)

Moreover, the variation of the generation of the h-th RESU can be expressed as

$$\Delta W_h=-\Delta W_{-,h}+\Delta W_{err,h}h=1,\dots ,N_W$$

(42)

Therefore, the balance of the power equation

$$\sum _{j=1}^{N_G}\Delta P_j+\sum _{h=1}^{N_W}\Delta W_h=\sum _{d=1}^{N_D}\Delta L_{err,d}$$

(43)

becomes

$$\sum _{j=1}^{N_G}\left(\Delta P_{+,j}-\Delta P_{-,j}+\Delta P_{R,j}\right)+\sum _{h=1}^{N_W}\left(-\Delta W_{-,h}+\Delta W_{err,h}\right)=\sum _{d=1}^{N_D}\Delta L_{err,d}$$

(44)

Equation (44) should be transformed taking into account the mentioned automatic frequency responses of DGUs, given by the following relation:

$$\Delta P_{R,j}=-{\zeta }_j\sum _{h=1}^{N_W}\Delta W_{err,h}-\sum _{d=1}^{N_D}\Delta L_{err,d}j=1,\dots ,N_G$$

(45)

with ${\zeta }_1+\dots +{\zeta }_{N_G}=1$. By substituting Equation (45) in Equation (44), it is possible to write:

$$\sum _{j=1}^{N_G}(\Delta P_{+,j}-\Delta P_{-,j})-\sum _{h=1}^{N_W}\Delta W_{-,h}=0$$

(46)

Active Power Flow Equality Constraints

Starting from the initial power flow $T_{0,ex,b_r}$ defined for each branch, after redispatching, the following relation should be satisfied for all $b_r=1,\dots ,N_{BR}$:

$$\begin{array}{cc}\hfill T_{0,ex,b_r}& +\sum _{j=1}^{N_G}\left(S_{b_{r}j}\Delta P_{+,j}-S_{br,j}\Delta P_{-,j}\right)-\sum _{h=1}^{N_W}{S}_{W,b_{r}h}\Delta W_{-,h}+\sum _{y=1}^{N_{PST}}{\alpha }_y{F}_{b_r}^{{\alpha }_y}+{\epsilon }_{-,b_r}\hfill \\ \hfill =T_{ex,b_r}& +{\epsilon }_{+,b_r}{b}_r=1,\dots ,N_{BR}\hfill \end{array}$$

(47)

where $S_{b_r,j}$ is the $(b_r,j)$-th element of the PTDF matrix S, $S_{W,b_{r}h}$ is the $(b_r,h)$-th element for the PTDF matrix $S_W$, while $F_{b_r}^{{\alpha }_y}$ is the PSDF related to the $b_r$-th branch and to the y-th PST. Notice that such an active power flow variation constraint is relaxed by introducing the relaxing variables ${\epsilon }_{-,l}$ and ${\epsilon }_{+,l}$.

By using relation (45), it is possible to obtain that the variation $d{T}_b{}_r$ of the active power flow on the $b_r$-th branch after redispatching is given by the following equation:

$$d{T}_{b_r}=\sum _{h=1}^{N_W}{S}_{W,b_{r}h}\Delta W_{err,h}-\sum _{d=1}^{N_D}{S}_{D,b_{r}d}\Delta L_{err,d}{b}_r=1,\dots ,N_{BR}$$

(48)

Active Power Flow Equality Constraints in case of Contingencies

In case the k-th contingency has occurred, starting from the initial power flow $T_{0,ex,b_r}^{\left(k\right)}$ defined for each branch, after redispatching the following relation should be satisfied for all $b_r=1,\dots ,N_{BR}$:

$$\begin{array}{cc}\hfill T_{0,ex,b_r}^{\left(k\right)}& +\sum _{j=1}^{N_G}\left(S_{b_{r}j}^{\left(k\right)}\Delta P_{+,j}-S_{br,j}^{\left(k\right)}\Delta P_{-,j}\right)-\sum _{h=1}^{N_W}{S}_{W,b_{r}h}^{\left(k\right)}\Delta W_{-,h}+\sum _{y=1}^{N_{PST}}{\alpha }_y{F}_{b_r}^{{\alpha }_y\left(k\right)}+{\epsilon }_{-,b_r}^{\left(k\right)}\hfill \\ \hfill =T_{ex,b_r}^{\left(k\right)}& +{\epsilon }_{+,b_r}^{\left(k\right)}{b}_r=1,\dots ,N_{BR},k=1,\dots ,N_{CTG}\hfill \end{array}$$

(49)

where $S_{b_r,j}^{\left(k\right)}$ is the $(b_r,j)$-th element for the PTDF matrix $S^{\left(k\right)}$, $S_{W,b_{r}h}^{\left(k\right)}$ is the $(b_r,h)$-th element for the PTDF matrix $S_W^{\left(k\right)}$, while $F_{b_r}^{{\alpha }_y}$ is the PSDF related to the $b_r$-th branch and to the y-th PST in case of the k-th contingency. Notice that such an active power flow variation constraint is relaxed by introducing the relaxing variables ${\epsilon }_{-,l}^{\left(k\right)}$ and ${\epsilon }_{+,l}^{\left(k\right)}$.

By using relation (45), it is possible to obtain that the variation $d{T}_b{{}_r}^{\left(k\right)}$ of the active power flow after redispatching on the $b_r$-th branch is given by the following equation:

$$d{T}_{b_r}^{\left(k\right)}=\sum _{h=1}^{N_W}{S}_{W,b_{r}h}^{\left(k\right)}\Delta W_{err,h}-\sum _{d=1}^{N_D}{S}_{D,b_{r}d}^{\left(k\right)}\Delta L_{err,d}{b}_r=1,\dots ,N_{BR},k=1,\dots ,N_{CTG}$$

(50)

where $S_{D,b_{r}d}^{\left(k\right)}$ is the $(b_r,d)$-th element of the PTDF matrix $S_D^{\left(k\right)}$.

Inequality Chance Constraints

The inequality constraints concern both the branch power flows, and the generator active power limits and are expressed in probabilistic terms by the following chance constraints:

• Active power flow limits at steady state on the $b_r$-th branch:

  $$\begin{array}{cc}\hfill \mathbb{P}\left(T_{b_r}<T_{lim,b_r}\right)& >1-{\epsilon }_T{b}_r=1,\dots ,N_{BR}\hfill \end{array}$$

  (51)

  $$\begin{array}{cc}\hfill \mathbb{P}\left(T_{b_r}>-T_{lim,b_r}\right)& >1-{\epsilon }_T{b}_r=1,\dots ,N_{BR}\hfill \end{array}$$

  (52)

  where ${\epsilon }_T$ is the relevant violation probability threshold. These inequalities impose that the probability (indicated with the notation $\mathbb{P}$) of violating the limits on the active power flows is greater than $1-{\epsilon }_T$.
• Active power flow limits after the k-th contingency on the $b_r$-th branch:

  $$\begin{array}{cc}\hfill \mathbb{P}\left(T_{b_r}^{\left(k\right)}<T_{lim,b_r}\right)& >1-{\epsilon }_{T_1}{b}_r=1,\dots ,N_{BR}\hfill \end{array}$$

  (53)

  $$\begin{array}{cc}\hfill \mathbb{P}\left(T_{b_r}^{\left(k\right)}>-T_{lim,b_r}\right)& >1-{\epsilon }_{T_1}{b}_r=1,\dots ,N_{BR}\hfill \end{array}$$

  (54)

  where ${\epsilon }_{T_1}$ is the violation probability threshold related to these inequality constraints. The meaning of these inequalities is the same expressed by Equations (51) and (52).
• Active power set-points of the j-th DGU:

  $$\begin{array}{cc}\hfill \mathbb{P}\left(\Delta P_j<P_{max,j}-P_{0,j}\right)& >1-{\epsilon }_{gen}j=1,\dots ,N_G\hfill \end{array}$$

  (55)

  $$\begin{array}{cc}\hfill \mathbb{P}\left(\Delta P_j>-P_{min,j}+P_{0,j}\right)& >1-{\epsilon }_{gen}j=1,\dots ,N_G\hfill \end{array}$$

  (56)

  where ${\epsilon }_{gen}$ is the violation probability threshold related to these constraints. Equations (55) and (56) force that the probability of violating the technical limits of the DGUs after the redispatch (i.e., after the application of the active power control) is greater than $1-{\epsilon }_{gen}$.
• Active power injection limit of the h-th RESU:

  $$\mathbb{P}\left(\Delta W_h+W_{0,h}>0\right)>1-{\epsilon }_{W}h=1,\dots ,N_W$$

  (57)

  where ${\epsilon }_W$ is the relevant violation probability threshold. This constraint is related to the maximum possible curtailment that the active power control can exploit. The meaning of Equation (57) is similar to the downward reserve for the DGUs.

Reformulation of the Chance Constraints

To write the proposed optimization procedure as a linear programming problem, the probabilistic Constraints (51)–(57) have to be reformulated as chance constraints.

As reported in Section 2.3, uncertainty variables $\Delta L_{err,d}$, $\Delta W_{err,h}$, $\Delta P_{R,j}$, $d{T}_{b_r}$, and $d{T}_{b_r}^{\left(k\right)}$ probability distributions are characterized as detailed in [16], using a combination of Point Estimate Method (PEM), Third order Polynomial Normal Transformation (TPNT) and Cornish–Fisher Expansion that allows taking into account their mutual correlation. As a result, the following left and right quantiles are defined:

• $q_{dT,b_r}$: left quantile of $d{T}_{b_r}$, given the violation probability ${\epsilon }_T$.
• ${\overline{q}}_{dT,b_r}$: right quantile of $d{T}_{b_r}$, given the violation probability ${\epsilon }_T$.
• $q_{dT,b_r}^{\left(k\right)}$: left quantile of $d{T}_{b_r}^{\left(k\right)}$, given the violation probability ${\epsilon }_{T1}$.
• ${\overline{q}}_{dT,b_r}^{\left(k\right)}$: right quantile of $d{T}_{b_r}^{\left(k\right)}$, given the violation probability ${\epsilon }_{T1}$.
• $q_j$: left quantile of $\Delta P_{R,j}$, given the violation probability ${\epsilon }_G$.
• ${\overline{q}}_j$: right quantile of $\Delta P_{R,j}$, given the violation probability ${\epsilon }_G$.
• $q_{W,h}$: left quantile of $\Delta W_{err,h}$, given the violation probability ${\epsilon }_W$.

Using these quantiles, chance Constraints (51)–(57) are transformed into the following linear inequalities (58)–(64):

• Active power flow limits at steady state on the $b_r$-th branch:

  $$\begin{array}{cc}\hfill T_{ex,b_r}& \le T_{lim,b_r}-{\overline{q}}_{dT,b_r}\forall b_r=1:N_{BR}\hfill \end{array}$$

  (58)

  $$\begin{array}{cc}\hfill T_{ex,b_r}& \ge {-T}_{lim,b_r}-q_{dT,b_r}\forall b_r=1:N_{BR}\hfill \end{array}$$

  (59)
• Active power flow limits after the k-th contingency on the $b_r$-th branch:

  $$\begin{array}{cc}\hfill T_{ex,b_r}^{\left(k\right)}& \le T_{lim,b_r}-{\overline{q}}_{dT,b_r}^{\left(k\right)}\forall k=1,\dots ,N_{CTG},\forall b_r=1:N_{BR}\hfill \end{array}$$

  (60)

  $$\begin{array}{cc}\hfill T_{ex,b_r}^{\left(k\right)}& \ge {-T}_{lim,b_r}-q_{dT,b_r}^{\left(k\right)}\forall k=1,\dots ,N_{CTG},\forall b_r=1:N_{BR}\hfill \end{array}$$

  (61)
• Active power set-points of the j-th DGU:

  $$\begin{array}{cc}\hfill \Delta P_{+,j}-\Delta P_{-,j}& \le P_{max,j}-P_{0,j}-{\overline{q}}_j\forall j=1,\dots ,N_G\hfill \end{array}$$

  (62)

  $$\begin{array}{cc}\hfill \Delta P_{+,j}-\Delta P_{-,j}& \ge P_{min,j}-P_{0,j}-q_j\forall j=1,\dots ,N_G\hfill \end{array}$$

  (63)
• Active power injection limit of the h-th RESU:

  $$\Delta W_{-,h}\le W_{0,h}+q_{W,h}\forall h=1,\dots ,N_W$$

  (64)

PST Inequality Constraints

Finally, Constraints (65) and (66) impose the limits on the angles related to the PSTs.

$$\begin{array}{cc}\hfill {\alpha }_y& \le {\alpha }_{max,y}\forall y=1,\dots ,N_{PST}\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill {\alpha }_y& \ge {\alpha }_{min,y}\forall y=1,\dots ,N_{PST}\hfill \end{array}$$

(66)

Summarizing, the active power control optimization algorithm is composed of the objective Function (40) and the Constraints (46), (47), (49), (58)–(66).

3.2. Volt/Var Control

This section describes the mathematical formulation of the optimization algorithm for the Volt/Var control. The scheme of the proposed two-stage strategy is to use the Volt/Var algorithm after the new working point provided by the active power control (described in Section 3.1.3) to manage the technical limits related to the node voltages, the reactive power of the available generators (DGUs and RESUs), and the reactive power flows. If these constraints are not satisfied, the algorithm described in this section attempts to solve this problem by adjusting the voltage set-points of the available generators. This formulation is based on the decoupled linearization of the AC power flow (see Section 2.1) that allows the exploitation of the proposed Volt/Var methodology in series to the stochastic optimization algorithm for the active power control.

After the first stage of the proposed procedure, the network model is updated with the output provided by the active power control (in terms of DGUs redispatching, RESUs curtailment, and the angle modifications of the PSTs). Finally, this new model is taken as input by the Volt/Var procedure that exploits the linear AC power flow previously described (see Section 2.1).

3.2.1. Problem Parameters

This section summarizes all the parameters for the Volt/Var algorithm (notice that the same notation presented for Section 2.1 and Section 3.1 is used, wherein $N_B$ is the number of buses, $n_1$ is the number of PQ nodes, $n_2$ is the number of PV nodes and $N_{GEN}=N_G+N_W$ is the number of generators):

• Network buses, branches, and topology parameters. In particular: initial voltage at network buses ($V_{0,g},g=1,\dots ,n_2$), initial reactive power at PQ nodes ($Q_{0,l},l=1,\dots ,n_1$), initial reactive power at PV nodes ($Q_{0,g},g=1,\dots ,n_2$), maximum and minimum voltage limits at all the grid nodes ($V_{max,b}$, $V_{min,b},b=1,\dots ,N_B$), and reactive power flow limits for all the branches ($Q_{flowlim,b{b}_1},b,b_1=1,\dots ,N_B$).
• Maximum and minimum reactive power limits ($Q_{max,z}$, $Q_{min,z},z=1,\dots ,N_{GEN}$) for all the generators (DGUs and RESUs).
• AC linear power flow parameters: matrices ${\left(M^{\prime }\right)}^{-1}$, $L^{\prime }$, $D^{\prime }$, $J^{\prime }$, $\beta$, $A^{\prime }$, $\gamma$, and vector W are evaluated according to the relations presented in Section 2.1, starting from the network model updated after the active power control.

Finally,

Table 3. Parameters for the optimal Volt/Var control algorithm.

Description	Symbol	Size	Unit
Initial voltage at PV nodes	$V_{0,g}$	$1\times n_2$	p.u.
Maximum voltage	$V_{max,b}$	$1\times N_B$	p.u.
Minimum voltage	$V_{min,b}$	$1\times N_B$	p.u.
Initial reactive power at PQ nodes	$Q_{0,l}$	$1\times n_1$	p.u.
Initial reactive power at PV node	$Q_{0,g}$	$1\times n_2$	p.u.
Maximum reactive power flow	$Q_{flowlim,b{b}_1,}$	$N_B\times N_B$	p.u.
Maximum reactive power for generator	$Q_{max,z}$	$1\times N_{GEN}$	p.u.
Minimum reactive power for generator	$Q_{min,z}$	$1\times N_{GEN}$	p.u.
Matrix defined in Equation (21)	${\left(M^{\prime }\right)}_{l{l}_1}^{-1}$	$n_1\times n_1$	−
Vector defined in Equation (19)	$W{}_l$	$1\times n_1$	−
Matrix defined in Equation (21)	$L_{lg}^{\prime }$	$n_1\times n_2$	−
Matrix defined in Equation (23)	$D_{gl}^{\prime }$	$n_2\times n_1$	−
Matrix defined in Equation (21)	$J_{g{g}_1}^{\prime }$	$n_2\times n_2$	−
Matrix defined in Equation (28)	${\beta }_{b{b}_1}$	$N_B\times N_B$	−
Matrix defined in Equation (18)	$A_{b{b}_1}^{\prime }$	$N_B\times N_B$	−
Matrix defined in Equation (29)	${\gamma }_{b{b}_1}$	$N_B\times N_B$	−

summarizes all the parameters for the Volt/Var optimization algorithm.

3.2.2. Optimization Variables

In this section, all the variables of the Volt/Var algorithm are summarized:

• Upward and downward voltage set-point variations at the b-th bus, indicated as $\Delta V_{+,b}$, $\Delta V_{-,b}$; both variables are supposed to be positive ($\Delta V_{+,b}$, $\Delta V_{+,b}\ge 0$).
• Voltage magnitude at the b-th node after the Volt/Var control, indicated with $V_b$. For the algorithm implementation also two sub-variables of $V={[V_1,\dots ,V_{N_B}]}^T={[V_l,V_g]}^T$ have been defined. This definition follows Equation (21) within the linearization of the AC power flow (see Section 2.1).
• Reactive power injection at the l-th PQ node ($Q_l$) and at the g-th PV node ($Q_g$) after the Volt/Var optimization.
• Reactive power related to the z-th generator (${\widehat{Q}}_z$) after the control application.
• Reactive power expected after redispatching on branch $(b,b1)$, represented with $Q_{flow,b{b}_1}$.
• Auxiliary variable for the l-th PQ nodes, indicated with $H_l^{\prime }$ (see Equation (24) of the AC linear power flow for its definition).

Finally,

Table 4. Variables for the optimal Volt/Var control algorithm.

Description	Symbol	Size	Unit
Node voltage	$V_b$	$1\times N_B$	p.u.
Voltage at PQ nodes	$V_l$	$1\times n_1$	p.u.
Voltage at PV nodes	$V_g$	$1\times n_2$	p.u.
Upward variation of the voltage set-point	$\Delta V_{+,b}$	$1\times N_B$	p.u.
Downward variation of the voltage set-point	$\Delta V_{-,b}$	$1\times N_B$	p.u.
Reactive power at PQ nodes	$Q_l$	$1\times n_1$	p.u.
Reactive power at PV nodes	$Q_g$	$1\times n_2$	p.u.
Generator reactive power	${\widehat{Q}}_z$	$1\times N_{GEN}$	p.u.
Reactive power flow	$Q_{flow,b{b}_1}$	$1\times 2N_B$	p.u.
Auxiliary variable defined in Equation (24)	$H_l^{\prime }$	$1\times n_1$	−

summarizes all the Volt/Var control algorithm variables with information about the size and the measurement unit.

3.2.3. Problem Formulation

The basic principle of this control stage, defined as a linear optimization problem, is to perform the minimum amount of voltage set-point adjustments to satisfy the constraints on voltages and reactive powers. More specifically, the problem minimizes the variations ($\Delta V_{+,b},\Delta V_{-,b}\ge 0$ for the b-th node) in the voltage set-points of the generators:

$$min\sum _{b=1}^{N_B}\left[\Delta V_{+,b}+\Delta V_{-,b}\right]$$

(67)

The sequel of the subsection describes the constraints of the problem. Equation (68) imposes the invariance of the reactive powers absorbed at PQ nodes.

$$Q_l=Q_{0,l}\forall l=1,\dots ,n_1$$

(68)

where $Q_l$ is the reactive power at the l-th PQ node, while $Q_{0,l}$ is the reactive power at the same node before the control application. Equation (69) defines $V_l$, i.e., the voltage at the l-th PQ node, as a function of $Q_l$ (see Equation (68)), and $V_g$ (voltage magnitude at the PV nodes) through the auxiliary Equation (70), according to the linear AC power flow formulation described in Section 2.1.

$$V_l=\sum _{l1=1}^{n_1}\left[{\left(M^{\prime }\right)}_{l{l}_1}^{-1}{Q}_l+{\left(M^{\prime }\right)}_{l{l}_1}^{-1}{H}_{l_1}^{\prime }\right]\forall l=1,\dots ,n_1$$

(69)

$$H_l^{\prime }=W_l-\sum _{g=1}^{n_2}{L}_{lg}^{\prime }{V}_g\forall l=1,\dots ,n_1$$

(70)

Constraint (71) evaluates the updated value of the voltage at PQ nodes adding the possible contribution given by $\Delta V_{+,g}$ and $\Delta V_{-,g}$.

$$V_g=V_{0,g}+\Delta V_{+,g}-\Delta V_{-,g}\forall g=1,\dots ,n_2$$

(71)

The following constraints indicate that the voltage magnitudes at PQ nodes can be affected by voltage set-point variations only at PV nodes:

$$\begin{array}{cc}\hfill \Delta V_{+,l}& =0\forall l=1,\dots ,n_1\hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill \Delta V_{-,l}& =0\forall l=1,\dots ,n_1\hfill \end{array}$$

(73)

The next equations represent the constraints of maximum ($V_{max,b}$ for the b-th node) and minimum ($V_{min,b}$ for the b-th node) voltage on the network buses:

$$\begin{array}{cc}\hfill V_b& \le V_{max,b}\forall b=1,\dots ,N_B\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill V_b& \ge V_{min,b}\forall b=1,\dots ,N_B\hfill \end{array}$$

(75)

The following constraint relates the reactive power at PV nodes ($Q_g$ for the g-th PV bus) with the reactive power at PQ nodes ($Q_l$) and the voltage at the PV nodes ($V_g$), following the linearization approach proposed in Section 2.1 (see Equation (23)).

$$Q_g=\left(\sum _{l=1}^{n_1}{D}_{gl}^{\prime }{Q}_l+D_{gl}^{\prime }{H}_l^{\prime }\right)+\left(\sum _{g_1=1}^{n_2}{J}_{g{g}_1}^{\prime }{V}_{g_1}\right)-W_g\forall g=1,\dots ,n_2$$

(76)

Equation (77) allows evaluating the actual reactive power contribution related only to the generators for all the PV nodes.

$$\sum _{z\in Se{t}_g}{\widehat{Q}}_z=Q_g-Q_{0,g}\forall g=1,\dots ,n_2$$

(77)

In Equation (77), $Se{t}_g$ indicates the set of the generators at the g-th PV node.

The following constraints force the maximum ($Q_{max,z}$) and minimum ($Q_{min,z}$) reactive power for the z-th generator.

$$\begin{array}{cc}\hfill {\widehat{Q}}_z& \le Q_{max,z}\forall z=1,\dots ,N_{GEN}\hfill \end{array}$$

(78)

$$\begin{array}{cc}\hfill {\widehat{Q}}_z& \ge Q_{min,z}\forall z=1,\dots ,N_{GEN}\hfill \end{array}$$

(79)

Equation (80) gives the expression of the reactive power flow of the line $(b,b_1)$ (still according to the AC linear power flow formulation, see Equation (27)).

$$Q_{flow,b{b}_1}={\beta }_{b{b}_1}{V}_b+A_{b{b}_1}^{\prime }{V}_{b_1}+{\gamma }_{b{b}_1}\forall b,b_1=1,\dots ,N_B$$

(80)

The next constraints impose the limits on the maximum Equation (81) and minimum Equation (82) reactive power flow for line $(b,b_1)$.

$$\begin{array}{cc}\hfill Q_{flow,b{b}_1}& \le Q_{flowlim,b{b}_1}\forall b,b_1=1,\dots ,N_B\hfill \end{array}$$

(81)

$$\begin{array}{cc}\hfill Q_{flow,b{b}_1}& \ge -Q_{flowlim,b{b}_1}\forall b,b_1=1,\dots ,N_B\hfill \end{array}$$

(82)

Summing up, the formulation proposed for the Volt/Var control is composed by the objective Function (67) and the Constraints defined in Equations (68)–(82).

4. Study Case

This section describes the considered test system and the simulation scenarios for the proposed control strategy.

4.1. Test System

The presented strategy has been tested and validated on the IEEE 118 bus system [23], where 10 out of 54 synchronous generators (at buses 46, 49, 61, 54, 59, 65, 91, 103, 105, 110) have been replaced by RESU injections (5 PV plants and 5 wind plants, see Figure 1) with the same active and reactive power rating of the original DGUs. Notice that the capability curves of all the generators have been supposed to be rectangular. This simplification allows the problem decoupling and a linear formulation for the proposed procedure. Moreover, since the branch ratings are not assigned in [23], they are set to $150\%$ of the power flows in state N, with a minimum value of 150 MVA, with the exception of line 68–81, with a rating equal to the $110\%$ of the initial power flow.

The feasible phase shift angles belong to the interval $[-{30}^{\circ },+{30}^{\circ }]$, while the minimum and the maximum bus voltage at all the buses have been set, respectively, equal to 0.94 p.u. and 1.06 p.u.

The redispatching costs considered for the DGUs have been set to $c_{+,j},c_{-,j}=10$ AMU/MW (Arbitrary Monetary Unit). A more precise selection of the costs can be performed considering the markets, the type of generator, the cost of the fuel, etc., while the RESU curtailment has been strongly penalized with $c_{curt,j}={10}^4$ AMU/MW. Finally, the costs related to the relaxing variables have been set very high, i.e., $G_{+,b_r},G_{-,b_r}={10}^5$ AMU/MW for all the lines of IEEE 118 bus system.

The set of contingencies considered in the simulations is composed of the N-1 outages related to all the lines (excluding radial connections) with an initial active power flow higher than a threshold (set to 150 MW in the present work).

In line with [25], the forecast errors for renewable injections are modeled as a beta distribution centered on the mean forecast error. Unless differently specified, the standard deviations for the forecast error distributions of aggregated renewable injections over a 24-h forecast horizon are equal to $5\%$ of the rated power of each injection.

As far as the correlation among forecast errors are concerned, a Toeplitz matrix is considered with a generating vector linearly decreasing from 1 to 0.1 with steps of width 0.14.

Parameters related to the linearization of the AC power flow, PTDF/PSDF matrices, initial active/reactive power flows, and active/reactive injections (see Table 1 and Table 3), both in state N and in N-1 scenario have been calculated using Matpower [26].

4.2. Simulation Scenarios

This section describes the simulation scenarios for validating the linearization techniques (concerning the AC power flow equations and the PST modeling) and for applying the proposed two-stage formulation.

Table 5. Simulation scenarios.

ID	Description	Goal
Validation of Linearization Techniques (VLT)	Comparison of the voltages and power flows obtained from the linearized models and the full AC models	Validation of the modeling techniques for AC power flow and PSTs
Active Power Control 1 (APC1)	Application of the active power control without the use of PSTs	Validation of the first stage of the proposed control
Active Power Control 2 (APC2)	Application of the active power control with the presence of a PST	Evaluation of the effectiveness of the first stage of the proposed control and the potential of the PSTs
Volt/Var Control (VVC)	Application of the Volt/Var control after the APC2 scenario	Validation of the second stage of the proposed control

summarizes the simulation scenarios reporting for each of them the denomination (ID), a brief description, and its main goals.

4.3. Validation of the Linearization Techniques

The validation of the linearized models integrated into the proposed control strategy is a fundamental step to guarantee the accuracy of the outcomes of the control itself. The AC power flow linearization described in Section 2.1 has been validated on the IEEE 118 bus system exploiting Matpower as a benchmark. In particular, the function runpf of Matpower has been used [26]. This function performs a full AC power flow.

The Key Performance Indicator (KPI) used in this comparison has been the Mean Absolute Error (MAE):

$$\mathrm{MAE}=\frac{1}{N}\sum _{i=1}^N|x_i-{\widehat{x}}_i|$$

(83)

where N is the number of samples, $x_i$ is the actual value and ${\widehat{x}}_i$ is its estimation. In this case, MAE has been evaluated to have information on the error committed on: angles, voltage magnitudes at PQ nodes, active/reactive power injections at PV nodes, and active/reactive power flows.

To validate the PST modeling described in Section 2.2, several simulations have been performed, including a PST in series with lines 30–17. These tests have been executed varying the phase shift ($\alpha ={20}^{\circ },{15}^{\circ },-5^{\circ },-{10}^{\circ }$). Again the benchmark exploited has been the DC load flow of Matpower (rundcpf function [26]), while the KPI has been the Mean Absolute Percentage Error (MAPE), defined as:

$$\mathrm{MAPE}=\frac{1}{N}\sum _{i=1}^N\frac{|x_i-{\widehat{x}}_i|}{x}$$

(84)

where all the quantities involved have been already defined for MAE in Equation (83). In this case, MAPE has been computed on the active power flows.

Table 6. KPI for the AC power flow linearization.

Quantity	MAE
Angles	1.56${}^{\circ }$
Voltage at PQ nodes	0.0006 p.u.
Active power at PV nodes	1.45 MW
Reactive power at PV nodes	2.59 MVar
Active power flows	2.62 MW
Reactive power flows	1.36 MVar

reports the errors in terms of MAE. Figure 2 shows the results related to the approximation of the voltage at PQ nodes, the reactive power at PV nodes, and the reactive power flows. In all these figures, the output of the proposed procedure is depicted in red, while the benchmark is represented in blue.

As seen from the KPIs evaluated in Table 6 and the plots reported in Figure 2, the results are very satisfactory, confirming that this linearization procedure can be exploited within the presented two-stage strategy.

Table 7. Results of the approximation of the active power flows with the presence of a PST.

Phase Shift	20°	15°	−5°	−10°
MAPE	1.56	2.07	1.57	3.21

summarizes the results in terms of MAPE for the approximation of the active power flows with the presence of a PST and for different values of the phase shift.

Figure 3 reports the results obtained with a phase shift equal to 20° (in red) and the comparison with the benchmark (in blue).

The results reported in Table 7 and Figure 3 confirm the accuracy of this model (already proposed in [19,20]), which can be employed to evaluate the possible contribution of the PST in managing the active power flows.

4.4. Active Power Control

Two scenarios have been introduced to evaluate the first stage (active power control) of the proposed procedure. In the first scenario (that will be indicated as Scenario APC1), no PSTs have been considered, while in the second (Scenario APC2), a PST is located in series with line 30–17. Therefore, the comparison between the values of the objective function obtained with these two scenarios will help to evaluate the possible contributions of the PSTs. This subsection presents the analysis of the results related to both scenarios.

4.4.1. Scenario APC1

The value of the objective function for this scenario is 2298.93 AMU, while the computational time is 0.117 s (obtained with a processor Intel i7, 16 GB RAM).

Figure 4 reports the results in terms of upward and downward set-point variation of DGUs (in blue and red, respectively) and the curtailment of RESUs (in yellow). It is worth noting that there is no curtailment since it is strongly penalized in the objective function (see Equation (40)) of the active power control.

4.4.2. Scenario APC2

This scenario aims to analyze the possible benefits related to the exploitation of PSTs in operating and economic terms for the optimal redispatching.

The value of the objective function passes from 2298.93 AMU (Scenario A) to 2266.33 AMU, proving that the presence of just one PST already brings an interesting advantage both in economic terms (decrease in the objective function imposing a phase shift equal to −3.20°), and in operational terms (all the constraints are satisfied with smaller variations in the set-points of the DGUs).

Figure 5 shows the results obtained in this scenario regarding variations of the active power set-points for DGUs and RESUs.

As can be seen from Figure 5, the presence of a PST allows a lower redispatching of the DGUs (see Figure 4 for Scenario A without PSTs).

Figure 6 illustrates the details of the DGUs. In particular, this figure reports the final set-points of DGUs after the application of the control (in red), the initial set-points of the DGUs (blue line), and the maximum/minimum active power limits (dotted black lines).

As demonstrated by Figure 5 and Figure 7, there is no curtailment for the RESUs (in fact, red and blue curves overlap in Figure 7) for the same motivations described in Scenario APC1. It is worth noting that in Figure 7, the black dotted line indicates the maximum level of curtailment (downward set-point variation of a RESU) imposed by the probabilistic Constraint (64).

Finally, Figure 8 reports the active power flows on network branches before the application of the control (in state N ($T_{0,ex,b_r}$) and in contingency operation ($T_{0,ex,b_r}^{\left(k\right)}$, $k=1,\dots ,N_{CTG}$)), while Figure 9 shows the same quantities after the active power control (expected power flows in state N ($T_{ex,b_r}$), and in case of contingencies ($T_{ex,b_r}^{\left(k\right)},k=1,\dots ,N_{CTG}$)). As can be seen from Figure 8, there are some branches with a flow out of the technical limit (see orange arrows), while after the application of the control (see Figure 9), these violations are solved and all the constraints are satisfied.

4.5. Volt/Var Control

The second stage of the proposed control strategy (Volt/Var control) has been tested after the optimal redispatching presented by the active power control on Scenario APC2. Before applying the Volt/Var control, the working point of the IEEE 188 bus system has been updated according to the output presented in Section 4.4.2 (see scenario VVC in Table 5).

The value of the objective function representing the sum of the modifications in p.u. of the voltage set-points is equal to 0.03 per unit. This means that the network does not need particular interventions to satisfy the constraints related to voltage/reactive power. The calculation time has been 0.079 s (obtained with a processor Intel i7-16 GB RAM). This confirms that the linear formulation allows having a low computational time and that this approach could be interesting with large electrical networks.

Figure 10 shows the results regarding the voltage profile. It is possible to observe the voltage magnitude before the Volt/Var control (in blue), the voltage profile obtained by the optimization (in red), and the voltage limits [0.94–1.06 p.u.] (in dotted black lines).

Figure 11 shows the results related to the reactive power of the generators. Moreover, in this case, the output of the optimization is depicted in red, the state before the application of the control is reported in blue, and the technical limits are reported in black. From this figure, it is worth noting that the proposed procedure can be useful in some cases to fulfill the constraints related to the maximum and minimum reactive power. Finally, in this scenario, there are no violations regarding the reactive power flow limits on the network branches.

5. Conclusions

This paper has proposed a probabilistic preventive control aimed at assuring N-1 security in presence of correlated uncertainties on renewable production and load demand. Exploiting an AC power flow linearization, the procedure is decoupled into two successive linear optimization problems. The first stage concerns active power control. Its goal is to find the optimal conventional generation redispatching, the renewable curtailment, and the PST settings at the minimum cost to assure the fulfillment of N and N-1 security constraints on branch active power flows under forecast uncertainties. The second stage is composed of another linear optimization algorithm, which aims to minimize the adjustments to generators’ set-point voltages to avoid violations at node voltages, reactive power flows, and generator reactive power limits. The approach has been tested by applications to the IEEE 118 bus system.

Future developments may address performance evaluation of the control methodology in the case of large electrical networks, with possible application of decomposition techniques (e.g., Benders decomposition) to further improve the computational speed.