New Distributed Optimization Method for TSO–DSO Coordinated Grid Operation Preserving Power System Operator Sovereignty

Abstract

Electrical power system operators (SOs) are free to realize grid operations according to their own strategies. However, because resulting power flows also depend on the actions of neighboring SOs, appropriate coordination is needed to improve the resulting system states from an overall perspective and from an individual SO perspective. In this paper, a new method is presented that preserves the data integrity of the SOs and their independent operation of their grids. This method is compared with a non-coordinated local control and another sequential method that has been identified as the most promising distributed optimization method in previous research. The time series simulations use transformer tap positioning as well as generation unit voltage setpoints and reactive power injections as flexibilities. The methods are tested on a multi-voltage, multi-SO, realistic benchmark grid with different objective combinations of the SOs. In conclusion, the results of the new method are much closer to the theoretical optimum represented by central optimization than those of the other two methods. Furthermore, the introduced method integrates a sophisticated procedure to provide fairness between SOs that is missing in other methods.

1. Introduction

Except for island grids, power systems are operated by multiple system operators (SOs) that may follow independent goals and strategies. The electrical system, however, knows no borders between SOs, only physical laws. Therefore, improved coordination between SOs offers the potential to improve the efficiency and security in the overall grid operation across all SOs. As shown in principle in Figure 1a, this applies both between transmission system operators (TSOs) as well as between TSOs and distribution system operators (DSOs). The latter has special relevance due to the energy system transition because the generation capacity and corresponding ancillary service of wind and solar plants are more often connected to DSOs than those of conventional plants.

Maximum coordination of variables at SO boundaries is provided by grid operation approaches that centralize the data and control of the SO. These, however, do not respect the sovereignty and data integrity of SOs, which operate independently and are cautious of sharing data on their critical infrastructure. On the contrary, grid operation methods should limit the data that must be shared and let SOs control the flexibilities of their grids. In addition, near-optimal system states are desired, while the SOs operate with limited knowledge about the overall system (Figure 1b).

Among different types of coordinated decentralized grid operation methods, this paper addresses this challenge by proposing a new method in the group of distributed optimization methods.

1.1. Distributed Optimization Methods Coordinating Multiple System Operators

In the literature, several distributed optimization methods have been proposed to provide a well-suited trade-off between sovereignty-preserving but suboptimally performing non-coordinated methods and a central optimization [1]. Distributed optimization methods apply mathematical optimizations which means optimal power flows (OPFs) calculations in this context. As depicted in Figure 2, three different categories exist: First, different algorithms decompose the central optimization problem to optimize and coordinate the problem parts iteratively. Second, grid equivalents with iterative updates represent the power systems of neighboring SOs within optimizations. Third, methods with fixed procedures without iterations are called sequential.

Several mathematical decomposition methods have been proposed and are still being further developed. The Lagrangian relaxation (LR) technique has been proposed and proven several times with electrical power systems [2,3,4,5,6,7,8,9]. Therein, LR is often applied to optimize voltage control and reactive power dispatch [6,7,8,9]. In contrast, other decomposition techniques include variables at the system operator borders to the optimization problems of both sides. For these copy variables, equality constraints must be added and considered. This second category includes the auxiliary problem principle (APP) [10], the predictor–corrector proximal multiplier method (PCPM) [19], the alternating direction method of multipliers (ADMM) [20,21], and the augmented-Lagrangian-based alternating direction inexact newton method (ALADIN) [13]. Applications to electrical power systems and comparisons are available in [11,12].

Similarly, distributed optimization methods using grid equivalents have applied reactive power dispatch. Simple and obvious grid equivalents at SO borders are $PQ$ loads [1,14,15,17] and slack elements to represent neighboring higher voltage SOs [1,14]. Additionally, $PV$ and $PQ\left(V\right)$ elements, as well as more complex grid equivalents, such as Extended Wards and REI equivalents, have been investigated with distributed optimization methods [1,17,22]. In [16], an equivalent susceptance matrix and an optimal reactive power injection vector were passed for coordination between SOs. Besides grid equivalent types, which can be easily shared via the CIM/CGMES protocol [23], updating procedures of the equivalents and further adaptions may differ in method.

In [18], TSO and DSO exchange operational limits (voltage magnitude or reactive power exchange) and setpoints at the border were introduced. This method, called the DSO–TSO–DSO chain, was analyzed with an adaption for reduced data availability of external grids and TSO–TSO interfaces in [1].

1.2. Contribution of the Paper and Structure

In [1], a review of most of the methods introduced in Section 1.1 is provided. The discussion determines the sequential method of the DSO–TSO–DSO chain as the most promising method for implementation in real grids, especially due to its reliability. However, the low improvement over the non-coordinated grid operation benchmark and an appropriate balance between interests of multiple SOs were identified as issues for future research. To close this gap, this paper introduces a new sequential, distributed optimization method (in green in Figure 2). It can handle different objective functions of the SOs fairly through implementing polynomial equivalent functions while preserving the sovereignty and data integrity of the TSOs and DSOs.

The paper is structured as follows. In Section 2, the objectives, control variables, and constraints of power system grid operation are defined. Section 3 introduces how the performance of the new method is measured and compared to existing methods. In Section 4, the new method is proposed. The time series simulation results are presented and analyzed in Section 5. In Section 6, concluding statements finalize the paper.

2. Grid Operation: Analyzed Constraints, Control Variables and Objectives

