Novel Harmonic Distortion Prediction Methods for Meshed Transmission Grids with Large Amount of Underground Cables ^†

Abstract

The tremendous and fast green transition in Denmark has initiated the large-scale grid-integration of renewable energy sources, electrification of energy consumption, and establishment of PtX and Energy Islands, setting goals for transmission grid development—such as the establishment of new connections—and for grid reconstruction—such as the extensive substitution of overhead lines (OHLs) with underground cables (UGCs). The share of UGCs in the Danish transmission grid is increasing. Presence of UGC has resulted in that resonances of the harmonic impedance characteristics of the transmission grid are brought within the harmonic order range coinciding with the harmonic emission sources and causing systemwide increase of the harmonic voltage distortion in the 400 kV transmission grid. The transformation of the 400 kV transmission grid has given rise to the need to predict harmonic voltage distortion using simulation models to secure an adequate power quality and support investment decisions and harmonic mitigation for the grid stage, which has not yet been established and which differs from the present grid. This paper presents the experiences of Energinet, the Transmission System Operator (TSO) in Denmark, with harmonic distortion in the Danish transmission grid due to the establishment of 400 kV UGCs, and the development of measurement-validated methods for harmonic distortion simulation and prediction. The paper also presents ongoing developments within, and research addressing, the prediction of harmonic distortion in meshed grids; for example, it explores where and how an analytical approach can replace observational studies with many numerical simulations. The methods shall make it possible to predict whether, where in the transmission grid, and for which harmonic orders connections that have not yet been commissioned may cause the violation of the planning levels, and which mitigations are necessary for bringing the harmonic distortion below the planning levels with respect to a given margin.

1. Introduction

The green transition of energy sectors in Denmark leads to a greater level of renewable energy utilization within electric energy supply, consumption, transportation, heating, and storage, instead of fossil energy sources. It also leads to the expansion and reconstruction of the transmission grid and a stronger coupling with neighboring systems. ‘Expansion’ means the commissioning of new transmission lines while ‘reconstruction’ is the utilization of underground cables (UGCs) instead of overhead lines (OHLs) in an existing transmission grid. Considering the power quality, the green transition may introduce, on the one hand, the evolution of harmonic emission sources, including more converter-interfaced units within supply, storage, consumption, and transportation [1,2,3,4,5], and, on the other hand, the development, reinforcement, and reconstruction of electricity infrastructure, altering the harmonic impedance characteristics and resonances of the transmission grid due to the introduction of more UGCs [6,7,8,9,10,11].

In Denmark, the first cases of severe harmonic amplification following the energization of the HVAC UGC have already been reported, such as in the Anholt case [12] and the Vejle-Ådal 400 kV UGC [13]. Therefore, securing adequate power quality, and mitigating excessive harmonic voltage distortion, are necessary for a successful green transition and for the development of the transmission grid. Being able to predict the harmonic voltage distortion in a future grid with a larger share of UGCs using measurement-validated simulation models, and to propose harmonic mitigation, are important goals for grid analysis and planning in many countries [6,7,8,9,10,11,12,13].

For learning more about potential barriers and solutions related to the establishment and operation of UGCs in meshed transmission grids, Energinet has initiated the DANPAC 2020 (DANish Power systems with Ac Cables) research project, with the development of a harmonic simulation approach and the mitigation for the meshed transmission grid with UGCs numbering among the goals. In 2020, this work resulted in the development of a novel, deterministic, measurement-validated method for the direct simulation and prediction of harmonic voltage distortion in the meshed transmission grid [14,15].

The harmonic assessment of transmission and distribution grids with converter-interfaced units can be viewed from two basic approaches. The first approach includes [1,2,3,4,5]:

• Detailed representation of the converter-interfaced units and their control systems, allowing EMT or frequency domain assessment of the power quality characteristics of the units;
• Radial type connections of the units to a bulk transmission grid model;
• The use of a frequency dependent network equivalent (FDNE), i.e., a Thevenin equivalent with a single harmonic voltage source behind a frequency-dependent impedance.

From a TSO perspective, the first approach does not show how the harmonic voltage magnitudes in the different substations within the physical meshed transmission grid will be influenced, how the new units will interact with the existing units and affect the harmonic distortion in different parts of the meshed transmission grid, or how the efficiency of the existing harmonic filters in the meshed grid will be affected by the connection of new units through a radial connection. The concerns that have arisen are real; the possibility that the establishment of radial connections with harmonic filters to one substation of the meshed transmission grid may reduce the efficiency of the existing harmonic filters in other substations of the grid has already been reported [1].

Further, a transmission grid assessment will normally include many combinations of the operation regimes of the units (active harmonic emission) and the (n − L), L ≥ 0, operation regimes (passive grid), setting up a requirement on fast execution of many simulations. The application of overly detailed EMT-type models for harmonic assessment will contradict with the abovementioned requirement.

The second approach includes detailed models of radial or meshed transmission grids, specifically:

• Exemplified grid representations, which can be sufficient for general findings on harmonic propagation in radial or meshed transmission grids [9] but are often presented without relation to a physical transmission grid or the possibility for validation using harmonic voltage measurements;
• The notification of challenges, such as difficulties with the modelling of harmonic impedance and emission of the LCC HVDC Converter Stations [7] but without giving sufficient details for the ready application of such unit models for the harmonic assessment of large, meshed transmission grids with several LCC HVDC Converter Stations;
• The application of impedance characteristics and their gain factors for the estimation of background harmonic distortion in a future grid using the harmonic voltage measurements in the present grid [8], which seems a conservative methodology and may lead to proposals concerning the extensive need for new harmonic filters.

From a TSO perspective, the goal is to apply a measurement-validated simulation model for the assessment (in the present grid stage) and prediction (in a future grid development stage) of harmonic voltage distortion in the physical, meshed transmission grid. The simulation model must be accurate in the prediction of both magnitudes and change tendencies of the harmonic voltage distortion, and must also be characterized by fast execution time.

This paper presents the experience with harmonic distortion and describes a novel method for the direct simulation of harmonic voltage distortion in the physical, meshed transmission grid of Denmark. Specifically, this method proposes for how a harmonic assessment with several LCC HVDC Converter Stations can be efficiently conducted, accounting for the different operation regimes of the stations and the grid itself. The described deterministic method utilizes a well-known technique of harmonic load flow, with the injection of either harmonic current or harmonic voltage into the specific nodes of the grid, for the simulation of the background harmonic voltage distortion. The novice of the described method includes the following characteristics:

• The described method works with a hypothesis that certain characteristics of the harmonic current or voltage vectors are identifiable as invariants of the operation conditions of the passive meshed transmission grid.
• These invariant characteristics are assigned using numerical tuning of the harmonic current and voltage vectors for different operation regimes of the meshed transmission grid.
• The result of numerical tuning should satisfy the measured harmonic voltage magnitudes in all substations of the meshed transmission grid simultaneously. In this study, it was conducted within the Western Danish 400 kV grid for all operation regimes included in the tuning.
• Once tuned, these invariant characteristics of the harmonic current and voltage vectors are locked, i.e., they are not allowed to be changed, and the simulation model can be applied for the harmonic assessment of other operation regimes of the present meshed grid or for the expansion of the transmission grid through means such as the establishment of new transmission lines (i.e., modifying the passive part of the meshed grid).
• A successful validation is a must for the simulation model of the meshed transmission grid with the tuned and locked characteristics of the harmonic current and voltage vectors, using measured harmonic voltage magnitudes in all substations with available measurements, that factors in available measurements.

Energinet currently applies this deterministic method for the prediction and mitigation of excessive harmonic distortion during the planning and projecting of new lines and the reconstruction of the transmission grid in Denmark. Denmark consists, in fact, of two asynchronous HVAC transmission grids (using the same nominal frequency of 50 Hz, but out of phase with each other):

• The first is Eastern Denmark, which includes the main island of Zealand with the capital Copenhagen and the islands of Lolland, Falster, and Moen. The Eastern Danish transmission grid is HVAC-connected to Sweden and belongs to the Nordic synchronous system [15];
• The second is Western Denmark, which is HVAC-connected to Germany and belongs to the Continental European synchronous system.

Our conference paper “Harmonic distortion prediction method for a meshed transmission grid with distributed harmonic emission sources—Eastern Danish transmission grid case study” [15] presented the validation cases for the Eastern Danish transmission grid, while this extended paper focuses on the Western Danish grid and provides more details regarding the theoretical background and numerical preparation of the harmonic source models. The application, and successful validation, of the presented deterministic method for both transmission grids of Denmark will strengthen confidence in the method.

The disposition of this paper is as follows. Section 2 gives an overview of the practical experiences, challenges, and findings so far regarding harmonic voltage distortion in the Danish 400 kV transmission grid. Section 3 and Section 4 describes the deterministic method used by Energinet and the preparation of the simulation model for harmonic assessment. Section 5 provides the validation examples of the simulation model and Section 6 presents a case involving the prediction of harmonic voltage distortion following the modification of a 400 kV transmission line within the meshed transmission grid. Section 7 discusses the uncertainties of the presented method and how these uncertainties should be handled. Section 8 presents a statistical assessment confirming a harmonic propagation pattern in the Danish 400 kV transmission grid which was previously observed in simulations using the deterministic method. Section 9 and Section 10 provide the method-specific and general conclusions of the presented work, respectively.

2. Grid Development and Experience with Harmonic Distortion

Specifically for Denmark, principles for developing the electricity infrastructure were defined in the political agreement in the Danish Parliament of October 2020:

• New 400 kV transmission lines must be established using UGCs to the extent that is technically feasible. Beyond what is technically feasible to establish using UGCs, OHLs must be used for 400 kV lines;
• Existing 132–150 kV OHLs are replaced with UGCs when the need for extensive reinvestment in the lines arises. Furthermore, existing 132–150 kV OHLs must be replaced with UGCs if the lines are in the vicinity of new 400 kV OHLs;
• New 132–150 kV transmission lines must be established using UGCs.

The share of UGCs in the meshed transmission grid increases, bringing harmonic resonances down within the harmonic order range coinciding with the background harmonic distortion. For the Danish 400 kV transmission grid, the most critical harmonic orders are the 3rd, 5th, 7th, 11th, 13th, 23rd, and 25th harmonic orders [12,13,14,15,16].

For the evaluation and definition of excessive magnitudes of background harmonic distortion in the Danish transmission grid, Energinet applies the IEC planning levels according to [17]:

• The application of the 95th weekly percentiles of the ten-minute average harmonic voltage magnitudes when verifying the compliance of the measured harmonic distortion to the planning levels in [17];
• The application of the maxima of the three phase-to-ground harmonic voltage magnitudes, with reference to the nominal-frequency voltage, for comparison to the planning levels in [17].

The coordination of power quality compliance with the neighboring transmission system operators, such as in Germany (due to the interconnectedness of the Western Danish and German 400 kV transmission systems), follows the European Commission Regulations [18,19]. The German TSOs also apply the IEC planning levels [17]. However, the three-phase root-mean-square harmonic voltage magnitudes are applied for comparison to the planning levels in [17].

Considering the power quality issues arising due to the increasing utilization of UGCs, the current Danish experience is briefly presented in the section below. In Figure 1, each project is marked A, and each of the locations with changed harmonic voltage distortion as a result of the projects is marked B.

Western Denmark includes the peninsula of Jutland, which borders Germany, and the island of Funen, and these locations are listed in the upper-right corner of Figure 1.

