Modelling of non-uniform lines using rational approximation and mode revealing transformation

The authors present a methodology to improve the rational modelling of non-uniform lines (NuLs). First, they formulate a segmented chain matrix representation of the NuL, which is converted into a nodal admittance formulation with respect to the two line ends. The admittance matrix is fitted with rational functions while utilising the so-called chain matrix in combination with a mode revealing transformation matrix, to improve the identification of poorly observable modes and poles. The procedure is demonstrated for a case of a very wide river crossing where the transmission towers are over 300 m high. The results compare favourably with the ones obtained using a numerical Laplace transform.


Introduction
A large number of methods have been developed for transmission line modelling (overhead lines and cables) with inclusion of frequency-dependent and distributed parameter effects. Most of this work has focused on modelling uniform lines [1-4], i.e. lines where the electrical per-unit-length parameters can be assumed constant along the line. The uniformity permits analytical solution of the voltage/current relations along the line, which enables efficient solution by the travelling wave method. There are however several situations where the uniformity assumption is not applicable, for instance in the case of river and fjord crossings where the conductor height varies drastically along the span such as the ones in the Amazon forest [5,6]. These are usually referred as non-uniform lines (NuLs), and there is in general no analytical solution for the voltage along the line in the frequency domain.
One of the early approaches to solve NuLs was based on the use of finite differences, which implies spatial discretisation along the line [7]. Later, it was proposed to represent a NuL by a cascade of uniform lines and then use the chain matrix, i.e. the matrix transfer function of voltage and current from one segment end to the other, to obtain an equivalent nodal admittance matrix [8]. More recently, it was proposed to represent the line directly by a matrix transfer function [9] with the time-domain responses obtained via numerical integration [10]. However, it was found that spurious oscillations might occur, e.g. see Fig. 11 in [10] or [11]. Thus, there is still a need to accurately represent a NuL in a time-domain simulation framework. Another possibility to obtain an approximate analytical solution is to assume that both the series impedance Z and shunt admittance Y have a known exponential dependency as function of position x along the line, Z = Z(x), Y = Y(x). This procedure is commonly referred to as the exponential line and it has been applied to lossless lines [12], lossy lines [13] and frequency-dependent tower models [14,15].
In this work, we proceed with the approach in [8], i.e. of cascading segments of uniform lines into a NuL. Each segment is represented by the line terminal admittance matrix, which is converted into its equivalent chain matrix. The chain matrices are combined into a single chain matrix, which is finally converted back into the equivalent terminal admittance matrix with respect to the ends of the total NuL. One difficulty encountered with the approach is that the small admittance elements associated with charging currents at low frequencies tend to be lost in the (large) short-circuit currents. We therefore introduce a mode-revealing transformation (MRT) [16] to make the small eigenvalues more observable in the admittance matrix. In order to obtain a model suitable for time-domain simulation, the transformed admittance matrix is subjected to rational function model extraction using vector fitting (VF) [17][18][19], followed by passivity removal by perturbation [20]. Finally, the model is transformed back into the physical domain using the inverse transformation. A river crossing, where the span between towers is relatively small, is considered as a test case to illustrate the approach. The same configuration was used in [9,10].

Difficulties with explicit representation of segments
The electrical characteristics of a non-uniform transmission line are governed by the relations where both per unit length impedance, Z(x,s) and shunt admittance Y(x,s) are functions of both frequency s = jω and the position x along the line. With the cascading approach, the line is divided into a number of short lines where each line is assumed uniform. This allows representing each line segment by the travelling wave method, also known as method of characteristics. However, for time-domain simulations the small travel times associated with the short line lengths imply a very small time step, making the whole simulation unnecessarily time consuming and demanding in memory from the computational point of view. This is one of the reasons that a NuL remains one of the main systems for the evaluation of transients based on numerical Laplace transform (NLT). Similar difficulties result with the application of finite-difference time-domain method.

Terminal admittance calculation using chain matrix
To overcome the difficulties related to the time step length, one may instead use a lumped-parameter (delay-less) representation of each line segment and combine all segments into a single element using matrix manipulation. The procedure is implemented using nodal analysis where the nodal admittance matrix for the segments is combined into a single matrix, followed by matrix reduction.
For situations with many segments, a much more efficient approach is the use of the chain matrix representation of each segment as one avoids the need for inverting a large (global) admittance matrix. To illustrate the procedure, consider the asymmetrical line span depicted in Fig. 1. Assume that this particular line span, which is in fact a NuL given the height variation the span, can be approximated by three uniform lines, namely L 1 , L 2 and L 3 (as shown in Fig. 1). Further, assume that each uniform line has an admittance matrix Y i 21 can be defined in the phase domain by [21] where Y c i is the characteristic admittance matrix and H i is the propagation function matrix for current of line segment L i , and I is the identity matrix. The nodal admittance of each uniform line (2) is converted into a matrix transfer function relating the voltage and current vector between the two ends v i + 1 This transfer matrix, also known as the chain matrix [8], is given below for the line segment L i Repeating the process for all line segments, we can obtain the transfer matrix for the entire NuL from the matrix product From the transfer matrix in (6), it is possible to obtain an equivalent nodal admittance matrix, Y eq , as shown below As the nodal admittance above is symmetric also for NuLs, only three of the four block matrices need to be calculated since we have This procedure can even consider the effect of distributed source along the NuL as shown in [22].