The area of grid operation is restricted by power flow (Equations (1) and (2)), generation unit limits (Equations (3) and (4)), voltage magnitude limits (Equation (5)), tap ratio magnitude limits (Equation (6)), and branch loading limits (Equations (7) and (8)). The definitions of limits make use of different sets of elements (buses: $\mathcal{B}$, generation units: $\mathcal{G}$, loads, respectively, demand: $\mathcal{D}$, junctions (lines and transformers $\mathcal{T}$): $\mathcal{J}$).

$$\begin{array}{cccc}& \sum _{i\in {\mathcal{G}}_b}{p}_{g_i}-\sum _{i\in {\mathcal{D}}_b}{p}_{d_i}=\sum _{\begin{array}{c}(i,f,t)\in \mathcal{J},\\ f=b\end{array}}{p}_{(i,f,t)}+\sum _{\begin{array}{c}(i,f,t)\in \mathcal{J},\\ t=b\end{array}}{p}_{(i,t,f)}\hfill & & \forall b\in \mathcal{B}\hfill \end{array}$$

(1)

$$\begin{array}{cccc}& \sum _{i\in {\mathcal{G}}_b}{q}_{g_i}-\sum _{i\in {\mathcal{D}}_b}{q}_{d_i}=\sum _{\begin{array}{c}(i,f,t)\in \mathcal{J},\\ f=b\end{array}}{q}_{(i,f,t)}+\sum _{\begin{array}{c}(i,f,t)\in \mathcal{J},\\ t=b\end{array}}{q}_{(i,t,f)}\hfill & & \forall b\in \mathcal{B}\hfill \end{array}$$

(2)

$$\begin{array}{cccc}& p_g^{\min }\le p_g\le p_g^{\max }\hfill & & \forall g\in \mathcal{G}\hfill \end{array}$$

(3)

$$\begin{array}{cccc}& q_g^{\min }\le q_g\le q_g^{\max }\hfill & & \forall g\in \mathcal{G}\hfill \end{array}$$

(4)

$$\begin{array}{cccc}& v_{m_b}^{\min }\le v_{m_b}\le v_{m_b}^{\max }\hfill & & \forall b\in \mathcal{B}\hfill \end{array}$$

(5)

$$\begin{array}{cccc}& {\tau }_i^{\min }\le {\tau }_i\le {\tau }_i^{\max }\hfill & & \forall (i,f,t)\in \mathcal{T}\hfill \end{array}$$

(6)

$$\begin{array}{cccc}& {\check{s}}_{(i,f,t)}^2\le 1\hfill & & \forall (i,f,t)\in \mathcal{J}\hfill \end{array}$$

(7)

$$\begin{array}{cccc}& {\check{s}}_{(i,t,f)}^2\le 1\hfill & & \forall (i,f,t)\in \mathcal{J}\hfill \end{array}$$

(8)

The power flows (used in Equations (1) and (2)) and loadings (used in Equations (7), (8) and (18)) at bus f into branch i towards bus t and vice versa ($i,t,f$) are defined for balanced three-phase power systems by Equations (9)–(12) and Equations (13)–(16), respectively. Similar to Equation (16), if the nominal voltage $V_{n_{i,f}}$ of the branch i at the side of bus f equals the nominal voltage of that bus $V_{n_f}$. Equation (15) can be simplified.

$$\begin{array}{cc}\hfill p_{(i,f,t)}=& (g_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.·g_{i,\mathrm{ground}})v_{m_f}^2{\tau }_i^2\hfill \\ \hfill & -g_i{v}_{m_f}{v}_{m_t}{\tau }_{i}cos({\theta }_f-{\theta }_t+{\theta }_{i,\mathrm{shift}})\hfill \\ \hfill & -b_i{v}_{m_f}{v}_{m_t}{\tau }_{i}sin({\theta }_f-{\theta }_t+{\theta }_{i,\mathrm{shift}})\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill p_{(i,t,f)}=& (g_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.·g_{i,\mathrm{ground}})v_{m_t}^2\hfill \\ \hfill & -g_i{v}_{m_f}{v}_{m_t}{\tau }_{i}cos({\theta }_f-{\theta }_t+{\theta }_{i,\mathrm{shift}})\hfill \\ \hfill & +b_i{v}_{m_f}{v}_{m_t}{\tau }_{i}sin({\theta }_f-{\theta }_t+{\theta }_{i,\mathrm{shift}})\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill q_{(i,f,t)}=& -(b_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.·b_{i,\mathrm{ground}})v_{m_f}^2{\tau }_i^2\hfill \\ \hfill & -g_i{v}_{m_f}{v}_{m_t}{\tau }_{i}sin({\theta }_f-{\theta }_t+{\theta }_{i,\mathrm{shift}})\hfill \\ \hfill & +b_i{v}_{m_f}{v}_{m_t}{\tau }_{i}cos({\theta }_f-{\theta }_t+{\theta }_{i,\mathrm{shift}})\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill q_{(i,t,f)}=& -(b_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.·b_{i,\mathrm{ground}})v_{m_t}^2\hfill \\ \hfill & +g_i{v}_{m_f}{v}_{m_t}{\tau }_{i}sin({\theta }_f-{\theta }_t+{\theta }_{i,\mathrm{shift}})\hfill \\ \hfill & +b_i{v}_{m_f}{v}_{m_t}{\tau }_{i}cos({\theta }_f-{\theta }_t+{\theta }_{i,\mathrm{shift}})\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\check{s}}_{(i,f,t)}^2=& \frac{p_{(i,f,t)}^2+q_{(i,f,t)}^2}{s_{i,\mathrm{rated}}^2}\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\check{s}}_{(i,t,f)}^2=& \frac{p_{(i,t,f)}^2+q_{(i,t,f)}^2}{s_{i,\mathrm{rated}}^2}\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\check{i}}_{(i,f,t)}=& \frac{{\check{s}}_{(i,f,t)}}{v_{m_f}}\frac{V_{n_{i,f}}}{V_{n_f}}\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill {\check{i}}_{(i,t,f)}=& \frac{{\check{s}}_{(i,t,f)}}{v_{m_t}}\frac{V_{n_{i,t}}}{V_{n_t}}\hfill \end{array}$$

