**Using a Low Complexity Numeric Routine for Solving Electromagnetic Transient Simulations**

Rafael Cuerda Monzani, Afonso José do Prado, Sérgio Kurokawa, Luiz Fernando Bovolato and José Pissolato Filho

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/48507

## **1. Introduction**

462 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

16, 2006, Salvador. Anais... Salvador: [s.n.], 2006. v. 1 p. 268-273,

[s.n.], 2003, New Orleans. Anais New Orleans: [S.l.], 2003. 5p.

v.20, n.2, p. 894-903, 2005.

n.2, p. 200-203, 2003.

1110, 1988.

2, 792p.

1293-1297, 1989.

Kurokawa, S.; Yamanaka, F. N. R.; Prado, A. J. Representação de linhas de transmissão por meio de variáveis de estado levando em consideração o efeito da freqüência sobre os parâmetros

Kurokawa, S.; Yamanaka, F. N. R.; Prado, A. J.; Pissolato, J. Representação de linhas de transmissão por meio de variáveis de estado considerando o efeito da frequência sobre os parâmetros longitudinais. In: CONGRESSO BRASILEIRO DE AUTOMÁTICA- CBA,

Mácias, J. A. R.; Expósito A. G.; Soler, A. B. A Comparison of techniques for state- space transient analysis of transmission lines. IEEE Transactions on Power Delivery, [S.l.],

Mamis, M. S. Computation of electromagnetic transients on transmission lines with nonlinear components. IEE. Proc. General Transmission and Distribution, [S.l.], v.150,

Mamis, M. S. State-space transient analysis of single-phase transmission lines with corona. In: INTERNATIONAL CONFERENCE ON POWER SYSTEMS TRANSIENTS- IPST,

Mamis, M. S.; Nacaroglu, A. Transient voltage and current distributions on transmission lines. IEE. Proc. General Transmission and Distribution., [S.l.], v.149, n. 6, p. 705-712, 2003. Martí, L. Simulation of transients in underground cables with frequency-dependent modal transformation matrices. IEEE Transactions on Power Delivery, [S.l.] , v. 3, n.3, p.1099-

Nelms, R. M.; Sheble, G. B.; NEWTON, S. R.; GRIGSBY, L. L.; Using a personal computer to teach power system transients. IEEE Transactions on Power Systems, [S.l.] v. 4, n. 3, p.

Swokowski, E.W. Cálculo com geometria analítica. São Paulo: Ed. Makron do Brasil, 1994. v.

Tavares, M. C.; Pissolato, J.; Portela, M. C. Quasi-modes multiphase transmission line model,

Yamanaka, F. N. R.; Kurokawa, S.; Prado, A. J.; Pissolato, J.; Bovolato, L. F. Analysis of longitudinal and temporal distribution of electromagnetic waves in transmission lines by using state-variable techniques. In: SIXTH LATIN-AMERICAN CONGRESS ON ELECTRICITY GENERATION AND TRANSMISSION, 16, 2005, Mar del Plata.

Electric Power Systems Research, [S.l.],v. 49, n. 3, p. 159-167, 1999.

Proceeding... Mar del Plata: [s.n.], 2005. p. 1-7.

longitudinais. Sba Controle& Automação, Campinas, v.18, n.3, p.337- 346, 2007.

Different methods can be applied to accomplish the analysis of transmission lines. Many mathematical tools can be used, the main tools used are: circuits analysis with the use of Laplace or Fourier Transform, State Variables and Differential Equations. These tools can be included in a numeric routine in order to obtain voltage and current values in simulation of electromagnetic transients, at any point of the circuit.

The EMTP (ElectroMagnetic Transient Program) [1] is the main kind of this software. The prototype was developed in 60's by professionals of power system area led by Dr. Hermann Dommel (University of British Columbia, in Vancouver, B.C., Canada), and Dr. Scott Meyer (Bonneville Power Administration in Portland, Oregon, U.S.A.). Currently, EMTP is the basis of electromagnetic transients simulations in power systems.

With EMTP type programs, the following analysis can be done: simulation of switching and lightning surges, transient and temporary overvoltage, electrical machines, resonance phenomena, harmonics, power quality and power electronics applications. The most known programs of EMTP type are:


© 2012 Monzani et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 Monzani et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

$$
\dot{\mathfrak{x}} = A\mathfrak{x} + Bu \tag{1}
$$

$$\left(\boldsymbol{\chi}[k+1]-\boldsymbol{\chi}[k]\right) = \frac{T}{2}\left(A\boldsymbol{\chi}[k+1]+Bu[k+1]+A\boldsymbol{\chi}[k]+Bu[k]\right) \tag{2}$$

$$\mathbf{x}[k+1] = \mathbf{x}[k] + \frac{T}{2} \left( A\mathbf{x}[k+1] + Bu[k+1] + Ax[k] + Bu[k] \right) \tag{3}$$

$$
\left[I - \frac{T}{2}A\right] \mathbf{x}[k+1] = \left[I + \frac{T}{2}A\right] \mathbf{x}[k] + \frac{T}{2}B \left[u[k] + u[k+1]\right] \tag{4}
$$

$$A'x[k+1] = A'"{x[k]} + B'[u[k] + u[k+1]]\tag{5}$$

$$A' = \left[I - \frac{T}{2}A\right]$$

$$A'^{\prime} = A^{\prime} \cdot \left[I + \frac{T}{2}A\right] \tag{6}$$

$$B^{\prime} = A^{\prime} \cdot \frac{T}{2}B$$

$$A = \begin{bmatrix} -\frac{G}{C} & -\frac{2}{C} & 0 & \dots & \dots & 0\\ \frac{1}{L} & -\frac{R}{L} & -\frac{1}{C} & 0 & \ddots & \vdots\\ 0 & \frac{1}{C} & -\frac{G}{C} & -\frac{1}{L} & 0 & \vdots\\ \vdots & 0 & \ddots & \ddots & \ddots & 0\\ \vdots & \ddots & 0 & \frac{1}{L} & -\frac{R}{L} & -\frac{1}{C}\\ 0 & \dots & \dots & 0 & \frac{2}{C} & -\frac{G}{C} \end{bmatrix} \tag{7}$$

