This repository has been archived by the owner on Aug 24, 2023. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 1
/
crossmatch_underdampedInput.m
102 lines (102 loc) · 4.08 KB
/
crossmatch_underdampedInput.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
%0.1.1 Generalized method of computing transfer function of an equivalent
%0.1.2 transformer RLC ladder network.
%1.0 Defining experimentally controlled variables
%1.1 Number of RLC ladders (physically speaking, coils in the TF)
npotential=7;
nactual=7;
%1.2.1 Partial discharge in the form of a current source Is
%1.2.2 (to be defined later) Is=?
%1.3.1 Location of node of partial discharge development (i):
%1.3.2 Right now only ALONG the winding
ipotential=7;
iactual=3;
%***********************************************************************
%2.0 Defining ladder nw constants:
Rs=1.33;Cs=0.6;Cg=0.933;Ls=0.4310;
M=[0.2392,0.1435,0.0947,0.0496,zeros(1,npotential-5)];
%***********************************************************************
%3.0 Constructing matrices useful for state space representation:
%3.1 Inductance (L) Matrix
L=zeros(npotential,npotential);
for r= 1:npotential
for c= 1:npotential
if r == c
L(r,c)=Ls;
else
L(r,c)=M(abs(r-c));
end
end
end
%3.2 Resistance (R) Matrix
R=Rs*eye(npotential);
%3.3.1 (S,SL1,SL2,SN1,SN2) Matrices
%3.3.2 L derived from when Live (node) shorted (to GND)
%3.3.3 N derived from when Neutral (node) shorted (to GND)
S=[-eye(npotential),zeros(npotential,1)]+[zeros(npotential,1),eye(npotential)];
SL1=S;SL1(:,1)=[];
SL2=zeros(npotential,1);SL2(ipotential-1,:)=1;
SN1=S;SN1(:,end)=[];
SN2=zeros(npotential,1);SN2(ipotential,:)=1;
%3.4 Capacitance (C) Matrix with spinoffs (CL1,CN1)
tmp=eye(npotential+1);tmp(1,1)=1/2;tmp(end,end)=1/2;
C=(Cg+2*Cs)*tmp+(-Cs)*[zeros(npotential,1),eye(npotential);zeros(1,npotential+1)]+(-Cs)*[zeros(1,npotential+1);eye(npotential),zeros(npotential,1)];
CL1=C;CL1(1,:)=[];CL1(:,1)=[];
CN1=C;CN1(end,:)=[];CN1(:,end)=[];
%***********************************************************************
%4.0 Define (A),(B),(C),(D) matrices (State Space Representation):
%But first, define (M),(-G),(T1):
ML=[L,zeros(npotential,npotential);zeros(npotential,npotential),CL1];
GL=-[-R,-SL1;SL1.',zeros(npotential,npotential)];
T1L=[zeros(npotential,1);SL2];
MN=[L,zeros(npotential,npotential);zeros(npotential,npotential),CN1];
GN=-[-R,-SN1;SN1.',zeros(npotential,npotential)];
T1N=[zeros(npotential,1);SN2];
AL=-ML\GL;
BL=ML\T1L;
AN=-MN\GN;
BN=MN\T1N;
tmpL=zeros(1,2*npotential);tmpL(:,1)=1;
tmpN=zeros(1,2*npotential);tmpN(:,npotential)=1;
CL=-Cs*AL(npotential+1,:)+tmpL;
DL=-Cs*BL(npotential+1,:);
CN=Cs*AN(2*npotential,:)+tmpN;
DN=Cs*BN(2*npotential,:);
%***********************************************************************
%5.0 Convert State Space Representation (SSR) into Transfer Function (TF)
[TFLb,TFLa]=ss2tf(AL,BL,CL,DL);
sysL=tf(TFLb,TFLa)
[TFNb,TFNa]=ss2tf(AN,BN,CN,DN);
sysN=tf(TFNb,TFNa)
%***********************************************************************
%6.0 Plot obtained TF responses, etc.
dt=1e-3;
t = 0:dt:1;
impulse= t==0;
testFilenameL=['fyL_' num2str(nactual) '_' num2str(iactual) '_underdampedInput.mat'];
fyLload=load(testFilenameL);
testFilenameN=['fyN_' num2str(nactual) '_' num2str(iactual) '__underdampedInput.mat'];
fyNload=load(testFilenameN);
fyL=[fyLload.fyL];
fyN=[fyNload.fyN];
xLimp=(fft(lsim(sysL,impulse,t)));
xNimp=(fft(lsim(sysN,impulse,t)));
len=length(xLimp); %to take the frequency axis of the harmonics.
q=-(len-1)/2:(len-1)/2; %divide the frequency compone
fxLimp=sqrt(xLimp.*conj(xLimp));
fxNimp=sqrt(xNimp.*conj(xNimp));
xL=fyL./fxLimp;
xN=fyN./fxNimp;
plot(q,xN,'-b')
hold on
plot(q,xL,':or')
hold off
title(['Fourier Spectrum Plots for n = ' num2str(npotential) ' and i = ' num2str(ipotential) '. Actual values being n = ' num2str(nactual) ' and i = ' num2str(iactual)]);
axis([0 2000 0 15]);
xlabel('Frequency');
ylabel('Amplitude');
ax = gca;
ax.FontSize = 13;
legend({'y = xN','y = xL'},'Location','northwest')
filename=['xN_vs_xL_actual_' num2str(nactual) '_' num2str(iactual) '_potential_' num2str(npotential) '_' num2str(ipotential) '_underdampedInput'];
saveas(gcf,filename,'png')
corrcoef(xN,xL)