**Appendix**

#### CODE EXECUTION

The following code provides the line model described in the paper and it is embedded into an application example. It simulates the simultaneous energizing of a 150 km long aerial line. At the source side the three voltage sources have a 600 Thevenin impedance. The program asks for the type of source (unit step or three phase sinusoids). At the load end the line is open. Figure 14 shows the geometry of the simulated line. Figure 15 shows the sending and receiving voltages for the unit step source, while Figure 16 shows the sending and receiving voltages for the sinusoidal source.

Note at Figure 15 that waveforms for phases A and C are equal and their plots are superposed. This is because the symmetry of the line and the excitation.

#### Main program

288 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

Proper design and operation of present-day power systems and apparatuses each time require accurate simulations of their electromagnetic transient performance. An important aspect of these simulations is the realistic representation of transmission lines by digital computer models. The ULM is the most general line model available today, mostly with EMTP–type programs. By being of relatively recent creation, this model still is a subject for substantial improvements in accuracy, stability and computational efficiency. It has been postulated in this work that, both, researchers and power system analysts will benefit considerably from the full understanding of the theoretical basis of the ULM, as well as from counting with a ULM–type code that is easy to understand and modify. It has been contended also that the best way to carry out ULM research and development is by providing a model version in an interpretive environment and Matlab has been the platform chosen for this. This chapter provides a comprehensive description of the theoretical basis of ULM, phase domain line model. In addition to this, full description of a ULM prototype in Matlab has been provided here, along with executable code and typical application

The following code provides the line model described in the paper and it is embedded into an application example. It simulates the simultaneous energizing of a 150 km long aerial line. At the source side the three voltage sources have a 600 Thevenin impedance. The program asks for the type of source (unit step or three phase sinusoids). At the load end the line is open. Figure 14 shows the geometry of the simulated line. Figure 15 shows the sending and receiving voltages for the unit step source, while Figure 16 shows the sending

Note at Figure 15 that waveforms for phases A and C are equal and their plots are

superposed. This is because the symmetry of the line and the excitation.

**8. Conclusions** 

examples.

**Author details** 

Jose Luis Naredo

**Appendix** 

CODE EXECUTION

Octavio Ramos-Leaños

*Cinvestav-Guadalajara, Mexico* 

Jose Alberto Gutierrez-Robles *University of Guadalajara, Mexico* 

*École Polytechnique de Montréal, Canada* 

and receiving voltages for the sinusoidal source.

clear all clc % m file to set line data LineData % Per unit length parameters [Zg,Zt,Zc,Yg,ZL,YL] = LineParameters(Mu,Eo,Rsu,Geom,Ncon,Ns,w); % Modal Parameters for k=1:Ns [Yc(:,:,k),Vm(k,:),Hm(k,:)] = ABYZLM(ZL(:,:,k),YL(:,:,k),lenght,w(k)); end % Characteristic Admittance Fitting [YcPoles,YcResidues,YcConstant,YcProportional] = YcFit(Yc,f,Ns,Ncon); % Hk fit [HkPoles,HkResidues,HkConstant,HkProportional,md] = HkFitTrace(Hm,Vm,ZL,YL,f,Ns,lenght,Ncon); % m file to execute simulation loop. SimulationLoop

#### Code to Load LineData

```
% Line Geometry
% column 1—conductor number
% column 2-- x position of each cond in m
% column 3-- y position of each cod in m
% column 4-- radii of each conductor
% column 5-- number of conductor in bundle
% column 6-- distance between conductors in bundle
% column 7—conductor resistivity 
% column 8—conductor relative permitivity
% column 9-- line lenght in m
Geom=[1 0 20 0.0153 3 0.4 2.826e-8 1e3 150e3
 2 10 20 0.0153 3 0.4 2.826e-8 1e3 150e3
 3 20 20 0.0153 3 0.4 2.826e-8 1e3 150e3]; 
 lenght = Geom(1,9); % Line lenght
Ncon = Geom(max(Geom(:,1)),1); % # of cond
Rsu = 100; % Earth resistivity Ohm-m
Mu = 4*pi*1E-7; % Henry's/meters
Eo = (1/(36*pi))*1E-9; % Farads/meters
Rhi = 9.09E-7; % Ohm-m resistivity of the iron.
Ral = 2.61E-8; % Ohm-m res of the aluminum.
Rhg = 2.71E-7; % Ohm-m res of the sky wires.
Ns = 500; % Number of samples
f = logspace(-2, 6, Ns); % Vector of log spaced Frequencies
w = 2*pi*f; % Vector of freqs in radian/sec.
```
#### Function LineParameters