(16)

Various control variables are available for SOs. Tap ratio magnitudes $\tau$ and phase shift angles ${\theta }_{\mathrm{shift}}$—both affected by the tap positions of on-load tap changers (OLTCs)—are important for multi-voltage level use cases. Furthermore, the reactive power provision $q_g$ and voltage magnitude setpoints of generation units $v_m$ are considered in this paper. The active power generation $p_g$, an effective variable to adapt power flows, is not considered as a control variable because it is determined by the economic dispatch [24] and redispatch [25]. The loads ($p_d$ and $q_d$) are assumed as fixed, as well.

As defined in Articles 31 and 40 of Directive (EU) 2019/944, SOs are responsible for operating under economic conditions a secure, reliable, and efficient electricity system. For that purpose, two objective functions $f\left(x\right)$ are derived for minimization, where x includes all independent $(\theta ,v_m,p_g,q_g,\tau ,{\theta }_{\mathrm{shift}})$ and dependent system variables. First, the objective function of Equation (19) achieves a maximum distance to critical grid states caused by voltage limit violations (Equation (17)) or line and transformer overloadings (Equation (18)). Therefore, the objective is appropriate to maximize the security and reliability of grid operation and to reduce corresponding congestion costs. This paper assumes $v_{m_{\mathrm{target}}}=1.03$, $a_1=250$, and $a_2=10$. Secondly, the objective of Equation (20) is to describe the system losses. This objective describes one other economic aspect of grid operation. It is widely used and, in addition to the first objective, well suited to test methods handling SOs with different objective properties.