Modelling using rational functions
The equivalent nodal admittance matrix Y eq (9) can now be calculated at a set of discrete frequencies (s k , Y eq (s k )) and fitted with a stable and passive rational model on pole-residue form (9). The resulting model can be included in EMT simulation programs using an equivalent circuit or via recursive convolution [1]

Rational modelling using MRT
The extraction of the model (9) is made difficult by the fact that Y eq has at low frequencies a combination of very large and very small eigenvalues. To prevent the small eigenvalues to be lost in the fitting process, we make use of the MRT introduced in [16]. MRT aims at improving the observability of this small eigenvalues by choosing a suitable transformation matrix T that preserves the physical properties of symmetry, realness, stability, causality and passivity. It can be understood as an alternative to the modal VF presented in [23,24]. The procedure is as follows. First we calculate the eigenvalues of Y eq and determine the frequency with the largest ratio between the largest and smallest eigenvalue in Y eq . For instance, if λ(Y eq ) is the eigenvalue matrix of Y eq , we calculate κ(s 0 ) as κ(s 0 ) = max s [ max (abs(λ(Y eq (s))))/ min (abs(λ(Y eq (s))))] (10) Then, the eigenvector matrix T at s 0 is rotated to minimise the imaginary part in the least squares sense and approximated by T 0 = Re(T). The MRT matrix is then obtained from the nearest orthogonal approximation of T 0 using SVD, i.e.
using Q we obtain a modified nodal admittance matrix and this new matrix is subjected to the rational approximation thus The model (13) is enforced to be passive by residue perturbation [20]. Finally, the inverse transformation is applied to the rational model to obtain the rational approximation of Y

Example: simple river crossing
To illustrate this procedure, consider the case of a simple river case described in [9,10]. The line crosses a river by a 600 m span as shown in Fig. 2 where the height of and the conductors varies between 28 m on one side to 230.4 m on the other side. It is a horizontal circuit with adjacent conductors 10 m apart. The phase conductors have a 2.54 cm radius, and internal radius 0.3645 cm and DC resistance is 6.1142·10 −5 Ω/m. The ground wires are 3/8′ extra high strength (EHS). The water resistivity is assumed to be 10 Ω·m.
To obtain the equivalent nodal admittance matrix of the 600 m span, we adopt a length of 20 m for the homogeneous (uniform) line sections. Given the significant height difference between these where q is the 'specific weight' of the conductor [9]. Thus, phase conductors present distinct sags when compared with ground wires. The constant height of conductors to be used in each segment Q i is obtained by the integration of (15), leading to where x 1 and x 0 are the x-coordinates that defines the length of the uniform line approximation. After obtaining the chain matrix for each segment, a simple multiplication allows to obtain the equivalent transfer matrix for the whole system, which is then converted to an equivalent nodal admittance matrix, as mentioned in Section 2. It is important mentioning that both Z(s) and Y(s) are calculated considering the same approach as in [25].

Direct fitting
A rational approximation of this configuration was presented in [26] and for sake of clarity is summarised here. First, we must obtain the characteristic admittance from each terminal of the NuL at infinite frequency to force the representation to be asymptotically correct, i.e.
where Y c i is the characteristic admittance using the line segments closest to each terminal, i.e. with Y c 1 as the characteristic admittance of the line segment closest to terminal #1 and Y c 2 being the one related to terminal #2. In practice, we choose a very high frequency (e.g. 100 MHz) to obtain those values. Given the distinct height of these line segments Y c 1 ≠ Y c 2 . We subject the equivalent nodal admittance Y eq to a rational approximation, Y fit , such as For the fitting, we assume the approximation to be strictly proper and a total of 699 samples from 1 Hz up to 2 MHz are considered. The samples are a combination of logarithmic and linear distribution. For the rational model, we considered 50 poles using inverse matrix norm as a weighting function. Twenty-three passivity violations were found and they were corrected by the procedure in [20]. Fig. 3a depicts the rational approximation of the equivalent nodal admittance matrix. Although a very good match is attained, the large ratio between eigenvalues caused a poor accuracy at the lower frequencies. Fig. 3b depicts the eigenvalues of the equivalent nodal admittance matrix and its rational approximation as a function of frequency. The smallest eigenvalues are identified correctly only above 10 kHz.

