×
INTELLIGENT WORK FORUMS
FOR ENGINEERING PROFESSIONALS

Contact US

Log In

Come Join Us!

Are you an
Engineering professional?
Join Eng-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!

*Eng-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Students Click Here

Matlab coding help

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ž

RE: Matlab coding help

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.

RE: Matlab coding help

(OP)
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

RE: Matlab coding help

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)

RE: Matlab coding help

(OP)
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.

RE: Matlab coding help

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)

RE: Matlab coding help

(OP)
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.

RE: Matlab coding help

(OP)
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 = @(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

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.

RE: Matlab coding help

(OP)
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.

RE: Matlab coding help

Not a problem, glad my inelegant solution helped out.

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Eng-Tips Forums free from inappropriate posts.
The Eng-Tips staff will check this out and take appropriate action.

Reply To This Thread

Posting in the Eng-Tips forums is a member-only feature.

Click Here to join Eng-Tips and talk with other members! Already a Member? Login


Resources

Low-Volume Rapid Injection Molding With 3D Printed Molds
Learn methods and guidelines for using stereolithography (SLA) 3D printed molds in the injection molding process to lower costs and lead time. Discover how this hybrid manufacturing process enables on-demand mold fabrication to quickly produce small batches of thermoplastic parts. Download Now
Design for Additive Manufacturing (DfAM)
Examine how the principles of DfAM upend many of the long-standing rules around manufacturability - allowing engineers and designers to place a part’s function at the center of their design considerations. Download Now
Taking Control of Engineering Documents
This ebook covers tips for creating and managing workflows, security best practices and protection of intellectual property, Cloud vs. on-premise software solutions, CAD file management, compliance, and more. Download Now

Close Box

Join Eng-Tips® Today!

Join your peers on the Internet's largest technical engineering professional community.
It's easy to join and it's free.

Here's Why Members Love Eng-Tips Forums:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close