Harmonic Propagation Model for Analyses in Meshed Power Systems

This paper presents a novel harmonic propagation model for the analysis of harmonic propagation in meshed transmission networks. This method is based on electromagnetic waves and uses scattering parameters known from microwave theory. The physical understanding of harmonic distribution in a power system is achieved. Changes may be quantified, and the method is validated against both physical harmonic distortion measurements and simulations. This method differs from traditional analysis methods in that the propagation pathways, as well as the contribution of individual sources to the propagated harmonic voltage in the entire network, can be determined unambiguously, allowing the influence of grid changes on the harmonic voltage to be predicted.


I. INTRODUCTION
P OWER systems around the world undergo significant changes these years due to the green transition and integration of renewable energy sources.Traditionally, power transmission networks have mainly been built using overhead lines (OHL)s.However, transmission system operators (TSO)s have increasingly used underground cables (UGC)s for network extensions due to public opposition, as well as political and environmental pressure on TSOs against the construction of additional OHL transmission infrastructure.
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TPWRD.2023.3334389.
Digital Object Identifier 10.1109/TPWRD.2023.3334389 Harmonic studies in meshed transmission systems are typically harmonic load flow calculations and frequency sweep calculations related to specific network nodes.We can study the difference in results for various study cases and conclude whether the harmonic distortion has increased or decreased, or if resonance frequencies have "moved" towards higher or lower frequencies [7], [8], [9], [10], [11], [12].It is not obvious with these analysis techniques to determine if a change in harmonic voltage is due to amplification phenomena or a change in the propagation throughout the meshed system, that is, a change in the distribution of harmonic voltages throughout the grid.
Harmonic propagation can be explained using standing wave theory, according to research [12], [15], [16].In [15], it is demonstrated how a simple change in line termination or the transition between an OHL and UGC line type affects harmonic propagation due to the changed terminal conditions.Harmonic penetration analyses are used in [16] on a simple model of a meshed system with both OHL and UGC lines.Standing wave phenomena were observed here, and the significance of reflection coefficients was discussed.Reference [16] discovers that it is difficult to retrieve reflection coefficients using an analytical approach in a large, meshed grid using the traditional wave theory approach.