function [Zg,Zt,Zc,Yg,ZT,YT]=LineParameters (Mu,Eo,Rsu,Geom,Ncon,Ns,w) % Function to compute the distances between conductor [Dij,dij,hij]=Height(Geom);

```
Zg = zeros(Ncon,Ncon,Ns); 
Zt = zeros(Ncon,Ncon,Ns); 
Zc = zeros(Ncon,Ncon,Ns); 
Yg = zeros(Ncon,Ncon,Ns); 
Zcd = zeros(Ncon,Ns); 
Zaf = zeros(Ncon,Ns); 
P = (1./sqrt(j*w*Mu/Rsu)); % Complex depth
Pmatrix = log(Dij./dij); % Potential Coeff. Matrix 
Pinv = inv(Pmatrix); % Inverse of Pmatrix
% Loop to compute matrices at all frequencies 
for kl = 1:Ns 
 % Geometric impedance
 Zg(:,:,kl) = (j*w(kl)*Mu/(2*pi))*Pmatrix; 
 % Earth impedance
 for km = 1:Ncon 
 for kn = 1:Ncon 
 if km == kn 
 Zt(km,km,kl) = (j*w(kl)*Mu/(2*pi))* 
 log(1+P(kl)./(0.5*hij(km,km))); 
 else
 num = hij(km,kn)^2 + 4*P(kl)*hij(km,kn) + 
 4*P(kl)^2 + dij(km,kn)^2; 
 den = hij(km,kn)^2 + dij(km,kn)^2; 
 Zt(km,kn,kl) = (j*w(kl)*Mu/(4*pi))* 
 log(num/den); 
 end
 end
 end
 % Geometric admittance
 Yg(:,:,kl) = (j*w(kl)*2*pi*Eo)*Pinv; 
end
% Conductor impedance
for kd = 1:Ncon; 
 Rcon = Geom(kd,4); % conductor radii in m.
 Nhaz = Geom(kd,5); % # of conductor in bundle
 Rpha = Geom(kd,7); % Resistivity in Ohm-m.
 Zcd(kd,:) = (1/Nhaz)*Rpha./(pi.*Rcon.^2); 
 Zaf(kd,:) = (1/Nhaz)*(1+j).*(1./(2.*pi.*Rcon)) .* 
 sqrt(0.5.*w.*Mu.*Rpha); 
 Zc(kd,kd,:) = sqrt(Zcd(kd,:).^2 + Zaf(kd,:).^2); 
end
% Outputs
ZT = Zg + Zt + Zc ; % Total impedance
YT = Yg ; % Total admittance
```

```
Function ABYZLM
```
function [Yc,Vm,Hmo] = ABYZLM(Z,Y,Lo,w) [M, Lm] = eig(Y\*Z); % Eigenvalues of YZ Minv = inv(M); % Inverse of eigenvectors matrix Yc = inv(Z)\*(M\*sqrt(Lm)\*Minv); % Characteristic Admittance

```
Gamma = sqrt(diag(Lm)); % Propagation constants.
Vm = w./imag(Gamma); % Modal Velocities
Hmo = diag(expm(-sqrtm(Lm)*Lo)); % Modal propag. Matrix H
```
#### Function YcFit

290 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

Zg = zeros(Ncon,Ncon,Ns); Zt = zeros(Ncon,Ncon,Ns); Zc = zeros(Ncon,Ncon,Ns); Yg = zeros(Ncon,Ncon,Ns); Zcd = zeros(Ncon,Ns); Zaf = zeros(Ncon,Ns);

for kl = 1:Ns

else

 end end end

end

end % Outputs

