Presently, many ideas have been proposed for constructing zonal models of the European transmission network. Most of them are based on transmission network reduction and bus clustering, which requires the ex-ante knowledge on the details (impedance, interconnections, etc.) of the network’s elements. Generally, the information is not publicly disclosed. This paper describes a method for serving the purpose, which only requires the ex-post (i.e., published) inter-zonal power flow data as input. It performs inverse modeling on the ex-post data to reveal the virtual reactance of the interconnectors in the zonal model. The reactance values are then used for computing the zonal Power Transfer Distribution Factors (PTDFs) of the model, which is the best representation of the loop flows. A large number of such zonal PTDFs could exist in one zonal model because of its dependency on the Generation Shift Keys (GSKs). Theoretically, the GSKs could exist as an infinite set of unique combinations but only a finite number of them applied in real system operations. As a solution, the authors propose a heuristic approach for grouping the ex-post power flow cases into clusters, and find the zonal PTDF for each of them; in order to effectively reduce the domain of the zonal PTDFs and enhance the tractability of the solution. The identified zonal PTDFs meaningfully represent the zonal model. Additionally, they can also be used for studying the best and worst case scenarios of cross-border security-constrained economic dispatch. Keywords: Zonal model, transmission network modeling, inverse problem, test case clustering, clustering heuristics I Introduction The last decade has seen several exciting developments in the European single market in electricity [1]. Among others, the development of a common market design through system couplings remain active [2]. As a result, the Flow-Based Market Coupling (FBMC) mechanism is established in the Central-Western Europe (CWE) in May 2015 [3]. Regulators in other regions too (for example, the Nordic Regional Security Coordinator) had begun to study the feasibility of FBMC in their respective region. Eventually, the model is expected to be extended to the entire Europe [3]. The FBMC model entails many contributions from the stakeholders (mainly the Transmission Systems Operators (TSOs)) to establish a common set of rules for the sharing of generation and transfer capacity between the coupled market zones. Its most important design is the use of a multi- zonal common grid model [4] (henceforth called as the zonal model) to properly model the cross- border loop flow of electric power. This is important for the market and operational planning involving Security Constrained Economic Dispatch (SCED). 2 For example, the 9-bus transmission network model shown in Figure 1 (a) can be represented by a 3-bus zonal model (Figure 1 (b)). The SCED could be calculated on the zonal models to derive the least-cost dispatch for serving the demand of the coupled markets. The calculation requires the knowledge of the aggregated Power Transfer Distribution Factors (PTDFs), 𝛽֕⟶֕ɱӴֈ͗ which is the magnitude of change in power transferred at a cross-border interconnection m (comprise of all lines that connects buses which belongs to two market zones) contributed by the transfer of one unit of power at a zone, (zone z) to a reference zone (zone 𝑧Ј); and the allocated cross-border transfer capacity. The transfer capacity can be allocated more precisely compared to the existing market coupling (the existing method is based on the available transfer capacity approach), because the loop flow pattern in the entire zonal model is well defined by the 𝛽֕⟶֕ɱӴֈ͗. (a) (b) Figure 1. (a) A 9-bus transmission network model, (b) the 3-bus zonal model of the transmission network II The Formulations The greatest advantage of FBMC is its ability to remove the impediments to competition by means of promoting the non-discriminatory access to the external power generation and transmission services through proper modeling of the cross-border loop flows. This is the key to the running of a successful regional electricity market [5]. Nonetheless, its implementation is not straightforward. Especially, an infinite number of unique zonal PTDF matrix Γ͗ (a two-dimension matrix containing all 𝛽֕⟶֕ɱӴֈ͗) could exist in the zonal model. Each of them is associated with a Generation Shift Keys (GSKs) matrix, Π֕ which comprise the commitment of the generators in zone z to the production of one unit of power. The Γ͗ is derived based on the nodal (bus-level) PTDFs matrix, Γ comprises all 𝛽ք⟶քɱӴֆ (the magnitude of change in power transfer at transmission line k corresponding to the transfer of a one unit of power at a bus (bus i) to a reference bus (bus 𝑖Ј)). It is calculated as follows [9] Γ=BցB֍−φ (1) 3 where, B is a matrix defined based on the susceptance of the transmission lines, 𝑏քօ= 1𝑥քօ⁄ (where 𝑥քօ is the reactance of a transmission line k which connects bus i and j), such that its non-diagonal elements are, Bքօ= Bօք= −𝑏քօ; ∀𝑖≠𝑗 and its diagonal elements are, Bքք=∑𝑏քօN օ≠ք. The B֍−φ is the inverse of the reduced B matrix, B֍ (the row and column in the B matrix associated with an arbitrarily selected reference bus is removed to form the B֍). The Bց is a two-dimension matrix where each row is associated with transmission line k; in each row, the value is 𝑏քօ at column i, and −𝑏քօ at column j. The column associated with the reference bus is as well removed. The derived Γ is a two-dimension matrix where the row index correspond to transmission line k and each column corresponds to bus i. The Γ͗ cannot be directly calculated by (1), since the reactance of its cross-border interconnectors (interconnector m), 𝑥ֈ is not known. Also, the 𝑥ֈ is not a simple sum of 𝑥քօ of all of its constituting transmission lines. Instead, it has to be derived using Γ and the GSKs of generators in each zone, Π֕=ृ𝜋φӴ֕,⋯,𝜋GӴ֕ॄ, where ∑𝜋ւӴ֕G ւ=φ= 1. For example, Πφ={0,0.4,0.5,0.1 } indicates the first generator in zone 1 is not involved in the dispatch of the next unit of power, the second generator will dispatch 40% of it, the third generator will dispatch 50% of it and the forth generator will dispatch 10% of it. The 𝛽֕⟶֕ɱӴֈA is calculated as follows 𝛽֕⟶֕ɱӴֈA=ం ঁం𝜋ւӴ֕𝛽ք∗⟶քɱӴֆֆং,∀𝑘∈𝐾ֈ,𝑧∈{𝑧Ј,𝑧} ւ,𝑧≠𝑧Ј (2) where, 𝑖∗ is the index of the bus which is connected with generator g, 𝑧Ј is the index of the reference zone, and 𝐾ֈ is the set of indices of all the constituting transmission lines of interconnector m. Table I shows a sample calculation of 𝛽ϵ⟶φӴφA of the zonal model in Figure 1 (b), which is the changes in power flow of interconnector I1 associated to the supply of 1 unit of marginal power in zone A2 to the reference zone, A1 (the 𝜋ւӴ֕ of generators in zone 1 is of negative sign). The result shows a 1 MW change in A2 will produce 0.6825 MW power flow on I1 in reverse direction (indicated by the negative sign). Table I. An example of calculation of 𝛽ϵ⟶φӴφA of the zonal model of Figure 1 (b) Zone g Bus 𝜋ւӴ֕ Line 𝛽ք∗⟶քɱӴֆ 𝜋ւӴ֕𝛽ք∗⟶քɱӴֆ A1 1 B1 0 L5 0 0 A1 1 B1 0 L6 0 0 A1 2 B2 -0.4 L5 -0.0090 0.0036 A1 2 B2 -0.4 L6 -0.0381 0.0152 A1 3 B3 -0.5 L5 -0.0186 0.0093 A1 3 B3 -0.5 L6 -0.0785 0.0393 A1 4 B4 -0.1 L5 -0.0312 0.0031 A1 4 B4 -0.1 L6 -0.1321 0.0132 A2 1 B8 0.7 L5 -0.5612 -0.3928 A2 1 B8 0.7 L6 -0.2590 -0.1813 A2 2 B9 0.3 L5 -0.1224 -0.0367 A2 2 B9 0.3 L6 -0.5180 -0.1554 𝛽ϵ⟶φӴφA -0.6825 4 The entire ΓA could be fully computed by repeating the calculation for all of the zones and interconnectors. Once the ΓA is computed, the power flow at the interconnector m, 𝐹ֈ can be computed as follows 𝐹ֈ=ం𝑃֕A𝛽֕⟶֕ɱӴֈAZ ֕=φ (3) where, 𝑃֕A is the sum of power generation in zone z. As shown in (2), the ΓA is a function of Π֕ because the 𝛽ք∗⟶քɱӴֆ is constant in value. This is a known shortcoming of FBMC [6] [7]. Generally, there could exist an infinite number of possible Π֕( by changing the 𝜋ւӴ֕ )which yield an infinite set of ΓA. As a different combination of 𝜋ւӴ֕ values (in a unique Π֕ )would yield a completely different power flow at the interconnectors, the domain of ΓA of a zonal model is non-linear. In fact, a change in power dispatched by a generator located at a different position relative to an interconnector would result in a different power flow (in magnitude and direction) on it. Table II shows an example on the relationship between the 𝛽ϵ⟶φӴφA and Πϵ of the 3-bus model (Figure 1 (b)). The value of 𝛽ϵ⟶φӴφA is different for a different Πϵ. Table II. A example of the non-linear relationship between ΓA and Πϵ Bus 𝜋ւӴ֕φ 𝜋ւӴ֕ϵ 𝜋ւӴ֕ϯ B8 0.7 0.9 0.5 B9 0.3 0.1 0.5 𝛽ϵ⟶φӴφբ -0.68 -0.72 -0.65 Table III. Three cases of Π֕ of the 3-bus zonal model (Figure 1 (b)) Zone Bus 𝜋ւӴ֕ref 𝜋ւӴ֕1 Δ𝜋ւӴ֕φ 𝜋ւӴ֕2 Δ𝜋ւӴ֕ϵ A1 B1 0 0 - 0 - A1 B2 -0.4 -0.400 - -0.400 - A1 B3 -0.5 -0.525 5 % -0.525 5 % A1 B4 -0.1 -0.075 -25 % -0.075 -25 % A3 B5 0.5 0.500 - 0.500 - A3 B6 0.4 0.420 5 % 0.420 5 % A3 B7 0.1 0.080 -20 % 0.080 -20 % A2 B8 0.7 0.700 - 0.100 -86 % A2 B9 0.3 0.300 - 0.900 200 % III The non-linearity of the zonal PTDFs Nonetheless, the changes in the ΓA is less sensitive to minor changes in Π֕ .An example is shown in Table III, which comprises a reference case with 𝜋ւӴ֕ref values, a case (case 1) with 𝜋ւӴ֕φ being slightly different from the 𝜋ւӴ֕ref (up to 25 %), and a case (case 2) with 𝜋ւӴ֕ϵ being substantially different than 𝜋ւӴ֕ref (up to 200%). Table IV shows the difference between 𝛽֕⟶֕ɱӴֈA of case 1 and 2 and the reference case. As shown, case 1 only entails a 0.61 % maximum shift in 𝛽֕⟶֕ɱӴֈA. Meanwhile, case 5 2 entails a greater shift in 𝛽֕⟶֕ɱӴֈA for up to 33.49 %. Thus, the ΓA of cases with low variance in Π֕ are approximately similar. For example, the ΓA of the reference case could be applied to approximate the 𝐹ֈ without a great compromise in its accuracy. Table IV. The 𝛽֕⟶֕ɱӴֈA value for each case given in Table III Zone Interface 𝛽֕⟶֕ɱӴֈAӴ ref 𝛽֕⟶֕ɱӴֈAӴ 1 Δ𝛽֕⟶֕ɱӴֈAӴ 1 𝛽֕⟶֕ɱӴֈAӴ 2 Δ𝛽֕⟶֕ɱӴֈAӴϵ A1 I1 0 0 0 0 0 A1 I2 0 0 0 0 0 A1 I3 0 0 0 0 0 A2 I1 -0.6825 -0.6842 -0.25 % -0.5763 15.56% A2 I2 0.3174 0.3157 -0.54 % 0.4237 33.49% A2 I3 -0.3174 -0.3157 0.54 % -0.4237 -33.49% A3 I1 -0.2118 -0.2105 0.61 % -0.2105 0.61% A3 I2 -0.2118 -0.2105 0.61 % -0.2105 0.61% A3 I3 -0.78881 -0.7894 -0.07 % -0.7894 -0.16% Hence, a few clusters of cases which share an approximately identical ΓA exist in a zonal model. The cases can be identified and grouped based on their proximity in Π֕ .Each of them posseses an approximate ΓA,∗, that could be used for approximating their actual 𝐹ֈ with minor errors. This greatly reduces the (otherwise infinite) number of ΓA in the zonal model, making it (and the FBMC) a more tractable problem. The concept also resembles the actual network operations where there are preferably less diversity of Π֕ in used. IV The inverse modeling of transmission network The grouping of the cases into clusters would effectively reduce the non-linear domain of zonal model. Nevertheless, the zonal modeling still could not be done as the details of the transmission network (i.e., the 𝑥ֆ and Π֕ )are not known. The authors applied the inverse modeling [4] to construct the zonal model and identify all of its ΓAӴ∗ based on the ex-post inter-zonal power flow data. 4.1 The method The inverse modeling of a transmission network is introduced in [4]. It is derived based on the linear power flow theory [8]. As shown in (1), knowing the 𝑥ֆ is the most important step for the computation of the Γ. Similarly, knowing the 𝑥ֈ allows us to compute the ΓA. The authors applied the inverse modeling and the clustering approach to compute the list of all ΓAӴ∗ in a zonal model. The linear power flow theory relates the voltage angle to the net power generated at the buses of a transmission network as follows P= BΘ (4) where, P = [𝑃φ,…,𝑃λ]⊺ contains the net active power at each bus, and Θ = [𝜃φ,…,𝜃N]⊺ contains the voltage angle at each bus. The knowledge of Θ is needed for calculating the power flow on each transmission line, which is computed as follows 6 Θ = B−φP (5) Then, the power flow on the transmission line k, 𝐹ֆ which connects buses i and j can be calculated as follows 𝐹ֆ=Bքօि𝜃ք−𝜃օी (6) The inverse modeling intends to find an approximated value of all 𝑥ֆ in the network, or X∗= [𝑥φ∗,…,𝑥K∗], that will meet the following objective: 𝜙=ం|𝐹ֆ∗−𝐹ֆ|≈0έ ֆ=φ (7) The 𝐹ֆ is the ex-post power flow, and 𝐹ֆ∗ is the power flow calculated using the X∗ based on (5) and (6). The approximation is successful when an X∗ that yield 𝐹ֆ∗ value closest possible to the known 𝐹ֆ for all of the transmission lines is found. There are two stages of calculation involved. Firstly, the approximated 𝑥ֆ has to be determined and the B matrix (in (5)) has to be calculated. Secondly, the Θ and 𝐹ֆ∗ has to be calculated. The problem has to be solved iteratively. The calculation of the approximated Γ, Γ∗ is performed using (1) after an optimal X∗ is found. The approach is directly applicable to zonal modeling, by designating the zones as buses, the cross-border interconnectors as the transmission lines, and the cross-border power flow as the 𝐹ֆ. The X∗ derived is used to compute the approximated ΓA, ΓAӴ∗. The authors applied the Genetic Algorithm (GA) as the solver. The details of the GA’s implementation is described in [4]. During implementation, it is possible to find several unique X∗ that would meet the 𝜙 (in (7)). The solution space of X∗ is thus unbounded. It is difficult to ensure the uniqueness of the solution by only inversing one set of the ex-post power flow data. In such case, the GA finds different sub-optimal X∗ in different repetitions, while only one of them is globally optimal. This problem could be solved by adopting multiple ex-post power flow data (henceforth denoted by index s) in 𝜙 to evaluate the quality of X∗: 𝜙=ంం|𝐹ֆ∗Ӵ֎−𝐹ֆ֎|≈0έ ֆ=φϢ ֎ (8) Generally, a globally optimal X∗ could be found by adopting 3 cases in (8). Interestingly, the concept of the cross-border flow clustering matches well with (8) as the clusters comprises a group of ex-post zonal flow cases with approximately identical X∗ and ΓAӴ∗. The main advantage of the inverse modeling is its capability in identifying all ΓAӴ∗ in the zonal network model without requiring the ex-ante knowledge of Π֕ .The method is especially useful for correctly identifying all of the ΓAӴ∗ of a zonal model. For transmission network model reduction (where the transmission network details are fully known), it would be sufficient to use (2) to calculate the ΓAӴ∗ of the zonal model. Then, a zonal model based SCED could be meaningfully implemented on the ΓAӴ∗. 4.2 The clustering of ex-post power flow cases The clustering of the cases in the ex-post data of a zonal model could be done using a simple heuristic method. The cases (denoted by notion s) in a cluster (that has the same ΓAӴ∗) would exhibit the following characteristics: a) Each of the zones in the cases will have an approximately identical value in power generation ratio, 𝑤֕P compared to their counterpart in the other cases of the same cluster. For each case, 7 the 𝑤֕P is calculated by dividing the aggregated total generation at zone z, 𝑃֕A by a weighting factor, W which is the largest absolute value of the 𝑃֕A in the case. For a case denoted with s, the 𝑊֎ is calculated as follows 𝑊֎=max(|𝑃֕AӴ֎|); ∀𝑧∈{1,…,Z} (9) Hence, the 𝑤֕PӴ֎ is defined as: 𝑤֕PӴ֎=𝑃֕AӴ֎ 𝑊֎ (10) The cases could potentially be grouped in the same cluster if they meet the following criteria: 𝑤֕PӴ֎≈𝑤֕PӴ֎ref; ∀𝑠∈{1,⋯,S}; ∀𝑧∈{1,⋯,Z} (11) where, 𝑠ref is an arbitrarily selected case. b) The cases which met the criteria of (11) could be grouped in the same cluster if they have an approximately identical value in their interconnector power flow ratio, 𝑤ֈΑӴ֎ compared with their counterpart in the other cases. The ratio is calculated by dividing the 𝐹ֈ֎ by the 𝑊֎ in (9): 𝑤ֈFӴ֎=𝐹ֈ֎ 𝑊֎ (12) The cases from the same cluster (possessing the same ΓբӴ∗) would meet the following criteria: 𝑤ֈFӴ֎≈𝑤ֈFӴ֎ref; ∀𝑠∈{1,…, S}; ∀𝑧∈{1,…, Z} (13) For instance, the 𝑤ֈF of I1, 𝑤φF of Figure 1 (b) of two arbitrary cases are formulated as follows 𝑤φFӴφ =𝐹φφ 𝑊φ (14) =𝑃ϵAӴφ𝛽ϵ⟶φӴφAӴφ+𝑃ϯAӴφ𝛽ϯ⟶φӴφAӴφ 𝑊φ =𝑤ϵPӴφ𝛽ϵ⟶φӴφAӴφ+ 𝑤ϯPӴφ𝛽ϯ⟶φӴφAӴφ and, 𝑤φFӴϵ =𝐹φϵ 𝑊ϵ (15) =𝑃ϵAӴϵ𝛽ϵ⟶φӴφAӴϵ+𝑃ϯAӴϵ𝛽ϯ⟶φӴφAӴϵ 𝑊φ =𝑤ϵPӴϵ𝛽ϵ⟶φӴφAӴϵ+ 𝑤ϯPӴϵ𝛽ϯ⟶φӴφAӴϵ In (14) and (15), if 𝑤φFӴφ≈𝑤φFӴϵ and [𝑤ϵPӴφ,𝑤ϯPӴφ ]⊺≈ [𝑤ϵPӴϵ,𝑤ϯPӴϵ]⊺, then [𝛽ϵ⟶φӴφAӴφ,𝛽ϯ⟶φӴφAӴφ]⊺= [𝛽ϵ⟶φӴφAӴϵ,𝛽ϯ⟶φӴφAӴϵ]⊺. Both case 1 and 2 would have the same ΓAӴ∗. Table V shows an example of three ex-post power flow cases of the same cluster, where three of them have approximately identical value of 𝑤֕PӴ֎ and 𝑤ֈFӴ֎. A Γ͗Ӵ∗ could be found for the cluster. Table VI, and VII respectively shows the ΓAӴ∗ found by GA (with 𝜙= 0.0013 ) and the precision of the resulting power flow at the interconnectors computed using the ΓAӴ∗. Basically, this is a hard clustering approach because the cases in each cluster must strictly fulfil (11) and (13). In reality, the ΓAӴ∗ of some of the distinct clusters are identical. The cases are identified as the members of different cluster because they have different pattern of 𝑃֕A. This issue can be handled by testing a new case with ΓAӴ∗ of identified clusters. The cases belong to the same cluster if the identified ΓAӴ∗ could reproduce the 𝐹ֈ of all the interconnectors in the new case. 8 Table V. An example of three ex-post power flow cases Zone 𝑃֕AӴφ 𝑤֕PӴφ 𝑃֕AӴϵ 𝑤֕PӴϵ 𝑃֕AӴϯ 𝑤֕PӴϯ A1 12 1.00 5.0 1.00 45 1.00 A2 -3 -0.25 -1.3 -0.26 -11 0.24 A3 -9 -0.75 -3.7 -0.74 -34 -0.76 Interconnector 𝐹ֈφ 𝑤ֈFӴφ 𝐹ֈϵ 𝑤ֈFӴϵ 𝐹ֈϯ 𝑤ֈFӴϯ I1 3.9537 0.3295 1.6709 0.3342 14.7087 0.3269 I2 0.9540 0.0795 0.3710 0.0742 3.7098 0.0824 I3 8.0451 0.6704 3.3286 0.6657 30.2868 0.6730 Table VI. The ΓAӴ∗ of the test cases (in Table V) Zone Interconnector 𝛽֕⟶֕ɱӴֈբ A1 I1 0 A1 I2 0 A1 I3 0 A2 I1 -0.684672 A2 I2 0.315328 A2 I3 -0.315328 A3 I1 -0.211141 A3 I2 -0.211141 A3 I3 -0.788859 Table VII. The precision of 𝐹ֈ computed using the ΓAӴ∗ (in Table VI) Case Interconnector 𝐹ֈ∗ 𝐹ֈ |𝐹ֈ∗−𝐹ֈ| 1 I1 3.9539 3.9537 0.0002 1 I2 0.9542 0.9540 0.0002 1 I3 8.0449 8.0451 0.0002 2 I1 1.6711 1.6709 0.0002 2 I2 0.3713 0.3710 0.0003 2 I3 3.3284 3.3286 0.0002 3 I1 14.7807 14.7087 0 3 I2 3.7098 3.7098 0 3 I3 30.2868 30.2868 0 𝜙 0.0013 V Discussions and conclusions The paper described a method to construct the zonal network model based on ex-post power flow data. It described the concept of ΓA, which is a function of Π֕ and 𝛽ք∗⟶քɱӴֆ. It also pointed out that the domain of ΓA is large and non-linear, because of its dependency of Π֕ which exist in unique combinations. Then it demonstrates that the ΓA is insensitive to reasonably small changes in Π֕ . 9 This allows the inverse modeling of cases with close proximity in ΓA to derive the ΓAӴ∗ through the iterative search of the X∗. The authors introduced a simple heuristic approach to identify cases with identical ΓAӴ∗. The clustering effectively reduce the domain of the inverse modeling problem, thus enhancing its tractability. The method is developed for serving market coupling studies under the FBMC initiatives. Its major strength is its independence on the knowledge of the physical transmission network. It only requires the TSO’s historical cross-border power flow data as input in its implementation. Besides the FBMC studies, the method can also be used for long-term (multi-stage) simulation of SCED (i.e., hydroelectric power scheduling). It has substantially better representation of the physical flow, compare to the point-to-point simulation approach. The multiple ΓAӴ∗ of the zonal network model derived could be used for studying the best and worst case grid operation scenarios.