Function Poles

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

 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));

legend('Sending end phase A' , 'Sending end phase B' , 'Sending end phase C' , 'Receiving end phase A' ,

for m = 1:Ng

for m = 1:NpYc

 IfarA(:,2:h3) = IfarA(:,1:h2); IfarB(:,2:h3) = IfarB(:,1:h2);

VO = Gys\*(Isr(:,k)+HistO);

vt = (0:Dt:length(Vi(1,:))\*Dt-(t0+4)\*Dt)';

figure(1),plot(vt,Vi(:,a1:a2),':',vt,Vf(:,a1:a2))

'Receiving end phase B' , 'Receiving end phase C')

Req(nc) = (Geom(nc,4).\*Geom(nc,5).\*k4(nc).^

ylabel('Amplitude in volts') xlabel('Time in seconds')

function[Dij,dij,hij]=Height(Geom) Ls = Geom(max(Geom(:,1)),1);

 VL = Gyr\*HistL; IO = Gy\*VO - HistO; IL = Gy\*VL - HistL; Vi(:,k) = VO; Vf(:,k) = VL; Ii(:,k) = IO; If(:,k) = IL;

 ZA(:,m) = Ai(m)\*ZA(:,m) + Di(:,:,m)\*VO; ZB(:,m) = Ai(m)\*ZB(:,m) + Di(:,:,m)\*VL;

 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);

 HistO = - sum(ZA(:,:),2) + sum(sum(YA(:,:,:),3),2); HistL = - sum(ZB(:,:),2) + sum(sum(YB(:,:,:),3),2);

end

 end for l = 1:Ng for m = 1:NpH

 end end

end Iri = Yi\*Vi; Irf = Yr\*Vf;

N = length(vt); a1 = t1+1; a2 = Nt+t1-3;

Function Height

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

```
function [A]=Poles(Fs,s,Pi,Ns,Ka); 
Np = length(Pi); % Length of vector containing starting poles
CPX = imag(Pi)~=0; % 0 for real pole and 1 for complex pole
```

```
rp = 0; % Initialize the index for real poles
cp = 0; % Initialize the index for complex poles
RePole = []; % Initialize the vector of real poles
CxPole=[];%Initialize the vector of complex poles
% Loop to separate real poles and complex poles
for k = 1:Np 
 if CPX(k) == 0 % Real Pole
 rp = rp + 1; 
 RePole(rp) = Pi(k); 
 elseif CPX(k) == 1 % Complex pole
 cp = cp + 1; 
 CxPole(cp) = Pi(k); 
 end
end
Lambda = Pi.'; 
RePole = sort(RePole); % Sort real poles
CxPole = sort(CxPole); % Sort complex poles
Lambda = [RePole CxPole]; % Concatenate poles
I = diag(ones(1,Np)); % Unit matrix
A = []; % Poles
B = ones(Ns,1); % the weight factor
C = []; % Residues
D = zeros(1); % Constant term
E = zeros(1); % Proportional term
KQA = ones(Ns,1); 
cpx = imag(Lambda)~=0; % 0 if pole is real and 1 if pole is 
% complex.
dix = zeros(1,Np); % Initializes vector of pole types
if cpx(1)~=0 % If the first pole is complex
 dix(1)=1; % real part
 dix(2)=2; % imag part
 k=3; % continue dix for third position
else
 k=2; % If the first pole is real continue dix for the second position
end
% complete the classification of the poles
for m=k:Np 
 if cpx(m)~=0 % If the pole is complex
 if dix(m-1)==1 
 dix(m)=2; % If the previous position has the real part put 2 
% to identifies the imag part
 else
 dix(m)=1; % 1 for the real part of a complex pole
 end
 end
end
% Creates matriz A divided in four parts A = [A1 A2 A3 A4]
% A1 = Dk
% A2 = B.*ones(Ns,1)
```