% Geometric impedance

 % Earth impedance for km = 1:Ncon for kn = 1:Ncon if km == kn

log(num/den);

% Geometric admittance

% Conductor impedance for kd = 1:Ncon;

sqrt(0.5.\*w.\*Mu.\*Rpha);

Function ABYZLM

Yg(:,:,kl) = (j\*w(kl)\*2\*pi\*Eo)\*Pinv;

 Rcon = Geom(kd,4); % conductor radii in m. Nhaz = Geom(kd,5); % # of conductor in bundle Rpha = Geom(kd,7); % Resistivity in Ohm-m. Zcd(kd,:) = (1/Nhaz)\*Rpha./(pi.\*Rcon.^2); Zaf(kd,:) = (1/Nhaz)\*(1+j).\*(1./(2.\*pi.\*Rcon)) .\*

Zc(kd,kd,:) = sqrt(Zcd(kd,:).^2 + Zaf(kd,:).^2);

function [Yc,Vm,Hmo] = ABYZLM(Z,Y,Lo,w) [M, Lm] = eig(Y\*Z); % Eigenvalues of YZ Minv = inv(M); % Inverse of eigenvectors matrix

Yc = inv(Z)\*(M\*sqrt(Lm)\*Minv); % Characteristic Admittance

ZT = Zg + Zt + Zc ; % Total impedance YT = Yg ; % Total admittance

P = (1./sqrt(j\*w\*Mu/Rsu)); % Complex depth Pmatrix = log(Dij./dij); % Potential Coeff. Matrix Pinv = inv(Pmatrix); % Inverse of Pmatrix % Loop to compute matrices at all frequencies

Zg(:,:,kl) = (j\*w(kl)\*Mu/(2\*pi))\*Pmatrix;

 Zt(km,km,kl) = (j\*w(kl)\*Mu/(2\*pi))\* log(1+P(kl)./(0.5\*hij(km,km)));

 den = hij(km,kn)^2 + dij(km,kn)^2; Zt(km,kn,kl) = (j\*w(kl)\*Mu/(4\*pi))\*

4\*P(kl)^2 + dij(km,kn)^2;

num = hij(km,kn)^2 + 4\*P(kl)\*hij(km,kn) +

```
function [YcPoles,YcResidues,YcConstant, YcProportional]=YcFit(Yc,f,Ns,Ncon) 
% Trace of characteristic admittance matrix
for k = 1:Ns 
 Ytrace(k,1) = trace(Yc(:,:,k)); 
end
Npol = 6; % Number of poles
[Ps] = InitialPoles(f,Npol); % Set initial poles
s = j*2*pi*f.'; % Vector of values of variable "s"
Ka=2; % 1.-Strictly proper, 2.-Proper, 3.-Improper
for khg=1:20 
 % Fit the trace of Yc (Poles)
 [YcPoles]=Poles(Ytrace.',s,Ps,Ns,Ka); 
 Ps=YcPoles; 
end
% Residues and constant term for Yc from poles of trace of Yc
for k = 1:Ncon 
 for l = 1:Ncon 
 Hs(:,1) = Yc(k,l,:); % k-l term of admittance
 [C,D,E]=Residue(Hs.',s,YcPoles,Ns,Ka); 
 YcResidues(k,l,:) = C; % k-l residues term 
 YcConstant(k,l) = D; % k-l constant term
 YcProportional(k,l)=E; %k-l proportional term
 end
end
```
#### Function HkFitTrace

```
function [HkPoles,HkResidues,HkConstant, HkProportional,md]=HkFit(Hm,Vm,ZL,YL,f,Ns, 
lenght,Ncon); 
% Minimum phase of each mode
md = ModeDelay(Hm.',f,lenght,Vm.',Ns,Ncon); 
% Computing Idempotents 
for k=1:Ns 
 % Function to calculate Idempotents of Y*Z
 [Hk] = HmIdem(ZL(:,:,k),YL(:,:,k),lenght,f(k), md,Hm(k,:)); 
 HkIdem(:,:,:,k) = Hk; % Idempotents
end
for m = 1:3 
 for k=1:Ns 
 TraceHk(m,k) = trace(HkIdem(:,:,m,k)); 
 end
end
s = j*2*pi*f.'; % Vector of the variable "s"
Ka =1;%1.-Strictly proper, 2.-Proper, 3.-Improper
Npol = 5; % Number of poles
[Ps] = InitialPoles(f,Npol); % Set the initial poles
```