A. Motivation
This paper presents a novel propagation model based on standing wave techniques and tools known from radio and microwave theory.The method uses the properties of standing harmonic voltage waves to determine the propagation of harmonic waves in a meshed transmission system and to quantify the preferred propagation paths through the grid.Thus, the purpose of this method is not equated with the need to determine harmonic distortion levels at substations.
By quantifying network propagation patterns, the influence of different network designs on these propagation patterns may be compared.Thus, several network topologies and, in particular, the effect of cable lengths on propagation paths can be analysed.Identification of these principal propagation paths can aid in the identification of sites in the meshed grid where mitigating equipment should preferably be implemented, or alternate network designs can be utilised to eliminate the harmonic problem.
Because transmission networks will be greatly expanded in the future with UGCs and a greater degree of interconnection, it is critical for a TSO to understand the physical impact of these grid changes on harmonic voltage propagation.It is possible to quantify the change in standing harmonic waves caused by natural reflections created by the transition between OHL and UGC or between a line and a substation with this method.
In the following, we first explain the relation of the proposed method to traditional methods such as harmonic penetration analysis and impedance analysis.Next, the merits of the method are outlined.
1) Harmonic Voltage and Current: Most commercial harmonics simulation packages include harmonic penetration solution methods based on a fundamental frequency load flow.These methods are frequently referred to as "harmonic load flow," but usually they do not perform a true load flow solution at harmonic orders, where "true load flow" refers to the ordinary iterative load flow performed for each harmonic frequency targeted in the analysis.The methods are as follows: r Harmonic penetration analysis: [10], [11] makes use of injected current calculated as a function of the fundamental bus voltage and source parameters.Iteratively determined injected current based on calculated harmonic voltages is used in true load flow techniques [10], [11].The sum of harmonic voltages at nodes and the sum of harmonic current flows are typical results for both techniques.
r The propagation model: makes use of the harmonic voltage at the source bus/node, which is determined by power system measurements or simulations.The output voltage is either the harmonic voltage propagating from each individual source or the sum of harmonic voltages at nodes/buses.Line currents can also be shown.
2) Harmonic Network Impedance: Most commercial harmonics simulation packages include frequency scan methods for displaying the network's frequency response as seen from a specific bus or node.Another admittance-based method for analysing the impact of resonance in a system is modal resonance analysis.The techniques are as follows: r Frequency scan: [11] recalculates the network admittance matrix at each frequency step and applies current injection with fixed magnitude to obtain the corresponding bus voltages.The frequency response of the network as seen from a specific bus or node is the output.
r Modal resonance analysis: [17], [18] uses frequency scan results as input for modal analysis to determine the frequency at which harmonic resonance excitation or observation is most likely identified as seen from a bus.The results are based on evaluation of eigenvectors.
r The propagation model: can determine the frequency re- sponse of the network at a specific bus or node using incident and reflected voltage waves, reflection coefficients, and characteristic impedances.
3) Harmonic Propagation: Harmonic voltages are the primary concern for TSOs using the IEC 61000 standard [19] for harmonic assessment because harmonic distortion limits based on voltage for HV and EHV power systems have been established in [19], and harmonic voltages are typically measured in the grid.When comparing study cases, harmonic penetration or load flow analyses can provide observations of bus voltages, but the distribution patterns throughout the grid are difficult to quantify.The following are the differences between the methods: r Harmonic penetration analysis: or harmonic load flow analysis can provide harmonic voltage magnitudes at buses that can be compared between study cases.
r Transfer gain methods: are used for propagation analyses.
References [20], [21] discusses the advantages, disadvantages, and applicability of transfer gain methods for various network topologies.It concludes that while their use yields accurate results in radial network topologies, the results in meshed power systems are overestimated due to multiple propagation pathways.
r The propagation model: provides propagated voltage throughout the grid due to reflections caused by network impedance transitions.The magnitude of reflections and propagated voltage can be quantified, and specific pathways can be identified.Harmonic voltage paths from each individual source to specific points in the grid can be identified and quantified individually.It follows from the above that the propagation model can provide the same types of results as traditional methods but using other calculation techniques.Because the methodology will be used for the analysis of propagation pathways in the network and not for determining the sum of voltages at substations, the following evaluations will be demonstrated in this paper: r Assessment of reflection and transmission coefficients as an element to explain the physical circumstances of OHL and UGC that affect the propagation and pinpoint the exact location(s) causing harmonic problems.
r Assessment of propagated waves throughout the network for a specific real-world case study to identify changes in propagation patterns caused by partial cable laying in a transmission line.
r Assessment of changes in propagation patterns in the grid from specific sources.

B. Paper Outline
In Section II, the concept of harmonic propagation is explained and in Sections III and IV, the propagation model and its elements are described.In Section V, the developed model is validated against real life measurements and a simulation model developed by [21], and in Sections VI and VII, study cases demonstrating the value of the novel propagation model are presented.Sections VIII and IX discuss the applicability and conclusions of the method.

II. CONCEPT OF HARMONIC PROPAGATION
The proposed model is based on electromagnetic waves and for a general transmission line, including loss effects, the propagation constant γ is complex valued and frequency dependant.The voltage at any point x along a line is [22], [23]: where u + is a forward (incident) travelling wave and u − is a backward (reflected) travelling wave.The voltage propagation throughout a complete meshed network varies for each study case in consideration, and because we want to quantify Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
case-specific changes to the propagation in comparison to a base case, we will define a propagation coefficient at node i as: Here, ζ i is the ratio between the sum of incident and the sum of reflected voltage waves at node i.The benefit of this parameter is explained in Section VI.

III. PROPAGATION MODEL DESIGN CONCEPT
This section presents the design of the propagation model.The process begins with a review of some basic microwave theory techniques [23], [24], followed by interpretation in relation to a power system.The individual component models are also presented here.

