Matlab coding help 1

[powerpi]

Hello

I am running a Matlab script and I would need a tip (or help) in coding on how to get things done. I am trying to simulate the CT saturation in a protection core CT. To do this I need to integrate the area of the secondary current (lower side of the plot) and compare it to a predefined value. Once the integration value is higher than this, the secondary current should drop to zero. It should stay at zero value until the zero crossing appears in the oposite direction.

I would appreciate any help. It is probably an easy task for a skilled programmer, but a paramount for me.

Thanks

Blaž

 

[Stick.]

My first thought is along these lines: Generate a vector of zeros of the appropriate size for the integral. Then inside of a for loop, solve/numerically approximate the integral (depending on how your data is defined) at the current index position in the data vector and compare it against your predefined value. If it's below the value, change the value of the position in the integral vector for the loop index you are on to the calculated value. If it's above, leave it as zero. You'll have to figure out how to handle if it's equal since I know nothing about your application. You could probably also do this with a while loop (it might be faster, so I would look in to it if you care). I'm a pretty bad programmer so there's probably a better way to do this, but I think that should work.

 

[powerpi]

I am not sure if I understand you correct. I have added a piece of code in my reply so you can se what I am doing. First part is of no importance as it is only the definition of constants, to be used later in the code. After the "for" loop it is when it all starts. I have calculated the primary current Ip, convert it with a CT ratio to secondary one "is" and the in begins. I need to integrate the area below "is" and in the moment it becomes higher than a PCT, the secondary current ought to fall to zero until the next zero crossing. Then all ought to repeat. The "is" should be partially the scaled value of ip, after the saturatuon it should be zero. Iss is an absolute value of scaled ip current.

I am also not a coder so this ought to be done much more elegant...
-----------------------------------------------------

%Model tokovnega transformatorja

%Naloga:

%Spreminjajte: amplitudo in casovno konstanto enosmerne komponente primarnega toka
% in koleno magnetilne krivulje
function [ALFd, KtFmax]=CTtest(Ik,CT)
%Opazujte: poteke primernega in sekundarnega toka ter fluksa tokovnega
% merilnega transformatorja med normalnim obratovanjem in v
% nasicenju

if nargin <2
CT=600;
end

if nargin <1
Ik=15000;
end

w=100*pi; %omega
Rb=2;
Rct=7;
Pn=20;
ALF=15; %primarni, skundarni ovoji in bremenski R
T2=0.5;
T1=0.1; % časovna konstanta mreže
Rmagnas=700; %zacetna vrednost fluksa
dt=0.001;
tk=0.07;
c=0;
m=ceil(tk/dt); %korak, koncni cas in stevilo korakov simulacije
ALFd=ALF*((Rct+Pn/(1*1))/(Rct+Rb))
ICT=ALF*(Rb+Rct);
Ktfmax=ALFd/(Ik/CT);

ip=zeros(m,1); %polje za primarni tok
is=zeros(m,1);
iss=zeros(m,1);
temp=zeros(m,1); %polje za sekundarni tok
flux=zeros(m,1); %polje za fluks
t=zeros(m,1);
u=zeros(m,1);
Ktf=zeros(m,1);
CTsat=ones(m,1);
PCT=ICT;
index1=0;
index2=0;

for k=1:1:m %simulacijska zanka
t(k)=k*dt;
ip(k) = Ik*(exp(-t(k)/T1)-cos(w*t(k))); %primarni tok
Ktf(k)=1+((w*T1*T2)/(T2-T1))*(exp(-t(k)/T2)-exp(-t(k)/T1)); %tranzientni faktor
is(k) =(ip(k)/CT); %nepopacen sek tok
iss(k)=abs(is(k)); %absolutna vrednost toka
temp(k)=trapz(iss);
...
end

 

[Stick.]

I was having trouble sorting out your code since there is a lot of extraneous stuff and also because I'm not at all familiar with your application. So I wrote my own code to calculate the definite integral for a sine function until the integral value reaches a cutoff value, then the integral drops to zero. I believe this should be analogous to your situation. Note here that because the sine function I used is differentiable, I evaluated the integral exactly. If your function isn't, then you will have to use your favorite numerical integration method. I'm almost certain you could do this faster with a while loop instead of a for loop, but I don't have the time to write that right now. And just to reiterate, this is far from elegant code, so hopefully someone comes along with a better method.

clear, clc, close all;

dt = 0.1; %time step size
t = 0:dt:10; %time vector
y = sin(t)+1; %arbitrary periodic function, shifted up by the plus 1 to have all positive values
cutoff = 8; %arbitrary cutoff value

integral = zeros(length(t)); %initializing the integral vector with all zeros

for i = 1:length(t) %start of for loop, running from index of i from 1 to the length of the time vector
TempIntegral = (t(i)-cos(t(i)))-(t(1)-cos(t(1))); %evaluating the definite integral of the function y from t(1) to t(i), which is done exactly here since the function for y is both known and differentiable. Other cases may require numerical integration such as the trapezoidal rule or other methods

if TempIntegral < cutoff
integral(i) = TempIntegral;
elseif
integral(i) = 0;
end