$$\frac{d\,i\_{10}}{dt} = \frac{i\_{10}}{L\_0} \left(-\sum\_{j=1}^{m} R\_j\right) + \frac{1}{L\_0} \left(\sum\_{j=1}^{m} R\_j \dot{i}\_{1j}\right) + \frac{1}{L\_0} u(t) - \frac{1}{L\_0} v\_1(t)$$

$$\frac{d\dot{i}\_{11}}{dt} = \frac{R\_1}{L\_1} \dot{i}\_{10} - \frac{R\_1}{L\_1} \dot{i}\_{11}$$

$$\frac{d\dot{i}\_{12}}{dt} = \frac{R\_2}{L\_2} \dot{i}\_{10} - \frac{R\_2}{L\_2} \dot{i}\_{12} \tag{8}$$

$$\frac{d\dot{i}\_{1m}}{dt} = \frac{R\_m}{L\_m} \dot{i}\_{10} - \frac{R\_m}{L\_m} \dot{i}\_{1m}$$

$$\frac{d\upsilon\_1(t)}{dt} = \frac{2}{C} \dot{i}\_{10} - \frac{G}{C} \upsilon\_1(t)$$

$$A = \begin{bmatrix} [A\_{11}] & [A\_{12}] & \cdots & [A\_{1n}] \\ [A\_{21}] & [A\_{22}] & \cdots & [A\_{2n}] \\ \vdots & \vdots & \ddots & \vdots \\ [A\_{n1}] & [A\_{n2}] & \cdots & [A\_{nn}] \end{bmatrix} \tag{9}$$

$$A' = \begin{bmatrix} -\frac{\sum\_{j=0}^{m} R\_j}{L\_0} & \frac{R\_1}{L\_0} & \frac{R\_2}{L\_0} & \cdots & \frac{R\_m}{L\_0} & -\frac{1}{L\_0} \\ \frac{R\_1}{L\_1} & -\frac{R\_1}{L\_1} & 0 & \cdots & 0 & 0 \\ \frac{R\_2}{L\_2} & 0 & -\frac{R\_2}{L\_2} & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & 0 & 0 \\ \frac{R\_m}{L\_m} & 0 & 0 & \cdots & -\frac{R\_m}{L\_m} & 0 \\ \frac{2}{C} & 0 & 0 & \cdots & 0 & -\frac{G}{C} \end{bmatrix} \tag{10}$$

$$\mathbf{x}^T = [\mathbf{x}\_1 \cdots \mathbf{x}\_n] \tag{11}$$

$$\boldsymbol{x}\_{k}^{T} = \begin{bmatrix} \boldsymbol{t}\_{k0} \ \boldsymbol{t}\_{k1} \ \cdots \ \boldsymbol{t}\_{km} \ \boldsymbol{v}\_{kn} \end{bmatrix} \tag{12}$$