A. Scattering Parameters
In power system analyses, matrix descriptions such as the impedance matrix Z or admittance matrix Y are commonly used to relate voltages and currents for two-port network models.Scattering parameters are another type of matrix parameter that is used for practical characterisation of microwave network structures in the kilo-hertz to giga-hertz frequency range [23], [24], [25].Scattering parameters, or S-parameters, apply to all frequencies, but they are most used in radio frequency or microwave frequency networks, where signal power and energy considerations are easier to quantify than currents and voltages.
Typically, S-parameters are calculated in terms of other parameters such as Z , Y , or transmission parameters (ABCDparameters) using the characteristic impedance, propagation constant, and, in the case of transmission lines, the line length.S-parameters, on the other hand, have a very different interpretation because they contain information about reflections carried by a propagating wave in an interconnected circuit.S-parameters are bidirectional, like Z , Y , and ABCD-parameters, but they include information about reflections that Z , Y , and ABCDparameters do not.Because standing wave patterns are formed by wave reflections, S-parameters are used in this paper for the propagation model.S-parameters are not discussed further because they are well described in the literature, and the reader is encouraged to study [23], [24], [25], [26] or similar.In this case, we will use the common expression: where − → a is the excitation of a network S , characterised by a square matrix of complex S-parameters and − → b is the network response when exited by − → a .Equation (3) can be written as: The first number in the subscript of S ij refers to the responding port, while the second number refers to the incident port.Thus S 21 means the response at port 2 due to a signal at port 1.The port reflection coefficients {S 11 , S 22 , . . ., S NN } refers to the ratio of the amplitude of the voltage that reflects from the port to the voltage incident to the port.The port reflection is due to impedance discontinuity in the transmission medium.The voltage gain parameters {S i1 , S i2 , . . ., S iN }, or the transmission coefficients, describe the transferred voltage from one port to another.The focus of this paper is on steady-state harmonics and therefore, we will assume that the system is an LTI-system, which implies linearity, causality, and time invariance.If non-linear components are to be considered, their models must be linearized to a specific working point or, if applicable, to a range of operation to be evaluated.

B. Interpretation in Relation to Power Systems
In connection with a power system, we may define scattering as a way in which the travelling currents and voltages are affected when encountering a discontinuity caused by the insertion of another network element in the transmission line or due to the transmission line terminal conditions.
Our definition corresponds then to a wave along a line impinging an impedance, which is different from the characteristic impedance of the line itself.This approach differs from the microwave technology and instead of considering a power system as a 'black box' and measuring power waves at the ports, we will model the individual transmission lines with a common reference impedance and connect them to nodes to determine the propagation of voltage waves throughout the entire meshed system.In this case, the voltage − → b is used as a measure of propagation, because − → b is the network response given by (4).Elements for transmission lines, sources, loads, and junctions are developed for propagation models.These components are shown in Appendix A.

IV. CONSTRUCTING A PROPAGATION MODEL
An example of cascaded elements in a network is shown in Fig. 1.Arrows indicate the direction of waves, where direction away from an element is the response b i and towards an element is the excitation a i , where i is the node index.
A source S on the far left has the response b 2 at the output port.The two-port network A on the left immediately after the source Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.The output b 9 , for instance, is the product of the wave b 2 travelling from the source through the networks A, B and C taking the reflections at all ports into consideration.The same type of relationship is true for the output wave b 11 at port 2 of network D. In this case, D is terminated in a load L causing a reflection b 12 which is travelling back to the source through networks D, B and A. Some of this wave is also scattered to network C.
Several methods exists for calculating the scattering matrix of an arbitrary microwave network and common for these methods is the idea of concatenating the individual scattering matrices of the elements in the network [27], [28], [29].A similar approach will be used in this case, with the exception that the methods developed by [27], [28] are designed to solve (3) as a multiport microwave circuit element.For this approach, the propagation of − → b throughout a network of individual connected transmission lines will be determined for a single value of a at a time, say, a 1 , in order to determine a single source contribution to the overall propagation.Therefore, a 2:N needs to be determined before solving (3).The approach is described in the following.

A. Model Description 1) Connection Matrix:
We will begin by building a connection matrix C that contains information about how the blocks that represent the lines/terminations are interconnected in the network.This can be expressed mathematically as: C describes how − → a and − → b is related as shown in Fig. 1.As there are 13 ports in the cascaded network on the figure, the connection matrix must be a 13x13 matrix and the principle is as follows: The connection between incident and reflected waves is simply designated with unity.For instance, the output b 2 is connected to the input a 3 as unity in the 3 rd row and 2nd column of C. The connections related to elements S, A and B in Fig. 1 are marked in ( 6) with a circle.