```
% A3 = B.*s
% A4 = -Dk*Fs
Dk=zeros(Ns,Np); % Initialize matrix with zeros
for m=1:Np % Iterative cycle for all poles
 if dix(m)== 0 % For a real pole
 Dk(:,m) = B./(s-Lambda(m)); 
 elseif dix(m)== 1 % For the real part
 Dk(:,m)=B./(s-Lambda(m)) + 
 B./(s-Lambda(m)'); 
 elseif dix(m)== 2 % For the imag part
 Dk(:,m) = i.*B./(s-Lambda(m-1)) - 
 i.*B./(s-Lambda(m-1)'); 
 end
end
% Creates work space for matrix A
A1 = Dk; 
A2 = B.*ones(Ns,1); 
A3 = B.*s; 
for col = 1:Np 
 A4(:,col) = -(Dk(:,col).*Fs.'); 
end 
% Asigns values to A
if Ka == 1 
 A = [A1 A4]; % Strictly proper rational fitting
elseif Ka == 2 
 A = [A1 A2 A4]; % Proper rational fitting
elseif Ka == 3 
 A = [A1 A2 A3 A4]; % Improper rational fitting
else
 disp('Ka need to be 1, 2 or 3') 
end
% Creates matrix b = B*Fs
b = B.*Fs.'; 
% Separating real and imaginary part
Are = real(A); % Real part of matrix A
Aim = imag(A); % Imaginary part of matrix A
bre = real(b); % Real part of matrix b
bim = imag(b); % Imaginary part of matrix b
An = [Are; Aim]; % Real and imaginary part of A
bn = [bre; bim]; % Real and imaginary part of b
% Routine to applies the Euclidian norm to An
[Xmax Ymax] = size(An); 
for col=1:Ymax 
 Euclidian(col)=norm(An(:,col),2); 
 An(:,col)=An(:,col)./Euclidian(col); 
end
% Solving system
Xn = An\bn; 
Xn = Xn./Euclidian.'; 
% Put the residues into matrix C
```
296 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

rp = 0; % Initialize the index for real poles cp = 0; % Initialize the index for complex poles RePole = []; % Initialize the vector of real poles CxPole=[];%Initialize the vector of complex poles % Loop to separate real poles and complex poles

for k = 1:Np

 rp = rp + 1; RePole(rp) = Pi(k);

 cp = cp + 1; CxPole(cp) = Pi(k);

Lambda = Pi.';

 end end

if CPX(k) == 0 % Real Pole

elseif CPX(k) == 1 % Complex pole

RePole = sort(RePole); % Sort real poles CxPole = sort(CxPole); % Sort complex poles Lambda = [RePole CxPole]; % Concatenate poles

I = diag(ones(1,Np)); % Unit matrix

B = ones(Ns,1); % the weight factor C = []; % Residues D = zeros(1); % Constant term E = zeros(1); % Proportional term

cpx = imag(Lambda)~=0; % 0 if pole is real and 1 if pole is

k=2; % If the first pole is real continue dix for the second position

dix(m)=2; % If the previous position has the real part put 2

dix(m)=1; % 1 for the real part of a complex pole

% Creates matriz A divided in four parts A = [A1 A2 A3 A4]

dix = zeros(1,Np); % Initializes vector of pole types

if cpx(1)~=0 % If the first pole is complex

k=3; % continue dix for third position

% complete the classification of the poles

if cpx(m)~=0 % If the pole is complex

A = []; % Poles

KQA = ones(Ns,1);

 dix(1)=1; % real part dix(2)=2; % imag part

% complex.

else

end

for m=k:Np

else

 end end end

% A1 = Dk

% A2 = B.\*ones(Ns,1)

if dix(m-1)==1

% to identifies the imag part

