Robust Placement and Control of Phase-Shifting Transformers Considering Redispatch Measures

Abstract

Flexible AC transmission systems (FACTSs) can maximize capacity utilization under time-varying grid usage patterns by actively controlling the power flow of the transmission lines, e.g., with phase-shifting transformers (PST). In this paper, we propose an algorithm to determine the minimum number of PSTs and their location such that the grid can operate robustly for any realization of the (active) power set points from a known, continuous uncertainty set. As we show in our experiments, only considering a few extreme grid scenarios cannot provide this guarantee. The proposed algorithm considers the trade-offs between PST placement and operational decisions, such as PST control and redispatch. By minimizing the worst-case redispatch cost, it yields two affine linear control policies for these as a byproduct. Power flow is modeled as a constrained linear system, and the control design and actuator minimization tasks are formulated as a mixed-integer linear program (MILP). We also design a greedy algorithm, whose optimal value differs less than 20% from the MILP solution while being one to two orders of magnitude faster to compute. The proposed algorithm is evaluated for a small demonstrative 3-bus example and the IEEE 39 bus test system.

1. Introduction

Variable renewable energies from wind and sun are important to reduce the global carbon footprint of energy systems. At the same time, their variability creates novel challenges for both the design and the operation of transmission grids. Especially when large renewables capacities are localized in few regions with good natural resources, their power generation may lead to line overloads in the grid. Redispatch is a common measure to counteract these congestions, but it may be expensive. Grid capacity extensions are a last resort but typically require a long time to implement and are subject to societal resistances. Thus, utilizing active power flow control technologies in the grid is an attractive and often cost-effective option for grid operators. Flexible alternating current transmission systems (FACTSs) utilize power electronics to enhance the grid’s controllability [1]. Examples include static var compensators (SVC), thyristor-controlled series capacitors (TCSCs), and thyristor-controlled phase-shifting transformers (PST), which are part of various grid development plans, e.g., in Germany (www.netzentwicklungsplan.de (accessed on 2 February 2023)).

An important question that naturally arises is where to place such FACTS devices in the power grid and how to control them. Maximizing the loadability of the system for one demand scenario considering various FACTS devices can be solved by genetic algorithms [2] or by the alternating direction method of multiplier (ADMM) [3]. Mixed integer linear programming (MILP) can be used for PST [4] or SVC [5] placement in the same one-scenario setting. Several scenarios are considered in approaches based on two-stage stochastic programming [6,7,8], which allows to account for the expected scenario-dependent generation costs as well. All of these approaches do not guarantee system feasibility or the (re-)dispatch costs for situations outside the (typically few) considered scenarios. In our experimental section, we show that relying only on extreme generation situations may be misleading about the grid’s feasibility for intermediate in-feeds.

In this paper, we propose an algorithm to find the minimal number and location of PSTs that can guarantee system feasibility and limit redispatch costs for all scenarios in a continuous uncertainty set. Considered uncertainty factors are renewable, and load fluctuations as well as the corresponding dispatch decisions are determined by markets or load–frequency control, see Figure 1. During PST placement, we take into account the grid operators’ situation-dependent control options, such as specifying PST settings and redispatch orders to the attached power plants. Our algorithm thus produces, as a byproduct of the PST placement task, control policies to dynamically determine the PST settings and the redispatch orders. These policies are assumed to be affine linear in the uncertain factors in our work. They are chosen optimally in the sense that they minimize the worst-case redispatch cost in this class of controls. Since we consider the often strongly non-linear dispatch decisions by markets and frequency control as part of the uncertainty set, the linearity assumption for small to medium-sized redispatch corrections is plausible.

In our proposed method, we model the power flow in the grid as a set of linear equations, similar to the approach used in [9,10]. However, we employ a robust optimization framework instead of relying on chance constraints to model the uncertainty in power injections. Doing so allows us to write the continuous part of the placement and control problem as a linear program, which is easier to solve than the second-order cone programs derived from chance-constrained optimization problems. The authors of [11] also use a robust optimization approach in the context of reserve scheduling in power systems with highly uncertain power injections. In our work, we extend this approach by also considering the optimal location of PSTs. The resulting robust optimization problem can be solved with MILP. It is worth noting that fuzzy logic represents a widely used alternative to handle uncertainty in power systems [12,13,14]. However, fuzzy methods may encounter scalability issues when dealing with the large-scale problems often encountered in power systems. Consequently, heuristics may be necessary to reduce solving times [15].

The uncertainty set considered in the robust optimization should be large enough to cover at least a selection of plausible scenarios and all mixtures of these. At the same time, the uncertainty set should not be too large in order to avoid excessive conservatism of the solution. We use principal component analysis (PCA) [16] based on a set of given scenarios to define the uncertainty set as a rotated hyperrectangle as in [17]. While the solution might be conservative, it is guaranteed to yield feasible system states for this continuous set.