2) S-Matrix:
Matrix S is a partitioned matrix or a block diagonal matrix.The network model of Fig. 1 is: (7) and S is also a 13x13 matrix with the component S-parameters at the diagonal in the same order as in (6).
3) Solving the Equations: If ( 3) is inserted into (5), with S given as in (7), we can write: From ( 8), it is evident that the input incident vector − → a depends not only on itself, but also of its distribution in the network, that is, its output is related by CS.Equation ( 8) can be rewritten as: where I is the identity matrix and M is the network matrix.
In Fig. 1, a source is connected at block A and assuming the incident wave a 1 is known, the remaining incident waves a 2:N can be solved after elimination of the first row in M , giving the following expression: Finally, the reflected waves b 1:N can be obtained from (3).In case of multiple sources, (10) is solved for each source and the resulting incident and reflected waves are calculated by superposition.

B. Modelling Summary
Fig. 2 shows the modelling approach.The propagation model is designed based on a given system layout.First, a network model consisting of two-port transmission lines, sources and loads as well as N-port junctions is constructed and ports are labelled according to their position in the matrices.The connection matrix C describing the interconnection of the network blocks is designed as well as the partitioned component matrix S .The next step is to calculate a network matrix M , incident waves − → a and reflected waves − → b .The results are stored and in case of multiple sources, M is re-ordered and new incident and reflected waves are calculated.Lastly, the resulting waves are determined by means of superposition.The node voltages can be determined as the sum of forward travelling waves and reflected travelling waves at a node by means of (1).

V. PROPAGATION MODEL VALIDATION
The propagation model method is validated in this section against a real-world study case utilising harmonic voltage measurements and simulation results provided by Energinet, the Danish Transmission System Operator [21].Although the suggested method is not designed for penetration analysis, harmonic Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.voltages computed at substations using the propagation model and harmonic penetration analysis done using a simulation model will be used for validation.400 kV lines and substations in the Denmark-West area as well as three HVDC converter stations are shown.Increased 11th (and 13th) harmonic distortion were measured at TRI and FGD due to the replacement of an approximately 7 km OHL section along the line LAG-MAL with two parallel UGC of 7 km, at Vejle-Ådal (VLA), approximately 80-90 km away from TRI and FGD, respectively.The work by [21] uses this real-life case for validation of harmonic assessment studies.
A model of Denmark-West 400 kV grid is constructed according to the flow chart in Fig. 2 and the propagation model is shown in Fig. 4, where solid lines is OHL type, dotted lines UGC type and substation names are according to Fig. 3. Node numbers for each element are shown.A parallel UGC section with two junctions and two UGC line sections is represented by a special symbol, made up of ten nodes.The grid is divided into three regions, with region I containing nodes 1-19 and substations FGD, KIN, LAG and KAS, region II with nodes 20-109 and substations ASR, REV and TJE.Region III contains nodes 110-178 and substations MAL, TRI, FER, VHA and NVV.Line data is given in appendix B.
A reference impedance is needed for determining the Sparameters [25].In this paper, Z 0,ref (11h) = 298.11Ω is used and this value is the real part of a complex characteristic impedance of an arbitrary OHL section in the grid and therefore, it is possible to evaluate the change in harmonic propagation with respect to the transition between OHLs and UGCs.
In [21], the HVDC converter stations located in VHA, FGD and TJE (according to Fig. 3) are the identified harmonic sources Fig. 3. 400 kV transmission grid of Denmark-West.Legend: solid -400 kV transmission lines and substations, dashed -HVDC and converter stations, dotted -UGC section at Vejle-Ådal (VLA).Substation names and abbreviations are shown.[21] together with the 400 kV interconnectors to Germany at substation KAS.Table I shows the applied source values for the 11th harmonic and these values are also used in this paper for the model evaluation.The method in [21] uses an extensive and difficult tuning process of both amplitudes and angles for both sources and loads, in order to achieve coincidence with the measured voltage amplitudes.This means that the source angles are provided empirically.

TABLE I HARMONIC SOURCE VALUES USED IN THE SIMULATION MODEL OF
A grid model of Denmark-West is also implemented in a commercial power system simulation software, in order to evaluate the presented method in Section IV against a simulation tool.This grid model is referred to as the evaluation model in this paper, and for simplicity, it does not take geometric asymmetries, transformers, shunt reactors, or loads into account, as was the case in [21].These components can however, easily be integrated as two-port models.As a result, asymmetry and damping caused by loads and reactive units are not represented, and the results cannot be expected to be completely comparable with the measurements.The evaluation model in the commercial software uses distributed parameters and data in appendix B. Equal parameters are used for parallel UGC sections.