```
if Ka == 1 
 C = Xn(Np+1:Ymax); % Strictly proper fitting
elseif Ka == 2 
 C = Xn(Np+2:Ymax); % Proper rational fitting
elseif Ka == 3 
 C = Xn(Np+3:Ymax);% Improper rational fitting
else
 disp('Ka need to be 1, 2 or 3') 
end
% C complex when the residues are complex
for m=1:Np 
 if dix(m)==1 
 alfa = C(m); % real part of a complex pole
 betta = C(m+1); % imag part of a complex pole
 C(m) = alfa + i*betta; % the complex pole
 C(m+1) = alfa - i*betta; % the conjugate
 end
end
% Now calculate the zeros for sigma
BDA = zeros(Np); 
KQA = ones(Np,1); 
% Loop to calculate the zeros of sigma which are the new poles
for km = 1:Np 
 if dix(km)== 0 % For a real pole
 BDA(km,km) = Lambda(km); 
 elseif dix(km)== 1 % For a cp with - imag part
 BDA(km,km) = real(Lambda(km)); 
 BDA(km,km+1) = imag(Lambda(km)); 
 KQA(km) = 2; 
 Aux = C(km); 
 C(km) = real(Aux); 
 elseif dix(km)== 2 % For a cp with + imag part
 BDA(km,km) = real(Lambda(km)); 
 BDA(km,km-1) = imag(Lambda(km)); 
 KQA(km) = 0; 
 C(km) = imag(Aux); 
 end
end
ZEROS = BDA - KQA*C.'; 
POLS = eig(ZEROS).'; 
%Forcing (flipping) unstable poles to make them stable
uns = real(POLS)>0; 
POLS(uns) = POLS(uns)-2*real(POLS(uns)); 
% Sort poles in ascending order. First real poles and then complex poles
CPX = imag(POLS)~=0; % Set to 0 for a real pole and to1 for a 
%complex pole
rp = 0; % Initialize index for real poles
cp = 0; % Initialize index for complex poles
RePole = []; % Initialize the vector of real poles
CxPole = []; % Initialize the vector of cp
```