Due to the computational expense of exact methods to determine the optimal placement of FACTS, various studies have suggested heuristics as a solution to this problem [18,19,20]. In this study, we also propose a greedy algorithm that yields near-optimal solutions while being much faster to solve than the MILP. It iteratively checks for improvement in the objective function when adding a PST to a transmission line in a hill-climbing fashion. In each iteration, the algorithm solves the linear programming relaxation of the proposed MILP and projects the fractional solution onto the feasible set of the original problem. While the maximum number of iterations grows linearly with the number of transmission lines in the grid, a parametrized stopping condition can drastically reduce the algorithm’s search space. We demonstrate that for a realistically sized transmission grid, the greedy algorithm is more than 40 times faster to compute than the MILP and finds a solution with an objective value that is 17% greater. Despite a relatively big optimality gap for a planning problem, the solution of the greedy algorithm is still useful, e.g., in preliminary grid analyses.

Two examples are evaluated numerically. The first example with three buses stylizes a typical grid situation in Germany, where centralized wind power leads to loop flows through neighboring countries. The computed affine linear policies are compared against optimizing the redispatch separately for each scenario. We can show that considering only extreme scenarios with no or maximal wind generation are not sufficient to guarantee grid feasibility in all situations. The second example uses the IEEE 39 bus test case and real-world time series.

The remainder of the paper is structured as follows. In Section 2, we define the linear power flow model employed throughout this work. We give a formal problem definition in Section 3. Section 4 presents our proposed MILP formulation assuming that the uncertainty set is a (convex) polytope. In Section 5, we show how to construct the polytopic uncertainty set given a collection of renewable and load scenarios and their corresponding dispatch results. The numeric evaluation is described in Section 7. Conclusions are drawn in Section 8.

2. Power Flow Model

We use the common DC approximation [21] to linearly model the power flow in a transmission grid with N buses and L lines. Power line flows are denoted as ${\mathbf{p}}^F\in {\mathbb{R}}^L$, voltage phase angles as $\mathsf{\theta }\in {\mathbb{R}}^N$, and the phase shifts potentially added by the PSTs as $\mathsf{\varphi }\in {\mathbb{R}}^L$. We then have

$$\begin{array}{c}\hfill {\mathbf{p}}^F(\mathsf{\theta },\mathsf{\varphi })=\mathbf{B}\left(\mathbf{A}\mathsf{\theta }+\mathsf{\varphi }\right),\end{array}$$

(1)

where $\mathbf{A}\in {\mathbb{R}}^{L\times N}$ is the grid’s incidence matrix and $\mathbf{B}\in {\mathbb{R}}^{L\times L}$ is a diagonal matrix with the line susceptances.

Let $\mathbf{x}\in \mathcal{X}\subset {\mathbb{R}}^P$ be the uncertain vector of power set points of the P generators/loads of the grid, where $\mathcal{X}$ is the corresponding uncertainty set. Moreover, let $\mathbf{y}\in {\mathbb{R}}^P$ denote the vector of power adjustments caused by redispatch actions. The power injections at the buses of the grid ${\mathbf{p}}^I(\mathbf{x},\mathbf{y})\in {\mathbb{R}}^N$ are then given by

$$\begin{array}{c}\hfill {\mathbf{p}}^I(\mathbf{x},\mathbf{y})=\mathbf{C}(\mathbf{x}+\mathbf{y}),\end{array}$$

(2)

where $\mathbf{C}\in {\mathbb{R}}^{N\times P}$ is the matrix that maps generators/loads to the buses, at which they are connected to the grid. According to Kirchhoff’s first law, the nodal power injections match the sum of the power flows on the connected lines, i.e., ${\mathbf{p}}^I(\theta ,\mathsf{\varphi })={\mathbf{A}}^T{\mathbf{p}}^F(\theta ,\mathsf{\varphi })$. By substituting (1) into the previous relationship and equalizing it with (2), it is possible to express the voltage phase angles of the buses $\mathsf{\theta }$ as

$$\begin{array}{c}\hfill \mathsf{\theta }(\mathbf{x},\mathbf{y},\mathsf{\varphi })={\left({\mathbf{A}}^T\mathbf{B}\mathbf{A}\right)}^{+}\left(\mathbf{C}\left(\mathbf{x}+\mathbf{y}\right)-{\mathbf{A}}^T\mathbf{B}\mathsf{\varphi }\right).\end{array}$$

(3)

Since a constant shift of the phase angles does not change the physical situation, ${\mathbf{A}}^T\mathbf{B}\mathbf{A}$ is not invertible, and we use the pseudo-inverse to obtain a minimum norm solution. Substituting (3) into (1), the power line flows ${\mathbf{p}}^F$ can be written as a linear function of the power set points $\mathbf{x}$, the redispatch $\mathbf{y}$ and the phase shifts added by the PSTs $\mathsf{\varphi }$ as

$$\begin{array}{c}\hfill {\mathbf{p}}^F(\mathbf{x},\mathbf{y},\mathsf{\varphi })=\underset{\tilde{\mathbf{C}}}{\underbrace{\mathbf{E}\mathbf{C}}}(\mathbf{x}+\mathbf{y})+\underset{\tilde{\mathbf{B}}}{\underbrace{\left(\mathbb{I}-\mathbf{E}{\mathbf{A}}^T\right)\mathbf{B}}}\mathsf{\varphi },\end{array}$$

(4)

with $\mathbf{E}=\mathbf{B}\mathbf{A}{\left({\mathbf{A}}^T\mathbf{B}\mathbf{A}\right)}^{+}$. Throughout this paper, $\mathbb{I}$ denotes the identity matrix with appropriate dimensions. This model allows us to express line capacity constraints as a linear function of the modeled uncertainties and the control decisions.