TABLE II MEASURED AND CALCULATED HARMONIC DISTORTIONS
The measured harmonic distortion as well as calculated values for substation TRI are shown in Table II.A value of 1.47 % is measured and a value of 1.42 % is simulated by Energinet [21].The results for both the validation model and the proposed propagation model are approx 1.55 % and this corresponds to 9.1 % difference between [21] and the propagation model results, and a 5.4 % deviation from the actual measured harmonic distortion at substation TRI.The results obtained with respectively the evaluation model and the propagation model for other nodes, excluding source nodes, also appear in the table.There is good agreement between the obtained results (less than 5 % deviation).Results for these nodes are not given in [21].There is sufficient agreement between the individual results, taking the model differences into account.

VI. EVALUATION OF PROPAGATION MODEL PARAMETERS
In this section, the propagation model parameters are evaluated and discussed.These parameters clarify the physical relationships for stationary waves at a line section transition.For this evaluation, three study cases given in Table III are used.
The first case 1 A is related to an incident wave travelling through line sections of the same type, whereas cases 1B and 1 C are related to an incident wave travelling through the boundary between two different characteristic impedances.For all cases,  Thévenin source types with zero impedance and 1∠0 • per unit is used as well as a matched load.The 550 Hz impedance and admittance according to Appendix B are used for these study cases.Propagation models according to Fig. 2 are set up.Fig. 5 shows the three study case models.
For each element, the element type and harmonic impedance and admittance are given.All UGC sections use equal parameters.

A. Evaluating S-Parameters
Table IV shows the S-parameters for the individual line sections of Fig. 5.The nodes given in brackets are for line sections according to Fig. 5(c).
Parameters S 12 and S 21 gives the transferred voltage along the line sections and it is seen that the UGC sections reduces the transferred voltage to approximately 83 % of the input voltage.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE V PROPAGATED WAVE AND PROPAGATION COEFFICIENT CASE 1A-1 C
It is also seen that the UGC section implies a larger angular displacement than the OHL section (−36.5 • − (−4.7 • ) = −31.8• ).The reflection coefficients S 11 and S 22 determines the reflection due to the impedance mismatch between line sections.For OHLs, minor reflections of 0.0003 − 0.0021 are seen due to attenuation effects and this has little practical influence on the harmonic propagation in power systems.For the transition from OHL to UGC sections, the reflection is significant as it is 0.5507, that is, approximately 55 % of the wave is reflected.

B. Evaluating Wave Propagation
The amplitude and angle of the propagated waves b i and propagation coefficient ζ 5 for node 5 are shown in Table V.
For study case 1 A it can be seen that |b| is close to zero (0.002-0.003) and ∠b ≈ 74 • at nodes 3, 5 and 7 along the line sections, indicating very low reflections.This is because the line sections are all OHL of same type and therefore, almost no reflection occurs as observed in table IV.In other words, a continuous stationary voltage wave is present from source to and all harmonic voltage propagates from source to load, only limited by the attenuation along the lines.At the load (node 9), b 9 = 0, indicating a matched load.
For study case 1B, |b 5 | is high (0.857) compared to study case 1 A (0.002) and ∠b 5 is also larger (−104.8• ) than for study case 1 A (−73.9 • ).This is because of the major reflection at node 5, where S 11 = 0.5507∠ − 126.6 • as seen from table IV.Due to the major reflection, the voltage at node 5 increases above the applied source value: and this increased voltage propagates back towards the source.The same conditions apply at node 6 and forward towards the load.b 5 and ζ 5 are given in Table V.
Study case 1 C is a case, where two UGCs are in parallel between nodes 5 and 14.There is a 3-port junction at these nodes.Evaluating the results with same approach as for study case 1B (11), the voltage at node 5 is 1.69∠ − 16.1 • pu, which is even higher (0.44 per unit) than for study case 1B.
The voltage level at a node is determined not only by the amplitude of incident and reflected wave, but also by the angle difference between the waves.Evaluating the phase difference between a and b for different study cases gives an indication of whether the voltage is likely to increase or decrease at a specific node as the voltage is a sum of two complex vectors.This can be evaluated by means of the propagation coefficient in (2).As the S-parameter angles is determined by means of γl, they are highly dependent of the line length as γ is a constant (at a given frequency) and therefore, the specific line length can indicate the angular displacement and hence the change in node voltage.