$$R = \frac{R' \cdot d}{n} \tag{13}$$

$$L = \frac{L' \cdot d}{n}$$

$$C = \frac{C' \cdot d}{n}$$

$$G = \frac{G' \cdot d}{n}$$

In order of results, both routines, with and without frequency influence, will be shown and discussed step by step. The objective of these routines is to show how to use MatLabTM by using simple functions or commands to present the wave propagation of current or voltage to undergraduate students. The functions and commands used into routine will be explained; also the loops and conditions used into it will also be explained below the routine.

Using a Low Complexity Numeric Routine for Solving Electromagnetic Transient Simulations 471

u(2\*e(cont)-1,1) = fun;

 fun = inline(fun); u(2\*e(cont),1) = fun;

e(cont)= input ('\n Indicate the pi circuit= ');

fun = input('\n Enter the function:');

 end if (i==2)

end

end %A matrix for j=1:(2\*P-1) h = rem(j,2); if h==1

else

 end end

%A (2\*P,2\*P) A(2\*P,2\*P)=-(G/C);

%B matrix B(2\*P,1)=0; for j=1:numel(u) h = rem(j,2); if h==1 B(j,j)=1/L;

else

 end end

B(j,j)=1/C;

%Constant terms

A2=A1\*A2; B1=(T/2)\*B; B2=A1\*B1; x(2\*P,1)=0;

A1=inv(eye(2\*P)-(T/2)\*A); A2=(eye(2\*P)+(T/2)\*A);

%Iterations to solve the linear system

cont=cont+1;

 A(j,j)= -(R/L); A(j,j+1)= -(1/L); A(j+1,j)= 1/C;

 A(j,j)= -(G/C); A(j,j+1)= -(1/C); A(j+1,j)= 1/L;

The first routine doesn't consider the frequency influence.


Using a Low Complexity Numeric Routine for Solving Electromagnetic Transient Simulations 471

```
 u(2*e(cont)-1,1) = fun; 
 end 
 if (i==2) 
 e(cont)= input ('\n Indicate the pi circuit= '); 
 fun = input('\n Enter the function:'); 
 fun = inline(fun); 
 u(2*e(cont),1) = fun; 
 end 
 cont=cont+1; 
end 
%A matrix 
for j=1:(2*P-1) 
 h = rem(j,2); 
 if h==1 
 A(j,j)= -(R/L); 
 A(j,j+1)= -(1/L); 
 A(j+1,j)= 1/C; 
 else 
 A(j,j)= -(G/C); 
 A(j,j+1)= -(1/C); 
 A(j+1,j)= 1/L; 
 end 
end 
%A (2*P,2*P) 
A(2*P,2*P)=-(G/C); 
%B matrix 
B(2*P,1)=0; 
for j=1:numel(u) 
 h = rem(j,2); 
 if h==1 
 B(j,j)=1/L; 
 else 
 B(j,j)=1/C; 
 end 
end 
%Constant terms 
A1=inv(eye(2*P)-(T/2)*A); 
A2=(eye(2*P)+(T/2)*A); 
A2=A1*A2; 
B1=(T/2)*B; 
B2=A1*B1; 
x(2*P,1)=0; 
%Iterations to solve the linear system
```
470 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

The first routine doesn't consider the frequency influence.

%Clear variables and command window

%Definitions of Entrance Elements disp('Transmission Line Analyze')

T = input ('Set time step [s]: ');

%Distributed Elements

R=R\*D/P; L=L\*D/P; G=G\*D/P; C=C\*D/P;

cont=1; while (i ~= 0)

clc;

P = input ('Enter the Number of Pi Circuits = '); D = input ('Enter the line length (km) = ');

t = input ('Enter the simulation total time [s]: ');

%e = Pi circuits which will receive entries

e(cont)= input ('\n Indicate the pi circuit = ');

fun = input('\n Enter the function:');

disp('Indicate the input kind:');

 disp('1 - Voltage'); disp('2 - Current'); disp('0 - Exit'); i = input(''); if (i==1)

fun = inline(fun);

R = input ('Enter the value of distributed resistance (ohm/km) = '); L = input ('Enter the value of distributed inductance (H/km) = '); C = input ('Enter the value of distributed capacitance (F/km) = '); G = input ('Enter the value of distributed conductance (S/km) = ');

routine.

clear; clc;

Close all;

**MatLab code** 

%Close windows

In order of results, both routines, with and without frequency influence, will be shown and discussed step by step. The objective of these routines is to show how to use MatLabTM by using simple functions or commands to present the wave propagation of current or voltage to undergraduate students. The functions and commands used into routine will be explained; also the loops and conditions used into it will also be explained below the

Using a Low Complexity Numeric Routine for Solving Electromagnetic Transient Simulations 473

For mounting the A matrix a loop **for** is used because the number of interactions is known, inside the loop the command **if** is used to construct the mainly, upper and lower diagonals

Finally, a loop for is used to solve the trapezoidal rule considering the time step and the total time elapsed. The variable *y* retains in the first column the value of steps of time, the second column retains the value of voltage, and the third column retains the value of

Fig. 4 shows the previous version of the routine, in that routine the command **syms** was used. This command creates a symbol as a variable to solve this. The command **subs** must be used, in order to substitute the variable in the function with the value requested. In this case, the function entry of current or voltage source didn't need to be inserted as a string, but the user always had to insert the function by using the variable specified, what was not always done. The second purpose used the command **inline** which gets a string and converts it into a function; with this command the user can use any variable, to solve this in the trapezoidal rule. It's used the command **feval**, which gets the function and substitutes

The second routine considers the frequency influence. Almost all the steps used in this routine were the same as shown in the first one. Only the different steps done here will be

 1 – Voltage 2 - Current 0 - Exit **1** 

 Indicate the pi circuit = **1**  Enter the function: **'1'** 

the values with the time step.

%Close all graphic windows

disp('Transmission Line Analysis');

P = input ('Enter the Number of Pi Circuits = '); D = input ('Enter the line length (km) = ');

%It's defined R(1) and L(1) which represents R0 and L0 respectively,

R(1) = input ('Enter the value of distributed resistance (ohm/km) = ');

%values of R1, L1, ... Rm, Lm, will be with a plus number.

%Clear variables and command window

described.

clear all; clc;

close all;

sprintf('\n');

**MatLab code** 

as shown in Eq. 7. The same is done in the B matrix.

current in the terminal of the line transmission.


In the beginning of the routine some commands are used like:


The command **disp** shows a message in the command window for users to know what happens in the routine, that way no one thinks the routine is not working, because sometimes, it does take a lot of time to finish all the procedures, thus, a message is important to enlighten that.

The command **input** asks the user to insert a value for a variable, instead of defining a constant value inside the routine, this command makes it more interactive, so the user can put any value wanted and analyze the answer for different values inserted.

The loop **while** is used with the purpose to insert voltage or current sources as many as the user wants. The user will specify if the source is of voltage (1) or current (2), after specifying the kind of source, it shall specify the ߨ circuit that will receive the source, and finally specify the function that represents the source. A loop while is used, to mount a vector of entries. The user will get out of the loop inserting the value (0). The function of entry (voltage or current) must be inserted as a string, as explained above. At least one source for the routine works correctly is necessary, as an example, the user can insert a step function in the first circuit. Entering the following:

Indicate the input kind:

 1 – Voltage 2 - Current 0 - Exit **1**  Indicate the pi circuit = **1**  Enter the function: **'1'** 

472 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

j=1;

for q= 0: T: t u1=feval(u,q); u2=feval(u,q+T); x=A2\*x+B2\*(u1+u2);

 y(j,1)=q; y(j,2)=x(2\*P,1); y(j,3)=x(2\*P-1,1);

 j=j+1; end

pause

plot(y(:,1),y(:,2),'b')

plot(y(:,1),y(:,3),'r')

ylabel('Current [A]'); xlabel('Time [s]');

important to enlighten that.

ylabel('Voltage [kV]'); xlabel('Time [s]');

title('Voltage in the end of transmission line');

title('Current in the end of transmission line');


the first circuit. Entering the following:

Indicate the input kind:

In the beginning of the routine some commands are used like:

The command **disp** shows a message in the command window for users to know what happens in the routine, that way no one thinks the routine is not working, because sometimes, it does take a lot of time to finish all the procedures, thus, a message is

The command **input** asks the user to insert a value for a variable, instead of defining a constant value inside the routine, this command makes it more interactive, so the user can

The loop **while** is used with the purpose to insert voltage or current sources as many as the user wants. The user will specify if the source is of voltage (1) or current (2), after specifying the kind of source, it shall specify the ߨ circuit that will receive the source, and finally specify the function that represents the source. A loop while is used, to mount a vector of entries. The user will get out of the loop inserting the value (0). The function of entry (voltage or current) must be inserted as a string, as explained above. At least one source for the routine works correctly is necessary, as an example, the user can insert a step function in

put any value wanted and analyze the answer for different values inserted.


For mounting the A matrix a loop **for** is used because the number of interactions is known, inside the loop the command **if** is used to construct the mainly, upper and lower diagonals as shown in Eq. 7. The same is done in the B matrix.

Finally, a loop for is used to solve the trapezoidal rule considering the time step and the total time elapsed. The variable *y* retains in the first column the value of steps of time, the second column retains the value of voltage, and the third column retains the value of current in the terminal of the line transmission.

Fig. 4 shows the previous version of the routine, in that routine the command **syms** was used. This command creates a symbol as a variable to solve this. The command **subs** must be used, in order to substitute the variable in the function with the value requested. In this case, the function entry of current or voltage source didn't need to be inserted as a string, but the user always had to insert the function by using the variable specified, what was not always done. The second purpose used the command **inline** which gets a string and converts it into a function; with this command the user can use any variable, to solve this in the trapezoidal rule. It's used the command **feval**, which gets the function and substitutes the values with the time step.

The second routine considers the frequency influence. Almost all the steps used in this routine were the same as shown in the first one. Only the different steps done here will be described.

## **MatLab code**  %Clear variables and command window clear all; clc; %Close all graphic windows close all; disp('Transmission Line Analysis'); sprintf('\n'); P = input ('Enter the Number of Pi Circuits = '); D = input ('Enter the line length (km) = '); %It's defined R(1) and L(1) which represents R0 and L0 respectively, %values of R1, L1, ... Rm, Lm, will be with a plus number. R(1) = input ('Enter the value of distributed resistance (ohm/km) = ');

Using a Low Complexity Numeric Routine for Solving Electromagnetic Transient Simulations 475

%Mainly matrix

%Putting the matrix in all the positions

AP(j\*(m+2)+1:(j+1)\*(m+2),j\*(m+2)+1:(j+1)\*(m+2))=A;

%Putting the lower and upper matrices in their places

AP=A;

end

for j=1:(P-1)

for j=1:(P-1)

A(m+2,1)=2/C;

end

cont=1; while (i ~= 0)

clc;

 end if (i==2)

 AP(j\*(m+2),j\*(m+2)+1)=-1/C; AP(j\*(m+2)+1,j\*(m+2))=1/L(1);

%Element of the first matrix

%e = Pi circuits which will receive entries

disp('Indicate the input kind:');

%cont sets which pi circuit will receive the entry

e(cont)= input ('\n Indicate the pi circuit = ');

e(cont)= input ('\n Indicate the pi circuit = ');

fun = input('\n Enter the function:');

u(e(cont)\*(m+2)-(m+1),1) = fun;

%Element of last matrix A(P\*(m+2),1)=2/C;

 disp('1 - Voltage'); disp('2 - Current'); disp('0 - Exit'); i = input(''); if (i==1)

fun = inline(fun);

%Lower matrix AI(1,m+2)=1/L(1); AI(m+2,m+2)=0; %Upper matrix AS(m+2,1)=-1/C; AS(m+2,m+2)=0;

```
L(1) = input ('Enter the value of distributed inductance (H/km) = '); 
C = input ('Enter the value of distributed capacitance (F/km) = '); 
G = input ('Enter the value of distributed conductance (S/km) = '); 
m = input ('Enter the number of coupling circuits = '); 
for i=2:m+1 
 R(i) = input (['Enter the value of resistance of the branch ' int2str(i-1) ' [ohm/km] = ']); 
 %Distributed Resistance 
 R(i) = R(i)*D/P; 
 L(i) = input (['Enter the value of inductance of the branch ' int2str(i-1) ' [H/km] = ']); 
 %Distributed Inductance 
 L(i) = L(i)*D/P; 
end 
%Time step and total time 
T = input ('Set time step [s]: '); 
t = input ('Enter the simulation total time [s]: '); 
%Distributed Elements 
R(1)=R(1)*D/P; 
L(1)=L(1)*D/P; 
G=G*D/P; 
C=C*D/P; 
%A matrix 
A=zeros(size(P*(m+2))); 
c=0; 
for i=1:m+1 
 A(i,i)=-R(i)/L(i); 
 A(1,i)=R(i)/L(1); 
 A(i,1)=-A(i,i); 
end 
%First term as positive 
A(1,1) = -A(1,1); 
for i=1:m+1 
 %c variable will get the values to put in A(1,1), uses recall 
 c=A(1,i)+c; 
end 
A(1,1)=-c; 
%Terminal elements of matrix 
A(1,m+2)=-1/L(1); 
A(m+2,1)=1/C; 
A(m+2,m+2)=-G/C;
```

```
%Mainly matrix 
AP=A; 
%Putting the matrix in all the positions 
for j=1:(P-1) 
 AP(j*(m+2)+1:(j+1)*(m+2),j*(m+2)+1:(j+1)*(m+2))=A; 
end 
%Lower matrix 
AI(1,m+2)=1/L(1); 
AI(m+2,m+2)=0; 
%Upper matrix 
AS(m+2,1)=-1/C; 
AS(m+2,m+2)=0; 
%Putting the lower and upper matrices in their places 
for j=1:(P-1) 
 AP(j*(m+2),j*(m+2)+1)=-1/C; 
 AP(j*(m+2)+1,j*(m+2))=1/L(1); 
end 
%Element of the first matrix 
A(m+2,1)=2/C; 
%Element of last matrix 
A(P*(m+2),1)=2/C; 
%e = Pi circuits which will receive entries 
%cont sets which pi circuit will receive the entry 
cont=1; 
while (i ~= 0) 
 clc; 
 disp('Indicate the input kind:'); 
 disp('1 - Voltage'); 
 disp('2 - Current'); 
 disp('0 - Exit'); 
 i = input(''); 
 if (i==1) 
 e(cont)= input ('\n Indicate the pi circuit = '); 
 fun = input('\n Enter the function:'); 
 fun = inline(fun); 
 u(e(cont)*(m+2)-(m+1),1) = fun; 
 end 
 if (i==2) 
 e(cont)= input ('\n Indicate the pi circuit = ');
```
R(i) = input (['Enter the value of resistance of the branch ' int2str(i-1) ' [ohm/km] = ']);

L(i) = input (['Enter the value of inductance of the branch ' int2str(i-1) ' [H/km] = ']);

L(1) = input ('Enter the value of distributed inductance (H/km) = '); C = input ('Enter the value of distributed capacitance (F/km) = '); G = input ('Enter the value of distributed conductance (S/km) = ');

m = input ('Enter the number of coupling circuits = ');

t = input ('Enter the simulation total time [s]: ');

%c variable will get the values to put in A(1,1), uses recall

for i=2:m+1

R(i) = R(i)\*D/P;

L(i) = L(i)\*D/P;

end

c=0;

end

end A(1,1)=-c;

for i=1:m+1 A(i,i)=-R(i)/L(i); A(1,i)=R(i)/L(1); A(i,1)=-A(i,i);

%Distributed Resistance

%Distributed Inductance

%Time step and total time T = input ('Set time step [s]: ');

%Distributed Elements

A=zeros(size(P\*(m+2)));

%First term as positive

%Terminal elements of matrix

A(1,1) = -A(1,1); for i=1:m+1

c=A(1,i)+c;

A(1,m+2)=-1/L(1); A(m+2,1)=1/C; A(m+2,m+2)=-G/C;

R(1)=R(1)\*D/P; L(1)=L(1)\*D/P; G=G\*D/P; C=C\*D/P; %A matrix

Using a Low Complexity Numeric Routine for Solving Electromagnetic Transient Simulations 477

z(kk,2)=x(kk\*(m+2),1);

z(kk,3)=x(kk\*(m+2),1);

z(kk,4)=x(kk\*(m+2),1);

z(kk,5)=x(kk\*(m+2),1);

z(kk,6)=x(kk\*(m+2),1);

z(kk,7)=x(kk\*(m+2),1);

z(kk,8)=x(kk\*(m+2),1);

z(kk,9)=x(kk\*(m+2),1);

 end end

 end end

 end end

 end end

 end end

 end end

 end end

 end end

 if (q==60\*T) for kk = 1:P

 if (q==90\*T) for kk = 1:P

 if (q==120\*T) for kk = 1:P

 if (q==150\*T) for kk = 1:P

 if (q==180\*T) for kk = 1:P

 if (q==210\*T) for kk = 1:P

 if (q==240\*T) for kk = 1:P

if (q==260\*T)

```
 fun = input('\n Enter the function:'); 
 fun = inline(fun); 
 u(e(cont)*(m+2),1) = fun; 
 end 
 cont=cont+1; 
end 
%Clear the command window 
clc; 
disp('PROCESSING...'); 
%B matrix 
B(P*(m+2),1)=0; 
for j=1:m+2:numel(u) 
 h = rem(j,2); 
 if h==1 
 B(j,j)=1/L(1); 
 else 
 B(j,j)=1/C; 
 end 
end 
A1=inv(eye(P*(m+2))-(T/2)*AP); 
A2=(eye(P*(m+2))+(T/2)*AP); 
A2=A1*A2; 
B1=(T/2)*B; 
B2=A1*B1; 
%x matrix 
x(P*(m+2),1)=0; 
j=1; 
for q = 0: T: t 
 u1=feval(u,q); 
 u2=feval(u,q+T); 
 x=A2*x+B2*(u1+u2); 
 y(j,1)=q; 
 y(j,2)=x(P*(m+2),1); 
 y(j,4)=x(P*(m+2)/2,1); %middle of line 
 y(j,3)=x(P*(m+2)-(m+1),1); 
 y(j,5)=x((P*(m+2)/2)-(m+1),1); %middle of line 
 if (q==30*T) 
 for kk = 1:P 
 z(kk,1)=kk*D/P;
```
Using a Low Complexity Numeric Routine for Solving Electromagnetic Transient Simulations 477

```
 z(kk,2)=x(kk*(m+2),1); 
 end 
 end 
 if (q==60*T) 
 for kk = 1:P 
 z(kk,3)=x(kk*(m+2),1); 
 end 
 end 
 if (q==90*T) 
 for kk = 1:P 
 z(kk,4)=x(kk*(m+2),1); 
 end 
 end 
 if (q==120*T) 
 for kk = 1:P 
 z(kk,5)=x(kk*(m+2),1); 
 end 
 end 
 if (q==150*T) 
 for kk = 1:P 
 z(kk,6)=x(kk*(m+2),1); 
 end 
 end 
 if (q==180*T) 
 for kk = 1:P 
 z(kk,7)=x(kk*(m+2),1); 
 end 
 end 
 if (q==210*T) 
 for kk = 1:P 
 z(kk,8)=x(kk*(m+2),1); 
 end 
 end 
 if (q==240*T) 
 for kk = 1:P 
 z(kk,9)=x(kk*(m+2),1); 
 end 
 end 
 if (q==260*T)
```
476 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

fun = input('\n Enter the function:');

fun = inline(fun);

end

end

clc;

%B matrix B(P\*(m+2),1)=0; for j=1:m+2:numel(u) h = rem(j,2); if h==1

B(j,j)=1/L(1);

B(j,j)=1/C;

A2=A1\*A2; B1=(T/2)\*B; B2=A1\*B1; %x matrix x(P\*(m+2),1)=0;

for q = 0: T: t u1=feval(u,q); u2=feval(u,q+T); x=A2\*x+B2\*(u1+u2);

y(j,1)=q;

 if (q==30\*T) for kk = 1:P z(kk,1)=kk\*D/P;

y(j,2)=x(P\*(m+2),1);

y(j,3)=x(P\*(m+2)-(m+1),1);

y(j,4)=x(P\*(m+2)/2,1); %middle of line

y(j,5)=x((P\*(m+2)/2)-(m+1),1); %middle of line

j=1;

else

 end end

cont=cont+1;

u(e(cont)\*(m+2),1) = fun;

%Clear the command window

A1=inv(eye(P\*(m+2))-(T/2)\*AP); A2=(eye(P\*(m+2))+(T/2)\*AP);

disp('PROCESSING...');

Using a Low Complexity Numeric Routine for Solving Electromagnetic Transient Simulations 479

This routine uses the loop **for** in the beginning to obtain all the values series branches. The user specifies the amount of series branches and enters the resistance and inductance values for each branch. In the last loop **for**, a sequence of **if** is used to obtain the wave propagation in different time instants. This is a very important part of the routine, because it shows how the voltage or current waves propagates into the line. This representation shows to undergraduate students that a signal put in the beginning of line will not appear instantaneously in the end of line. It takes a little time, like milliseconds to arrive at the end, because the length of the line is considered, differently from a bipole, the signal put in a terminal is at the other terminal at the same time. Another observation is that, with the EMTP type programs, the user can only analyze a specific point of the circuit and in this case in a time range. The routine shows how the wave propagates into the line for different

For all simulations the following values were used: the number of π circuits was 100 and the transmission line has 10 kilometers. The resistance value was ͲǤͲͷߗȀܭ݉. The inductance value was ͳ݉ܪȀܭ݉. The capacitance values was ͳͳǤͳͳ݊ܨȀܭ݉. The conductance value was

Resistors (Ω) Inductors (mH) R0 0,026 L0 2,209 R1 1,470 L1 0,74 R2 2,354 L2 0,12 R3 20,149 L3 0,10 R4 111,111 L4 0,05 **Table 1.** Longitudinal parameters values using the routine considering the frequency influence

A simulation was made for voltage input unitary step signal. It was a step function with a unitary step after ݐ ൌ Ͳݏ. The result of this simulation can be seen in Fig. 5. the voltage output at the receipting line end is shown. By using the routine without frequency influence, from the results of Fig. 5, it is observed that there is a period time related to the propagation time of the signal through the line. So, it represents a time delay between the input signal and the output signal. After the time delay, there are oscillations associated to wave reflections on the sending and receipting line ends that compose the shown voltage output.

ͲǤͷͷߤܵȀܭ݉. The time step was ͷͲ݊ݏ and the period of simulation was ͲͲݏߤ.

xlabel('Time [s]');

title('Voltage in line path'); ylabel('Voltage [kV]'); xlabel('Length (km)');

line points in the time instant.

**4. Obtained results** 

plot(z(:,1),z(:,2),'b',z(:,1),z(:,5),'r',z(:,1),z(:,11),'g',z(:,1),z(:,15),'k')

pause

```
 for kk = 1:P 
 z(kk,10)=x(kk*(m+2),1); 
 end 
 end 
 if (q==290*T) 
 for kk = 1:P 
 z(kk,11)=x(kk*(m+2),1); 
 end 
 end 
 if (q==310*T) 
 for kk = 1:P 
 z(kk,12)=x(kk*(m+2),1); 
 end 
 end 
 if (q==330*T) 
 for kk = 1:P 
 z(kk,13)=x(kk*(m+2),1); 
 end 
 end 
 if (q==360*T) 
 for kk = 1:P 
 z(kk,14)=x(kk*(m+2),1); 
 end 
 end 
 if (q==390*T) 
 for kk = 1:P 
 z(kk,15)=x(kk*(m+2),1); 
 end 
 end 
 j=j+1; 
end 
plot(y(:,1),y(:,2),'b') 
title('Voltage in the end of transmission line'); 
ylabel('Voltage [kV]'); 
xlabel('Time [s]'); 
pause 
plot(y(:,1),y(:,3),'r') 
title('Current in the end of transmission line'); 
ylabel('Current [A]');
```
xlabel('Time [s]'); pause plot(z(:,1),z(:,2),'b',z(:,1),z(:,5),'r',z(:,1),z(:,11),'g',z(:,1),z(:,15),'k') title('Voltage in line path'); ylabel('Voltage [kV]'); xlabel('Length (km)');

This routine uses the loop **for** in the beginning to obtain all the values series branches. The user specifies the amount of series branches and enters the resistance and inductance values for each branch. In the last loop **for**, a sequence of **if** is used to obtain the wave propagation in different time instants. This is a very important part of the routine, because it shows how the voltage or current waves propagates into the line. This representation shows to undergraduate students that a signal put in the beginning of line will not appear instantaneously in the end of line. It takes a little time, like milliseconds to arrive at the end, because the length of the line is considered, differently from a bipole, the signal put in a terminal is at the other terminal at the same time. Another observation is that, with the EMTP type programs, the user can only analyze a specific point of the circuit and in this case in a time range. The routine shows how the wave propagates into the line for different line points in the time instant.

## **4. Obtained results**

478 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

for kk = 1:P

 if (q==290\*T) for kk = 1:P

 if (q==310\*T) for kk = 1:P

 if (q==330\*T) for kk = 1:P

 if (q==360\*T) for kk = 1:P

 if (q==390\*T) for kk = 1:P

plot(y(:,1),y(:,2),'b')

plot(y(:,1),y(:,3),'r')

ylabel('Current [A]');

ylabel('Voltage [kV]'); xlabel('Time [s]');

 end end

 end end

 end end

 end end

 end end

 end end j=j+1; end

pause

z(kk,10)=x(kk\*(m+2),1);

z(kk,11)=x(kk\*(m+2),1);

z(kk,12)=x(kk\*(m+2),1);

z(kk,13)=x(kk\*(m+2),1);

z(kk,14)=x(kk\*(m+2),1);

z(kk,15)=x(kk\*(m+2),1);

title('Voltage in the end of transmission line');

title('Current in the end of transmission line');

For all simulations the following values were used: the number of π circuits was 100 and the transmission line has 10 kilometers. The resistance value was ͲǤͲͷߗȀܭ݉. The inductance value was ͳ݉ܪȀܭ݉. The capacitance values was ͳͳǤͳͳ݊ܨȀܭ݉. The conductance value was ͲǤͷͷߤܵȀܭ݉. The time step was ͷͲ݊ݏ and the period of simulation was ͲͲݏߤ.


**Table 1.** Longitudinal parameters values using the routine considering the frequency influence

A simulation was made for voltage input unitary step signal. It was a step function with a unitary step after ݐ ൌ Ͳݏ. The result of this simulation can be seen in Fig. 5. the voltage output at the receipting line end is shown. By using the routine without frequency influence, from the results of Fig. 5, it is observed that there is a period time related to the propagation time of the signal through the line. So, it represents a time delay between the input signal and the output signal. After the time delay, there are oscillations associated to wave reflections on the sending and receipting line ends that compose the shown voltage output. By using the routine with frequency influence at longitudinal parameters, it is obtained Fig. 6. In this figure, it is noticed that are not so many transients as in Fig. 5, because in this case the series branches dampen the signal.

Using a Low Complexity Numeric Routine for Solving Electromagnetic Transient Simulations 481

In Figs. 7 and 8, it's possible to observe the influence of frequency in the current. Without frequency influence the signal is very disturbed. On the other hand, when it is used the routine that considers frequency influence, it is noticed that in the end of the line it should not have any current, supposing an end line opened, but, there are some transients because of the last branch. As in voltage signal, the current signal is damped because of the series branches in circuit. By using the routine of frequency influence, it can also be obtained how

**Figure 8.** Current in the end of transmission line by using routine with frequency influence.

Voltage in line path

0 1 2 3 4 5 6 7 8 9 10

Length (km)

the voltage signal goes to the path of line for the time in Fig. 9.

**Figure 9.** Voltage in line path.

0

0.2

0.4

0.6

Voltage [kV]

0.8

1

1.2

1.4

**Figure 5.** Voltage in the end of transmission line by using routine without frequency influence.

**Figure 6.** Voltage in the end of transmission line by using routine with frequency influence

**Figure 7.** Current in the end of transmission line by using routine without frequency influence.

In Figs. 7 and 8, it's possible to observe the influence of frequency in the current. Without frequency influence the signal is very disturbed. On the other hand, when it is used the routine that considers frequency influence, it is noticed that in the end of the line it should not have any current, supposing an end line opened, but, there are some transients because of the last branch. As in voltage signal, the current signal is damped because of the series branches in circuit. By using the routine of frequency influence, it can also be obtained how the voltage signal goes to the path of line for the time in Fig. 9.

**Figure 8.** Current in the end of transmission line by using routine with frequency influence.

**Figure 9.** Voltage in line path.

480 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

the series branches dampen the signal.

By using the routine with frequency influence at longitudinal parameters, it is obtained Fig. 6. In this figure, it is noticed that are not so many transients as in Fig. 5, because in this case

**Figure 5.** Voltage in the end of transmission line by using routine without frequency influence.

**Figure 6.** Voltage in the end of transmission line by using routine with frequency influence

**Figure 7.** Current in the end of transmission line by using routine without frequency influence.

Using a Low Complexity Numeric Routine for Solving Electromagnetic Transient Simulations 483

*Department of Electrical Engineering at FEIS/UNESP – The University of São Paulo State,* 

*Department of Electrical Engineering, DSCE/UNICAMP – The State University of Campinas,* 

[1] Microtran Power System Analysis Corporation, Transients Analysis Program Reference

[2] R. M. Nelms, Steven, R. Newton, G. B. Sheble, L. L. Grigsby. "Simulation of

[3] R. M. Nelms, Steven, R. Newton, G. B. Sheble, L. L. Grigsby. "Using a personal Computer to teach power system Transients", IEEE Transaction on Power System, Vol.

[4] S. R. Newton, B. B. Reid, G. B. Sheble, R. M. Nelms, L. L. GRIGSBY. "Electromagnetic

[5] H. W. Dommel, EMTP Theory Book, Department of Electrical Engineering, University

[6] J. R. Marti, "Accurate modeling of frequency-dependent transmission lines in electromagnetic transients simulations", IEEE Trans. on PAS, vol. 101, pp. 147–155,

[8] S. Kurokawa, F. N. R. Yamanaka, A. J. Prado, L. F. Bovolato, J. Pissolato, "Representação de Linhas de Transmissão por meio de variáveis de Estado Levando em Consideração o Efeito da Freqüência Sobre os Parâmetros Longitudinais", SBA. Sociedade Brasileira de Automática, Controle & Automação, v. 18, p. 337-346,

[9] S. Kurokawa, F. N. R. Yamanaka, A. J. Prado, J. Pissolato, L. F. Bovolato, "Utilização de variáveis de estado no desenvolvimento de modelos de linhas de transmissão: Inclusão do efeito da freqüência nas matrizes de estado", XIX Seminário nacional de produção e transmissão de energia elétrica (XIX SNPTEE), p. 1-8, Rio de Janeiro.

[10] S. Kurokawa, F. N. R. Yamanaka, A. J. Prado, J. Pissolato, "Representação de linhas de transmissão por meio de variáveis de estado considerando o efeito da freqüência sobre os parâmetros longitudinais", XVI Congresso Brasileiro de Automática, v. 1, p. 268-273,

[7] E. Clarke, Circuit Analysis of AC Power Systems, vol. I, Wiley, New York, 1950.

transient simulator for large-scale spacecraft power systems". IEEE.

**Author details** 

José Pissolato Filho

**6. References** 

*Brazil* 

*Brazil* 

Rafael Cuerda Monzani, Afonso José do Prado, Sérgio Kurokawa and Luiz Fernando Bovolato

Manual, Vancouver, Canada, 1992.

Transmission Line Transients", IEEE.

of British Columbia, Vancouver, 1996.

4, No. 3, August 1989.

January 1982.

2007.

2007.

Salvador, 2006.

## **5. Conclusions**

By using a mono-phase circuit representation of transmission lines associated to the state variables and the trapezoidal rule, ߨ circuits are applied for electromagnetic transient simulations. The ߨ circuits represent the mono-phase circuit of transmission lines and considering no frequency influence, they are included in a numeric routine through a sparse matrix where only three diagonal lines have non-null elements. Considering the frequency dependent longitudinal line parameters, the parameters are included in a matrix, by using the first line, the first column and the main diagonal of this matrix. This matrix is introduced in another large matrix, rounded by other matrices with one non-null element each. In the upper matrix, the parameter is introducedȂ ͳȀܥ. In the lower matrix, the parameter is introducedͳȀܮ. The obtained linear system, which is based on state variables, it is solved by using trapezoidal rule integration method. The numeric solution uses the trapezoidal rule method, which solves the numeric integration by the addition of infinitesimal trapeziums areas. This method increases the accuracy of the calculation, comparing with the Euler method for the same time step. The numeric routine is applied to a mathematical matricial program and, because of this, it is possible to simulate the propagation of electromagnetic transients on transmission lines. Using MatLabTM software, it is possible to describe the input signal through mathematical functions that are included as an initial datum of the numeric routine. It is not included in the structure of the routine, but during the routine running. By using the mentioned routine, examples of a unitary step were applied as voltage input signal and some electromagnetic transients generated from these signals are shown in this chapter. The numeric routine permits the inclusion of any function describing the voltage or current source that represents the simulated transient. The included function can be located on any point of the represented transmission line.

The shown numeric routine is simple and because of this, it can be used by undergraduate students for their first contact with electromagnetic transient phenomena, by traveling wave propagation, transmission line analyses. For the first contact with wave propagation simulations for undergraduate students, the proposed routine is an excellent tool. It is simple and its use is easy. So, by using basic concepts of traveling wave, linear systems and computing, it is possible to manipulate the numeric routine, observing the wave propagation characteristics, such as time delays, wave reflection and refraction, numeric oscillations (Gibbs' oscillations), transient oscillations. For undergraduate students, it is possible to compare the results and the modeling of the electromagnetic transients in transmission lines considering or not frequency dependent line parameters. This is possible by accessing the proposed routine numeric code. In courses related to the transmission line area, the manipulation of this code can make the course more interesting for the undergraduate students. These advantages are added to the other ones that are shown in this chapter and related to the application of the MatLabTM software. The use of this software makes the implementation of the mentioned routine extremely easy.

## **Author details**

482 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

By using a mono-phase circuit representation of transmission lines associated to the state variables and the trapezoidal rule, ߨ circuits are applied for electromagnetic transient simulations. The ߨ circuits represent the mono-phase circuit of transmission lines and considering no frequency influence, they are included in a numeric routine through a sparse matrix where only three diagonal lines have non-null elements. Considering the frequency dependent longitudinal line parameters, the parameters are included in a matrix, by using the first line, the first column and the main diagonal of this matrix. This matrix is introduced in another large matrix, rounded by other matrices with one non-null element each. In the upper matrix, the parameter is introducedȂ ͳȀܥ. In the lower matrix, the parameter is introducedͳȀܮ. The obtained linear system, which is based on state variables, it is solved by using trapezoidal rule integration method. The numeric solution uses the trapezoidal rule method, which solves the numeric integration by the addition of infinitesimal trapeziums areas. This method increases the accuracy of the calculation, comparing with the Euler method for the same time step. The numeric routine is applied to a mathematical matricial program and, because of this, it is possible to simulate the propagation of electromagnetic transients on transmission lines. Using MatLabTM software, it is possible to describe the input signal through mathematical functions that are included as an initial datum of the numeric routine. It is not included in the structure of the routine, but during the routine running. By using the mentioned routine, examples of a unitary step were applied as voltage input signal and some electromagnetic transients generated from these signals are shown in this chapter. The numeric routine permits the inclusion of any function describing the voltage or current source that represents the simulated transient. The included function can be located on any point of the represented

The shown numeric routine is simple and because of this, it can be used by undergraduate students for their first contact with electromagnetic transient phenomena, by traveling wave propagation, transmission line analyses. For the first contact with wave propagation simulations for undergraduate students, the proposed routine is an excellent tool. It is simple and its use is easy. So, by using basic concepts of traveling wave, linear systems and computing, it is possible to manipulate the numeric routine, observing the wave propagation characteristics, such as time delays, wave reflection and refraction, numeric oscillations (Gibbs' oscillations), transient oscillations. For undergraduate students, it is possible to compare the results and the modeling of the electromagnetic transients in transmission lines considering or not frequency dependent line parameters. This is possible by accessing the proposed routine numeric code. In courses related to the transmission line area, the manipulation of this code can make the course more interesting for the undergraduate students. These advantages are added to the other ones that are shown in this chapter and related to the application of the MatLabTM software. The use of this

software makes the implementation of the mentioned routine extremely easy.

**5. Conclusions** 

transmission line.

Rafael Cuerda Monzani, Afonso José do Prado, Sérgio Kurokawa and Luiz Fernando Bovolato *Department of Electrical Engineering at FEIS/UNESP – The University of São Paulo State, Brazil* 

José Pissolato Filho

*Department of Electrical Engineering, DSCE/UNICAMP – The State University of Campinas, Brazil* 

## **6. References**

	- [11] F. N. R. Yamanaka, S. Kurokawa, A. J. Prado, J. Pissolato, L. F. Bovolato, "Analysis of longitudinal and temporal distribution of electromagnetic waves in transmission lines by using state-variable techniques", Sixty Latin-American Congress: Electricity Generation and Transmission – VI CLAGTEE, Mar Del Plata, 2005.

Generation and Transmission – VI CLAGTEE, Mar Del Plata, 2005.

[11] F. N. R. Yamanaka, S. Kurokawa, A. J. Prado, J. Pissolato, L. F. Bovolato, "Analysis of longitudinal and temporal distribution of electromagnetic waves in transmission lines by using state-variable techniques", Sixty Latin-American Congress: Electricity

## *Edited by Vasilios N. Katsikis*

This excellent book represents the final part of three-volumes regarding MATLABbased applications in almost every branch of science. The book consists of 19 excellent, insightful articles and the readers will find the results very useful to their work. In particular, the book consists of three parts, the first one is devoted to mathematical methods in the applied sciences by using MATLAB, the second is devoted to MATLAB applications of general interest and the third one discusses MATLAB for educational purposes. This collection of high quality articles, refers to a large range of professional fields and can be used for science as well as for various educational purposes.

MATLAB - A Fundamental Tool for Scientific Computing and Engineering Applications -

Volume 3

MATLAB

A Fundamental Tool for Scientific Computing

and Engineering Applications - Volume 3

*Edited by Vasilios N. Katsikis*

Photo by koya79 / iStock