3. Formal Problem Statement

The optimization of the number and placement of PSTs needs to consider the control actions to be performed with them as well as the effect of redispatch measures that are important additional operational measures available to the grid operators. In this paper, we assume that the control policies for setting the phase shifts $\mathsf{\varphi }$ and the power adjustments $\mathbf{y}$ are affine linear functions of the system state, i.e., they depend linearly on the uncertain power set points $\mathbf{x}$. We, thus, can write the policies as

$$\mathsf{\varphi }(\mathbf{x})=\mathbf{S}\mathbf{x}+\mathbf{w},$$

(5)

$$\mathbf{y}(\mathbf{x})=\mathbf{T}\mathbf{x}+\mathbf{q},$$

(6)

where $\mathsf{\varphi }(\mathbf{x})$ is parametrized by $\mathbf{S}\in {\mathbb{R}}^{L\times P}$ and $\mathbf{w}\in {\mathbb{R}}^L$, and $\mathbf{y}(\mathbf{x})$ by $\mathbf{T}\in {\mathbb{R}}^{P\times P}$ and $\mathbf{q}\in {\mathbb{R}}^P$.

We now can formally define the optimization task we aim to solve in this work. We aim at determining the minimal set of PSTs and their location along with the control policies (5) and (6) such that the total worst-case redispatch costs and the PST installations costs are minimized while ensuring that the grid operates within its feasible region for any realization of the uncertain elements $\mathbf{x}$. With $\mathbf{u}\in {\left\{0,1\right\}}^L$ indicating the placement of a PST at a transmission line and $\gamma \in \mathbb{R}$ being the worst-case redispatch cost, the optimization problem can then be written as

$$\begin{array}{cc}\hfill \underset{\mathbf{u},\gamma ,\mathbf{S},\mathbf{w},\mathbf{T},\mathbf{q}}{min}& \mu {\mathbf{1}}^T\mathbf{u}+\gamma \hfill \\ \hfill \mathrm{s}.\mathrm{t}.& |{\mathbf{p}}^F(\mathbf{x},\mathbf{y},\mathsf{\varphi })|\le \mathbf{c},\forall \mathbf{x}\in \mathcal{X},\hfill \\ \hfill & \left|\mathsf{\varphi }(\mathbf{x})\right|\le \overline{\varphi }\mathbf{u},\forall \mathbf{x}\in \mathcal{X},\hfill \\ \hfill & {\mathbf{1}}^T\mathbf{y}(\mathbf{x})=0,\forall \mathbf{x}\in \mathcal{X},\hfill \\ \hfill & \underline{\mathbf{h}}(\mathbf{x})\le \mathbf{y}(\mathbf{x})\le \overline{\mathbf{h}}(\mathbf{x}),\forall \mathbf{x}\in \mathcal{X},\hfill \\ \hfill & {\mathsf{\kappa }}^T\mathbf{y}(\mathbf{x})\le \gamma ,\forall \mathbf{x}\in \mathcal{X},\hfill \\ \hfill & (4),(5),(6).\hfill \end{array}$$

(7)

Here, $\mathbf{1}$ is a vector of ones with appropriate dimensions, and $\mu \in {\mathbb{R}}_{\ge 0}$ is a weighting factor between the cost for placing the PSTs in the grid and the worst-case redispatch cost. It should be chosen depending on the assumed probability of the worst-case redispatch situation. Vector $\mathbf{c}\in {\mathbb{R}}_{\ge 0}^L$ denotes the maximum transport capacity of the lines, and $\overline{\varphi }\in {\mathbb{R}}_{\ge 0}$ is the PST technical constraint on the maximum allowed phase shift. Redispatch power adjustments $\mathbf{y}$ sum to zero and should not violate the situation-dependent physical limitations $\overline{\mathbf{h}}(\mathbf{x})$, $\underline{\mathbf{h}}(\mathbf{x})$ of generation or consumption for each element of $\mathbf{x}$. An example of such bounds is $\underline{h}(\mathbf{x})=-\mathbf{x}$ and $\overline{\mathbf{h}}(\mathbf{x})=0$ for wind and solar power plants. We generally assume a linear relation, i.e., $\overline{\mathbf{h}}(\mathbf{x})=\overline{\mathbf{h}}+\overline{\mathbf{H}}\mathbf{x}$ and $\underline{\mathbf{h}}(\mathbf{x})=\underline{\mathbf{h}}+\underline{\mathbf{H}}\mathbf{x}$, where $\overline{\mathbf{H}}$ and $\underline{\mathbf{H}}$ are diagonal matrices. Moreover, $\mathsf{\kappa }\in {\mathbb{R}}_{\ge 0}^P$ represents the cost of increasing the power output of a generator/load by one unit.

4. Robust Optimization with Polytopic Uncertainty

In this section, we formulate a MILP model to solve the optimization task defined in the last section. By assuming a polytopic uncertainty set $\mathcal{X}$, we are able to express the infinite number of robustness constraints as a finite set of linear constraints using duality properties.

To this end, we first use expression (4) and the parametrizations (5) and (6) to rewrite the constraints of (7) in compact matricial form, i.e.,