VII. EVALUATION OF PROPAGATION MODEL DENMARK-WEST
In this section, the propagation of voltage in the Denmark-West model of Fig. 4 is evaluated and discussed.For this evaluation, two study cases are used, given in Table VI.Study case 2 A is a case with VLA as an OHL section and 2B with VLA as two parallel UGC sections, both cases with a line length of 7 km for VLA.Thévenin source types with zero impedance and 1∠0 • per unit is used.Case 2 A is the reference or base case.
The study from Energinet revealed, among other things, the following [21]: 1) According to [21], the transfer gain predicts an amplification of the 11th harmonic distortion from substation TJE towards substation TRI.However, statistical analyses of harmonic voltage measurements have found no correlation between harmonic distortion levels in substations TJE and TRI [3], [30], indicating that using gain factors for prediction of harmonic propagation may lead to erroneous Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Fig. 6.Real value of voltage at nodes as consequence of propagated harmonics from surrounding substations.Σb x is propagated voltage for study cases 2 A and 2B.ΔΣb is the difference in propagated voltage at each node.ΔΣu is the difference in harmonic voltage at each node.Results are shown for all nodes in Fig. 4 according to the coloured regions.
results for meshed grids.In Section VII-B, it will be shown that the propagation model proposed by the authors confirms the measurement analyses.
2) The 11th (and 13th) harmonic voltage distortion in 400 kV substations FER, KAS, LAG, TJE, and VHA remains almost unchanged after commissioning of the 400 kV cable system at Vejle-Ådal.
3) The 11th harmonic voltage in FGD increased significantly (≈ 2.5 times) after commissioning of the Vejle-Ådal 400 kV cable system.The following assessment will address the above conclusions from [21] and explain the observed behavior by means of the propagation model.The focus is on case-specific changes to the propagation and the direct influence caused by the change in the network, which is not immediately obvious using the traditional methods of harmonic analysis, whereas the previous section discussed the physical aspects of voltage propagation generally determined by the method.

A. Propagation Model Results
Fig. 6 shows the voltage (Σb) at nodes as a consequence of propagated harmonics from surrounding substations for case 2 A and 2B, where 2 A is the base case and 2B is the augmented case.The difference in propagated voltage (ΔΣb) = (Σb 2A − Σb 2B ) and the difference in node voltages (ΔΣu) = (Σu 2A − Σu 2B ) is also shown.Δ indicates the difference between the two cases and Σ indicates it is the sum of contributions from all sources at each node.Evaluating the real part of the complex per unit voltage gives the possibility to determine the changes in propagation as real numbers instead of an angular displacement determined by the imaginary values and therefore, the real part of the per unit voltage is shown on the y-axis.The corresponding node according to Fig. 4 is shown along the x-axis.The difference in propagated voltage (ΔΣb) between nodes 20 and 108 (Region II) is zero and this is also the case for (ΔΣu).In the intervals {1, . . ., 19} (Region I) and {109, . . ., 178} (Region III), the difference in propagation and as a consequence of this, a difference in node voltage occur.As (ΔΣb) is zero in Region II, the propagation along the lines KAS-REV-ASR-TJE are unchanged for both cases which is consistent with item 2 in the list above.

B. Voltage Propagation in Region I+III
Fig. 7 shows the difference in propagated voltage from the individual sources at specific nodes {1, 11, 8, 150, 152, 114, 132} (Regions I+III).The smaller the value of (Δb) the less is the difference in propagation from that specific source for the case under evaluation.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.It is evident that the change in propagation for the source located in substation TJE has minor influence on the overall change in propagation because this curve (with diamond legend) is close to zero at all nodes, which was also concluded by [3], [21], [30] (item 1).In addition, at substation TRI, Δb is seen to be -0.03,confirming that the transfer gain amplification at TRI caused by the source in TJE does not exist.
At node 132 (VHA), (Δb) is also minor for all source contributions (−0.30, . . ., 0.13), indicating minor changes in propagation and hence a minor change in the harmonic voltage level at substation VHA, which was concluded by [21] (item 2).The opposite is seen at the source at node 1 (FGD).Here, (Δb) is major for the contributions from VHA (0.67), KAS (−0.57), and FGD (−1.46, the source at the node itself) and at nodes close to this source node, and this indicates major changes in harmonic voltage when the OHL section changes to UGC at VLA (between nodes 138 and 149), which was also concluded by [21] (item 3).
The results presented in [21] were achieved through brute force simulations, a comprehensive tuning of harmonic sources, as well as engineering experience in-depth knowledge of the network.With the propagation model presented in this paper, the same conclusions can be obtained by a simple consideration of the contribution of individual sources to the harmonic distortion.Furthermore, knowledge is gained about harmonic propagation throughout the network and how individual network changes result in a changed distribution of harmonic voltages, as determined precisely by electromagnetic wave scattering.As a result, the effect of different grid changes on harmonic propagation as a screening method can be compared with less effort than comprehensive brute force simulations.