```
for m = 1:3 
 Hk = TraceHk(m,:); 
 for khg=1:10 
 [HkPol]=Poles(Hk,s,Ps,Ns,Ka); 
 Ps=HkPol; 
 end
 HkPoles(m,:)=Ps; 
end
% Residues for Idempotent matrices of
% Hm from the poles of each trace.
for m = 1:3 
 for k = 1:Ncon 
 for l = 1:Ncon 
 Hs(:,1) = HkIdem(k,l,m,:); % k-l term
 [C,D,E]=Residue(Hs.',s,HkPoles(m,:),Ns,Ka); 
 HkResidues(k,l,m,:) = C; % k-l-m term
 HkConstant(k,l,m) = D; % k-l-m constant
 HkProportional(k,l,m) = E; % k-l-m prop
 end
 end
end
```
#### SimulationLoop program

Ts = 0.016; % Observation time Nt = fix(Ts/Dt); % Number of time steps t1 = fix(md./Dt); % Delay of H in samples t0 = fix(max(md)./Dt); % Maximum time delay as expressed in % number of samples t = (0:Dt:(Nt+t1)\*Dt-Dt); % Vector of time Ks = menu('CHOOSE THE TYPE OF INPUT SOURCE' , '1 -unit step' , '2 -sinusoidal'); if Ks == 1 % unit step source Isr = ones(Ncon,Nt+t0); elseif Ks ==2 % sinusoidal source Isr(1,:) = sin(337\*t); Isr(2,:) = sin(337\*t+2\*pi/3); Isr(3,:) = sin(337\*t+4\*pi/3); end NpYc = length(YcPoles); % Number of poles of Yc NpH = length(HkPoles); % Number of poles for the first % Idempotent matrix Ng = 3; %Number of groups % Initialize the states for both nodes ZA = zeros(Ncon,NpYc); % State variables ZB = zeros(Ncon,NpYc); % State variables YA = zeros(Ncon,NpH,Ng); % State variables YB = zeros(Ncon,NpH,Ng); % State variables IfarA = zeros(Ncon,t0+3); % Current at node A IfarB = zeros(Ncon,t0+3); % Current at node B VO = zeros(Ncon,1); % Voltage at node A Vi = zeros(Ncon,Nt+t0); % Voltage at node A

292 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

for m = 1:3

Ps=HkPol;

for m = 1:3 for k = 1:Ncon for l = 1:Ncon

 end end end

HkPoles(m,:)=Ps;

end

end

 Hk = TraceHk(m,:); for khg=1:10

[HkPol]=Poles(Hk,s,Ps,Ns,Ka);