Fitting using MRT
To apply the MRT to (18), we assume the rational approximation to be strictly proper, i.e.
For the passivity enforcement, we consider the following rational approximation: To obtain a stable model using 50 poles it was necessary to use a high number of inner-loop iterations in the routine RPdriver.m [27] to remove the 18 passivity violations that were present in the original model. The passivity enforcement via residue perturbation was successful in obtaining a stable model; the results are depicted in Fig. 4a. There are some noticeable differences when compared with the direct fitting approach. It is observed that the fitting errors are higher in the elements of Y (Fig. 4a versus Fig. 3a), although the accuracy is acceptable over the full frequency range. On the other hand, the small eigenvalues of Y are much more accurately represented in the relative sense when modelling via MRT (Fig. 4b  vs. Fig. 3b). With MRT, the small eigenvalues are accurate down to about 100 Hz whereas the model by direct fitting gives large errors below 1 kHz for the small eigenvalues.

Time-domain simulation
The impact of the corrupted low-frequency eigenvalue can be easily observed in an open-circuit response test. Consider the circuit shown in Fig. 5 where a unit step voltage is applied at node 1 with all other terminals remaining open. For comparison, we also analyse the same configuration using the NLT [28][29][30]. It is important mentioning that the NLT response can be considered an accurate result since it is the time response of the modelled system without any rational approximation. The time-domain responses were obtained using an EMTP-like simulation built in MATLAB using recursive convolutions for the realisation of the rational model [31,32]. The results for node 4 are shown in Fig. 6. A very good match is obtained between all the approaches. For the induced voltages, the scenario is rather different. Figs. 7a and b depict the voltages at node 5, while Figs. 7c and d at note 6. For the initial time, all the approaches provided similar results. However, the direct approach starts to deviate from the NLT results as time increases. This behaviour is not found when the rational model is obtained using the MRT.
We next replace the voltage excitation in Fig. 7 with a unit step current excitation. Figs. 8 and 9 show, respectively, the voltage at the node 4, and the voltage at nodes 5 and 6. It is observed that in this case, the direct method gives substantial errors also in node 4, in addition to the induced voltages (nodes 5 and 6). The reason for the error in the node 4 voltage is mainly the inability of the direct approach to represent the line-to-earth capacitances associated with the small eigenvalues. For the equivalent rational modelling considering the MRT, such erroneous behaviour is not observed.

Discussion
NuL modelling remains a challenge for time-domain simulation with the typical approach of cascading small line segments modelled by travelling wave-type models (modal domain or phase domain). The main difficulty is that the travelling wave method requires the time step length Δt to be smaller than the line travel time τ on the considered segment, in practice no more than onefifth of τ. For the 600 m river crossing example, we used a spatial discretisation of 20 m which amounts to a travel time of τ = 67 ns with a propagation speed 300 m/µs. A travelling wave model representation would therefore dictate a time step length not exceeding Δt = τ/5 = 13.4 ns. The use of such small time step length in an EMTP-type simulation would easily result in excessive run times and memory requirements. In contrast, the usage of the presented rational model allowed us to simulate both cases using a time step of 400 ns. Even larger time steps can be used if one is not interested in high-frequency transients.
The direct fitting approach with inverse magnitude as weighting function cannot accurately capture the behaviour of small eigenvalues at lower frequencies as they are only weakly observable in the admittance matrix elements. The case of the simple river crossing illustrated well this behaviour. Despite the  The limitations associated with the direct fitting were overcome by the usage of MRT. One issue related to the use of MRT is to avoid fitting elements that are numerically close to zero, i.e. below 10 −15 or less. This is achieved by limiting the value of the (inverse) weighting function when solving the associated least squares problem as reported in [16].
The models (direct fitting and MRT) presented some small oscillations associated with the Gibbs oscillations due to the finite upper frequency limit in the rational fitting, barely observable in Fig. 8. It is important to note that the spurious oscillations associated with the time-domain modelling of a NuL as shown in [8] are of a different nature and do not arise with the rational modelling approach proposed here.
For the test case considered here, the MRT resulted in a model that provided faster passivity enforcement, i.e. it required a lesser number of interactions to be stable. Further research is needed to evaluate whether this was a particular behaviour for the case considered here or it comes from the higher accuracy that the MRT can provide.

Conclusions
A new approach has been presented for modelling NuLs which offers a combination of high accuracy and high efficiency in timedomain simulations. The method is based on efficient calculation of the terminal admittance matrix in the frequency domain by cascading line sections using the chain matrix approach. The problem of error magnification caused by small eigenvalues of the admittance matrix is overcome by the usage of a MRT prior to fitting with rational functions and subsequent passivity enforcement.
Application to a river crossing highlighted the accuracy advantages of using MRT in the fitting process. It was shown that MRT is able to mitigate the loss of accuracy of the small eigenvalues, which was shown to lead to large error magnifications in time-domain simulations with some terminal conditions. Compared with a traditional modelling approach of cascading travelling wave models, the proposed method offers higher computational efficiency by permitting the use of a larger time step length.