$$\begin{array}{cc}\hfill & f_{\mathrm{profile}}\left(x\right)=\sum _{b\in \mathcal{B}}{\left(v_{m_b}-v_{m_{\mathrm{target}}}\right)}^2\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & f_{\mathrm{loadings}}\left(x\right)=\sum _{(i,f,t)\in \mathcal{J}}\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left({\check{i}}_{(i,f,t)}^2+{\check{i}}_{(i,t,f)}^2\right)\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & f_{\mathrm{profile},\mathrm{loadings}}\left(x\right)=a_1·f_{\mathrm{profile}}\left(x\right)+a_2·f_{\mathrm{loadings}}\left(x\right)\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & f_{\mathrm{losses}}\left(x\right)=\sum _{(i,f,t)\in \mathcal{J}}\left(p_{(i,f,t)}+p_{(i,t,f)}\right)\hfill \\ \hfill & \phantom{f_{\mathrm{losses}}\left(x\right)}=\sum _{(i,f,t)\in \mathcal{J}}\left((g_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.g_{i,\mathrm{ground}})v_{m_f}^2{\tau }_i^2\right.\hfill \\ \hfill & \phantom{f_{\mathrm{losses}}\left(x\right)\sum _{(i,f,t)\in \mathcal{J}}(}\left.+(g_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.g_{i,\mathrm{ground}})v_{m_t}^2-2g_i{v}_{m_f}{v}_{m_t}{\tau }_{i}cos({\theta }_f-{\theta }_t)\phantom{v_{m_t}^2}\right)\hfill \end{array}$$

(20)

Since transformer tap positions are discrete values, the OPF, resulting from minimizing Equations (19) or (20) subjected to Equations (1)–(16), is a mixed-integer nonlinear programming (MINLP) problem. From a SOs perspective, the integration of appropriate grid equivalents for neighboring SOs must be added to this task (Section 3.2 and Section 4). From the perspective of central OPFs crossing SO borders, the objectives of the SOs needs to be concluded in one overall objective function (Section 3.3).

3. Previous Methods and Quantification of Methods’ Performances

To measure the performance of the new method proposed in Section 4, a representative of non-coordinated operation (Section 3.1) and the DSO–TSO–DSO chain (Section 3.2) are introduced for comparison. Furthermore, providing the procedure of the DSO–TSO–DSO chain at this point of the paper enhances the motivation and introduction the new method. In Section 3.3, the presented algorithm considers objective functions of multiple SOs fairly. The resulting overall objective function is appropriate as a key performance indicator and can be applied for central optimizations which represent the theoretical optimum when ignoring SO sovereignty.

3.1. Non-Coordinated Operation with Local Control

As assumed in [1], distributed energy resources (DERs) operate with a discrete local control for OLTC tap positions and different local control procedures, depending on the type of DER. Offshore wind applies $Q\left(v_m\right)$ control, half of PV plants are controlled by a $cos\phi \left(p\right)$ curve relation, and all other DER types deploy $Q\left(v_m\right)$ control for one half and $cos\phi \left(p\right)$ control for the other. Even though these control curves and the voltage-controlling generation units in the TSOs are appropriately set for the objective Equation (19), this non-coordinated benchmark should be outperformed by distributed optimization methods.

3.2. DSO–TSO–DSO Chain Method

The DSO–TSO–DSO chain is presented in [18] and extended to the problem formulation of Section 2 in [1]. As the name suggests, the DSO–TSO–DSO chain contains three steps: First, the DSOs determine limits for the reactive power flow at the interfaces and the voltage magnitudes at the boundary buses to provide these limits to the TSOs. Second, OPFs of the TSOs with the TSOs’ objectives define voltage magnitude setpoints at the boundary buses, $v_{m_{\otimes },\mathrm{set}}$. Finally, OPFs of the DSOs determine their control variables by minimizing deviations from the voltage magnitude setpoints.

TSOs consider neighboring TSOs as fixed and DSOs as variable $PQ$ injections. DSOs replace a neighboring TSO with one slack element and $PV$ injections at further boundaries.

The DSO–TSO–DSO chain has a clear procedure and has proved reliable in simulations. However, two essential drawbacks exist: First, power flow mismatches at the interfaces occur due to different values of the OPFs of connected SOs. Second, there is no appropriate balancing of interests; voltage magnitude setpoints at TSO–DSO interfaces are determined considering only the TSOs interests. A weighting of the deviations of the setpoints and the DSO objectives may allow the DSOs to incorporate their interests into the resulting grid conditions, but then the TSO control of Step 2 is no longer optimal.

3.3. Overall Objective Function as a Performance Indicator

As a multi-objective optimization algorithm to select a fair solution from the Pareto front, Equations (21)–(26) are applied, as introduced in [26]. According to this, each SO z determines its optimum $x^{{*}_z}$ and optimal objective value $f_z\left(x^{{*}_z}\right)$ (Equation (21)), where c and d represent the constraints introduced in Equations (1)–(8) and $\mathcal{Z}$ defines the set of SOs involved.

$$\begin{array}{cc}\hfill & min{f}_z\left(x\right)\mathrm{subject}\mathrm{to}c\left(x\right)\le 0\wedge d\left(x\right)=0\forall z\in \mathcal{Z}\hfill \end{array}$$

(21)

From the optima $x^{*}$ determined by Equation (21), an overall objective $f_{\mathrm{oo}}$ can be defined (Equation (22)) to minimize the distance to the SO individual optima on normalized axes. The normalization factors $\varsigma$ and $\chi$ allow a fair combination of different SO objectives. ${\varsigma }_z$ provides the value ranges of the objectives $f_z$ for the optima of the other SOs (Equation (23)). ${\chi }_z$ indicates how much the objectives of the other SOs are impaired by the optimum of z (Equation (24)). Consequently, $\varsigma$ and $\chi$ act to avoid “non-cooperative” objectives and objectives with high values being considered with too much weight in the overall objective.

$$\begin{array}{cc}\hfill & min{f}_{\mathrm{oo}}=min\sum _{z\in \mathcal{Z}}{w}_z{\left(\frac{f_z\left(x\right)-f_z\left(x^{{*}_z}\right)}{{\varsigma }_z·{\chi }_z}\right)}^2\hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & {\varsigma }_z=\sum _{j\in \mathcal{Z}}\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\left|\mathcal{Z}\right|$}\right.\left(f_z\left(x^{{*}_j}\right)-f_z\left(x^{{*}_z}\right)\right)\forall z\in \mathcal{Z}\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & {\chi }_z=\sum _{j\in \mathcal{Z}}\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\varsigma }_j$}\right.\left(f_j\left(x^{{*}_z}\right)-f_j\left(x^{{*}_j}\right)\right)\forall z\in \mathcal{Z}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & c\left(x\right)\le 0\forall z\in \mathcal{Z}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & d\left(x\right)=0\forall z\in \mathcal{Z}\hfill \end{array}$$

(26)

This algorithm with normalized costs (Equations (21)–(26)) is an a posteriori approach. In total, the number of $\left|\mathcal{Z}\right|+1$ optimizations need to be calculated. No predefined parameters are needed and all SOs are treated equally. However, from an overall perspective, it is reasonable to give larger SOs larger weights in the overall objective. Therefore, weighting factors w were added to Equation (22). This paper assumes weights that are independent from variable grid operation states but defined by line length sums l in km as well as last year’s sums of active power supply of connected consumers and lower-voltage SOs E in GWh:

$$\begin{array}{c}\hfill w_z=\frac{\left|\mathcal{Z}\right|}{2}(\frac{l_z}{{\sum }_{z\in \mathcal{Z}}{l}_z}+\frac{E_z}{{\sum }_{z\in \mathcal{Z}}{E}_z})\forall z\in \mathcal{Z}\end{array}$$

(27)

4. New Method with Equivalent Functions

In this section, a new sequential method is proposed for voltage control and reactive power dispatch. This new method, hereafter referred to as the equivalent function method, builds on the strengths of the DSO–TSO–DSO chain and aims at overcoming the methodological drawbacks. To reduce the power flow mismatch, the equivalent function method considers setpoints for voltage magnitudes of the boundary buses and the reactive power flow at the interface to each neighboring SO, see Section 4.1. Furthermore, to fairly weigh the interests of all players, a multi-objective algorithm with equivalent functions is introduced in Section 4.2.

4.1. Overall Procedure