$$\begin{array}{c}\hfill \tilde{\mathbf{A}}(\mathbf{S},\mathbf{T})\mathbf{x}\le \tilde{\mathbf{c}}(\mathbf{u},\mathbf{w},\mathbf{q},\gamma ),\forall \mathbf{x}\in \mathcal{X},\end{array}$$

(8)

with $\tilde{\mathbf{A}}(\mathbf{S},\mathbf{T})$ and $\tilde{\mathbf{c}}(\mathbf{u},\mathbf{w},\mathbf{q},\gamma )$ being defined in (9) and (10), respectively, where · represents zero entries:

$$\begin{array}{c}\hfill \tilde{\mathbf{A}}(\mathbf{S},\mathbf{T})=\left[\left[\begin{array}{r}\tilde{\mathbf{C}}\\ -\tilde{\mathbf{C}}\\ ·\\ ·\\ ·\\ ·\\ -\overline{\mathbf{H}}\\ \underline{\mathbf{H}}\\ ·\end{array}\right]+\left[\begin{array}{r}\tilde{\mathbf{B}}\\ -\tilde{\mathbf{B}}\\ \mathbb{I}\\ -\mathbb{I}\\ ·\\ ·\\ ·\\ ·\\ ·\end{array}\right]\mathbf{S}+\left[\begin{array}{r}\tilde{\mathbf{C}}\\ -\tilde{\mathbf{C}}\\ ·\\ ·\\ {\mathbf{1}}^T\\ -{\mathbf{1}}^T\\ \mathbb{I}\\ -\mathbb{I}\\ {\mathsf{\kappa }}^T\end{array}\right]\mathbf{T}\right],\end{array}$$

(9)

$$\begin{array}{c}\hfill \tilde{\mathbf{c}}(\mathbf{u},\mathbf{w},\mathbf{q},\gamma )=\left[\begin{array}{c}\mathbf{c}-\tilde{\mathbf{B}}\mathbf{w}-\tilde{\mathbf{C}}\mathbf{q}\\ \mathbf{c}+\tilde{\mathbf{B}}\mathbf{w}+\tilde{\mathbf{C}}\mathbf{q}\\ \overline{\varphi }\mathbf{u}-\mathbf{w}\\ \overline{\varphi }\mathbf{u}+\mathbf{w}\\ -{\mathbf{1}}^T\mathbf{q}\\ {\mathbf{1}}^T\mathbf{q}\\ \overline{\mathbf{h}}-\mathbf{q}\\ -\underline{\mathbf{h}}+\mathbf{q}\\ \gamma -{\mathsf{\kappa }}^T\mathbf{q}\end{array}\right].\end{array}$$

(10)

Next, we assume that the uncertainty set $\mathcal{X}$ of the power set points $\mathbf{x}$ is a polytope defined as the intersection of M halfspaces. That is, given the exterior representation of $\mathcal{X}$

$$\begin{array}{c}\hfill \mathcal{X}=\left\{\mathbf{x}\in {\mathbb{R}}^P|\mathbf{D}\mathbf{x}\le \mathbf{b}\right\},\end{array}$$

(11)

with $\mathbf{D}\in {\mathbb{R}}^{M\times P}$ and $\mathbf{b}\in {\mathbb{R}}^M$, we assume that $\exists \upsilon \in {\mathbb{R}}_{>0}$ such that $\Vert \mathbf{x}\Vert \le \upsilon$, $\forall \mathbf{x}\in \mathcal{X}$. In our experiments, the uncertainty set is guaranteed to be a polytope by construction as shown in Section 5.

Using (8) and (11), our optimization task is equivalent to the following min-max optimization problem:

$$\begin{array}{cc}\hfill \eta =& \underset{\mathbf{u},\gamma ,\mathbf{S},\mathbf{w},\mathbf{T},\mathbf{q}}{min}\mu {\mathbf{1}}^T\mathbf{u}+\gamma \hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}\begin{array}{cc}\hfill \underset{\mathbf{x}}{max}& {\tilde{\mathbf{A}}}_j(\mathbf{S},\mathbf{T})\mathbf{x}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{D}\mathbf{x}\le \mathbf{b}\hfill \end{array}\end{array}\right\}\le {\tilde{\mathbf{c}}}_j(\mathbf{u},\mathbf{w},\mathbf{q},\gamma ),\hfill \\ \hfill & \forall j\in \left[K\right],\hfill \end{array}$$

(12)

where K is the number of rows of ${\tilde{\mathbf{A}}}_j(\mathbf{S},\mathbf{T})$, $\left[K\right]$ is the set $\left\{1,\cdots ,K\right\}$, ${\tilde{\mathbf{A}}}_j(\mathbf{S},\mathbf{T})$ is the j-th row of $\tilde{\mathbf{A}}(\mathbf{S},\mathbf{T})$, and ${\tilde{\mathbf{c}}}_j(\mathbf{u},\mathbf{w},\mathbf{q},\gamma )$ is the j-th element of $\tilde{\mathbf{c}}(\mathbf{u},\mathbf{w},\mathbf{q},\gamma )$.

We transform the min-max problem (12) into a single-stage linear program using duality theory. As $\mathcal{X}$ is nonempty and bounded, the inner optimization problem is always feasible, and hence its primal and dual problems have the same objective value by strong duality. We can thus rewrite (12) as

