# Circuit based implementation of the Universal Line Model

Amjad Mouhaidali, Martin Henriksen, Leo Dalmar

► To cite this version:

Amjad Mouhaidali, Martin Henriksen, Leo Dalmar. Circuit based implementation of the Universal Line Model. Jicable, Jun 2023, Lyon, France. hal-04163571

HAL Id: hal-04163571

https://hal.science/hal-04163571

Submitted on 17 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Circuit based implementation of the Universal Line Model

Amjad Mouhaidali, Martin Henriksen, Leo Dalmar; SuperGrid Institute (Villeurbanne, France), amjad.mouhaidali@supergrid-institute.com, martin.henriksen@supergrid-institute.com, leo.dalmar@supergrid-institute.com

☑ Young Researcher (Proved full-time engineering and science university researchers and Ph.D. students under 35 YO)

## ABSTRACT

This paper presents an improved version of the Universal Line Model (ULM) integrated into a graphical programming environment GPE (Simulink for example) for time domain transient simulations. The new model is based on an electrical circuit representation of the shunt line admittance and has improved efficiency using the Rational Krylov technique. The accuracy and performance of the new model are assessed through application examples and compared to results from electromagnetic transient software.

## KEYWORDS

Electromagnetic transient, frequency dependency, transmission line model.

## INTRODUCTION

The increasing integration of renewable energy sources into the transmission grid has led to a significant transformation in the grid's infrastructure, resulting in a growing number of high voltage direct current (HVDC) links that use power. The use of these technologies is crucial to achieving the goal of climate neutrality by 2050 [1]. Manufacturers of power converters face the challenge of developing solutions for multiterminal and multivendor HVDC and medium voltage DC (MVDC) systems, which require interoperability and protection strategies for reliable operation. Therefore, accurate and reliable simulation tools and models are essential for studying electrical systems.

Various commercial software tools and simulation environments are used to simulate electrical systems, such as EMTP, PSCAD, or MATLAB/Simulink. While EMTP and PSCAD provide comprehensive cable and power grid models [2], MATLAB/Simulink is widely used by power electronics manufacturers to model converters because of its low-level control. However, MATLAB/Simulink has less reliable and detailed models of cables and power grids compared to EMTP and PSCAD.

To address this issue, this article presents an improved version of the Universal Line Model (ULM) [3], also known as the Wide-Band (WB) model, that has been integrated into the MATLAB/Simulink environment for time-domain transient simulations of electrical networks based on power electronics. The ULM is based on an electrical circuit representation of the shunt line admittance and has been enhanced using the Rational Krylov technique to achieve improved efficiency with lower-order approximations [4]. For multiconductor modeling, the characteristic admittance is represented as an equivalent two-port Y parameter system.

This article reviews the frequency-dependent multiconductor transmission line modeling in the frequency

and time domains, describes the new implementation step-by-step to reach the equivalent electrical circuit, and assesses the accuracy and numerical performance of the proposed implementation using representative application examples, comparing it to EMTP. The improved ULM combines the capability of reliable transient studies with similar performances to traditional grid simulation software.

## UNIVERSAL LINE MODEL

### Traveling wave formulation

Let us suppose a system of $n$ parallel lines of length $l$ above an infinite, perfectly conducting ground plane. Thus, the $n$ line voltages and $n$ line current in the frequency domain, in matrix form, are obtained as:

$$ \frac{dV(s)}{dx} = Z(s) \cdot I(s) \qquad (1) $$

$$ \frac{dI(s)}{dx} = Y(s) \cdot V(s) \qquad (2) $$

The per unit length series impedance $Z(s) \in \mathbb{C}^{N,N}$, and admittance $Y(s) \in \mathbb{C}^{N,N}$, matrices are given by: $Z(s) = R(s) + sL(s)$ and $Y(s) = G(s) + sC(s)$. And contains the per-unit length resistance $R(s)$, inductance $L(s)$, conductance $G(s)$ and capacitance.

<figure><img src="image_1.png" /><figcaption>Figure 1: Illustration of the current and voltage waves for a multiconductor</figcaption></figure>

Taking into consideration the terminal conditions at both ends, the currents and voltages can be expressed as forward and backward-traveling waves as shown in Figure 1. Subscripts $f$ and $b$ indicate the forward and backward traveling waves. The currents and voltages at both ends can be expressed as:

$$ I_0(s) = Y_c(s)V_0(s) - H(s)[Y_c(s)V_L(s) + I_L(s)] \qquad (3) $$

$$ I_L(s) = Y_c(s)V_L(s) - H(s)[Y_c(s)V_0(s) + I_0(s)] \qquad (4) $$

Where $Y_c(s) \in \mathbb{C}^{N,N}$ and $H(s) \in \mathbb{C}^{N,N}$ are matrices representing the characteristic admittance and the propagation function, respectively. $Y_c$ and $H$ being obtained by:
$$Y_c(s) = Z(s)^{-1} \sqrt{Z(s)Y(s)} \quad (5)$$

$$H = e^{-\sqrt{Z(s)Y(s)}l} \quad (6)$$

The backward traveling can be expressed as:

$$I_{b,0}(s) = H(s)I_{f,L}(s) = H(s)[Y_c(s)V_L(s) + I_L(s)] \quad (7)$$

$$I_L(s) = H(s)I_{f,0}(s) = H(s)[Y_c(s)V_0(s) + I_0(s)] \quad (8)$$

## Rational approximation in the ULM with RKFIT

In the ULM, $Y_c$ and $H$ are approximated using VFIT [3]. The rational approximations are used for time domain implementation in EMTP-type simulation tools [2].

A new rational approximation method has emerged recently, and it is based on the Rational Kylov fitting RKFIT [4]. This method is based on Rational Arnoldi Decomposition (RAD) for poles relocation. The used orthonormal RADs deliver a well-conditioned basis for the numerical linear algebra problems, and a more robust and faster convergent algorithm. The main difference between VFIT and RKFIT is the poles relocation method. Both methods solve a Least Square (LS) problem at each iteration, however VFIT uses a projection space represented in the partial fraction basis and RKFIT uses a projection space represented in discrete-orthogonal rational functions.

In this paper, a fitting procedure for $Y_c$ and $H$ is presented

### 1) Fitting the characteristic admittance:

The characteristic admittance $Y_c$ is easily fitted since it doesn't have resonance peaks. All the elements of $Y_c$ are fitted with the same set of poles expressed in the following form:

$$Y(s) \approx Y_{c\_fitted} = \sum_{n=1}^{N} \frac{R_n}{s - a_n} + D \quad (9)$$

Where $s = j\omega$ in rad/s, the model order $N$; the poles $a_n$ and residues $R_n$ are either real or come in complex conjugate pairs; $D$ is a constant matrix.

### 2) Fitting the propagation function:

The approximation of the propagation function is more difficult since its elements contain propagation delays. The modal contributions of $H$ are calculated and a proper delay is assigned to each mode. The delay is calculated using the optimization technique proposed in [5]. The criterion proposed in [5] is used to group repetitive or close delays. Time delays are removed from the modal components and then fitted using RKFIT. The propagation function is then fitted using the poles from the modal approximation. The only remaining unknowns are the residues, and $H$ is approximated as:

$$H(s) \approx H_{c\_fitted} = \sum_{i}^{K} \left( \sum_{m=1}^{N_i} \frac{R_{m,i}}{s - a_{m,i}} \right) e^{(-s\tau_i)} \quad (10)$$

where $K$ is the number of grouped modes, $\tau_i$ are the collapsed time delay, $R_{m,i}$ is the matrix of residues calculated for the corresponding poles $a_{m,i}$ of the grouped modes and $H_{c\_fitted}$ stands for the rational approximation of $H$ considering the grouped modes.

## Passivity

For a multi-port electrical network, the passivity is guaranteed if and only if its transfer admittance matrix $Y_n(s)$ has a purely positive real part through the whole frequency range [6]. The transfer function in terms of the cable model parameters is given by [6]. $Y_n(s)$ is positive definite if and only if its Hermitian part $Y_H(s)$ has only positive eigenvalues:

$$eig(Y_H(s)) = eig(Y_n(s) + Y_n^T(s)) \quad (11)$$

Where $Y_n^T(s)$ is the conjugate of $Y_n$.

## Advantages of the new fitting approach

The proposed fitting approach has several advantages over the fitting approach implemented in the EMT software:

* The proposed approach delivers accurate fitting results with low RMS error and guarantees convergence of the algorithm.

* The model order is automatically calculated by the RKFIT algorithm. The user only defines the RMS fitting error tolerance.

* The user requires no a priori knowledge of the pole location. Poles can be equal to the iee arithmetic representation for positive infinity, represented in Matlab by the value inf.

* An automated order reduction procedure is implemented, and it takes place if an RMS fitting error below the tolerance is achieved [7].

## TIME DOMAIN IMPLEMENTATION

Let's consider a system in the frequency domain, the input to the system is denoted $U(s)$, the output $Y(s)$ and $F(s)$ the transfer function, we can write

$$Y(s) = F(s)U(s) \quad (12)$$

In the time domain, the response is obtained via the convolution integral:

$$y(t) = \int_{-\infty}^{+\infty} f(t - \tau)u(\tau)d\tau = f(t) * u(t) \quad (13)$$

Assuming that the system is causal and the excitation $f(t) = 0$ for $t < 0$, the convolution is reduced to:

$$y(t) = \int_{0}^{t} f(t - \tau)u(\tau)d\tau \quad (14)$$

It is possible to derive a particularly efficient implementation of the convolution in equation (13) that provides the output $y(t)$ in terms of the input $u(t)$. Let's consider the transfer function $F(s)$ in the following form:

$$F(s) = \frac{R}{s - a} \quad (15)$$

Where $R$ is the residue, $a$ a pole and $s$ the angular frequency. The impulse response of the system is:

$$f(t) = L^{-1}F(s) = Re^{at}\theta(t) \quad (16)$$

So that the convolution in equation (16) simplifies to:

$$y(t) = Rx(t) \text{ with } x(t) = \int_{0}^{t} e^{a(t-\tau)}u(\tau)d\tau \quad (17)$$

Where $x(t)$ is an auxiliary state variable. Now we will consider an uniform discretization of the time axis with a fixed time step $\Delta t$, leading to discrete time point $t_k = k\Delta t$
for $k = 0, \dots, K$. We calculate $x(t)$ for $t = t_k$ and we perform the following derivation:

$$ 
\begin{aligned}
x(t_k) &= \int_{0}^{t_k} e^{a(t_k-\tau)} u(\tau) d\tau \\
&= e^{a\Delta t} x(t_{k-1}) \int_{t_{k-1}}^{t_k} e^{a(t_k-\tau)} u(\tau) d\tau
\end{aligned}
\quad (18) $$

This expression shows that the state time $t_k$ can be expressed in a recursive form, as a function of the same state at time $t_{k-1}$, plus a term that accumulates the contribution from the input $u(\tau)$, restricted to the interval $[t_k - 1, t_k]$, and suitably weighted by the term $e^{a(t_k-\tau)}$. We still must compute the integral in equation (18). Semlyen and Dabuleanu proposed in [8] a linear approximation on the input. Assuming that the input signal has slow variations in $[t_k - 1, t_k]$:

$$ u(\tau) \approx \frac{u(t_{k-1})(t_k - \tau) + u(t_k)(\tau - t_{k-1})}{\Delta t} \quad (19) $$

Substituting this approximation in equation (18):

$$ 
\begin{aligned}
x_k &= e^{a\Delta t} x_{k-1} \\
&\quad + \frac{[1 + (a\Delta t - 1)e^{a\Delta t}]u_{k-1} + (-1 - a\Delta t + e^{a\Delta t})u_k}{a^2\Delta t}
\end{aligned}
\quad (20) $$

This formulation is usually denoted as recursive convolution [8]. There are other various possible approximations depending on the discretization rule. One of the other main methods is the trapezoidal rule. All the discretization schemes can be collected in the compact form:

$$ 
\begin{cases} x_k = \alpha x_{k-1} + \beta_0 u_k + \beta_1 u_{k-1} \\ y_k = R x_k \end{cases} \quad (21) $$

Where $\alpha, \beta_{0,1}$ depend on the adopted discretization strategy:

<table><caption>Table 1: Integration coefficient formulas</caption>
  <thead>
    <tr>
      <th>Recursive convolution</th>
      <th>Trapezoidal rule</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>α = e<sup>aΔt</sup></td>
      <td>α = <sup>1 + <sup>aΔt</sup>&frasl;<sub>2</sub></sup>&frasl;<sub>1 - <sup>aΔt</sup>&frasl;<sub>2</sub></sub></td>
    </tr>
    <tr>
      <td>β<sub>0</sub> = -<sup>1</sup>&frasl;<sub>a</sub>(1 + <sup>1 - α</sup>&frasl;<sub>aΔt</sub>)</td>
      <td>β<sub>0</sub> = <sup>Δt</sup>&frasl;<sub>2</sub>(1 - <sup>aΔt</sup>&frasl;<sub>2</sub>)</td>
    </tr>
    <tr>
      <td>β<sub>1</sub> = -<sup>1</sup>&frasl;<sub>a</sub>(<sup>α - 1</sup>&frasl;<sub>aΔt</sub> - α)</td>
      <td>β<sub>0</sub> = <sup>Δt</sup>&frasl;<sub>2</sub>(1 - <sup>aΔt</sup>&frasl;<sub>2</sub>)</td>
    </tr>
  </tbody>
</table>

## Transmission line time domain solution

Let's consider the lossy transmission line with frequency-dependent characteristic admittance $Y_c(s)$ and propagation function $H(s)$, denoting their (exact) inverse Laplace transform as $y_c(t)$ and $h(t)$, respectively. Equations (7) and (8) are casted in the time domain through the following convolutions:

$$ i_0 = y_c * v_0 - h * (i_L + y_c * v_L) \quad (22) $$

$$ i_L = y_c * v_L - h * (i_0 + y_c * v_0) \quad (23) $$

It is possible to rewrite $i_0$ as:

$$ i_0 = i_{sh,0} - i_{b,0} \quad (24) $$

Where:

$$ i_{sh,0} = y_c * v_0 \quad (25) $$

$$ i_{b,0} = h * i_{f,L} = h * (i_L + y_c * v_L) \quad (26) $$

By exchanging the subscripts 0 and L in equations (15) and (16), the corresponding set of equations for $i_L$ can be obtained. The time domain equations are interfaced with Electromagnetic Transient (EMT) simulation software using the Norton equivalent in Figure 2.

<figure><img src="image_2.png" /><figcaption>Figure 2: Norton equivalent of a transmission system Time domain approach</figcaption></figure>

For a un transmission line with a core and screen conductors, the current and voltages at one end are expressed as follows

$$ \begin{bmatrix} I_{0core} \\ I_{0screen} \end{bmatrix} = \begin{bmatrix} Y_{c_{1,1}} & Y_{c_{1,2}} \\ Y_{c_{1,2}} & Y_{c_{2,2}} \end{bmatrix} \begin{bmatrix} V_{0core} \\ V_{0screen} \end{bmatrix} + \begin{bmatrix} I_{b,0core} \\ I_{b,0screen} \end{bmatrix} \quad (27) $$

Where

$$ I_{sh0core} = Y_{c_{1,1}} V_{0core} + Y_{c_{1,2}} V_{0screen} \quad (28) $$

$$ I_{sh0screen} = Y_{c_{1,2}} V_{0core} + Y_{c_{2,2}} V_{0screen} \quad (29) $$

Equation (27) is the matrix representation of a two port network, the equivalent electrical circuit is shown in Figure 3.

<figure><img src="image_3.png" /><figcaption>Figure 3: Two port representation of the transmission line sending end</figcaption></figure>

If a passive two-port can be expressed with admittance parameters and the admittance matrices are symmetrical by definition, then that two port is equivalent to a PI representation as shown in Figure 4.
<figure><img src="image_4.png" /><figcaption>Figure 4: PI representation of the transmission line sending end</figcaption></figure>

Where:

$$Y_x = \frac{1}{(Y_{c1,1} + Y_{c1,2})} \qquad (30)$$

$$Y_x = \frac{1}{(Y_{c2,2} + Y_{c1,2})} \qquad (31)$$

$$Y_Z = \frac{1}{Y_{c1,2}} \qquad (32)$$

The full circuit of a transmission line with a core and screen conductors is represented in Figure 5.

<figure><img src="image_5.png" /><figcaption>Figure 5: Transmission line implementation in Simulink</figcaption></figure>

## NUMERICAL RESULTS