C. Evaluation of Propagation at Substation LAG
The changes in propagation are due to replacement of the OHL section between nodes 138 and 149 with two parallel UGCs and we have seen how this has affected the magnitude of (Δb).We can also evaluate changes of the propagation coefficient ζ and in this case, we will evaluate |Δζ i:j | at the terminals of substation LAG for each source contribution.As an example, Fig. 8 shows a and b for both cases 2 A and 2B at nodes 6 and 9 for the source contribution from KAS.The magnitudes of |Δζ 6:9 | are given in Table VII.The contribution from KAS is changed at nodes 6 and 9.As |Δζ 6,KAS | < 1, then by the definition in (2), a 6 > b 6 , meaning that more voltage is transmitted through the junction towards nodes 7-9.Then |Δζ 9,KAS | > 1, meaning that b 9 > a 9 and less voltage is transmitted from node 9 towards nodes 6-8 and hence, more propagates back along the line away from 9. Δζ is: that is, the difference in propagation at node i due to the network change from case 2 A to case 2B.Using the propagation coefficient one can in general evaluate several cases against a base case, say, case 2 A, and determine the case specific changes to the propagated voltage at all nodes.This indicates how much the harmonic voltage changes and due to which sources.

VIII. DISCUSSION
The method presented in this paper can unequivocally determine the propagation of harmonic voltages along a single transmission line or in a meshed transmission network consisting of several transmission lines.It is simple to compare different study cases against each other, or against a base-case.With this method, it is possible to quantitatively determine changes in the harmonic voltages by given restructuring in the power system, as well as to determine the contribution of individual harmonic sources at nodes in the grid, evaluated broadly between individual study cases.Because of this, it is also possible to evaluate how the voltage from future harmonic sources will propagate in the meshed grid and thus draw attention to necessary mitigating measures already in an early screening process.
System component models are based on S-parameters known from microwave theory.General models of transmission lines, nodes, sources, and loads have been established.Expanding the model library with passive models for series-and parallelcoupled reactive components is simple, as these can be described by their characteristic impedance and reflection coefficient [23], [24].It will similarly be possible to describe active components, however, such model descriptions can be expected to be complex [10], [11], [12], [31], [32].To evaluate the propagation of harmonic voltage in the grid while taking active components into account, the harmonic impedance of the components for specific working points must be determined, followed by a set of reflection coefficients that can be used to evaluate the impact of the active components on propagation.It has no effect on how voltages propagate in the passive grid because the grid remains unchanged, but it does affect the magnitude of voltage at nodes, which is determined by the superposition of individual source contributions.
With the presented method, a general understanding of the physical relationship between harmonic propagation and amplification in a meshed system can be achieved.Changes in propagation are caused by changes in impedance conditions at terminal points, as well as changes in angular displacement along a line or when transitioning between two lines with different characteristic impedance.This causes altered incident and reflected voltage waves and more or less angular difference between these waves, which in turn causes a change in the voltage at the point in question.It has been shown that this change can be quantified by reference to a base-case using the propagation coefficient.
The methodology could be developed into a screening tool to quickly assess harmonic distortion when planning future projects in a power system and possibly contribute in a simpler way to better mitigation strategies.For the latter, the following steps are recommended: r Evaluate grid propagation patterns: the propagated voltage for a study case relative to a reference case is determined as |Δ b|.The preferred paths for voltage waves in the grid are thus determined.
r Evaluate source emission impact: the propagated voltage from each source is determined as |Δ b source |.The magnitude can be evaluated at each station, and problematic sources or points in the grid are determined.
r Evaluate the root cause: With the scatter parameters, a direct evaluation of the physical impact of line types and lengths is made, allowing exploration of network configurations and cable lengths.Since harmonic voltage and current in a power system can be described as stationary electromagnetic waves, it is precisely the propagation of the electromagnetic waves that is interesting and with the presented method, these waves can be determined.This also means that the method can be used in general to also determine reflection coefficients, incident and reflected waves as starting conditions in connection with, for example, transient studies.