$$\begin{array}{cc}\hfill \eta =& \underset{\mathbf{u},\gamma ,\mathbf{S},\mathbf{w},\mathbf{T},\mathbf{q}}{min}\mu {\mathbf{1}}^T\mathbf{u}+\gamma \hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}\begin{array}{cc}\hfill \underset{{\mathsf{\lambda }}_j}{min}& {\mathbf{b}}^T{\mathsf{\lambda }}_j\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {\mathbf{D}}^T{\mathsf{\lambda }}_j={\tilde{\mathbf{A}}}_j{(\mathbf{S},\mathbf{T})}^T\hfill \\ \hfill & {\mathsf{\lambda }}_j\ge 0\hfill \end{array}\end{array}\right\}\le {\tilde{\mathbf{c}}}_j(\mathbf{u},\mathbf{w},\mathbf{q},\gamma ),\hfill \\ \hfill & \forall j\in \left[K\right]\hfill \end{array}$$

(13)

where ${\mathsf{\lambda }}_j\in {\mathbb{R}}^M$ are the dual variables of the j-th inner problem. Optimization problem (13) can then be formulated as a single-level minimization problem:

$$\begin{array}{cc}\hfill \eta =\underset{\mathbf{u},\gamma ,\mathbf{S},\mathbf{w},\mathbf{T},\mathbf{q},\mathsf{\lambda }}{min}& \mu {\mathbf{1}}^T\mathbf{u}+\gamma \hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {\mathbf{b}}^T{\mathsf{\lambda }}_j\le {\tilde{\mathbf{c}}}_j(\mathbf{u},\mathbf{w},\mathbf{q},\gamma ),\forall j\in \left[K\right],\hfill \\ \hfill & {\mathbf{D}}^T{\mathsf{\lambda }}_j={\tilde{\mathbf{A}}}_j{(\mathbf{S},\mathbf{T})}^T,\forall j\in \left[K\right],\hfill \\ \hfill & {\mathsf{\lambda }}_j\ge 0,\forall j\in \left[K\right].\hfill \end{array}$$

(14)

The equivalence between (13) and (14) can be verified by noticing that an optimal solution for (14) is also a feasible solution for (13) and that both objective functions have the same value. On the other hand, an optimal solution for the outer problem of (13) implies that there exists a $\mathsf{\lambda }$ that satisfies the constraints of the inner problem, and, thus, the solution would also be feasible in (14), and both objective functions would have the same value.

5. Robust Uncertainty Set

The optimization formulated above enables us to solve the robust PST placement problem as a MILP model given a polytopic uncertainty set $\mathcal{X}$. As sketched in Figure 1, the influencing factors to be described via the uncertainty set are not only the renewable and load uncertainties but also the market and frequency control measures. These potentially large effects are determined via complex mechanisms, often showing highly non-linear behavior with respec to the loads and renewables. For example, the classic merit-order algorithm uses one power plant after the other in full, but it does not proportionally scale production levels. By including these factors in the uncertainty set, the linearity assumption for the remaining minor deviations from redispatch and active flow control can be justified.

To create realistic uncertainty sets, we propose to generate samples by creating a set of plausible renewable generation and load scenarios either from measurements or probabilistic models. For each scenario a market model, e.g., a simple merit-order dispatch, is used to compute the dispatch decisions. The full information is then used to generate a continuous uncertainty set $\mathcal{X}$ that covers the modeled scenarios but also the space in between, i.e., all (convex) combinations of them.

Given S scenarios of the uncertain vector of power set points ${\mathbf{x}}^{(k)}$, $k=\left\{1,\cdots ,S\right\}$, a first set enclosing all these points can be defined as

$$\begin{array}{c}\hfill \tilde{\mathcal{X}}={\times }_{i=1}^P\left[{\widetilde{\underline{\mathbf{x}}}}_i,{\widetilde{\overline{\mathbf{x}}}}_i\right],\end{array}$$

(15)

where ${\widetilde{\underline{\mathbf{x}}}}_i$ and ${\widetilde{\overline{\mathbf{x}}}}_i$ are defined as ${\widetilde{\underline{\mathbf{x}}}}_i=min({\mathbf{x}}_i^{(1)},\cdots ,{\mathbf{x}}_i^{(S)})$ and ${\widetilde{\overline{\mathbf{x}}}}_i=max({\mathbf{x}}_i^{(1)},\cdots ,{\mathbf{x}}_i^{(S)})$ for all dimensions $i=\left\{1,\cdots ,P\right\}$.

However, as shown in Figure 2, this axis-aligned hyperrectangle can be large and lead to overly conservative solutions in case the components of $\mathbf{x}$ are highly correlated, which is usually the case in power systems (consider, for instance, that the wind or photovoltaic production in close by locations will typically be similar).

To leverage these correlations between the different dimensions of the scenarios, we use PCA to compute a tighter, rotated uncertainty set as in [17].

