## 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