Regarding the equivalent function method, the TSO–TSO interfaces (Steps 1 and 2) are coordinated before the TSO–DSO interfaces (Steps 3 and 4). In both cases, setpoints for voltage magnitudes of the boundary buses $v_{m_{\otimes }}$ are determined first (Steps 1 and 3), followed by the setpoints for the reactive power flow at the interface to each neighboring SO $q_{\otimes }$ (Steps 2 and 4). For TSO–DSO interfaces (Step 4), the procedure has been simplified by concluding the sum of the reactive power flow at each interface as $q_{\otimes }$. In the final Step 5, all SOs operate to meet the coordinated setpoints and to adjust control variables according to their own objectives. This is applied by OPFs with objective functions as defined for SO z in Equation (28). Therein, the quadratic penalization of deviations alongside sufficiently large weighting factors ensures that negligible deviations occur. Furthermore, the factors weigh voltage deviations and reactive power deviations against each other. For the simulations in Section 5, $a_3=2.5$ and $a_4={10}^5$ were chosen in accordance with the value ranges of $f_z$, $v_{m_{\otimes }}$, and $q_{\otimes }$.

$$\begin{array}{c}\hfill f_{z,\mathrm{applied}}=f_z+a_3{\left(q_{\otimes }-q_{\otimes ,\mathrm{set}}\right)}^T\left(q_{\otimes }-q_{\otimes ,\mathrm{set}}\right)+a_4{\left(v_{m_{\otimes }}-v_{m_{\otimes ,\mathrm{set}}}\right)}^T\left(v_{m_{\otimes }}-v_{m_{\otimes ,\mathrm{set}}}\right)\end{array}$$

(28)

Figure 3 outlines the procedure, including substeps (a–e) introduced in Section 4.2. While all TSOs operate equally and all DSOs operate equally, Figure 3 depicts the equivalent function method from the viewpoint of TSO 1. The coordination with connected TSO k and connected DSO l represents the coordination that is performed with all connected SOs. Assuming that buses are connected to more than two TSOs, all coordination of setpoints is performed bilaterally. Horizontal arrows illustrate data exchange between SOs. Passing of data to Steps 1.d, 2.d, and 4.d implies exchange also, as discussed in Section 4.4.

While Steps 1, 2, and 4 consider multi-objective optimization problems, Step 3 follows a simplified approach to determine the voltage magnitude setpoints at TSO–DSO interfaces. Assuming that transformer tap shifting allows control over voltage magnitudes at the connection points as needed for the DSOs, the voltage magnitude setpoints are defined by the TSOs in Step 3.d’. The OPF of TSO 1 considers fixed reactive power flows to connected DSOs, named $q_{\otimes ,{\mathrm{set}}^{\prime }}$. Steps 3.a’ and 3.b’ determine $q_{\otimes ,{\mathrm{set}}^{\prime }}$ as an estimate for the final setpoints that are calculated in Step 4.d with the use of the voltage magnitude setpoints of Step 3.d’.

The method implies the boundary of neighboring SOs as illustrated in Figure 1. In the case of TSO–DSO interfaces, the boundary branches represent the transformers and the boundary buses are located at the TSO side. Neighboring TSOs are considered as variable $PV$ injections and also provide a slack element if needed. Neighboring DSOs are considered as $PQ$ injections, fixed in the first two steps and variable after that.

4.2. Multi-Objective Algorithm with Equivalent Functions

This subsection proposes substeps within Step 1 to 4 to determine setpoints for the boundary variables fairly concerning multiple, unpublished objectives and non-shared grid data. It applies the algorithm with normalized costs (Equations (21)–(26)) which, however, is by design an algorithm for centralized optimizations. It relies on known system variables x and subproblem objectives $f_z$. To get rid of the dependency on system variables of all SOs and to avoid the need for revealing the objectives of the SOs, polynomial equivalent functions $\widetilde{f_z}$ are defined which depend on the boundary variables $x_{\otimes }$ only. Instead of the real objectives of the SOs $f_z$, all existing system variables x, and constraints c and d, the equivalent functions $\widetilde{f_z}$, the boundary variables $x_{\otimes }$, and their constraints as defined in Equation (29) are applied to the algorithm with normalized costs (Equations (21)–(26)). Furthermore, Equation (23) is replaced by Equation (30). In this way, $\varsigma$, representing the value ranges of the objectives, is determined using not only the optima of the SOs but all available information. This information comes from additional samples to determine which lowest objective values are reachable for the SOs at defined boundary variables. The set of sample values $\mathcal{S}$ is explained in the next paragraphs alongside Figure 4.

It is worth noting that using the equivalent functions reduces the OPFs of the algorithm with normalized cost to small optimization problems that are much less complex to solve.

$$\begin{array}{cc}\hfill & x_{\otimes }^{\min }\le x_{\otimes }\le x_{\otimes }^{\max }|x_{\otimes }=v_{m_{\otimes }}\vee x_{\otimes }=q_{\otimes }\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & {\varsigma }_z=\sum _{s\in \mathcal{S}}\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\left|\mathcal{S}\right|$}\right.\left(f_z\left(x^s\right)-f_z\left(x^{{*}_z}\right)\right)\forall z\in \mathcal{Z}\hfill \end{array}$$

(30)

Including pre- and post-processing, the multi-objective algorithm with equivalent functions may contain up to five substeps, of which Substeps a and e can be partially omitted:

a

Each SO determines the limits of the boundary variables.

b

Each SO determines its optimum.

c

Each SO determines optimal objective values for the optimum variable values of its neighbors as well as for further sample values as input data to generate the equivalent functions (Substep d).

d