Let $\mathbf{X}={\left[\begin{array}{ccc}{\mathbf{x}}^{(1)}& \cdots & {\mathbf{x}}^{(S)}\end{array}\right]}^T\in {\mathbb{R}}^{S\times P}$ be the data matrix. We first center the data by subtracting them by their empirical mean ${\widehat{\mathbf{z}}}_j=\frac{1}{S}{\sum }_{i=1}^S{\mathbf{X}}_{ij}$ to obtain the matrix $\mathbf{Z}=\mathbf{X}-\widehat{\mathbf{z}}{\mathbf{1}}^T$. Next, we calculate the empirical covariance matrix $\mathsf{\Sigma }\in {\mathbb{R}}^{P\times P}$ as $\mathsf{\Sigma }=\frac{1}{S-1}{\mathbf{Z}}^T\mathbf{Z}$. Symmetric matrix $\mathsf{\Sigma }$ can be diagonalized, i.e.,  $\mathsf{\Sigma }=\mathbf{Q}\mathsf{\Lambda }{\mathbf{Q}}^{-1}$, where $\mathbf{Q}\in {\mathbb{R}}^{P\times P}$ has the eigenvectors of $\mathsf{\Sigma }$ in its columns and $\mathsf{\Lambda }\in {\mathbb{R}}^{P\times P}$ is a diagonal matrix with the corresponding eigenvalues on its diagonal. This allows us to construct a rotated, data-aligned hyperrectangle $\widehat{\mathcal{X}}$ as

$$\begin{array}{c}\hfill \widehat{\mathcal{X}}=\left\{\widehat{\mathbf{x}}\in {\mathbb{R}}^P|\mathbf{D}\widehat{\mathbf{x}}\le \mathbf{b}\right\},\end{array}$$

(16)

where $\mathbf{D}={\left[\begin{array}{cc}\mathbf{Q}& -\mathbf{Q}\end{array}\right]}^T$ and $\mathbf{b}={\left[\begin{array}{cc}\widehat{\overline{\mathbf{x}}}& -\widehat{\underline{\mathbf{x}}}\end{array}\right]}^T$, with ${\widehat{\overline{\mathbf{x}}}}_j={max}_i{\left(\mathbf{X}\mathbf{Q}\right)}_{ij}$ and ${\widehat{\underline{\mathbf{x}}}}_j={min}_i{\left(\mathbf{X}\mathbf{Q}\right)}_{ij}$ being the j-th elements of the row vectors $\widehat{\overline{\mathbf{x}}}$ and $\widehat{\underline{\mathbf{x}}}$ containing the maximum and minimum values among the S scenarios in the direction of the principal components, respectively.

With this approach, we are able to fit a possibly much tighter set to the scenarios by defining the uncertainty set $\mathcal{X}$ as the intersection between $\tilde{\mathcal{X}}$ and $\widehat{\mathcal{X}}$, i.e.,

$$\begin{array}{c}\hfill \mathcal{X}=\tilde{\mathcal{X}}\cap \widehat{\mathcal{X}},\end{array}$$

(17)

which is shown in Figure 2.

6. Greedy Algorithm for Solving Robust MILP

As the size of the grid increases, the runtime for solving (14) to optimality becomes too long as shown in our experiments (see Section 7). We, thus, propose a hill-climbing-like heuristic that iteratively checks whether adding a PST to a line decreases the value of the objective function. Each iteration relies on solving the linear programming relaxation of (14) and projecting the fractional solution onto the feasibility set of the original MILP. The algorithm performs, at most, L iterations, thus yielding polynomial time complexity. However, the solution of the proposed method is only guaranteed to be a local minimum on the solution space of the original MILP by nature of the hill-climbing procedure.

The algorithm first determines a lower bound for the objective function of (14) by solving its linear programming relaxation, as any solution of the (mixed-)integer program is a feasible solution of its linear relaxation. If the optimal solution of the relaxed problem has all values of $\mathbf{u}$ as 0 or 1, it will also be the optimal solution to the MILP. To obtain an upper bound, we can solve (14) with $\mathbf{u}$ fixed to the value of the solution to the relaxed problem rounded to the closest integer.

The algorithm proceeds in a hill-climbing fashion to improve the integer solution. With ${\mathbf{u}}^{\ast }$ being the value of $\mathbf{u}$ from the solution to the relaxed problem, the algorithm places a PST into the line corresponding to the largest element of ${\mathbf{u}}^{\ast }$, i.e., it fixes element $j=arg{max}_i{\mathbf{u}}_i^{\ast }$ of $\mathbf{u}$ to 1. If the line already has a PST, the algorithm adds a PST to the line corresponding to the next largest element of ${\mathbf{u}}^{\ast }$. With one of the elements of $\mathbf{u}$ fixed, the algorithm solves the relaxed problem and uses the rounded solution to update the upper bound of the objective function. The algorithm continues adding PSTs to the grid until the value of the objective function stops decreasing between iterations, or if the largest element of ${\mathbf{u}}^{\ast }$ is smaller than some parameter $ϵ\in \left[0,1\right]$.

The pseudocode for the proposed greedy algorithm is presented in Algorithm 1, where ⌀ is the empty set, $\mathcal{L}$ is the set of lines with PST, $\lfloor ·\rceil$ is the round operator, and $|·|$ represents the cardinality of a set. While the Main procedure implements the hill-climbing algorithm, the solve procedure solves the linear relaxation of the original MILP (14) with the variables corresponding to the lines that are already equipped with PST being fixed to one.

Algorithm 1 Greedy algorithm for solving robust MILP.