% Residues for Idempotent matrices of % Hm from the poles of each trace.

 Hs(:,1) = HkIdem(k,l,m,:); % k-l term [C,D,E]=Residue(Hs.',s,HkPoles(m,:),Ns,Ka); HkResidues(k,l,m,:) = C; % k-l-m term HkConstant(k,l,m) = D; % k-l-m constant HkProportional(k,l,m) = E; % k-l-m prop

SimulationLoop program

if Ks == 1 % unit step source Isr = ones(Ncon,Nt+t0); elseif Ks ==2 % sinusoidal source

% number of samples

 Isr(1,:) = sin(337\*t); Isr(2,:) = sin(337\*t+2\*pi/3); Isr(3,:) = sin(337\*t+4\*pi/3);

% Idempotent matrix

end

Ts = 0.016; % Observation time Nt = fix(Ts/Dt); % Number of time steps t1 = fix(md./Dt); % Delay of H in samples

t = (0:Dt:(Nt+t1)\*Dt-Dt); % Vector of time

NpYc = length(YcPoles); % Number of poles of Yc NpH = length(HkPoles); % Number of poles for the first

Ng = 3; %Number of groups % Initialize the states for both nodes ZA = zeros(Ncon,NpYc); % State variables ZB = zeros(Ncon,NpYc); % State variables YA = zeros(Ncon,NpH,Ng); % State variables YB = zeros(Ncon,NpH,Ng); % State variables IfarA = zeros(Ncon,t0+3); % Current at node A IfarB = zeros(Ncon,t0+3); % Current at node B VO = zeros(Ncon,1); % Voltage at node A Vi = zeros(Ncon,Nt+t0); % Voltage at node A

t0 = fix(max(md)./Dt); % Maximum time delay as expressed in

Ks = menu('CHOOSE THE TYPE OF INPUT SOURCE' , '1 -unit step' , '2 -sinusoidal');

```
VL = zeros(Ncon,1); % Voltage at node B
Vf = zeros(Ncon,Nt+t0); % Voltage at node B
IO = zeros(Ncon,1); % Current at node A
Ii = zeros(Ncon,Nt+t0); % Current at node A
IL = zeros(Ncon,1); % Current at node B
If = zeros(Ncon,Nt+t0); % Current at node B
Iri = zeros(Ncon,Nt+t0); % Current at Y source
Irf = zeros(Ncon,Nt+t0); % Current at Y charge
IfarAint = zeros(Ncon,Ng); % Current at node A
IfarBint = zeros(Ncon,Ng); % Current at node B
% Constants for the state ZA and ZB
Ai(:,1) = (1+(Dt/2)*YcPoles)./(1-(Dt/2)*YcPoles); 
Au(:,1) = ((Dt/2)./(1-(Dt/2)*YcPoles)); 
Bi(:,1) = (Ai+1).*Au; 
Gy = zeros(Ncon,Ncon); 
for nm = 1:NpYc 
 Di(:,:,nm) = YcResidues(:,:,nm)*Bi(nm); 
 Gy = Gy + YcResidues(:,:,nm)*Au(nm); 
end
% Constants for the states YA and YB
for k = 1:Ng 
 K1(:,k) = (1+(Dt/2)*HkPoles(:,k))./(1-(Dt/2)*HkPoles(:,k)); 
 Ka(:,k) = (((Dt/2))./(1-(Dt/2)*HkPoles(:,k))); 
 Ku(:,k) = (K1(:,k)+1).*Ka(:,k); 
end
for k = 1:Ng 
 for nm = 1:NpH 
 K2(:,:,nm,k) = HkResidues(:,:,nm,k).*Ka(nm,k); 
 K3(:,:,nm,k) = HkResidues(:,:,nm,k).*Ku(nm,k); 
 end
end
Gy = Gy + YcConstant; % Admitance of the Ish
Yi = diag(eye(3)*[1/600; 1/600; 1/600]); % Admittance of the source, connected at node A
Gys = inv(Gy + Yi); % Impedance to calculate VO
Yr =diag(eye(3)*[1/1e6; 1/1e6; 1/1e6]); % Admittance of load connected at node B
Gyr = inv(Gy + Yr); % Impedance to calculate VL
% Contants terms to perform the interpolation
tm =md - t1*Dt; % Time for the interpolation
% Linear interpolation constants
c1 = tm/Dt; 
c2 = 1-c1; 
c3 = ones(Ng,1); 
% Pointers for the interpolation and the buffer
h1 = t1+1; 
h2 = t1+2; 
h3 = t1+3; 
for k = t0+2:Nt+t0-3 
 IfarA(:,1) = IL + Gy*VL + sum(ZB(:,:),2); 
 IfarB(:,1) = IO + Gy*VO + sum(ZA(:,:),2); 
 % Linear interpolation
```

```
 for m = 1:Ng 
 IfarAint(:,m) = c2(m)*IfarA(:,t1(m)) + c3(m)*IfarA(:,h1(m)) + c1(m)*IfarA(:,h2(m)); 
 IfarBint(:,m) = c2(m)*IfarB(:,t1(m)) + c3(m)*IfarB(:,h1(m)) + c1(m)*IfarB(:,h2(m)); 
 end
 IfarA(:,2:h3) = IfarA(:,1:h2); 
 IfarB(:,2:h3) = IfarB(:,1:h2); 
 for m = 1:NpYc 
 ZA(:,m) = Ai(m)*ZA(:,m) + Di(:,:,m)*VO; 
 ZB(:,m) = Ai(m)*ZB(:,m) + Di(:,:,m)*VL; 
 end
 for l = 1:Ng 
 for m = 1:NpH 
 YA(:,m,l) = K1(m,l)*YA(:,m,l) + K2(:,:,m,l)*IfarAint(:,l); 
 YB(:,m,l) = K1(m,l)*YB(:,m,l) + K2(:,:,m,l)*IfarBint(:,l); 
 end
 end
 HistO = - sum(ZA(:,:),2) + sum(sum(YA(:,:,:),3),2); 
 HistL = - sum(ZB(:,:),2) + sum(sum(YB(:,:,:),3),2); 
 VO = Gys*(Isr(:,k)+HistO); 
 VL = Gyr*HistL; 
 IO = Gy*VO - HistO; 
 IL = Gy*VL - HistL; 
 Vi(:,k) = VO; 
 Vf(:,k) = VL; 
 Ii(:,k) = IO; 
 If(:,k) = IL; 
end
Iri = Yi*Vi; 
Irf = Yr*Vf; 
vt = (0:Dt:length(Vi(1,:))*Dt-(t0+4)*Dt)'; 
N = length(vt); 
a1 = t1+1; 
a2 = Nt+t1-3; 
figure(1),plot(vt,Vi(:,a1:a2),':',vt,Vf(:,a1:a2)) 
ylabel('Amplitude in volts') 
xlabel('Time in seconds') 
legend('Sending end phase A' , 'Sending end phase B' , 'Sending end phase C' , 'Receiving end phase A' , 
'Receiving end phase B' , 'Receiving end phase C')
```
#### Function Height

function[Dij,dij,hij]=Height(Geom) Ls = Geom(max(Geom(:,1)),1); Req = zeros(Ls,1); % Equivalent bundle radii k4 = sqrt(2\*(Geom(:,6)/2).^2); for nc = 1: Ls; if Geom(nc,5)==1 Req(nc) = Geom(nc,4); else Req(nc) = (Geom(nc,4).\*Geom(nc,5).\*k4(nc).^

```
(Geom(nc,5)-1)).^(1./Geom(nc,5)); 
 end
end
% Direct and image distances among conductors
for xl = 1:Ls; 
 for yl = 1:Ls; 
 if xl==yl 
 dij(xl,yl)=Req(xl); 
 y1=Geom(yl,3); 
 hij(xl,yl)=2*y1; 
 Dij(xl,yl)=hij(xl,yl); 
 else
 x=abs(Geom(yl,2)-Geom(xl,2)); 
 y=abs(Geom(yl,3)-Geom(xl,3)); 
 dij(xl,yl)=sqrt(x^2 + y^2); 
 y1=Geom(xl,3); 
 y2=Geom(yl,3); 
 hij(xl,yl)=y1+y2; 
 x=abs(Geom(yl,2)-Geom(xl,2)); 
 y=hij(xl,yl); 
 Dij(xl,yl)=sqrt(x^2 + y^2); 
 end
 end
end 
Function InitialPoles 
function [Ps]=InitialPoles(f,Npol) 
 even = fix(Npol/2); % # of complex initial poles
p_odd = Npol/2 - even; % Auxiliary variable to check if number 
% of initial poles is odd
disc = p_odd ~= 0; % 0 for even Nr of initial poles & 1 – for 
% odd Nr.
% Set a real pole in case of disc == 1
if disc == 0 % Even Nr of initial poles
 pols = []; 
else % Odd Nr of initial poles
 pols = [(max(f)-min(f))/2]; 
end
% Set the complex initial poles
bet = linspace(min(f),max(f),even); 
for n=1:length(bet) 
 alf=-bet(n)*1e-2; 
 pols=[pols (alf-j*bet(n)) (alf+j*bet(n)) ]; 
end
Ps = pols.'; % Column vector of initial poles
```