Each SO or a joint IT infrastructure derives one second-degree polynomial equivalent function $\widetilde{f_z}$ per SO from the values of Substep c via a regression analysis. Additionally, performing the algorithm with normalized costs and equivalent functions provides the desired multi-objective optimal setpoints.

e

If any SO cannot reach the setpoints determined, the SOs agree on the closest reachable set of boundary variables.

Figure 4 is presented to explain the selection of further sample values as input data to generate the equivalent functions $\widetilde{f_z}$. The left part of Figure 4 depicts the voltage setpoint coordination at the TSO1–TSO2 interface of Steps 1.c–1.d, whereas the right part visualizes the reactive power coordination of Steps 4.c–4.d at the TSO1–DSO3 interface.

m-variable polynomial functions with a degree of two are defined by a number of $\left(\genfrac{}{}{0pt}{}{m+2}{2}\right)=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.(m^2+3m+2)$ parameters. In Steps 1 and 2, $m=n_{\otimes }$ applies, where $n_{\otimes }$ is the number of boundary buses per SO interface. In Step 4, $m=1$ applies because the sum of exchanged reactive power is considered. The optima of two adjacent SOs at one interface provide only two points for the equivalent functions. Consequently, further sample values are needed for a least-squares approximation. For $m=2$, see Figure 4 (left), five further sample values are used. One is the midpoint between the two optima. The other four complete a circle in 60 degree increments around the midpoint. For $m=1$, three sample values are spread between the optima and outer limits, as plotted in Figure 4 (right).

4.3. Scalability and Realizability

The SOs can perform optimizations of the substeps in parallel.

Table 1. Number of optimizations performed consecutively within the equivalent function method (completely parallelized).

Step		1					2					3			4					5		Sum
Substep	b	c	d		a	b	c	d	e		a	b’	d’		a	b	c	d
OPFs	1	$6{}^{\alpha }$			2	1	$6{}^{\alpha }$		2		2	1	1		2	1	4			1		30
Less complex optimizations			2					2										2				6

^α To be adapted if the number of boundary buses at TSO–TSO interfaces is not two.

concludes the remaining number of optimizations performed successively. In Steps 1.c and 2.c, six OPFs are set for the investigated case of two boundary buses at the TSO1–TSO2 interface. It is the only number that needs to be increased if the grid configuration exhibits a higher number of boundary buses. That number should be $\ge \left(\genfrac{}{}{0pt}{}{m+2}{2}\right)$ if the system to determine the equivalent functions remains overdetermined. The number of boundary buses is usually small, even for large TSOs; thus, the increase in the number of optimizations and thus the computing time is acceptable. Otherwise, several options for reduction exist: considering derivatives at the sample points or historical data of previous time steps, combining electrical near boundary buses to determine the same setpoints for combined buses, or not determining selected polynomial coefficients.

The boxplot in Figure 5 illustrates the sum of the computing time of successively executed optimizations per time step executed on a common laptop. Even though the communication time between the SOs is added to real grid operation, more than a quarter-hourly operation of the equivalent function method is possible, even without further acceleration efforts.

4.4. Sovereignty and Data Exchange

The equivalent function method allows the SOs to exchange a small amount of information and to still operate the grid in accordance with their objectives. Exchanged data include operational limits (Substep a), optima and their objective values (Substep b), objective values of further sample values (Substep c), and adjusted setpoints (Substep e). No grid data or information on the actual objective function are shared.

Three risks are common to all methods that include coordination for better grid operations: first, non-cooperative manipulation of the exchanged data by SOs; second, manipulation by third parties; and third, outage of communication due to information and communication technology (ICT) failures.

To avoid manipulations by SOs, various options to certify the IT procedure are conceivable. One is that SOs must commit to reporting changes in their objective functions to the responsible regulator, such as the Bundesnetzagentur in Germany, and to backup relevant data of the grid conditions for a certain period. Within that period, the regulator has the opportunity to check the SOs to determine if the data exchanged to connected SOs are valid or show manipulations.

To prevent third-party manipulations, SOs need to protect the ICT infrastructure, as already necessary.

In case of ICT failure, similar to other methods, the SOs need to establish a backup procedure for such periods. That might be to estimate missing data, to agree on a characteristic curve of boundary variables depending on variables such as the active power exchange, or to reactivate non-coordinated grid operation.

5. Simulation Results

As introduced in Section 3, the time series simulation results of the equivalent function method were compared with the results of non-coordinated local control, the DSO–TSO–DSO chain, and the central OPF which represents the theoretical optimum when ignoring the sovereignty of SOs.

5.1. Case Study Setup

SimBench [27] provides an open, multi-voltage dataset that is well suited for balanced power flow analyses and classified as a benchmark grid [28]. Here, an excerpt of the SimBench dataset, “1-EHVHV-mixed-all-0-no_sw”, is used, which includes two TSOs and two DSOs located in northern Germany, see Figure 6. The first two days of the complete year of the time series data with a time resolution of 15 min are appropriate for validation and illustration because the power flow and in particular the DER injection vary significantly in that time period. Applying Equation (27) to this data, the weights to consider the different sizes of SOs within $f_{\mathrm{oo}}$ are shown in Equation (31).

Table 2. Analyzed combinations of SO objectives.

Combination	Objective	TSO1	TSO2	DSO3	DSO4
1	$f_{\mathrm{profile},\mathrm{loadings}}$	✓	✓	✓	✓
	$f_{\mathrm{losses}}$