1: procedure Main($ϵ$) 2:     ${\underline{\eta }}^{\ast }$, ${\mathbf{u}}^{\ast }\leftarrow$ Solve(·, ⌀) 3:     ${\overline{\eta }}^{\ast }$, $·\leftarrow$ Solve($\lfloor {\mathbf{u}}^{\ast }\rceil$, $\left[L\right]$) 4:     $\mathcal{L}\leftarrow ⌀$ 5:     while ${max}_{i\notin \mathcal{L}}{\mathbf{u}}_i^{\ast }>ϵ$ and ${\mathbf{u}}^{\ast }\notin {\left\{0,1\right\}}^L$ do 6:         $\mathcal{L}\leftarrow \mathcal{L}\cup \left\{arg{max}_{i\notin \mathcal{L}}{\mathbf{u}}_i^{\ast }\right\}$ 7:         ·, ${\widehat{\mathbf{u}}}^{\ast }\leftarrow$ Solve(${\left\{1\right\}}^{\left|\mathcal{L}\right|}$, $\mathcal{L}$) 8:         ${\widehat{\overline{\eta }}}^{\ast }$, $·\leftarrow$ Solve($\lfloor {\widehat{\mathbf{u}}}^{\ast }\rceil$, $\left[L\right]$) 9:         if ${\widehat{\overline{\eta }}}^{\ast }\ge {\overline{\eta }}^{\ast }$ then 10:            return $\lfloor {\mathbf{u}}^{\ast }\rceil$ 11:         end if 12:         ${\overline{\eta }}^{\ast }\leftarrow {\widehat{\overline{\eta }}}^{\ast }$, ${\mathbf{u}}^{\ast }\leftarrow {\widehat{\mathbf{u}}}^{\ast }$ 13:     end while 14:     return $\lfloor {\mathbf{u}}^{\ast }\rceil$ 15: end procedure 16:   17: procedure Solve(${\mathbf{u}}^{\ast }$, $\mathcal{Q}$) 18:     solve relaxation of MILP (14) with ${\mathbf{u}}_{i\in \mathcal{Q}}$ fixed to ${\mathbf{u}}^{\ast }$ 19:     return value of $\eta$, value of $\mathbf{u}$ 20: end procedure

7. Experiments

We demonstrate our proposed approach for two examples, a small 3-bus test grid and the IEEE 39 bus system. Our implementation uses the modeling language JuMP [22] and Gurobi [23] to solve the MILP model (14).

7.1. 3-Bus Test Grid

We first apply the proposed algorithm to the grid in Figure 3. The figure also contains the values of the loads, generator limits and unit costs, as well as line capacities and susceptances. The unit cost of wind production has a negative value, mirroring precedence for renewable power generation over conventional one, as well as standard subsidy policies.

Power is represented in the per unit system (p.u.) for some base power value and costs in a generic currency ¤. The example is a stylized sketch of the situation in Germany, with much wind and lignite generation in the East and the load center and many gas power plants in the Southwest [24]. Bus 3 would then represent neighboring countries.

In Figure 4, we examine the operation patterns of the test grid when varying the available wind power from zero to its maximum capacity. Both wind sites are scaled proportionally. For each wind in-feed, we determine the cost optimal dispatch using the merit order algorithm. Demands that cannot be covered from wind alone are supplied by the lignite power plant first and then by the gas power plant. For each scenario, we then obtain a cost-optimal redispatch by solving (7).

As can be seen in Figure 4, no redispatch is needed when the wind is minimal or maximal: the merit-order dispatch is directly feasible. With no wind, the generation pattern is 2 p.u. at bus 1 and 2 p.u. at bus 2. With maximum wind, the generation is 3 p.u. at bus 2 and 1 p.u. at bus 3. A grid analysis based only on these two scenarios could conclude that the grid can be feasibly operated for all possible wind outputs without additional measures, such as PST placement or redispatch. However, this would be the wrong conclusion, since for 2 p.u. of wind, the generation pattern is 3.5 p.u. at bus 2 and 0.5 p.u. at bus 3, which leads to a capacity violation on line (2,3). This example thus proves that relying on few scenarios, as is often done for two-stage stochastic programming approaches due to computational limitations, is not enough to ensure universal grid feasibility. In contrast, our robust approach considers a continuous uncertainty with all (convex) scenario mixtures included. While for this low-dimensional example, one can obviously include average scenarios into a scenario-based analysis, this would be much more difficult in more complex, high-dimensional setups, where the number of required scenarios could be huge. Note that our approach scales only in the dimension of the uncertainties but not in the number of scenarios, which can grow exponentially with the number of uncertainty dimensions if an even cover of the scenario space is desired.

In Figure 4a,b, we also show the redispatch policies computed by solving (14). For Figure 4a, we choose $\mathsf{\mu }$ such that no PST is installed, but higher redispatch costs are preferable, whereas in Figure 4b, the parameter $\mu$ is chosen to be small enough such that one PST on line (1,3) is optimal. Note that the derived values of the computed redispatch policy are not linear in the wind power production. However, they are affine linear in the dispatch, including the non-linear merit order results. Even for this low-dimensional example, the redispatch policies qualitatively resemble the optimal redispatch values computed for each scenario individually. For higher-dimensional examples, the linearity condition might even be less restrictive. For both cases, with and without PST, the implied worst-case cost by the computed affine linear policy matches the worst-case cost of the scenario-wise optimized redispatch schedule as shown in Figure 4c.

