Matlab coding help
Matlab coding help
(OP)
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ž
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ž
RE: Matlab coding help
RE: Matlab coding help
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
RE: Matlab coding help
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)
RE: Matlab coding help
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.
RE: Matlab coding help
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)
RE: Matlab coding help
RE: Matlab coding help
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 = @(y) find(y(:).*circshift(y(:), [-1 0]) <= 0); % Returns Zero-Crossing Indices Of Argument Vector
zx = zci(y)
subplot(2,1,1),plot(t,y,'gx-');
subplot(2,1,2),plot(t,integral,'b');
RE: Matlab coding help
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.
RE: Matlab coding help
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.
RE: Matlab coding help