```
% Loop to separate real and complex poles
```
for k = 1:Np if CPX(k) == 0 % Real Pole rp = rp + 1; RePole(rp) = POLS(k); elseif CPX(k) == 1 % Complex pole cp = cp + 1; CxPole(cp) = POLS(k); end end RePole = sort(RePole); % Sort real poles CxPole = sort(CxPole); % Sort complex poles % For conjugate pairs store first the one with positive imag part CxPole = (CxPole.')'; NewPol = [RePole CxPole]; A = NewPol.'; % Output

#### Function Residue

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

if Ka == 1

elseif Ka == 2

elseif Ka == 3

for m=1:Np if dix(m)==1

disp('Ka need to be 1, 2 or 3')

else

end

 end end

BDA = zeros(Np); KQA = ones(Np,1);

 KQA(km) = 2; Aux = C(km); C(km) = real(Aux);

 KQA(km) = 0; C(km) = imag(Aux);

ZEROS = BDA - KQA\*C.'; POLS = eig(ZEROS).';

uns = real(POLS)>0;

%complex pole

 end end

for km = 1:Np

C = Xn(Np+1:Ymax); % Strictly proper fitting

C = Xn(Np+2:Ymax); % Proper rational fitting

C = Xn(Np+3:Ymax);% Improper rational fitting

% C complex when the residues are complex

 alfa = C(m); % real part of a complex pole betta = C(m+1); % imag part of a complex pole C(m) = alfa + i\*betta; % the complex pole C(m+1) = alfa - i\*betta; % the conjugate

% Loop to calculate the zeros of sigma which are the new poles

% Now calculate the zeros for sigma

 if dix(km)== 0 % For a real pole BDA(km,km) = Lambda(km);

 elseif dix(km)== 1 % For a cp with - imag part BDA(km,km) = real(Lambda(km)); BDA(km,km+1) = imag(Lambda(km));

 elseif dix(km)== 2 % For a cp with + imag part BDA(km,km) = real(Lambda(km)); BDA(km,km-1) = imag(Lambda(km));

%Forcing (flipping) unstable poles to make them stable

% Sort poles in ascending order. First real poles and then complex poles CPX = imag(POLS)~=0; % Set to 0 for a real pole and to1 for a

POLS(uns) = POLS(uns)-2\*real(POLS(uns));

rp = 0; % Initialize index for real poles cp = 0; % Initialize index for complex poles RePole = []; % Initialize the vector of real poles CxPole = []; % Initialize the vector of cp

function [C,D,E]=Residue(Fs,s,Pi,Ns,Ka); Np = length(Pi); CPX = imag(Pi)~=0; % 0 for a rp and 1 for cp rp = 0; % Initialize the index for real poles cp = 0; % Initialize the index for complex poles RePole = []; % Initialize the vector of real poles CxPole=[]; %Initialize the vector of complex poles % Loop to separate real poles and complex poles for k = 1:Np if CPX(k) == 0 % Real Pole rp = rp + 1; RePole(rp) = Pi(k); elseif CPX(k) == 1 % Complex pole cp = cp + 1; CxPole(cp) = Pi(k); End end RePole = sort(RePole); % Sort real poles CxPole = sort(CxPole); % Sort complex poles CxPole = (CxPole.')'; Lambda = [RePole CxPole]; I = diag(ones(1,Np)); % Unit diagonal matrix A = []; % Poles B = ones(Ns,1); % weight factor C = []; % Residues D = zeros(1); % Constant term E = zeros(1); % Proportional term cpx = imag(Lambda)~=0; % 0 for rp and 1 for cp dix = zeros(1,Np); % Vto identifies poles if cpx(1)~=0 % If the first pole is complex dix(1)=1; % put 1 in dix(1) for the real part dix(2)=2; % put 2 in dix(2) for the imag part

 k=3; % continue dix for the third position else k=2; % If the first pole is real continue dix for the second % position end % complete classification of the poles for m=k:Np if cpx(m)~=0 % If the pole is complex if dix(m-1)==1 dix(m)=2; % If the previous position has the real part, set to % 2 to identify the imag part else dix(m)=1; % put 1 for the real part of a cp end end end % Output matrices: Dk=zeros(Ns,Np); for m=1:Np if dix(m)==0 % Real pole Dk(:,m) = B./(s-Lambda(m)); elseif dix(m)==1 % Complex pole, 1st part Dk(:,m) = B./(s-Lambda(m)) + B./(s-Lambda(m)'); elseif dix(m)==2 % Complex pole, 2nd part Dk(:,m) = i.\*B./(s-Lambda(m-1)) - i.\*B./(s-Lambda(m-1)'); end end % Creates work space for matrices A and b AA1=Dk; AA2=B.\*ones(Ns,1); AA3=B.\*s; if Ka == 1 AA = [AA1]; % Strictly proper rational fit elseif Ka == 2 AA = [AA1 AA2]; % Proper rational fit elseif Ka == 3 AA = [AA1 AA2 AA3]; % Improper fit else disp('Ka must be 1, 2 or 3') end bb = B.\*Fs.'; AAre = real(AA); % Real part of matrix A AAim = imag(AA); % Imaginary part of matrix A bbre = real(bb); % Real part of matrix b bbim = imag(bb); % Imaginary part of matrix b AAn = [AAre; AAim]; % Real and imag part of A bbn = [bbre; bbim]; % Real and imag part of b [Xmax Ymax] = size(AAn); for col=1:Ymax Eeuclidian(col)=norm(AAn(:,col),2); AAn(:,col)=AAn(:,col)./Eeuclidian(col);

```
end
% Solving system X
Xxn=AAn\bbn; 
X=Xxn./Eeuclidian.'; 
% Putting residues into matrix C
C=X(1:Np); 
% C is complex when the residues are complex
for m=1:Np 
 if dix(m)==1 
 alfa = C(m); % real part of a complex pole
 betta = C(m+1); % imag part of a complex pole
 C(m) = alfa + i*betta; % the complex pole
 C(m+1) = alfa - i*betta; % the conjugate
 end
end
% Outputs
if Ka == 1 
 A = Lambda.'; % Poles
 C = C; % Residues
 D = 0; % Constant term
 E = 0; % Proportional term
elseif Ka == 2 
 A = Lambda.'; % Poles
 C = C; % Residues
 D = X(Np+1); % Constant term
 E = 0; % Proportional term
elseif Ka == 3 
 A = Lambda.'; % Poles
 C = C; % Residues
 D = X(Np+1); % Constant term
 E = X(Np+2); % Proportional term
End
```
300 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

dix(m)=2; % If the previous position has the real part, set to % 2 to identify the imag part

k=3; % continue dix for the third position

% complete classification of the poles

if cpx(m)~=0 % If the pole is complex

dix(m)=1; % put 1 for the real part of a cp

elseif dix(m)==1 % Complex pole, 1st part

% Creates work space for matrices A and b

AA = [AA1]; % Strictly proper rational fit

AA = [AA1 AA2]; % Proper rational fit

AA = [AA1 AA2 AA3]; % Improper fit

AAre = real(AA); % Real part of matrix A AAim = imag(AA); % Imaginary part of matrix A bbre = real(bb); % Real part of matrix b bbim = imag(bb); % Imaginary part of matrix b AAn = [AAre; AAim]; % Real and imag part of A bbn = [bbre; bbim]; % Real and imag part of b

disp('Ka must be 1, 2 or 3')

[Xmax Ymax] = size(AAn);

 Eeuclidian(col)=norm(AAn(:,col),2); AAn(:,col)=AAn(:,col)./Eeuclidian(col);

 Dk(:,m) = B./(s-Lambda(m)) + B./(s-Lambda(m)'); elseif dix(m)==2 % Complex pole, 2nd part

Dk(:,m) = i.\*B./(s-Lambda(m-1)) - i.\*B./(s-Lambda(m-1)');

k=2; % If the first pole is real continue dix for the second

else

% position end

for m=k:Np

else

 end end end

 end end

AA1=Dk;

AA3=B.\*s; if Ka == 1

elseif Ka == 2

elseif Ka == 3

bb = B.\*Fs.';

for col=1:Ymax

else

end

AA2=B.\*ones(Ns,1);

if dix(m-1)==1

% Output matrices: Dk=zeros(Ns,Np); for m=1:Np

 if dix(m)==0 % Real pole Dk(:,m) = B./(s-Lambda(m));

**Figure 14.** Transversal geometry of aerial line in example.

**Figure 15.** Voltage responses at sending and receiving ends. Unit step excitation.

**Figure 16.** Voltage responses at sending and receiving ends. Sinusoidal excitation.

#### **9. References**



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

**Figure 15.** Voltage responses at sending and receiving ends. Unit step excitation.

100



Amplitude in volts

0

200

400

600

200

300

Amplitude in volts

400

500

600 700

<sup>0</sup> 0.005 0.01 0.015 0.02 <sup>0</sup>

Sending end phase A Sending end phase B Sending end phase C Receiving end phase A Receiving end phase B Receiving end phase C

Sending end phase A Sending end phase B Sending end phase C Receiving end phase A Receiving end phase B Receiving end phase C

Time in seconds

**Figure 16.** Voltage responses at sending and receiving ends. Sinusoidal excitation.

Bode, H. W. (1945). *Network Analysis and Feedback Amplifier Design*, D Van Nostrand

<sup>0</sup> 0.005 0.01 0.015 0.02 -600

Time in seconds

Brandao Faria, J. A. & Borges da Silva, J. F. (1986). Wave Propagation in Polyphase Transmission Lines a General Solution to Include Cases Where Ordinary Modal Theory Fails, *Power Delivery, IEEE Transactions on*, vol.1, no.2, (April 1986), pp.(182-

Dommel H. W. (1969). Digital computer solution of electromagnetic transients in single and multiphase networks. *IEEE Trans. On Power App. Syst.,* Vol. 88, (July 1969),

**9. References** 

189).

pp.(388).

Company, London, 1945.

	- Semlyen, A., Abdel-Rahman, M. (1982). A state variable approach for the calculation of switching transients on a power transmission line, *Circuits and Systems, IEEE Transactions on*, vol.29, no.9, (September 1982), pp. (624-633)
	- Wedepohl, L. M., (1965). Electrical characteristics of polyphase transmission systems with special reference to boundary-value calculations at power-line carrier frequencies, *Electrical Engineers, Proceedings of the Institution of*, vol.112, no.11,(November 1965), pp.(2103-2112)