The energization of the relatively short (with a line length of about 7 km) 400 kV double UGC system of Vejle-Ådal, between the substations Landerupgård and Malling in July 2017, has led to a significant amplification of the 11th harmonic voltage distortion in the 400 kV substations Trige and Fraugde [13,14] which are located approx. 90 km and 80 km from Vejle-Ådal, respectively. The harmonic distortions in the other 400 kV substations with harmonic voltage measurements and of the other harmonic orders have remained almost unchanged. The measured 11th harmonic voltage distortion is shown in Figure 2. The 11th harmonic voltage in the substation Trige has exceeded the relevant IEC planning level [17], while in the substation Fraugde, it came close to but did not exceed the relevant planning level.

The total line length of the combined OHL and UGC connections between the substations Landerupgård and Malling is approximately 80 km, implying that the Vejle-Ådal 400 kV UGC is a relatively short section causing a system-wide harmonic amplification.

The topic of harmonic amplification becomes a real concern for further 400 kV cabling projects in Denmark, defining a goal for the development of measurement-validated methods for the simulation and prediction of harmonic voltage distortion in the country’s meshed transmission grids.

The simulations for the Vejle-Ådal 400 kV UGC project are compared to the harmonic voltage measurements in Section 6 as part of the model validation and prediction demonstration.

The energization of a new 400 kV OHL connection between the substation Kassø in Western Denmark (Southern Jutland, marked in Figure 1) and the substation Handewitt in Germany in July 2020 was followed by a small increase of the 11th harmonic voltage distortion in Eastern Jutland in the 400 kV substation Landerupgård, located about 40 km from Kassø, and the substation Trige, located about 120 km from Kassø. The harmonic distortion in the substation Kassø remained almost unaffected. At the same time, a former 220 kV OHL connection to Kassø was disconnected. At that time, the deterministic model for harmonic assessment described in this paper had become available and reached the same conclusions as had been made in Eastern Jutland, resulting in a small increase of the 11th harmonic voltage distortion in Eastern Jutland after enabling the 400 kV OHL connection between Kassø and Handewitt and removing the former 220 kV connection.

The measured 5th harmonic voltage distortion in the 150 kV substation Stovstrup in Western Jutland (see marking in Figure 1) used to exceed the relevant IEC planning level [7]. The harmonic assessment by simulation for the 150 kV grid reconstruction, replacing the OHL with a UGC between the substations Kassø and Ribe in Southern Jutland, predicted the reduction of the excessive 5th harmonic voltage distortion in Stovstrup to below the relevant IEC planning level. The substation Stovstrup is located about 90 km from the assessed 150 kV grid reconstruction project. The 5th harmonic voltage reduction in Stovstrup was confirmed by measurements after the first 150 kV overhead lines between Kassø and Ribe were removed and replaced with a UGC beginning in November 2021 [15]. The 5th harmonic voltage distortion in Stovstrup, measured between Nov. 2021 and March 2023, is still being reduced, while the 150 kV grid reconstruction project is ongoing. The latter project shall be completed at the end of 2023, when new harmonic voltage measurements will be acquired for benchmarking the simulation model.

The Westcoast 400 kV connection from the substation Idomlund in Western Jutland to the German border via the substation Endrup shall be completed in 2025 and will deploy the combined OHL and UGC line sections. The harmonic assessment conducted by using simulations has shown a significant increase of the 11th and 13th harmonic voltage distortions in the 400 kV substations Landerupgård, Malling, and Trige in Eastern Jutland and in the 400 kV substations Kingstrup and Fraugde in Funen (see marking in Figure 1). The 11th harmonic distortion will become extremely high in the abovementioned substations because the predicted increase occurs already from the high magnitudes of the 11th harmonic voltage distortion in the present grid after the energization of the Vejle-Ådal 400 kV UGC in July 2017.

The assessment has also revealed that the establishment of the 400/220 kV transformers, in the substation Idomlund, for the grid-connection of the Thor offshore wind power plant (about 1 GW nominal power) will result in a violation of the 13th harmonic planning level in the 400 kV substation Stovstrup (yet to be commissioned as part of the 400 kV Westcoast connection) in certain (n − 1) operation conditions of the 400 kV grid. The explanation of the 13th harmonic voltage increase will be that the harmonic self-impedance will increase in the 400 kV substation Stovstrup when connecting the 400/220 kV transformers, creating a resonance of the harmonic self-impedance between the transformers and the UGC capacitance of the Westcoast 400 kV connection.

For the Westcoast 400 kV connection, the mitigation of excessive 11th and 13th harmonic voltage distortions will be needed. The proposed mitigation shall include a more extensive utilization of the existing harmonic filters of the HVDC Converter Stations in Fraugde and V. Hassing and the establishment of a new harmonic filter in the 400 kV substation Stovstrup. The total mitigation solution, including the new harmonic filter design and specification to the vendor, has been prepared and evaluated by simulations using the presented deterministic simulation model for harmonic assessment.

The harmonic assessment for the Revsing–Landerupgård 400 kV connection, which shall be established by 2027 with a line length of approx. 25 km, has shown that the full UGC connection will cause a small increase in the harmonic voltage distortion. The mitigation to be established before 2025 for the Westcoast 400 kV connection shall already be sufficient for the Revsing–Landerupgård 400 kV UGC.

When considering harmonic propagation, the Revsing–Landerupgård 400 kV UGC connection is fundamentally different from the Vejle-Ådal 400 kV UGC section. Both the engineering evaluation of the harmonic voltage measurements and the results of the deterministic simulation model show that the Western Danish 400 kV transmission grid is the two systems with uncorrelated harmonic distortion: Eastern Jutland with initially high 11th harmonic distortion, and Western Jutland with low harmonic distortion (see markings in Figure 1). This observation will be confirmed by a statistical approach, as described in Section 8.

The establishment of the Revsing–Landerupgård 400 kV UGC connection will result in harmonic propagation between these two systems of the Western Danish 400 kV grid, with the dampening of the 11th harmonic distortion in Eastern Jutland. Further analysis has shown that the harmonic impedance imposed on the Eastern 400 kV system by the 400 kV UGC results in the strong propagation and dampening of the excessive 11th harmonic distortion in Eastern Jutland.

The results presented above emphasize the necessity of conducting harmonic assessment before the projection of new 400 kV connections with UGCs, including both the direct simulation of the harmonic voltage distortion and the frequency sweeps of the harmonic impedance by applying measurement-validated simulation models. Thus, measurement-validated simulation models for the harmonic assessment of the relevant meshed transmission grids become a necessity in in the grid planning and project design.

The Danish experience confirms that the direct simulation of the harmonic voltage distortion in the meshed transmission grid is a complex task requiring the accurate representations of both the harmonic impedance characteristics and the harmonic emission sources. The harmonic impedance characteristics are simulated using the tried and proven method of frequency sweeps. These frequency sweeps will identify the harmonic resonances, or just harmonic orders with large impedance magnitudes, which represent the ability of the grid to increase the harmonic voltage distortion in coincidence with the harmonic emission sources. Amplification of the harmonic voltage distortion will occur if there is background harmonic distortion (harmonic emission sources) exciting the resonances. If there are no such background distortion (i.e., no sources of specific harmonic orders coinciding the resonances), then there is nothing to amplify and there will be no excessive harmonic voltages. Therefore, calculations of the harmonic impedance characteristics using the frequency sweeps will only emphasize a risk of harmonic amplification for specific harmonic orders.

The frequency sweeps cannot conclude whether, or by how much, the harmonic voltage magnitude violates the planning levels, since the method does not represent harmonic emission sources. However, the frequency sweeps are useful for providing explanations of the results of direct simulation of the harmonic voltage distortion, e.g., changed harmonic distortion due to harmonic impedance shifts as a result of the establishment of a new component.

Energinet requires the delivery of harmonic source models of plants connected to the transmission grid to receive grid-connection permits in Denmark. The harmonic source models include magnitudes of the harmonic current or voltage versus the power or current of the plants in the points-of-connection to the grid. Usually, such characteristics are nonlinear, and this has also been reported in [1,2,4,5,7].

Factual harmonic distortion in the meshed transmission grid will be the result of the superposition of the harmonic source vectors in the different substations and will depend on the passive grid operation conditions. Recognizing the significance of the superposition of the harmonic vectors for the harmonic distortion, both amplitudes and phase-angles shall be adequately represented in the harmonic source vectors for the harmonic assessment of the meshed transmission grid with multiple harmonic sources [14,15].

The passive transmission grid can be in different operation conditions, implying that some connections, transformers, shunts, and harmonic filters are out of service. Different operation conditions in the present grid, or the establishment of new lines, transformers, shunts, or harmonic filters, will change the harmonic self-impedance in the substations (points of evaluation—POEs) and the harmonic coupling impedance between the substations (i.e., between POEs and sources) in the meshed transmission grid. The harmonic voltage distortion is a vectoral superposition of the harmonic impedance and harmonic emission sources. The worst-case harmonic distortion in a point-of-evaluation (POE) within the meshed transmission grid will not necessarily occur at simultaneous maximal magnitudes of the harmonic emission sources. Although the harmonic emission magnitudes of the individual sources can be high, the resulting harmonic distortion in the POE becomes low if the background distortion propagates from the harmonic emission sources to the POE with opposite phase-angles and each reduces the contributions of the other. The propagation depends on the harmonic coupling impedance between the different substations, which depends on the operation conditions of the meshed transmission grid. In principle, the harmonic voltage distortion may increase in one part and reduce in another part of the meshed transmission grid due to changes of the grid operation conditions or the establishment of new components, although the harmonic emission sources do not change.

The goal is that the methods and simulation model shall accurately predict the harmonic voltage magnitudes in the different substations for different operation conditions of the present transmission grid, for being able to predict the harmonic voltage distortion, and evaluate whether and which mitigation is needed in a future grid development stage.

3. Deterministic Method of Harmonic Distortion Simulation

The analytical modelling methods are developed, described, and successfully validated for radial connections with a single or a distributed harmonic emission source in one end of, and a POE in the other end of, the radial connection [2,3,4,5,17]. However, the transmission grid is meshed and includes multiple harmonic emission sources. To our knowledge, there is no publicly available, proven-by-harmonic-voltage-measurements, analytical method for the simulation of harmonic voltage distortion in meshed transmission grids with multiple harmonic emission sources. This section presents a deterministic method used by Energinet for the preparation of the simulation model for the harmonic assessment of the meshed transmission grid in Denmark. This method is deterministic because it was empirically developed using the harmonic voltage measurements in the transmission grid.

The term “preparation of the simulation model” entails the definition of the harmonic emission sources, when the electro-geometrical and electrical data of the passive transmission grid are known, for being able to simulate the harmonic voltage distortion in the meshed transmission grid in different grid operation regimes.

First, a theoretical background of the method including the invariant characteristics of the harmonic vectors will be given with reference to the present grid development stage. Second, how the numerical tuning of the harmonic vectors is conducted in practice will be explained in the necessary level of detail using LCC HVDC representations as examples. Third, the application of the deterministic method for the prediction of harmonic voltage distortion in a future grid development stage, i.e., with modified passive grid, and the handling of the model uncertainties will be explained.

3.1. Theoretical Background

In the presented deterministic method, the harmonic emission sources are either Norton equivalents with harmonic current vectors, J_Nk, or Thevenin equivalents with harmonic voltage vectors, E_Nk, where N is the source identifier and k is the harmonic order. To maintain the simplicity of the presentation of the method, the harmonic emission sources are denoted as harmonic current vectors, J_Nk. In the presented method, the harmonic current vectors do not explicitly include impedances. The approach with no explicit impedance of the harmonic current vectors is applied due to the following reasons:

• Validated data of the harmonic impedance of the HVDC Converter Stations as well as of the impedance of many distributed loads and generation units are not readily available [7].
• Specifically, the harmonic impedance of the LCC HVDC Converter Stations with nonlinear dependencies on both operation regime and control of the Converter Station itself and the operation regime of the grid is extremely difficult to obtain [7].
• The application of some kind of generic-level, not-validated data instead of unavailable vendor-specific or equipment-specific data would introduce another kind of inaccuracy to the simulation model. Such introduced inaccuracy would have unpredictable, uncontrollable influence on the results and shall hence be avoided.

Instead, the harmonic current vectors, with no explicitly defined impedance, are established in parallel with the passive components under the 150/60 kV distribution transformers for representing distributed harmonic emission sources. In the applied setup, the distributed loads, generation units, and 60 kV network equivalents already include their passive harmonic characteristics [20]. Such characteristics serve implicitly as harmonic impedance of the current vectors; this is illustrated in Figure 3.

The harmonic current vectors of the HVDC Converter Stations are established in parallel with the harmonic filters of the HVDC Converter Stations. The harmonic current magnitudes are measured and provided by the vendors at the AC terminals of the physical Converter Stations, i.e., the current injected into the grid-connection substation with the already deducted contribution absorbed in the “internal” converter impedance. The HVDC Converter Stations are LCC-type and do not go in operation without a minimum required number of the harmonic filters that have been connected in advance. Thus, the harmonic current vectors of these relatively large, centralized, harmonic emission sources will always have some harmonic filters in operation which are connected in parallel. In the model, the harmonic filters will serve as impedance in parallel with the harmonic current injected by the HVDC Converter Station model and help against numerical overexcitation, such as in cases with high harmonic voltage magnitudes. The harmonic filters will implicitly serve as harmonic impedance of the equivalents representing the HVDC Converter Stations in the applied model setup, as illustrated in Figure 4.

Further, numerical overexcitation can and shall be prevented through careful and proper tuning of the harmonic current vectors under the considerations described above. The sufficiency of the applied numerical tuning will be shown in the validation cases with measured high magnitudes of the harmonic voltage distortion in Section 5.

The harmonic current vectors, J_Nk, propagate through the meshed grid represented by the harmonic impedance matrix [Z_NMk] and induce the harmonic voltages in the substations U_Mk where M is the substation identifier. The model equations become:

$$\left[\begin{array}{c}{U}_{1k}\\ \begin{array}{c}{U}_{2k}\\ \begin{array}{c}⋮\\ U_{Mk}\end{array}\end{array}\end{array}\right]=\left[\begin{array}{cc}\begin{array}{cc}{Z}_{11k}& Z_{12k}\\ Z_{21k}& Z_{22k}\end{array}& \begin{array}{cc}\dots & Z_{1Nk}\\ \cdots & Z_{2Nk}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ Z_{M1k}& Z_{M2k}\end{array}& \begin{array}{cc}\ddots & ⋮\\ \dots & Z_{MNk}\end{array}\end{array}\right]·\left[\begin{array}{c}{J}_{1k}\\ \begin{array}{c}{J}_{2k}\\ \begin{array}{c}⋮\\ J_{Nk}\end{array}\end{array}\end{array}\right],$$

(1)

where the kth harmonic impedances Z_MNk are in complex numbers, with the two equal indexes, i.e., M = N, representing the self-impedances in the substations M, and otherwise, i.e., M ≠ N, representing the coupling impedances between the substations M and N. The harmonic voltages U_Mk and harmonic currents J_Nk are vectors, i.e., complex numbers with magnitudes and phase-angles. The harmonic impedance matrix [Z_MNk] of Equation (1) is known from the passive part of the transmission grid model with the electro-geometrical data of the OHLs and the UGCs, the electrical data of the transformers, shunts and harmonic filters, and harmonic representations of passive components in the distribution grids such as loads. Equation (1) applies for each harmonic order, k, represented in the simulation model, meaning that there are k vectoral equations as Equation (1) to be prepared and solved.

The main properties of the harmonic emission sources, J_Nk, are:

• Their substations in the grid, denoted by the index N;
• The harmonic orders, denoted by the index k;
• The harmonic magnitudes and harmonic phase-angles;
• Asymmetry representing possible unbalance of the harmonic phase-current magnitudes.

The hypothesis of the presented deterministic method is that the harmonic current vectors can be represented using the relative components with reference to the fundamental magnitudes and phase-angles which are invariant characteristics of the harmonic current vectors with regard to the operation regimes of the passive meshed transmission grid. First, the relative components of the harmonic source vectors will be defined. Then, it will be explained how the numerical tuning of the relative components of the harmonic source vectors is conducted with reference to Equation (1). Asymmetry will be omitted in the mathematical presentation below.

3.2. Relative Components of Harmonic Source Vectors

The harmonic source vectors, J_N, in symmetrical components are [21]:

$$\begin{gathered} J_{NA1}=a_{N1}·e^{j·({\omega }_{1}t+{\phi }_{N1})}, J_{NB1}=a_{N1}·e^{j·({\omega }_{1}t-\frac{2}{3}\pi +{\phi }_{N1})}, J_{NC1}=a_{N1}·e^{j·({\omega }_{1}t+\frac{2}{3}\pi +{\phi }_{N1})}, \\ {J}_{NAk}=a_{Nk}·e^{j·(k·{\omega }_{1}t+{\phi }_{Nk})}, J_{NBk}=a_{Nk}·e^{j·(k·{\omega }_{1}t-\frac{2}{3}\pi +{\phi }_{Nk})}, J_{NCk}=a_{Nk}·e^{j·(k·{\omega }_{1}t+\frac{2}{3}\pi +{\phi }_{Nk})}, \end{gathered}$$

(2)

where ω₁ = 2·π∙f₁ is the fundamental angular frequency and f₁ is the fundamental frequency, the indexes A, B, and C refer to the three phases, the index 1 refers to the fundamental frequency, the index k means the kth harmonic order, a_NK represents the kth harmonic magnitudes, φ_Nk is the kth harmonic phase-angle of the harmonic source vectors in the substation N. The kth harmonic source vectors by Equation (2) can be rewritten thus:

$$\begin{gathered} J_{NAk}=b_{NAk}·\left(\frac{a_{Nk}}{a_R}\right)·e^{j·\left({\phi }_{Nk}-k·{\phi }_{N1}\right)}, J_{NBk}=b_{NBk}·\left(\frac{a_{Nk}}{a_R}\right)·e^{j·\left({\phi }_{Nk}+\frac{2·\left(k-1\right)}{3}\pi -k·{\phi }_{N1}\right)}, \\ {J}_{NCk}=b_{NCk}·\left(\frac{a_{Nk}}{a_R}\right)·e^{j·\left({\phi }_{Nk}-\frac{2·\left(k-1\right)}{3}\pi -k·{\phi }_{N1}\right)}, \end{gathered}$$

(3)

with (a_Nk/a_R) and (φ_Nk − k∙φ_N1) = α_Nak, (φ_Nk + ⅔π∙(k − 1) − k∙φ_N1) = α_NBk, and (φ_Nk − ⅔π∙(k − 1) − k∙φ_N1) = α_NCk expressing the relative magnitudes and relative phase-angles of the kth harmonic order with reference to the fundamental magnitudes and phase-angles in the substation N, and where a_R is the rated current.

In the presented method, b_Nak, b_NBk, and b_NCk, denote the components of the harmonic source vectors which are dependent from the grid operation conditions and so dependent from the fundamental load-flow solution. The remaining components in Equation (3) contain the relative magnitudes and relative phase-angles. In the presented method, the relative components (a_Nk/a_R), α_Nak, α_NBk, and α_NCk, are considered as being the harmonic characteristics of the sources and independent from the grid operation conditions.

Applying Equations (2) and (3), the harmonic source vectors in symmetrical components can be written thus:

$$\begin{gathered} \left[\begin{array}{c}{J}_{NAk}\\ J_{NBk}\\ J_{NCk}\end{array}\right]=\left[\begin{array}{ccc}{b}_{NAk}& 0& 0\\ 0& b_{NBk}& 0\\ 0& 0& b_{NCk}\end{array}\right]·\left(\frac{a_{Nk}}{a_R}\right)·\left[\begin{array}{c}{e}^{j·\left({\phi }_{Nk}-k·{\phi }_{N1}\right)}\\ e^{j·\left({\phi }_{Nk}-\frac{2}{3}\pi -k·{\phi }_{N1}\right)}\\ e^{j·\left({\phi }_{Nk}+\frac{2}{3}\pi -k·{\phi }_{N1}\right)}\end{array}\right]k=3, 6, 9, 12\dots \\ \left[\begin{array}{c}{J}_{NAk}\\ J_{NBk}\\ J_{NCk}\end{array}\right]=\left[\begin{array}{ccc}{b}_{NAk}& 0& 0\\ 0& b_{NBk}& 0\\ 0& 0& b_{NCk}\end{array}\right]·\left(\frac{a_{Nk}}{a_R}\right)·\left[\begin{array}{c}{e}^{j·\left({\phi }_{Nk}-k·{\phi }_{N1}\right)}\\ e^{j·\left({\phi }_{Nk}+\frac{2}{3}\pi -k·{\phi }_{N1}\right)}\\ e^{j·\left({\phi }_{Nk}-\frac{2}{3}\pi -k·{\phi }_{N1}\right)}\end{array}\right]k=2, 5, 8, 11\dots \\ \left[\begin{array}{c}{J}_{NAk}\\ J_{NBk}\\ J_{NCk}\end{array}\right]=\left[\begin{array}{ccc}{b}_{NAk}& 0& 0\\ 0& b_{NBk}& 0\\ 0& 0& b_{NCk}\end{array}\right]·\left(\frac{a_{Nk}}{a_R}\right)·\left[\begin{array}{c}{e}^{j·\left({\phi }_{Nk}-k·{\phi }_{N1}\right)}\\ e^{j·\left({\phi }_{Nk}-k·{\phi }_{N1}\right)}\\ e^{j·\left({\phi }_{Nk}-k·{\phi }_{N1}\right)}\end{array}\right]k=1, 4, 7, 10, 13\dots \end{gathered}$$

(4)

Considering the harmonic orders characteristic for the Danish 400 kV transmission grid, the 3rd harmonic vector is the zero-sequence rotation, the 5th, 11th, and 23rd harmonic vectors constitute the negative-sequence rotation, and the fundamental vector and the 7th, 13th, and 25th harmonic vectors constitute the positive-sequence rotation.

The relative phase-angles of the kth harmonic vectors in the substation N are:

$$\begin{gathered} {\alpha }_{NAk}=\left({\phi }_{Nk}-k·{\phi }_{N1}\right), \mathrm{for} \mathrm{all} \mathrm{harmonic} \mathrm{orders} k \\ {\alpha }_{NBk}={\alpha }_{NAk}+\{\begin{array}{c}-\frac{2}{3}\pi , k=3, 6, 9, 12\dots \\ \frac{2}{3}\pi , k=2, 5, 8, 11\dots \\ 0, k=4, 7, 10, 13\dots \end{array} \\ {\alpha }_{NCk}={\alpha }_{NAk}+\{\begin{array}{c}\frac{2}{3}\pi , k=3, 6, 9, 12\dots \\ -\frac{2}{3}\pi , k=2, 5, 8, 11\dots \\ 0, k=4, 7, 10, 13\dots \end{array} \end{gathered}$$

(5)

In the presented method, the relative magnitudes (a_Nk/a_R) and the relative phase-angles α_Nak, α_NBk, and α_NCk, presented above, are applied for the solving of Equation (1).

3.3. Numerical Solution of Harmonic Source Vectors

The goal of the numerical tuning is unambiguously to define the invariant characteristics of the harmonic current vectors so that:

• High magnitudes of the harmonic voltage distortion, i.e., those approaching and violating the planning levels, will be accurately simulated by the model;
• Low magnitudes, i.e., those significantly lower than the planning levels, will be accurately predicted as low but some larger relative discrepancy between the model and the harmonic voltage measurements can be accepted;
• The accuracy defined above shall be achieved for all harmonic orders included in the numerical tuning, all 400 kV substations of the Danish transmission grid, and all validation cases;
• High magnitudes imply the exceeding of 70% of the IEC planning level [17], the definition of which will be explained in Section 3.4.

The presented deterministic method works with the two types of harmonic emission sources:

• Vendor-specific sources;
• Distributed harmonic sources.

The available data of the vendor-specific sources, such as of the HVDC Converter Stations depicted in Figure 1, include specific grid-connection substations in the grid and worst-case magnitudes of the harmonic orders as nonlinear dependencies of the power transport of the HVDC Converter Station. The harmonic phase-angles are usually not available from the provided harmonic data of the vendors.

For the distributed harmonic sources, the harmonic magnitudes or phase-angles are not readily available. Hence the locations in the grid and harmonic orders are only identifiable using the harmonic voltage measurements. Thus, the harmonic magnitudes and phase-angles of the distributed harmonic sources J_1k…J_Nk in Equation (1) are not known in detail, which should be necessary for the conduction of straightforward harmonic assessment.

In the Danish transmission grid, the measured harmonic voltage magnitudes, |U_Mk|, are available for the majority (but not necessarily for all) of the 400 kV substations year-round as ten-minute average phase-to-ground voltage magnitudes (|U_Mak|, |U_MBk|, and |U_MCk|). The harmonic phase-angles are not included in the harmonic voltage measurements.

From the above discussion, it becomes obvious that the number of unknowns, meaning the unknown magnitudes and phase-angles of the harmonic sources and phase-angles of the harmonic voltages, exceeds the number of linear equations in Equation (1). The analytical unambiguous solving of Equation (1) for finding the harmonic source vectors J_1k…J_Nk does not seem possible.

Instead of seeking analytical solution of Equation (1), the deterministic method has been developed to find the numerical solution for the relative magnitudes and relative phase-angles of the harmonic source vectors. The target of the numerical tuning is that the insertion of the relative magnitudes (a_Nk/a_R) and relative phase-angles α_Nak, α_NBk and α_NCk in the models of the harmonic source vectors, as in Equation (5), shall result in the maximum of the simulated harmonic phase-to-ground voltage magnitudes reaching the maximum of the measured ten-minute average harmonic phase-to-ground voltage magnitudes, with smallest possible discrepancy, for all M substations with the available harmonic voltage measurements, for all (characteristic) harmonic orders k, for all grid operation regimes.

The numerical tuning includes several grid operation regimes. The grid operation regimes are combinations of various power transports and harmonic filters in service of the HVDC Converter Stations (and, generally, all harmonic source models utilizing vendor-specific data) and with various (n − L) conditions, L ≥ 0, referring to the transmission lines and power transformers that are out of service. The inclusion of the various grid operation regimes with various (n − L) conditions introduces variations in the harmonic impedance matrix [Z_MNk]. These variations in the [Z_MNk] matrix correspond to the different but specific measured magnitudes of the harmonic voltages in the substations, i.e., |U_Mak|, |U_MBk|, and |U_MCk|. The inclusion of the different power transports, supply, or consumption of the harmonic sources with vendor-specific data introduces the relative magnitudes and relative phase-angles of the harmonic sources in Equation (1), which are tightened up to suit these different conditions.

It is important to include the periods with several consecutive snapshots so that not only stationary harmonic voltage magnitudes, but also changes in the harmonic voltage magnitudes in response to changes of the passive grid operation conditions [Z_MNk]_p1 → [Z_MNk]_p2 are present in the numerical tuning of the harmonic source vectors and later in the validation of the simulation model as a whole.

The numerical tuning will define an (empirically developed) solution set of the magnitudes, (a_Nk/a_R), and phase-angles, α_Nak, α_NBk, and α_NCk, of the relative harmonic vectors so that the application of those magnitudes (a_Nk/a_R) and phase-angles α_Nak, α_NBk, and α_NCk results in the simulated harmonic voltage magnitudes converging with the measured harmonic voltage magnitudes, for all included in the tuning grid operation regimes p = 1…P, for all substations, M, and for all characteristic harmonic orders, k.

The various grid operation regimes are interpreted as the equation system, depicted thus:

$$\{\begin{array}{c}{\left[\begin{array}{c}{U}_{1k}\\ \begin{array}{c}{U}_{2k}\\ \begin{array}{c}⋮\\ U_{Mk}\end{array}\end{array}\end{array}\right]}_{p=1}={\left[\begin{array}{cc}\begin{array}{cc}{Z}_{11k}& Z_{12k}\\ Z_{21k}& Z_{22k}\end{array}& \begin{array}{cc}\dots & Z_{1Nk}\\ \cdots & Z_{2Nk}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ Z_{M1k}& Z_{M2k}\end{array}& \begin{array}{cc}\ddots & ⋮\\ \dots & Z_{MNk}\end{array}\end{array}\right]}_{p=1}·{\left[\begin{array}{c}{J}_{1k}\\ \begin{array}{c}{J}_{2k}\\ \begin{array}{c}⋮\\ J_{Nk}\end{array}\end{array}\end{array}\right]}_{p=1},\\ {\left[\begin{array}{c}{U}_{1k}\\ \begin{array}{c}{U}_{2k}\\ \begin{array}{c}⋮\\ U_{Mk}\end{array}\end{array}\end{array}\right]}_{p=2}={\left[\begin{array}{cc}\begin{array}{cc}{Z}_{11k}& Z_{12k}\\ Z_{21k}& Z_{22k}\end{array}& \begin{array}{cc}\dots & Z_{1Nk}\\ \cdots & Z_{2Nk}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ Z_{M1k}& Z_{M2k}\end{array}& \begin{array}{cc}\ddots & ⋮\\ \dots & Z_{MNk}\end{array}\end{array}\right]}_{p=2}·{\left[\begin{array}{c}{J}_{1k}\\ \begin{array}{c}{J}_{2k}\\ \begin{array}{c}⋮\\ J_{Nk}\end{array}\end{array}\end{array}\right]}_{p=2},\\ ⋮\\ {\left[\begin{array}{c}{U}_{1k}\\ \begin{array}{c}{U}_{2k}\\ \begin{array}{c}⋮\\ U_{Mk}\end{array}\end{array}\end{array}\right]}_{p=P}={\left[\begin{array}{cc}\begin{array}{cc}{Z}_{11k}& Z_{12k}\\ Z_{21k}& Z_{22k}\end{array}& \begin{array}{cc}\dots & Z_{1Nk}\\ \cdots & Z_{2Nk}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ Z_{M1k}& Z_{M2k}\end{array}& \begin{array}{cc}\ddots & ⋮\\ \dots & Z_{MNk}\end{array}\end{array}\right]}_{p=P}·{\left[\begin{array}{c}{J}_{1k}\\ \begin{array}{c}{J}_{2k}\\ \begin{array}{c}⋮\\ J_{Nk}\end{array}\end{array}\end{array}\right]}_{p=P},\end{array}$$

(6)

where the index p denotes the included grid operation regimes from 1 to P, the index k is the harmonic order, and the harmonic source vectors from 1 to N are presented in symmetrical components according to Equations (4) and (5):

$${\left[\begin{array}{c}{J}_{1k}\\ J_{2k}\\ ⋮\\ J_{Nk}\end{array}\right]}_p<=\left\{\begin{array}{c}{\left[\begin{array}{c}{J}_{1Ak}\\ J_{2Ak}\\ ⋮\\ J_{NAk}\end{array}\right]}_p={\left[\begin{array}{cccc}{b}_{1Ak}& 0& \cdots & 0\\ 0& b_{2Ak}& \cdots & ⋮\\ 0& ⋮& \ddots & 0\\ 0& 0& \cdots & b_{NAk}\end{array}\right]}_p·\left(\frac{a_{Nk}}{a_R}\right)·\left[\begin{array}{c}{e}^{j·{\alpha }_{1Ak}}\\ e^{j·{\alpha }_{2Ak}}\\ ⋮\\ e^{j·{\alpha }_{NAk}}\end{array}\right]\\ {\left[\begin{array}{c}1Bk\\ J_{2Bk}\\ ⋮\\ J_{NBk}\end{array}\right]}_p={\left[\begin{array}{cccc}{b}_{1Bk}& 0& \cdots & 0\\ 0& b_{2Bk}& \cdots & ⋮\\ 0& ⋮& \ddots & 0\\ 0& 0& \cdots & b_{NBk}\end{array}\right]}_p·\left(\frac{a_{Nk}}{a_R}\right)·\left[\begin{array}{c}{e}^{j·{\alpha }_{1Bk}}\\ e^{j·{\alpha }_{2Bk}}\\ ⋮\\ e^{j·{\alpha }_{NBk}}\end{array}\right],\\ {\left[\begin{array}{c}{J}_{1Ck}\\ J_{2Ck}\\ ⋮\\ J_{NCk}\end{array}\right]}_p={\left[\begin{array}{cccc}{b}_{1Ck}& 0& \cdots & 0\\ 0& b_{2Ck}& \cdots & ⋮\\ 0& ⋮& \ddots & 0\\ 0& 0& \cdots & b_{NCk}\end{array}\right]}_p·\left(\frac{a_{Nk}}{a_R}\right)·\left[\begin{array}{c}{e}^{j·{\alpha }_{1Ck}}\\ e^{j·{\alpha }_{2Ck}}\\ ⋮\\ e^{j·{\alpha }_{NCk}}\end{array}\right]\end{array}\right.$$

(7)

where the indices p = 1…P mark the parts of the numerical solution of the harmonic source vectors J_1k…J_Nk, being dependent from the grid operation regime p. This part of the solution includes the fundamental load-flow solution for the grid operation regime p, which is necessary to obtain before acquiring the harmonic solution and conducting the numerical tuning of the relative magnitudes and relative phase-angles.

The part of the numerical solution in Equation (7) without such indexes p includes the relative magnitudes (a_Nk/a_R) and the relative phase-angles α_Nak, α_NBk, and α_NCk (in symmetrical components), which are independent from the grid operation regimes and keep the same values throughout all operation regimes p = 1…P. This part of the solution will provide the input parameters to the models of the harmonic source vectors, J_1k…J_Nk, which does not depend on the fundamental load-flow solution.

Accepting the hypothesis that the relative magnitudes (a_Nk/a_R) and the relative phase-angles α_Nak, α_NBk, and α_NCk are invariants of the grid operation regimes, the number of unknowns, such as the relative magnitudes and phase-angles of the harmonic source vectors, becomes fixed and does not increase with the introduction of more grid operation regimes for numerically solving Equation (6).

On contrary, the harmonic source vectors, J_1k…J_Nk are not necessarily independent from the grid operation regimes. If the target was seeking the solution as the harmonic source vectors, J_1k…J_Nk, then the inclusion of more grid operation regimes into Equation (6) would also increase the number of unknowns. Thus, the inclusion of more grid operation regimes would not necessarily help in the numerical solving of the equation system of Equation (6) with the “whole” harmonic source vectors, J_1k…J_Nk, as unknowns.

The numerical tuning process, accounting for the number of equations and unknowns, is briefly described below, for the case of the HVDC Converter Stations.

Energinet applies the DIgSILENT PowerFactory^® software for conducting grid studies including harmonic assessment in means of the unbalanced harmonic load-flow. For the numerical tuning, the simulation model with the already-established and measurement-verified passive grid representation is applied. This means that the harmonic impedance matrix, according to Equation (1), of the interconnected HVAC grids of Western Denmark and Germany is established within the simulation model. The 400 kV part of the transmission grid model is illustrated in Figure 5. There will be a distinguishment between the passive grid operation regimes and those of the HVDC Converter Stations.

By bringing the passive components such as the 400 kV transmission lines and 400/150 kV transformers out-of-service, the (n − L) grid operation regimes are established in the simulation model. The different grid operation regimes correspond to modifications of the harmonic impedance matrix, and these operations are handled internally in the simulation model. In other words, the impedance matrix remains within the simulation software and all modifications of the impedance matrix are handled automatically.

The models of the HVDC Converter Stations and the Anholt Offshore Wind Power Plant (OWPP) are implemented using the AC current source models of the simulation software. The operation of the HVDC Converter Stations and Anholt OWPP is divided into the four regimes by 25% steps. This division with the steps by 25% is sufficient for the harmonic assessment conducted by Energinet. The exception is for the Skagerrak 1,2 HVDC Converter Station, which is represented with full transport. The Skagerrak 1,2 is connected to the 150 kV grid, imposing less influence on the harmonic voltage distortion in the 400 kV grid, and this justifies the simplification. For each operation regime, the HVDC Converter Station and Anholt OWPP will inject the required amount of the nominal-frequency active and reactive power. The HVDC Converter Stations will also be simulated with the number of the harmonic filters, which is measured at the time of the operation regime (for validation at the present grid stage) or required by the operation guideline (for assessment of a future grid development stage).

The harmonic characteristics are defined by the relative magnitudes (a_Nk/a_R) and relative phase-angles α_Nak, α_NBk, and α_NCk. The ratios (a_Nk/a_R) correspond to the harmonic current magnitudes in percent which are dependent from the power transport and found in the vendor-specific data for the HVDC Converter Stations. For the numerical tuning, the ratios (a_Nk/a_R) for the four operation regimes with 25% steps are applied, for each harmonic order and each HVDC Converter Station. The ratios (a_Nk/a_R) will not be adjusted during the numerical tuning because these ratios represent the nonlinearity of the harmonic emission at different operation regimes of specific HVDC Converter Stations. However, the current rating, a_R, for the harmonic emission will be part of the numerical tuning. The current ratings, a_R, for the harmonic emission:

• Can be different for the different HVDC Converter Stations;
• Shall be the same for the different operation regimes of the same HVDC Converter Station.

The conditions listed above are illustrated in Figure 5 by only assigning a single unknown parameter a_R for each HVDC Converter Station and Anholt OWPP.

The relative phase-angle α_Nak is an independent parameter, while the other two relative phase-angles α_NBk and α_NCk can be derived using α_Nak according to Equation (5). Thus, the relative phase-angle α_Nak is only the unknown to be found by using the numerical tuning.

The different operation regimes of the HVDC Converter Stations and Anholt OWPP will be characterized by the different relative phase-angles, α_Nak, which are denoted in Figure 5 in the list of unknowns.

In addition, the harmonic voltage source representing the measured harmonic distortion in the German 400 kV transmission grid will be assigned with its magnitude and phase-angle.

The numerical tuning procedure is as follows:

• The grid operation regimes (the passive part) and the operation regimes of the HVDC Converter Stations and Anholt OWPP are implemented in the simulation model of the Western Danish transmission grid.
• The AC load flow solutions are conducted for each operation regime.
• The parameters a_R and α_Nak are assigned to the harmonic source models (HVDC, Anholt, and the harmonic voltage source in Germany) according to their operation regimes.
• The harmonic load flow solutions are conducted for each operation regime, and the maximum of the three line-to-ground harmonic voltage magnitudes are recorded for each harmonic order of the tuning and for each substation.
• The simulated harmonic voltage magnitudes are compared to the measured harmonic voltage magnitudes (maximum of the three line-to-ground harmonic voltages) for the 400 kV substations with available harmonic measurements.
• If the discrepancy is not acceptable, the parameters a_R and α_Nak are adjusted and the procedure for solving the harmonic load flow and comparing the simulated and measured harmonic voltage magnitudes continues.
• When the discrepancy between the simulated and measured harmonic voltage magnitudes becomes small and acceptable, the tuning stops. The relative magnitudes and relative phase-angles for the representation of the harmonic emission from the HVDC and Anholt OWPP and the contribution of the German transmission grid in Denmark are fixed and locked.

The parameters b_Nk are calculated by the simulation software for each operation regime, but these parameters are not applied in the numerical tuning described above. The number of unknowns a_R and α_Nak for the HVDC Converter Stations, Anholt OWWP, and the harmonic voltage source in Germany, becomes 24 (as will be confirmed by an inspection of Figure 5).

The number of different operation regimes (see the snapshots referring to Section 4 for detailed explanation of how the data for the regimes were acquired) applied for the numerical tuning and validation is P = 31. Each operation regime, p, corresponds to one equation system as in Equation (1) when expressed in complex numbers or two equation systems when converted to real numbers. The unknowns a_R and α_Nak are real numbers. The number of equations applied for the numerical tuning of the harmonic emission sources of the HVDC Converter Stations, Anholt OWPP, and the harmonic voltage source in Germany exceeds the number of unknowns. In practice, the number of equations should be greater than the number of unknowns for the inclusion of as-many-as-possible various combinations of the power transports (different active harmonic emission regimes) and the passive grid operation regimes with outage of lines and power transformers, as well as inclusion of the periods with high, varying between high and low, and low magnitudes of the harmonic voltage distortion in the different substations of the meshed transmission grid. Consequently, application of the equation system as in Equation (6) with several grid operation regimes (and consecutive snapshots) increases the number of equations with an unchanged number of the unknowns, (a_Nk/a_R), and α_Nak, α_NBk, and α_NCk. Therefore, the described tuning process improves the possibility of finding an unambiguous solution—the (empirically defined) solution set of the relative magnitudes (a_Nk/a_R) and relative phase-angles α_Nak, α_NBk, and α_NCk—and so prepares the simulation model for the harmonic assessment of the transmission grid with multiple harmonic emission sources [6].

When the numerical tuning is completed, the values of the relative magnitudes (a_Nk/a_R) and relative phase-angles α_Nak, α_NBk, and α_NCk to be inputted in the representations of the harmonic source vectors J_1k…J_Nk are locked [15]. The relative magnitudes and relative phase-angles will not be changed during simulations of the harmonic voltage distortion in the present grid development stage. It is expected that the harmonic voltage distortion in new grid operation regimes in the existing meshed transmission grid, which have not been part of the numerical tuning, shall also be accurately simulated using the prepared model; see Section 5.3.

3.4. Prediction of Harmonic Distortion following Grid Expansion

When tuned properly, the solution set of the relative magnitudes (a_Nk/a_R) and relative phase-angles α_Nak, α_NBk, and α_NCk of the harmonic source vectors shall accurately simulate the harmonic voltage distortion in the various (n − L) operation regimes of the present grid development stage. The simulation model shall simulate the consecutive snapshots with both (n − 0) → (n − L) for the disconnection and (n − L) → (n − 0) for the reconnection of the transmission lines and transformers of the present grid stage. When simulating harmonic distortion, such operation regimes and snapshots correspond to the modification of the harmonic impedance matrix [Z_MNk], albeit without changing the matrix topology and at unchanged values of the solution set with the relative magnitudes (a_Nk/a_R) and relative phase-angles α_Nak, α_NBk, and α_NCk of the harmonic source vectors.

Grid expansion and grid reconstruction imply that L₁ connections are removed and L₂ new connections are added to the present transmission grid. In terms of the simulations, grid expansion and grid reconstruction correspond to modifications of the harmonic impedance matrix [Z_MNk] where the indexes M and N are not necessarily identical to those before the grid expansion and grid reconstruction [15]. The removal of L₁ connections and addition of L₂ connections can formally be expressed as (n − L₁ + L₂), which represents the (n − 0) state of the new grid development stage. However, the total numbers of connections n in the present and in the new grid development stages are not necessarily identical. For avoidance of confusion with the commonly applied (n − 0) and (n − L) terminologies for the existing connections in or out of service, the (n − L₁ + L₂) terminology will not be applied for the expanded or reconstructed grid with the permanent removal of some connections and the commissioning of other new connections. Instead, the term “present” will cover the present grid development stage and the term “next” will cover the new grid development stage with all reasonable-to-include (n − L), L ≥ 0, operation regimes in both grid stages, respectively.

Considering that the relative magnitudes (a_Nk/a_R) and relative phase-angles α_Nak, α_NBk, and α_NCk of the harmonic source vectors in the present grid development stage remain unchanged in the next grid development stage, the harmonic voltage distortion in the next grid development stage can be predicted by simulations thus:

$${\left[\begin{array}{c}{U}_{1k}\\ \begin{array}{c}{U}_{2k}\\ \begin{array}{c}⋮\\ U_{Mk}\end{array}\end{array}\end{array}\right]}_{NEXT}={\left[\begin{array}{cc}\begin{array}{cc}{Z}_{11k}& Z_{12k}\\ Z_{21k}& Z_{22k}\end{array}& \begin{array}{cc}\dots & Z_{1Nk}\\ \cdots & Z_{2Nk}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ Z_{M1k}& Z_{M2k}\end{array}& \begin{array}{cc}\ddots & ⋮\\ \dots & Z_{MNk}\end{array}\end{array}\right]}_{NEXT}·{\left[\begin{array}{c}{J}_{1k}\\ \begin{array}{c}{J}_{2k}\\ \begin{array}{c}⋮\\ J_{Nk}\end{array}\end{array}\end{array}\right]}_{NEXT},$$

(8)

where the harmonic source vectors in the next grid expansion stage apply the relative magnitudes and relative phase-angles determined for the present grid expansion stage:

$${\left[\begin{array}{c}{J}_{1k}\\ J_{2k}\\ ⋮\\ J_{Nk}\end{array}\right]}_p<=\left\{\begin{array}{c}{\left[\begin{array}{c}{J}_{1Ak}\\ \begin{array}{c}{J}_{2Ak}\\ \begin{array}{c}⋮\\ J_{NAk}\end{array}\end{array}\end{array}\right]}_{NEXT}={\left[\begin{array}{c}{b}_{1Ak}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}\begin{array}{c}0\\ \begin{array}{c}{b}_{2Ak}\\ ⋮\\ 0\end{array}\end{array}\begin{array}{c}\cdots \\ \begin{array}{c}\cdots \\ \ddots \\ \cdots \end{array}\end{array}\begin{array}{c}0\\ \begin{array}{c}⋮\\ 0\\ b_{NAk}\end{array}\end{array}\right]}_{NEXT}·{\left(\frac{a_{Nk}}{a_R}\right)}_{PRESENT}·{\left[\begin{array}{c}{e}^{j·{\alpha }_{1Ak}}\\ e^{j·{\alpha }_{2Ak}}\\ ⋮\\ e^{j·{\alpha }_{NAk}}\end{array}\right]}_{PRESENT}\\ {\left[\begin{array}{c}{J}_{1Bk}\\ \begin{array}{c}{J}_{2Bk}\\ \begin{array}{c}⋮\\ J_{NBk}\end{array}\end{array}\end{array}\right]}_{NEXT}={\left[\begin{array}{c}{b}_{1Bk}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}\begin{array}{c}0\\ \begin{array}{c}{b}_{2Bk}\\ ⋮\\ 0\end{array}\end{array}\begin{array}{c}\cdots \\ \begin{array}{c}\cdots \\ \ddots \\ \cdots \end{array}\end{array}\begin{array}{c}0\\ \begin{array}{c}⋮\\ 0\\ b_{NBk}\end{array}\end{array}\right]}_{NEXT}·{\left(\frac{a_{Nk}}{a_R}\right)}_{PRESENT}·{\left[\begin{array}{c}{e}^{j·{\alpha }_{1Bk}}\\ e^{j·{\alpha }_{2Bk}}\\ ⋮\\ e^{j·{\alpha }_{NBk}}\end{array}\right]}_{PRESENT}\\ {\left[\begin{array}{c}{J}_{1Ck}\\ \begin{array}{c}{J}_{2Ck}\\ \begin{array}{c}⋮\\ J_{NCk}\end{array}\end{array}\end{array}\right]}_{NEXT}={\left[\begin{array}{c}{b}_{1Ck}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}\begin{array}{c}0\\ \begin{array}{c}{b}_{2Ck}\\ ⋮\\ 0\end{array}\end{array}\begin{array}{c}\cdots \\ \begin{array}{c}\cdots \\ \ddots \\ \cdots \end{array}\end{array}\begin{array}{c}0\\ \begin{array}{c}⋮\\ 0\\ b_{NCk}\end{array}\end{array}\right]}_{NEXT}·{\left(\frac{a_{Nk}}{a_R}\right)}_{PRESENT}·{\left[\begin{array}{c}{e}^{j·{\alpha }_{1Ck}}\\ e^{j·{\alpha }_{2Ck}}\\ ⋮\\ e^{j·{\alpha }_{NCk}}\end{array}\right]}_{PRESENT}\end{array}\right.$$

(9)

The presented approach becomes useful for the direct simulation of the harmonic voltage distortion in a future grid development stage, where new connections will either replace obsolete connections (i.e., reconstruction) or be added (i.e., expansion) to the present grid. However, the presented method and simulation results will be subject to uncertainty and risks due to the following factors:

• Electro-geometrical and electrical data of new, not-yet-commissioned connections which are available at the time of the harmonic assessment may differ from the as-built data after commissioning;
• New harmonic emission sources with yet-unknown characteristics can be introduced in the transmission grid during or after the commissioning of the new connections;
• Initial inaccuracies in the numerical tuning of the relative magnitudes or relative phase-angles;
• Changes within foreign neighboring grids which are not properly addressed in the simulation model of the domestic transmission grid.

The first uncertainty relates to the risk of missing harmonic impedance peaks due to data inaccuracy and can be noticeably reduced by including parameter variations in the available (either typified or as-planned) data of the new connections when conducting the harmonic assessment. The parameter variations can be settled using previous experience from the commissioning of similar connection types. For example, the capacitances can be varied within a ±10% range of the typified data of not-yet-commissioned UGCs.

The second and third uncertainties are rather hard to mitigate by simple means. However, such uncertainties can be reduced by including many various operation regimes in the harmonic assessment, including unusual regimes with less-than-normally harmonic filters in-service, and looking after patterns causing excessive harmonic voltage magnitudes in the simulation results. Further, tendencies, changes, and magnitudes of the harmonic voltages will be included in evaluations of whether the harmonic voltage distortion will become excessive and need mitigation in the future grid development stage.

The negative consequences of the fourth uncertainty can be reduced through stronger cooperation with the foreign system operators, since, for example, the simulation model of the Western Danish 400 kV transmission grid includes part of the Northern German grid with electro-geometrical line data. Harmonic voltage measurements in geographically spread locations have also been received for the numerical tuning of the foreign contribution to the background harmonic distortion in Western Denmark in connection with the joint work on the establishment of the Westcoast 400 kV connection.

When assessing a future grid development stage, a model uncertainty margin to the planning levels can be applied. When the model uncertainty margin is 30% and the IEC planning levels [17] are applied, then the simulated harmonic distortion is concluded to be excessive, and it is concluded to likely require mitigation when it is exceeding 70% of the IEC planning level.

When the new grid expansion stage is brought into operation and the parameters of the new passive components are all measured and settled, the simulation model for harmonic assessment representing the new grid expansion stage will be revalidated and, if deemed necessary, recalibrated; i.e., the harmonic source vectors shall be revised and retuned.

The revalidation and recalibration process shall be conducted after the completion of a larger grid expansion or reconstruction project and before beginning the harmonic assessment of the following grid projects. For the simulation model of the Western Danish 400 kV transmission grid, revalidation (and recalibration) shall be scheduled to 2025, after the commissioning of the Westcoast 400 kV connection, and again to 2027, after the commissioning of the Revsing–Landerupgård 400 kV UGC connection.

4. Simulation Model Development Process

The development process of the model for the direct simulation of the harmonic voltage distortion is shown in Figure 6. The process includes three major steps.

The first step is the acquisition and processing of grid data and measurements. This step will result in time-synchronized data that combines the harmonic voltage measurements, the measured power transports—including transports of the HVDC Converter Stations with harmonic filters and converter-interfaced wind power plants (or other relevant large units in the grid)—and the grid operation conditions (n − L), L ≥ 0. The time-synchronized data are divided into the periods and snapshots to be used for the numerical tuning of the relative magnitudes and relative phase-angles of the harmonic source vectors and for the validation of the simulation model. Each period includes several consecutive snapshots so that not only the stationary magnitudes but also the changes and steps of the harmonic voltage distortion due to changes of the grid operation regimes are represented and validated in the simulation model. The number of all included snapshots is P as depicted in Equations (6) and (7).

The second step is the preparation of the passive-part grid model, including the electro-geometrical data of the transmission lines and the electrical data of the power transformers, shunts, and harmonic filters. These data are acquired using Energinet’s grid model database and the data received from the cable vendors and validated from the (fundamental) impedance measurements. The (n − L), L ≥ 0, operation conditions represent changes within the harmonic impedance matrix according to Equation (6). Thus, the included conditions, i.e., from 1 to P, define the area of the transmission grid where the developed model is applicable for the simulation of the harmonic voltage distortion in response to changes in the harmonic impedance matrix. Since the (n − L) operation conditions are for the Western Danish 400 kV transmission grid, the model developed and presented in this paper and the results for the direct simulation of the harmonic voltage distortion are applicable within the entire 400 kV transmission grid of Western Denmark (for the present grid development stage). The model includes the foreign contribution to the background harmonic distortion in Western Denmark, but the model is not readily prepared for the harmonic assessment of the foreign grid.

The third step is the definition and numerical tuning of the relative magnitudes and relative phase-angles of the harmonic source vectors according to Equations (6) and (7). The locations of the harmonic sources are their grid-connection substations in the domestic transmission grid, which for the main units are seen in Figure 1 and presented in

Table 1. Harmonic sources in the 400 kV transmission grid model of Western Denmark.

Unit	Substation	Main Harmonic Orders
Konti-Skan 1, 2 (740 MW) 1	V. Hassing 400 kV	11, 13, 23, 25
Skagerrak 3 (500 MW) 1	Tjele 400 kV	11, 13, 23, 25
Skagerrak 1, 2 (550 MW) 1	Tjele 150 kV	11, 13, 23, 25
Storebælt (600 MW) 1	Fraugde 400 kV	11, 13, 23, 25
Anholt (400 MW) 2	Anholt 220 kV 3	11, 13
Horns Rev C (400 MW) 2	Horns Rev C 220 kV 3	11, 13 4
Distributed sources	150 kV substations	2, 3, 5, 7

¹ LCC HVDC Converter Station. ² Offshore wind power plant. ³ Offshore platform. ⁴ Minor contribution.

.

This step shall result in the development of an empirically defined solution set, i.e., a set with the relative magnitudes and relative phase-angles, for the harmonic current and voltage vectors. The application in the simulation model of the tuned relative magnitudes and relative phase-angles of the harmonic current and voltage vectors shall result in good and sufficient agreement with the measured harmonic voltage distortion in the 400 kV substations, M, for all harmonic orders, k, and for all included snapshots, p.

When the validation is completed, the relative magnitudes and relative phase-angles are locked, i.e., they are not permitted to be changed in the assessment with modifications of the harmonic impedance matrix [Z_MNk]. Then, the empirically defined solution set of the harmonic emission sources becomes ready for the direct simulation of the harmonic voltage distortion in the present grid development stage and the prediction of the harmonic voltage distortion in future grid stages.

5. Validation

The major part of the model development and the numerical tuning of the harmonic source vectors was completed in 2020, with the data and measurements in the 400 kV transmission grid completed in the period 2017–2020. The Danish transmission grid is in continuous development as new connections and new converter-interfaced units have been brought into operation since its initial establishment. The major grid development has occurred in the 150 kV transmission grid, with only a minor influence on the harmonic distortion in the 400 kV grid. Therefore, the model validation cases are still relevant for the present grid stage as of 1 January 2023.

At this stage, it is important to emphasize that the numerical tuning of the harmonic source vectors is made for the grid development stage after the energization of the Vejle-Ådal 400 kV UGC, i.e., with amplified 11th harmonic voltage distortion in the 400 kV substations Fraugde and Trige. Why this emphasis is important will be explained in Section 6.

The simulated plots shown in the following subsections are in three phases, depicted using the following color marks:

• solid red—simulated phase A;
• solid green—simulated phase B;
• solid blue—simulated phase C;
• shaded red background—measured maximum of the three phases for the grid stage after the energization of the Vejle-Ådal 400 kV UGC.

5.1. Period of Varying 11th Harmonic Distortion

This validation period includes four consecutive snapshots, A, B, C, and D, where changed operation conditions of the Konti-Skan and Storebælt HVDC Converter Stations have resulted in varying magnitudes of the 11th harmonic voltages in Eastern Jutland (in the substations Landerupgård and Trige) and in Funen (in the substation Fraugde).

The operation conditions during the four snapshots are presented in Figure 7. In snapshot C, the Konti-Skan HVDC Converter Station is out of service, and thus this is considered as depicting an (n − 1) operation regime for the model validation. In snapshot D, both Konti-Skan HVDC Converter Station and the 400 kV line Ferslev–Nordjyllandsværket are out of service, and thus the snapshot is considered as depicting an (n − 2) operation regime.

The comparisons between the measured and simulated harmonic voltages in the 400 kV substations Trige, Fraugde, V. Hassing, and Landerupgård are shown in Figure 8, Figure 9, Figure 10, and Figure 11, respectively.

The model slightly overestimates the 11th harmonic voltage magnitudes in the 400 kV substations Fraugde and Trige, and slightly underestimates them in the 400 kV substation Landerupgård. For all characteristic harmonic orders, the model successfully captures both the magnitudes and tendencies for the harmonic voltages to change due to changed operation conditions in the Western Danish 400 kV transmission grid.

The IEC planning level of the 11th harmonic voltage distortion is 1.5% of the fundamental voltage [17]. The largest 11th harmonic distortion magnitudes in the 400 kV substations Fraugde, Trige, and Landerupgård occur when the Storebælt HVDC Converter Station is at 50% of the power transport and with a single harmonic filter in-service, and when the Konti-Skan HVDC Converter Station is at full power transport. In Trige and Landerupgård, the 11th harmonic voltage magnitudes are high, i.e., approaching or exceeding 70% of the IEC planning level, when the Konti-Skan HVDC Converter Station is at full power transport.

The measured and simulated harmonic voltage magnitudes are in good agreement in the presented validation period.

5.2. Period Containing High 11th Harmonic Voltage Magnitudes

This period includes four consecutive snapshots E, F, G, and H, where the operation conditions of the Konti-Skan and Storebælt HVDC Converter Stations result in the high magnitudes of the 11th harmonic voltages in the substations Trige and Landerupgård. The definition “high” means approaching or exceeding 70% of the IEC planning level.

The operation conditions during the four snapshots are presented in Figure 12. All shown snapshots are of (n − 0) operation regimes.

The comparison between the measured and simulated harmonic voltages in the 400 kV substations Trige and Landerupgård are shown in Figure 13 and Figure 14, respectively.

The simulation results may slightly differ from the measurements for the 11th harmonic voltage magnitudes in the shown 400 kV substations, and normally remain within a range of tolerances of the grid data, i.e., ± a few percent points.

The fact that the discrepancy is so small and consists of both signs for the same substations under different operation regimes implies that the numerical tuning is just right, i.e., it is neither overestimating nor underestimating the specific harmonic orders of the background distortion. Further attempts of numerical tuning may improve the simulation results in one group of substations but, at the same time, they may introduce a large discrepancy in another group of substations through various operation regimes applied for model development.

For all characteristic harmonic orders, the model captures both magnitudes and tendencies of how the harmonic voltages change due to changed operation conditions in the Western Danish 400 kV transmission grid.

Specifically, the model successfully predicts which harmonic orders and which substations are with high and low magnitudes. High magnitudes of the harmonic voltage distortion, especially when approaching and exceeding the IEC planning levels [17], may call for mitigation and investment decisions in the transmission grid. Therefore, the fact that the model accurately simulates such high magnitudes is a significant result of the presented deterministic method.

The model can also accurately predict when the harmonic voltage distortion will, for given operation conditions, have low magnitudes. However, the accuracy requirement to the simulated magnitudes can be lowered because there is no need for the mitigation or investment decisions for such low harmonic distortions. The lowered accuracy requirement means there is no need for exact matching between the simulated magnitudes and low measured magnitudes of the harmonic voltage distortion, and that rather, it is enough to just accurately predict the range of such low magnitudes.

5.3. Determining Harmonic Voltages of Not-Previously-Simulated Regimes

Recently, a harmonic assessment should be conducted for the grid-connection of a client under the 400 kV substation Endrup. For securing the accuracy of the harmonic assessment of a new operation regime after the grid-connection of the client, the simulation model has been benchmarked to the harmonic voltage measurements in Endrup for the present grid stage. For this purpose, the maximum 95th weekly percentiles measured over several weeks in 2022 have been compared to the simulations of different operation regimes, which should produce the highest magnitudes of the harmonic voltages in the 400 kV substation Endrup. The analysis of the measurements has revealed that the maximum 11th harmonic voltage magnitude in Endrup may occur when the Askær-Revsing 400 kV line is out-of-service.

Comparison of the measured 95th weekly percentiles and the simulation results for different operation regimes is given in Figure 15.

Simulated operation regime 4 is with the disconnected Askær-Revsing 400 kV line. As can be seen, the simulation model accurately simulates the maximum harmonic voltage distortion for the given case, confirming the statement in Section 3 that the model, once tuned and with locked characteristics of the harmonic sources, will be able to accurately predict harmonic voltage magnitudes in other operation regimes, which have not been included in the previously conducted model development.

6. Prediction of Harmonic Voltage Distortion

At present, the 400 kV connections which could significantly affect the harmonic voltage distortion in Western Denmark are either in the process of commissioning—for e.g., consider the Westcoast 400 kV connection to Germany—or in the design stage, as in the Revsing–Landerupgård 400 kV UGC. Therefore, the prediction of the harmonic voltage distortion in a changed transmission grid, i.e., how the method works, will be demonstrated using the Vejle-Ådal 400 kV UGC.

In the simulated, hypothetical, grid development case, the Vejle-Ådal 400 kV UGC will be replaced with a 400 kV OHL section at the Landerupgård–Malling transmission line, as seen in Figure 2. The idea of the hypothetical case is to simulate the grid development stage before the establishment of the Vejle-Ådal 400 kV UGC; i.e., before the significant amplification of the 11th harmonic voltage distortion in the 400 kV substations Trige and Fraugde occurred in July 2017, as shown in Figure 2. The other projects, which were completed in the 150 kV grid between July 2017 and the present, remain unrevoked in the hypothetical grid development case, as those projects are not considered to be significantly influencing the 11th harmonic distortion in the 400 kV transmission grid of either Eastern Jutland or Funen (see marking in Figure 1).

The operation regime snapshot B (depicted in Figure 7), which is shown to be resulting in excessive 11th harmonic voltage distortion, has been chosen for the demonstration of the prediction method. The relative magnitudes (a_Nk/a_R) and relative phase-angles of the harmonic source vectors α_Nak, α_NBk, and α_NCk in Equations (4) and (5) have been tuned for the present grid development stage, which is for the grid after establishment of the Vejle-Ådal 400 kV UGC. The tuned relative parameters are locked and must not be changed during simulation of the hypothetical grid development case which is considered as before the Vejle-Ådal 400 kV UGC.

For demonstration that the method works, the simulation of the hypothetical grid development case shall predict significantly reduced 11th harmonic voltage magnitudes in the 400 kV substations Trige and Fraugde (but not in the other substations) for snapshot B. The simulated phase-to-ground voltage magnitudes for snapshot B and the measured maximum of the three phase-to-ground 11th harmonic voltage magnitudes are shown in Figure 16, with the following color codes:

• solid blue—simulated magnitudes for the grid stage before the Vejle-Ådal UGC;
• solid green—simulated magnitudes for the grid stage after the Vejle-Ådal UGC;
• red shaded—measured maximum for the grid stage after the Vejle-Ådal UGC.

In the measured plots in Figure 16, the blue arrow marks the energization time of the Vejle-Ådal 400 kV UGC in the physical grid in July 2017.

Revoking the Vejle-Ådal 400 kV UGC in the passive grid model corresponds to the modification of the harmonic impedance matrix [Z_MNk]. Thus, the prediction method applies the locked characteristics of the harmonic vectors to simulate the change of the harmonic voltage distortion at a modified harmonic impedance of the grid, according to Equations (8) and (9). The measured and simulated harmonic voltages in Figure 16 are in good agreement as shown by the following observations:

• Significantly changed harmonic voltage magnitudes in Trige and Fraugde;
• High and unchanged magnitudes in Landerupgård;
• Low and unchanged magnitudes in V. Hassing.

The other interesting outcome is that the correlation between the 11th harmonic distortion in the 400 kV substations Trige and Landerupgård has changed after energization of the Vejle-Ådal 400 kV UGC. The correlations are presented in Figure 17 for the present grid stage and in Figure 18 for the grid stage before the establishment of the Vejle-Ådal 400 kV UGC.

From the plots in Figure 17 and Figure 18, the following can be seen:

• In the present grid, there is a strong positive correlation between the 11th harmonic voltage magnitudes in Trige and Landerupgård, i.e., there are simultaneously increasing and decreasing magnitudes;
• In the grid before the energization of the Vejle-Ådal 400 kV UGC, there is a possibility of weak negative correlation. Negative correlation means increasing the 11th harmonic voltage magnitudes in one substation while decreasing them in the other substation in specific regimes. The correlation is pronounced weak due to the round shape of the measured dependency, as depicted in Figure 18a.
• The correlation patterns have not been conclusively demonstrated by observational analysis using the simulation model, but they will be examined in more detail in Section 8.

The simulated tendencies fit with the measured positive correlation of the 11th harmonic voltages for the grid development stage after the energization of the Vejle-Ådal 400 kV UGC.

For the grid stage before Vejle-Ådal, the model shows no clear correlation as both negative and positive tendencies are found among different consecutive operation snapshots. Thus, the model shows no clear correlation though the conducted simple analysis of the harmonic voltage measurements points at the possibility of a weak negative correlation between the 11th harmonic voltages in the substations Trige and Landerupgård.

As demonstrated in this section, the method works well for the prediction of the harmonic voltage distortion in a new topology of the transmission grid, predicting significant change of the 11th harmonic magnitudes in the substations Trige and Fraugde as well as unchanged low or high magnitudes in the substations V. Hassing and Landerupgård, respectively.

7. Accuracy of Deterministic Model

The deterministically developed model accurately represents both the high and low magnitudes of the harmonic voltages in different operation conditions of the Western Danish 400 kV transmission grid of the present grid development stage, and can predict the harmonic voltage magnitudes in a changed grid. For most substations, the discrepancy between the simulated and measured harmonic orders with high magnitudes is within ± 10 percent of the measured magnitudes. However, there can be certain operation conditions where the tendency of the simulated magnitudes is correctly simulated, but the discrepancy is larger than 10 percent, such as in the snapshot G (shown in Figure 13).

The model accurately represents the harmonic orders with small magnitudes as well. However, the discrepancy in percentage form of the measured magnitudes can be higher than 10 percent because the comparison is made between two small numbers, i.e., a low simulated magnitude and a low measured magnitude. In practice, such discrepancy will be insignificant because small harmonic magnitudes do not require mitigation. However, the ability to represent either high or small magnitudes for specific harmonic orders in specific operation conditions of the meshed transmission grid is crucial for the overall accuracy of the presented simulation model and the method.

After analyzing the results of various validation cases, the simulation model is deemed to be 80–90% accurate; i.e., 80–90% of the operation regimes will be simulated with high accuracy and small discrepancy comparing to the harmonic voltage measurements, while 10–20% of the operation regimes may have larger discrepancy for certain harmonic orders in certain substations.

The application of the deterministically developed model will be done by using observational studies including a large number of simulation cases that combine different transports and different (n − 1) operation conditions for finding patterns for the harmonic voltage distortion and mitigating results with potentially larger discrepancies (if such results are deviant from the common patterns).

The next section introduces a statistical method to explain harmonic patterns in the 400 kV grid from harmonic voltage measurements, with the aim of identifying the most significant operating regimes and reducing the number of observational studies.

8. Analysis of Harmonic Patterns from Measurements

The measured harmonic voltages from different substations are analyzed in this section using Principal Component Analysis (PCA), a statistical method for exploring large data sets. The analysis in this paper is limited to demonstrating the relationship between the Eastern and Western 400 kV systems in Jutland (see marking in Figure 1), specifically the relationship between the measured 11th and 13th harmonic voltages and the number of filters in operation at HVDC stations.

8.1. Preparing a Dataset

A dataset with variables according to

Table 2. Measured variables from a set of substations in the 400 kV grid of Figure 1.

Variable	Description
xxx_11 (kV)	Measured 11th harmonic voltage; xxx indicates substation with measurement
xxx_13 (kV)	Measured 13th harmonic voltage; xxx indicates substation with measurement
VHA. CZF.1.5.SUM	Number of harmonic filters in operation at substation V. Hassing 400 kV (VHA) with HVDC system KS12
VHA. CZF2	AC filter (of other type) in operation at substation VHA
FGD. CZF.1.4.SUM	Number of harmonic filters in operation at substation Fraugde 400 kV (FGD) with HVDC system SB
I. CZF.1.4.SUM	Number of harmonic filters in operation at substation Tjele 400 kV (I) with HVDC system SK3

is prepared. Substations with measurements include Endrup (EDR), Ferslev (FER), Fraugde (FGD), Kassø (KAS), Landerupgård (LAG), Malling (MAL), Tjele (TJE), Trige (TRI), and V. Hassing (VHA). Substation names are abbreviated, and this abbreviation is used as a variable name.

This data set has a total of 20 columns, with each observation containing 1 hour of aggregated 10-minute variable values. In this case, the dataset contains 1464 observations, with the first 176 measured before Vejle-Ådal and the last 1288 measured after Vejle-Ådal. As a result, this limited dataset contains two types of data.

8.2. Principal Component Analysis

This paper introduces PCA briefly, but for a detailed description of the method, refer to relevant literature such as the references [22,23,24] or similar papers. The description that follows is based on this literature. PCA is a technique for analyzing large data sets with many variables per observation, increasing data interpretation while preserving the maximum amount of information and enabling multidimensional data visualization. PCA is a statistical technique for reducing the dimensionality of a data set by linearly transforming the data into a new coordinate system (referred to as principal component space or PC-space) where most of the variation in the data can be described with fewer dimensions than the original data. Many studies use this technique to plot data in two dimensions and visually identify clusters of closely related data points. The results of PCA depend on the scaling of the variables, and by scaling each variable with its standard deviation, one ends up with dimensionless variables with unit variation, therefore finding the variable correlations to be comparable [23,24].

The PCA model can be described mathematically thus:

$$\left[X\right]=\left[T\right]{[P]}^T+\left[E\right]$$

(10)

where [X] is the matrix with data that one wishes to explain, [T][P]^T is the explained part of the data, and [E] is the part of the data which is not explained by the analysis. Matrix [E] consists of residuals or noise. The matrix [P] is referred to as Loads, and it contains unit vectors of the principal component axes (PCs). These components are mutually orthogonal and oriented in the variable space to capture the direction of maximum variation of data points. Matrix [T] is referred to as Scores and it contains the coordinates of the projected datapoints on the PC-space formed by the orthogonal PCs. Loadings and scores are used to interpret PCA results. The more PCs there are, the less noise is left, and the more data are explained [22,23,24].

The loadings for both groups of data and the first two principal components (more exists, but they are not presented in this paper) are shown in Figure 19. Data before Vejle-Ådal are shown on the left (a) and data after Vejle-Ådal are shown on the right (b).

A loading plot shows the relationship between two PCs, and the percentage indicates the amount of variance explained by the component. It depicts how much each variable influences a PC and the relationship between the variables. The stronger the influence is, the closer the variable is to a PC axis. Loadings range from −1 to 1 and a high absolute value (towards 1 or −1) indicates that the variable has a significant variation. Values near 0 indicate that the variable has minor variation. The sign of a loading (+ or −) indicates whether a variable and a PC are correlated positively or negatively.

8.2.1. Before Vejle-Ådal

Loadings for this case are shown in Figure 19a. Approximately 46% of the variance in data is explained by PC1 (Comp 1—25.74%) and PC2 (Comp 2—20.31%). The biggest sources of variation are variables TRI_13, VHA_13, and FER_13 as they are farthest away from the origin along PC1 in the loadings plot. These variables seem to be positively correlated with each other: that is, they vary in the same direction. Along PC2, LAG_11, LAG_13, and FGD_11 are the biggest sources of variation. The first two variables seem to be positively correlated with each other. The third variable is likely to be negatively correlated with the other two.

Variables TRI_11 and LAG_11 is seen to have almost the same variation (0.25–0.30) but they are almost orthogonal to each other in the plane spanned by PC1 and PC2, and therefore, they are most likely not correlated. In Section 6, the interpretation was that these two variables were negatively correlated in the simple analysis of the harmonic voltage measurements and with no clear correlation in the simulation model, but from this data set, it is evident that they are most likely not correlated (this is in agreement with the model, see no correlation results in Figure 18b,c).

8.2.2. After Vejle-Ådal

The loading plot in Figure 19b is for the case after Vejle-Ådal, and it shows a different pattern. Here, TRI_13 is still the biggest source of variation along PC1 and now TJE_13, KAS_13, and EDR_13 is the biggest sources of variation along PC2. Approximately 55% of the variance in data is explained by PC1 and PC2 put together. It is seen that the variables TRI_11 and LAG_11 now most likely has a correlation as they are close in the PC-space; this was also the interpretation made in Section 6, see results in Figure 17.

We can also assess the variation caused by the number of filters in use at HVDC-stations. Filters for the Konti-Skan 1,2 HVDC station in V. Hassing (VHA.CZF.1.5.SUM) correlate positively with PC1, while filters for the Storebælt HVDC station in Fraugde (FGD.CZF.1.4.SUM) correlate negatively with PC1. This means that the number of filters in operation has the largest influence on harmonic level in the substations Ferslev and Trige, because these variables correlate with PC1, whereas the number of the harmonic filters in Fraugde will ‘pull in the opposite direction’. This means that more filters in Fraugde will cause the harmonic voltages in Ferslev and Trige to decrease.

This can readily be seen from a scores plot, as shown in Figure 20. Scores are colored based on the number of filters in use for Konti-Skan 1,2 in V. Hassing in (a) and in (b), scores are colored based on the number of filters in use for Storebælt in Fraugde.

Recall the location of variables FER_11 and TRI_13 in Figure 19b; then, according to Figure 20a, the higher the 11th harmonic voltage in Ferslev and the 13th harmonic voltage in Trige are, the more filters will be in use in V. Hassing, as the scores for many filters in operation are located on the far right in this plot.

The opposite is seen in Figure 20b, as a high number of filters in operation in Fraugde are situated on the far left in the scores plot. This means that if more filters are in operation in Fraugde, lower harmonic voltage will most likely be seen in Ferslev and Trige, because the harmonic filters in Fraugde and V. Hassing work against each other.

The variables EDR_11, EDR_13, KAS_13, I_13, and I.CZF.1.4.SUM correlate with PC2 (loading plot on Figure 19b) and thus are orthogonal to the variables along PC1. This strongly suggests that harmonic voltage in the substations Tjele, Kassø, and Endrup is correlated with that in the Skagerrak 3 HVDC station in Tjele and its filters, but has little or no correlation with the other two, i.e., the Konti-Skan 1,2 and Storebælt HVDC stations in V. Hassing and Fraugde, respectively.

From a practical perspective, it can be said that the electrical properties of the transmission lines Ferslev–Tjele and Kassø–Landerupgård now limit the propagation of harmonics between the Eastern and Western parts of the 400 kV network in Jutland due to the insertion of the 400 kV UGC at Vejle-Ådal; this is illustrated in Figure 21.

An assumption regarding uncorrelated harmonic distortion in these two systems of the 400 kV network was also made in Section 2 after conducting and analyzing an extensive number of simulation cases. However, the presented statistical analysis, including the operating conditions of the HVDC facilities by means of the harmonic filter status, allows one to reach the same conclusion using simpler efforts.

This section has demonstrated that advanced, statistical methods are powerful instruments for the validation of the simulation model. An overly simplified analysis of the measurements may lead to an erroneous conclusion.

9. Analysis of the Results Received

Comparison of the simulated and measured magnitudes of the harmonic voltage distortions has been presented for the operation conditions of the 400 kV transmission grid for the substations Fraugde, Landerupgård, and Trige, as explained below:

• Varying between low and high magnitudes of the 11th harmonic voltages, in the consecutive snapshots A through D as presented in Section 5.1;
• High magnitudes of the 11th harmonic voltages, as shown in the snapshots E through H in Section 5.2.

The simulation model accurately simulates both the magnitudes and the change tendencies, i.e., the variation of the magnitudes, of the harmonic voltage distortion following changes of the grid operation conditions.

For the model preparation, an acceptable range of discrepancies between the simulated magnitudes and the measured magnitudes was ±10%; the given range applies for the harmonic orders with high magnitudes, i.e., those approaching or exceeding 70% of the IEC planning levels [17]. For the harmonic orders with low magnitudes, the model accurately simulates that the magnitudes are low—for e.g., as found in the measurements—but the discrepancy in % can be higher than ±10%. Higher discrepancies have been accepted for the low harmonic voltage magnitudes due to the common difficulty faced by the numerical methods when achieving high accuracy for simulations of small numbers, and because low harmonic voltages will be accepted without requiring mitigation proposals.

The presented work is not an academic exercise targeting marginally perfect alignment between the measurements and simulations for some selected few cases. Instead, the presented deterministic method is for the accurate simulation and prediction of the harmonic voltage distortion in as-many-as-possible various grid operation conditions with an acceptable accuracy range of the simulation results when compared to the harmonic voltage measurements.

An analysis of both presented and conducted (but not presented in this paper) validation cases reveals that the accuracy of the simulation model is 80–90%. This figure indicates that the discrepancy between the simulated and measured magnitudes will be according to the acceptable range for 80–90% of all conducted cases. In the remaining 10–20% of the cases, the discrepancy can be larger than ±10%: for example, in Figure 13 for the snapshot G; the tendency seems correct, but the simulated magnitude has been overestimated when compared to the measurement.

The model can accurately simulate the cases with other operation conditions of the grid, which have not been included in the model preparation. This outcome was demonstrated in Section 5.3, where the maximum of the 95th weekly percentiles of the measured harmonic voltage in the 400 kV substation Endrup was successfully compared to the simulated magnitudes for the cases with operation conditions resulting in high harmonic magnitudes. The presented case for Endrup was conducted after the model preparation was completed.

Finally, the model accurately simulates the behavior of the harmonic voltage distortion in the 400 kV transmission grid of Western Denmark, with changed 11th harmonic voltage magnitudes in Fraugde and Trige but not in the other substations and not for the other harmonic orders, due to the replacement of an OHL section with a UGC on a 400 kV line: the Vejle-Ådal case presented in Section 6. The Vejle-Ådal case has not been part of the model preparation, the successful simulation of which can strengthen the confidence in the simulation model and the presented deterministic method.

10. Conclusions

The green transition in electric energy generation and consumption as well as the construction of the Energy Islands require the accelerated development and reconstruction of the transmission grid. In Denmark, the grid development will, for the most technically possible extent, utilize HVAC UGCs instead of OHLs. The increasing share of UGCs has brought the resonances of harmonic grid impedances down into the range of background harmonic distortion. Further grid development may increase the harmonic voltage distortion in some parts, and reduce the harmonic distortion in other parts, of the meshed transmission grid. Therefore, securing adequate levels of power quality in future grid development stages and developing prediction methods for harmonic voltage distortion by using simulations become essential for a successful green transition.

This paper has proposed a method for the direct simulation of the harmonic voltage distortion in the meshed transmission grid, which has many harmonic emission sources. In the method, different operation regimes of the passive transmission grid, (n − L), L ≥ 0, represent variations of the harmonic impedance matrix of the present grid stage. The harmonic emission sources are numerically adjusted for the best possible matching the measured harmonic voltage magnitudes in different substations for the different operation regimes. The numerical tuning converges to an empirically defined solution set of the relative magnitudes and relative phase-angles of the harmonic emission sources, so that the same solution set matches the measured harmonic voltage distortion for the most possible operation conditions in the meshed transmission grid.

The method interprets the grid development projects with new or modified connections as modifications in the harmonic impedance matrix and applies the solution set with locked relative magnitudes and phase-angles of the harmonic emission sources for the simulation and prediction of the harmonic voltage distortion in a future grid. The primary goal of the method has been to predict whether, how severely, and where in the transmission grid the harmonic voltage distortion will become excessive and need mitigation for grid development projects. The strength of the method is its ability to point out excessive harmonic distortion and propose mitigation before a project with new or modified connections is commissioned and brought into operation.

The method has been illustrated using a simulation model of the Western Danish 400 kV transmission grid. The method is deterministic and subject to uncertainties influencing the accuracy of the simulated harmonic voltage distortion. Such uncertainties include the tolerances of the passive grid data, measurements applied for the model setup, and representations of the harmonic emission sources. Therefore, the simulation model will be revalidated and, when deemed necessary, recalibrated after the establishment of new passive-grid components and new harmonic sources.

It is also demonstrated that, despite the lack of voltage angle measurements, harmonic voltage measurements support the developed method when analyzed using statistical methods: in this case, PCA-analysis. Historical measurements can thus be used to determine possible harmonic correlations between substations in the meshed grid, as well as their causes. Advanced, statistical methods are necessary because an overly simplified analysis of the measurements may lead to erroneous conclusions regarding the harmonic behavior within the meshed transmission grid, especially when validating the simulation model.