end

plot(t,integral)




%(If the program doesn't run, double check the syntax on the for loops and if statements since I originally wrote this in Octave and edited it to what I believe is the correct matlab syntax)

 

[powerpi]

Hello again

Thanks for the code. I have run it and I have found two issues.

First, I need to numericaly integrate the signal ("y" in your code, "is" in mine). So I would need to use one of the integral functions available. I thought of trapz, but I would have to tell it when to stop and when to reset and start again.

And that brings me to the second issue. Please check the time plot of "y" and "integral" signals. The integral drops to zero at the cutoff value correctly, but it would have to restart the calculation from zero at every zero crossing, after the cutoff value is reached.

Otherwise, thanks for the effort and help.

 

[Stick.]

Here is some updated code that fixes two errors I caught (using the zeros function incorrectly and using an elseif where I should have used an else). I also changed it so that when the integral drops to zero, it restarts integrating from that point, which is what I understood your second issue to be (if it's not, chalk it up to my lack of electrical knowledge). I also threw in a commented out line that uses the trapz function for the integral, but I only tested it once so definitely double check me. Using cumtrapz may be possible as well, not sure if that could make it faster if you care about that.

clear, clc, close all;

dt = 0.1; %time step size
t = 0:dt:10; %time vector
y = sin(t)+1; %arbitrary periodic function, shifted up by the plus 1 to have all positive values
cutoff = 4; %arbitrary cutoff value

integral = zeros(1,length(t)); %initializing the integral vector with all zeros
j = 1; %index for the start of the interval for the integral

for i = 1:length(t) %start of for loop, running from index of i from 1 to the length of the time vector
TempIntegral = (t(i)-cos(t(i)))-(t(j)-cos(t(j))); %evaluating the definite integral of the function y from t(j) to t(i), which is done exactly here since the function for y is both known and differentiable. Other cases may require numerical integration such as the trapezoidal rule or other methods
%TempIntegral = trapz(dt,y(j:i));

if TempIntegral < cutoff
integral(i) = TempIntegral;
else
integral(i) = 0;
j = i; %if the integral is reset to zero, then the index for the start of the interval to evaluate the definite integral is changed to the current loop index i
end

end

plot(t,integral)

 

[powerpi]

a, ok, now I see. I have tested it with the normal sine function and it works ok. now I only need to restart the integral with some delay, after the next zero crossing. Otherwise thanks. Much appreciated.

 

[powerpi]

I have added two lines for finding the zero crossing in the original signal "y". The "zx" returns an index of the points where zero crossing takes place. Perhaps this information can be used to further extend the delay of "integral" until this index is reached. Now, the "integral" function, as soon as it hits zero, starts to count back, but it ought to wait...

clear, clc, close all;

dt = 0.1; %time step size
t = 0:dt:10; %time vector
y = sin(t); %arbitrary periodic function, shifted up by the plus 1 to have all positive values
cutoff = 1; %arbitrary cutoff value

integral = zeros(1,length(t)); %initializing the integral vector with all zeros
j = 1; %index for the start of the interval for the integral

for i = 1:length(t) %start of for loop, running from index of i from 1 to the length of the time vector
%TempIntegral = (t(i)-cos(t(i)))-(t(j)-cos(t(j))); %evaluating the definite integral of the function y from t(j) to t(i), which is done exactly here since the function for y is both known and differentiable. Other cases may require numerical integration such as the trapezoidal rule or other methods
TempIntegral = trapz(dt,y(j:i));

if TempIntegral < cutoff
integral(i) = TempIntegral;
else
integral(i) = 0;
j = i; %if the integral is reset to zero, then the index for the start of the interval to evaluate the definite integral is changed to the current loop index i
end

end

zci = @ find(y).*circshift(y), [-1 0]) <= 0); % Returns Zero-Crossing Indices Of Argument Vector
zx = zci
subplot(2,1,1),plot(t,y,'gx-');
subplot(2,1,2),plot(t,integral,'b');

 

[Stick.]

Ah, I didn't realize that you needed a delay before restarting the integration. Having the locations of the restarts should help. Possibly you could replace the j=i statement with j=(somehow determine the correct index from zx). If you did that you would also have to deal with the case in the TempIntegral line where j is greater than i, so that may take an if statement and a bit more work. Take a crack at it and let me know how it goes.

If the integral doesn't exceed the cutoff by the time it gets to a zero crossing index, what happens then (or is that physically impossible?)? Will it still reset to zero, or will it keep growing until the cutoff? I'm wondering if it's possible to evaluate the integral between each index in zx (so multiple definite integrals in one vector) and then just change any value in that result that's above the cutoff to zero.

 

[powerpi]

Yes, now I have added additional delay to the "j" index and it works fine. Thank you very much for your help. I wouldn't have make it without you.

In general the negative an positive values of integral should compensate each other out, but when unevenly distributed zero crossings are present, one side (positive one) prevails. Yes, It would also be possible to make parital integrals I guess and then compare those with the cutoff value.

 

[Stick.]

Not a problem, glad my inelegant solution helped out.