The example also shows that the redispatch cost can be reduced by increasing the controllability of the power flow in the grid with the addition of a PST. Figure 4d shows the angle shift added by the PST following the affine control policy as well as the amount of power flow redirected from branch (2,3) to (2,1). These settings ensure a feasible grid power flow in all situations. It is worth noting that in the “no wind” situation, the PST also redirects energy from line (2,3) to line (1,2), which helps increase the distance to the feasibility border since without PST interaction the line (2,3) would be fully loaded.

7.2. IEEE 39 Bus Test Case

We now examine a larger example, namely the IEEE 39 test case, which contains 39 buses, 46 transmission lines, 10 generators, and 21 loads as shown in Figure 5. The properties of the transmission lines, generators, and loads are based on the MATPOWER test cases [25]. To generate diverse usage scenarios, we proceed as follows. We start with an hourly load profile curve aggregated for a single country. Here, we use the Danish load from 28 September 2020 and 27 September 2021 to generate 8735 scenarios [26]. To obtain nodal load time series that are both realistically correlated but are also partly independent of each other, we normalize the given load curve by its maximum value to yield $\ell (t)$, $t=\{1,\cdots ,8735\}$, and add randomness as follows:

$$\begin{array}{c}\hfill d_i(t)=d_i^N\left(\ell (t)+u_i(t)\right),\end{array}$$

(18)

where $d_i^N$ is the nominal value given in the IEEE test case for load i, and $u_i(t)$ is a sample of the uniform distribution ${\mathcal{U}}_i(-0.2,0.2)$. Merit order is used for the dispatch of the generators where the dispatch cost of a generator is taken proportional to the inverse of its maximum active power output.

In Figure 6a, we show the dispatch of the generators for the first 200 consecutive time steps. The system without redispatch or PSTs is infeasible for the generated scenarios in the sense that the power flow over some transmission lines is larger than their capacity. We show in Figure 5 the transmission lines with the largest violations in red.

In Figure 6b, we show the power output of the generators when following the computed redispatch policy without any PST. This is achieved for very high PST installation costs in the MILP (14) problem. The optimal redispatch policy makes the power flow feasible for all scenarios, and the worst-case redispatch cost is $3.30$ ¤. The cheap generator G1 produces less power than is economically optimal without considering grid constraints.

For Figure 6c, we solve the MILP (14), now lowering the PST installation costs significantly. With these settings, the algorithm proposes to install 1 PST between buses 1 and 39 as shown in Figure 5. Given this decision, it is possible to reduce the worst-case redispatch cost to $1.93$ ¤ by using more power from the cheaper generator G1.

Optimality versus computation time is examined in

Table 1. Performance evaluation of MILP and greedy algorithm $\mathcal{G}(ϵ)$.

	Cost [¤]	Time [$\mathbf{min}$]	PSTs
MILP	1.93	348	1
$\mathcal{G}$($0.2$)	3.30	3	0
$\mathcal{G}$($0.06$)	2.25	8	1

, comparing the MILP and the greedy algorithm proposed in Section 6 for the same setting. The greedy algorithm is run with two different values of $ϵ$. With $ϵ$ equal to $0.2$, the algorithm does not propose the installation of any PST. If the value of $ϵ$ is decreased to $0.06$, the algorithm finds a better solution that yields $2.25$ ¤ for the worst-case redispatch cost, an increase of $17\%$ compared to the optimal solution of the MILP. In this solution, a PST is installed between buses 2 and 3.

Using a local compute node with 22 cores, MILP optimality with a 10% gap was proven after 5 h and 48 min. While acceptable for grid-planning purposes, the computation time is relatively high. In the configuration with $ϵ$ equal to $0.2$, the computation time of the greedy algorithm is 116 times faster than solving the MILP, and it is $43.5$ times faster for $ϵ$ equal to $0.06$. One could further improve the computational time of both the MILP and the proposed algorithm by, for example, restricting the affine linear maps to depend on principal components of the in-feeds only, thereby reducing the dimensions of $\mathbf{S}$ and $\mathbf{T}$. Nevertheless, Table 1 shows that the proposed greedy algorithm outputs near-optimal objective values while being much faster to solve than the original MILP. The results show that the hill-climbing approach could be a valuable tool in preliminary analyses of PST placement, where sub-optimal solutions are acceptable. Then, once the final model is defined, the MILP can be run to find the optimal placement solution.

8. Conclusions

In this work, we presented an algorithm to find the minimal set of PSTs and the optimal phase shift and redispatch policies that enable the feasible operation of the transmission system for a continuous set of uncertain in-feed scenarios. We demonstrated our algorithm with two examples, proving that only considering a few extreme scenarios may not be sufficient to guarantee feasible grid states under all conditions and that the PST placement can significantly reduce the worst-case redispatch cost. The experiments also show that a proposed greedy algorithm can compute solutions much faster than the MILP, making it suitable for preliminary analyses on large grids. The proposed methodology can be easily extended to solve the optimal placement problem of different FACTS devices as long as a linear relationship model exists between their actuation and the grid power flow.