The proposed implementation is used with two underground transmission lines and the results are compared with the corresponding results from a commercial transient simulation software (EMTP).

In all simulations, the characteristic admittance and the propagation function matrices are calculated for the frequency range from 0,01Hz up to 1MHz with a sampling rate of logarithmically spaced 10 points/ decade.

Let's consider 2 HVDC cables installed as shown in Figure 6. The geometrical and materials proprieties of the cables are listed in Table 2.

The transient simulation is done in Matlab/Simulink. The short circuited cable end is presented in Figure 7, where the first core is energized with a unit step at t = 0.

<figure><img src="image_6.png" /><figcaption>Figure 6: Cross section of two underground cable</figcaption></figure>

<table><caption>Table 2: Geometrical and material parameters of the HVDC cables</caption><tr><td>Core</td><td>r = 29.95mm, ρ = 1.78e-8Ωm</td></tr></table>

<table><tr><td>Inner semiconductor</td><td>t = 1.5mm, ρ = 1Ωm</td></tr><tr><td>Insulation</td><td>t = 18.7mm, ϵ = 2.5</td></tr><tr><td>Outer semiconductor</td><td>t = 1.5mm, ρ = 1Ωm</td></tr><tr><td>Metallic sheath</td><td>t = 4.3mm, ρ = 2.84e-8Ωm</td></tr><tr><td>Jacket</td><td>t = 6mm, ϵ = 2.5</td></tr></table>

<figure><img src="image_7.png" /><figcaption>Figure 7: Short circuited cable end</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 8: Short circuit current after energization</figcaption></figure>
Figure 8 shows the short circuit current at the sending and receiving ends of the core and screen conductors of the fist cable. The results from the developed wide Band model are compared with the results from EMTP.

<figure><img src="image_9.png" /><figcaption>Figure 9: cable energization by step voltage with opened end</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 10: Unit step voltage response at the cable receiving end</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 11: Unit step voltage response at the cable receiving end</figcaption></figure>

Figure 9 shows the cable energization circuit with opened end with a step voltage of 1V. Figure 10 and 11 compares the results from the wide band model and EMTP.

## CONCLUSIONS

In conclusion, this paper has presented an efficient and accurate version of the Universal Line Model (ULM) integrated into a graphical programming environment. The new model is based on an electrical circuit representation

of the shunt line admittance and utilizes the Rational Krylov technique to improve its efficiency.

The application examples presented in this paper demonstrate the model's performance and its ability to provide results comparable to those obtained from electromagnetic transient software.

The improved ULM presented in this paper is a valuable contribution to the field of time domain transient simulations, providing an effective tool for engineers and researchers working in this area.

## REFERENCES

[1] Jovcic, D. (2019). *High voltage direct current transmission: converters, systems and DC grids*. John Wiley & Sons.

[2] Dommel, H. W. (1986). *Electromagnetic Transients Program Reference Manual:(EMTP) Theory Book*. Bonneville Power Administration.

[3] Morched, A., Gustavsen, B., & Tartibi, M. (1999). A universal model for accurate calculation of electromagnetic transients on overhead lines and underground cables. *IEEE Transactions on Power Delivery*, *14*(3), 1032-1038.

[4] Mouhaidali, A., Tromeur-Dervout, D., Chadebec, O., Guichon, J. M., & Silvant, S. (2020). Electromagnetic Transient Analysis of Transmission line based on rational Krylov approximation. *IEEE Transactions on Power Delivery*, *36*(5), 2913-2920.

[5] Gustavsen, B. (2004, May). Time delay identification for transmission line modeling. In *Proceedings. 8th IEEE Workshop on Signal Propagation on Interconnects* (pp. 103-106). IEEE.

[6] Grivet-Talocia, S., & Gustavsen, B. (2015). *Passive macromodeling: Theory and applications*. John Wiley & Sons.

[7] Berljafa, M., & Güttel, S. (2015). Generalized rational Krylov decompositions with an application to rational approximation. *SIAM Journal on Matrix Analysis and Applications*, *36*(2), 894-916.

[8] Semlyen, A., & Dabuleanu, A. (1975). Fast and accurate switching transient calculations on transmission lines with ground return using recursive convolutions. *IEEE Transactions on Power Apparatus and Systems*, *94*(2), 561-571.