2	$f_{\mathrm{profile},\mathrm{loadings}}$
	$f_{\mathrm{losses}}$	✓	✓	✓	✓
3	$f_{\mathrm{profile},\mathrm{loadings}}$			✓	✓
	$f_{\mathrm{losses}}$	✓	✓
4	$f_{\mathrm{profile},\mathrm{loadings}}$		✓		✓
	$f_{\mathrm{losses}}$	✓		✓

outlines four combinations of SO objectives that are analyzed to test the equivalent function method for different objectives.

$$\begin{array}{c}\hfill w={\left[\begin{array}{cccc}1.005& 1.790& 0.581& 0.624\end{array}\right]}^T\end{array}$$

(31)

An open-source tool to model and analyze power systems, pandapower [29], is applied to process grid data and perform power flow calculations within time series simulations. Additionally, optimizations are modeled via IEEOptTool [30], a flexible software written in AMPL [31] to solve OPF in the field of power system operation and planning. IPOPT from the COIN-OR project and KNITRO from Artelys are applied as solvers for nonlinear programming (NLP) and MINLP.

5.2. Summary of Time Series Results and Overall Objective Explanation

Figure 7 outlines the results of time series simulations. It depicts the mean values of the overall objective $f_{\mathrm{oo}}$ of the compared methods and the four objective combinations. The values of the equivalent function method are significantly lower than those of the assumed local control and the DSO–TSO–DSO chain. Hence, a satisfactory performance for different objective combinations is shown.

As with several other multi-objective functions, an intuitive understanding of $f_{\mathrm{oo}}$ is lacking. For that reason, exemplary suboptimalities of individual objectives of SOs are considered in the following to demonstrate the effects on $f_{\mathrm{oo}}$. According to Equation (32), $k_z$ describes the relative suboptimality of the objective of SO z. If $k_z=1$ applies, SO z does not contribute to $f_{\mathrm{oo}}$. The curves plotted in Figure 8 arise by only one SO contributing to $f_{\mathrm{oo}}$ at a time. The left end of the curves, for example, results from ${\phantom{k}}^{\mathrm{TSO}1}k$, where only the first value differs from 1, see Equation (33). To define the four curves, Equation (34) defines a relative suboptimality vector $\kappa$ for each curve. Finally, the $f_{\mathrm{oo}}$ values illustrated in Figure 8 can be calculated using the single SO objective optima and weighting factors $\varsigma$ and $\chi$ of the first time step (Equations (35)–(37)) and the predefined weights w (Equation (31)).

The three colored curves show the difference in $f_{\mathrm{oo}}$ for equal relative suboptimalities. In contrast, the black curve shows that all SOs contribute $1.98\times {10}^{-3}$ to $f_{\mathrm{oo}}$, although the relative suboptimalities range from $1.084$ to $1.446$. These differences result from the weighting factors w, $\varsigma$, and $\chi$.

The sum of these four equal contributions is $7.94\times {10}^{-3}$, which is achieved by the equivalent function method for objective combination 3, see Figure 7. This value also results if all SOs achieve their individual optima, except for DSO3, where $k_3=1.291$, see the point of the orange curve at DSO3 in Figure 8.

$$\begin{array}{cc}\hfill f_z\left(x\right)& =k_z·f_z\left(x^{{*}_z}\right)\forall z\in \mathcal{Z}|k_z\ge 1\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \left[\begin{array}{cccc}{\phantom{k}}^{\mathrm{TSO}1}k& {\phantom{k}}^{\mathrm{TSO}2}k& {\phantom{k}}^{\mathrm{DSO}3}k& {\phantom{k}}^{\mathrm{DSO}4}k\end{array}\right]& =\left[\begin{array}{cccc}\left(\begin{array}{c}{\kappa }_1\\ 1\\ 1\\ 1\end{array}\right)& \left(\begin{array}{c}1\\ {\kappa }_2\\ 1\\ 1\end{array}\right)& \left(\begin{array}{c}1\\ 1\\ {\kappa }_3\\ 1\end{array}\right)& \left(\begin{array}{c}1\\ 1\\ 1\\ {\kappa }_4\end{array}\right)\end{array}\right]\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill \left[\begin{array}{cccc}{\phantom{\kappa }}^{10\%}\kappa & {\phantom{\kappa }}^{29.1\%}\kappa & {\phantom{\kappa }}^{80\%}\kappa & {\phantom{\kappa }}^{\mathrm{variable}}\kappa \end{array}\right]& =\left[\begin{array}{cccc}\left(\begin{array}{c}1.1\\ 1.1\\ 1.1\\ 1.1\end{array}\right)& \left(\begin{array}{c}1.291\\ 1.291\\ 1.291\\ 1.291\end{array}\right)& \left(\begin{array}{c}1.8\\ 1.8\\ 1.8\\ 1.8\end{array}\right)& \left(\begin{array}{c}1.446\\ 1.084\\ 1.145\\ 1.194\end{array}\right)\end{array}\right]\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill f_z\left(x^{{*}_z}\right)& =\left[\begin{array}{cccc}21.16& 80.76& 108.33& 35.46\end{array}\right]\forall z\in \mathcal{Z}\hfill \end{array}$$

(35)

$$\begin{array}{cc}\varsigma \hfill & ={\left[\begin{array}{cccc}39.80& 65.00& 32.76& 48.52\end{array}\right]}^T\hfill \end{array}$$

(36)

$$\begin{array}{cc}\chi \hfill & ={\left[\begin{array}{cccc}4.81& 3.76& 5.64& 1.78\end{array}\right]}^T\hfill \end{array}$$