IX. CONCLUSION
In this paper, a novel propagation model for determining and understanding the propagation of harmonic voltages in a meshed transmission system has been presented.
The model is based on the propagation of electromagnetic waves in the grid, which precisely determine the harmonic voltages, because they are stationary or standing waves.Techniques known from microwave theory are used in the form of scattering parameters, which describe the physical characteristics of the elements in the model and their behaviour towards electromagnetic waves.
The developed model can simply describe the harmonic distribution in a meshed transmission network.The model is validated against physical measurements of the distortion levels in the Denmark-West 400 kV transmission network.It is possible to evaluate changes in harmonic distortion between different study cases or in relation to a base case and, thus, quantify changes in the harmonic distortion for given topological changes in the transmission network.
The approach is different from conventional harmonic assessments, which typically use brute force simulations to assess changes in the harmonic distortion for particular nodes.The propagation model gives a clear understanding of how each source contributes to harmonic distortion, as well as the significance of the propagation when a transmission line is changed or replaced with another type.It can be both OHL and UGC.A physical understanding of the phenomenon of propagation can be achieved with this model.
The perspectives are that the method can be developed into a screening tool for TSOs, for use in future changes in a transmission network and new network expansions, and to easily assess mitigation methods such as filters.

1) Transmission Line Model:
The transmission line model is a two-port model given as a 2x2 matrix [24]: According to (1), the complex propagation constant γ determines the magnitude and angular displacement of the voltage wave along the line and is thus directly affected by the line length.As a result, the two-port S-parameter network must reflect this displacement along the transmission line, and the choice of transmission line parameters, that is, lumped or distributed, affects the two-port S-parameter model.
2) Source Model: Thévenin source models are used because the focus is on steady-state harmonic voltage propagation.Therefore, non-linearities and cross-coupling effects are neglected.A Thévenin source is also considered as a change in impedance seen from an electromagnetic wave point of view.Therefore, a two-port model can be used as a generic source model with the following network matrix: where the Thévenin voltage is applied to port 1 (that is, a 1 ) and leaves the two-port at b 2 and the reflection at port 2, determined by the relationship between the network characteristic impedance and the Thévenin impedance are designated as the complex reflection coefficient Γ S .For a zero-impedance Thévenin source Γ S := −1.This propagation model is constructed using Thévenin voltage sources, but current sources can be used.Because the concept is based on voltage waves, Norton-type current sources must be transformed into an equivalent two-port voltage source.
3) Load Model: The concept of a load model is similar to a source model: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 1 .
Fig. 1.Generic network model used for derivation of the propagation model.Arrows indicate direction of waves.Node IDs for element connection points are shown.
has input at port 1, that is, a 3 and b 3 and output at port 2 given as a 4 and b 4 .The response b 4 from A is input to the three-port network B, that is, a 5 = b 4 and the reflected wave b 5 is input at port 2 of A, i.e., a 4 = b 5 .The transmitted waves b 6 and b 7 at the output ports of B are input to the rightmost networks C and D, respectively, that is, a 8 = b 6 and a 10 = b 7 .

Fig. 4 .
Fig. 4. Developed propagation model for Denmark-West transmission system.Solid lines are OHL type and dotted lines are UGC type and the node number for each element is shown.A parallel UGC section with two junctions and two UGC line sections is represented by a special symbol.Substation names according to Fig. 3.The grid is divided into three regions I, II and III.

Fig. 5 .
Fig. 5. Three study case models.(a) -Case 1 A with three OHL sections.(b) -Case 1B with OHL followed by UGC and OHL.(c) -Case 1 C with OHL followed by two parallel UGC and OHL.Element impedance and admittance values are given.

Fig. 8 .
Fig. 8. Incident and reflected waves at nodes 6 and 9 at substation LAG for the source contribution from KAS for both cases 2 A and 2B are shown.The propagation coefficient for nodes 6 and 9 are also shown.

TABLE IV S
-PARAMETERS FOR OHL AND UGC SECTIONS IN STUDY CASE 1A-1 C

TABLE VI STUDY
CASES FOR DENMARK-WEST EVALUATION