(37)

5.3. Time Series Results in Temporal Resolution

In Figure 9, the courses of the overall objectives, the objectives of the individual SOs, the voltage magnitudes, and the loadings of the lines of the multi-SO power system are presented for objective combination 3. As with the other three objective combinations, the equivalent function method consistently achieves better overall objective values than the method with local control or the DSO–TSO–DSO chain. Moreover, every individual SO benefits from better results which are similar to the central OPF. The resulting voltage profiles and line loadings similar to the central OPF corroborate this assertion. However, the maximum voltages of the equivalent function do not completely reach the upper limit of the allowed voltage area ($0.9\le v_m\le 1.1$). The reason is that these simulations apply a conservatively narrowed voltage area ($0.92\le v_m\le 1.08$) within the equivalent function method to ensure that the resulting overall system state satisfies the original voltage limits.

In contrast, the DSO–TSO–DSO chain shows no appropriate coordination of the different objectives of the SOs. In particular, suboptimalities arise with the objectives $f_{\mathrm{profiles},\mathrm{loadings},3}$ and $f_{\mathrm{losses},1}$. These are related to the overvoltages induced by the method applied to the configuration of the studied grid. Among the deviations of the resulting boundary variables from the boundary variables assumed in the OPFs, the deviation of the voltage magnitude at the northern boundary bus between TSO1 and TSO2 is the most significant. This deviation in the OPF from TSO1 carries over into the setpoints for DSO3. Consequently, the DSO–TSO–DSO chain cannot be applied in the studied form to multi-TSO, mutli-DSO power systems.

5.4. Limits of the Equivalent Function Method

To provide a realizable, SO-sovereignty-preserving method for grid operation, the equivalent function method uses simplifications and assumptions. These are responsible for the difference compared to the optimum represented by the central OPF. In general, for different power system characteristics and objectives, it is not guaranteed that the simplifications of the method remain appropriate.

Firstly, the individual objectives of the SOs should be sufficiently accurately reproduced by the polynomial functions created by sample values. Figure 4 exemplifies the approximation for Steps 1.d and 4.d, the first time step, TSO1, and objective combination 3. Furthermore, it illustrates the vertical distance between the equivalent functions and the sample values, the maximum distances of which are about 0.005 MW and 0.3 MW. Additionally,

Table 3. Maximum vertical distance between the equivalent function values and the sample values during time series analyses.

Step	1.d			2.d			4.d
System Operator	TSO1	TSO2		TSO1	TSO2		TSO1	TSO2	DSO3	DSO3
Objective combination 1	0.02	0.04		1.83	0.68		1.06	0.00	1.33	0.6
Objective combination 3	0.11 MW	0.09 MW		1.65 MW	17.45 MW		0.37 MW	0.02 MW	2.98 MW	0.05 MW
Objective combination 2	0.12 MW	0.09 MW		1.65 MW	17.32 MW		0.37 MW	0.02 MW	6.89	0.23
Objective combination 4	0.34 MW	0.10		1.11 MW	0.78		0.21 MW	0.00	0.15 MW	0.06

provides the maximum distances of the complete time series simulations for all steps, SOs, and objective combinations. Generally, no significant distances occur. This is supported by truncating the top 5% and bottom 5% of the range from minimum to maximum reactive power flows at the interfaces which might be available only at high cost and thus are not desirable anyway. Only in a few time steps (afternoon and evening of the second day), when TSO2 applies $f_{\mathrm{losses}}$, are the optima of TSO1 so suboptimal for TSO2 that the corresponding high objective values cannot be appropriately approximated by the equivalent function of Step 2.d. This, however, does not impair the performance of the method significantly.

Secondly, if a TSO has weak abilities to control the voltage magnitudes and reactive power exchanges at its interfaces, suboptimal adjustments of the setpoints are needed. In all investigated studies, this happened only to TSO1 in very few cases. It is possible to also add Substep e in Steps 1, 3, and 4 if necessary.

Thirdly, although the equivalent function method does not change any active power injection or load, the interaction of active power with controlled voltage magnitudes and controlled reactive power occurs if interconnections between multiple boundary buses allow circular power flows. If the equivalent function method is applied to other power systems where this interaction significantly impairs the performance, applied external grid representations (TSOs as $PV$ injections and DSOs as $PQ$ injections) should be replaced by grid equivalents that include an active power response, e.g., via impedances between the boundaries.

6. Conclusions

A new method with equivalent functions is proposed to coordinate reactive power, voltage, and transformer tap control in multi-SO grid operation. The preservation of the sovereignty of all SOs is integrated into the method. Similar to the sequential DSO–TSO–DSO chain method, the equivalent function method shows no divergence issues and the step-by-step procedure provides easy-to-understand setpoints at the SO interfaces. To the best of the authors’ knowledge, it is the first method that involves a fair balancing of the interests of all SOs. Time series simulations with four different combinations of objectives of the SOs show that the equivalent function method outperforms non-coordinated local control and DSO–TSO–DSO chain methods. Besides performance analyses, issues in realizing distributed optimization methods for real grid operations are discussed. For example, strategies to avoid non-cooperative manipulations are proposed and the time to compute the steps of the method is assessed as being sufficiently low. Applying the equivalent function method to different power system characteristics and SO objectives, the appropriateness of the described simplifications and assumptions must be validated. Therefore, future work should prove the method in field tests and with further grid characteristics and configurations.