Why for loop doesn't work in my code

sajad_ghassemi · July 30, 2021, 8:14pm

Hi dear friends,
I have problem with for loop. I don’t know my code doesn’t enter and repeat in loop.

using Plots
plotly()
Plots.PlotlyBackend()
currentcap=zeros(1,250001);
voltagecap=zeros(1,250001);
currentind=zeros(1,250001);
V= zeros(9,250004);
t_end = 250e-3;
dt = 1e-6;
t = 0:dt:t_end;
vp=20;
freq=50;
vs = vp*sin.(2*π*freq.* t')
L=1e-3;
R=1000;
C=1e-6;
Roff=10^10;
i=2;
c1=0;
c2=0;
for i=2:250001
         c1=c1+1;
         Y = [dt/L -dt/L 0 0 1 0 0 0 0
             -dt/L dt/L+1/Roff -c/Roff -a/Roff 0 1 0 0 -1
              0 -c/Roff (1/R)+(C/dt)+1/Roff -1/R-C/dt 0 -1 -1 0 0
              0 -a/Roff -1/R-C/dt 1/R+C/dt+1/Roff  0 0 0 0 1 1
              1 0 0 0 0 0 0 0 0
              0 a -a 0 0 b 0 0 0
              0 0 -c 0 0 0 e 0 0
              0 0 0 f 0 0 0 g 0
              0 -h 0 0 0 0 0 0 n];
              x1=V[5,i]-currentind[1,i-1]; x2=currentind[1,i-1]-V[6,i]+V[9,i]; x3=V[6,i]-currentcap[1,i-1]; x4=currentcap[1,i]-V[8,i]-V[9,i] ;
              I =[x1
                  x2
                  x3
                  x4
                  vp*sin.(2*π*freq.*t'[1,i])
                  0
                  0
                  0
                  0];
                x=inv(Y)*I;
                V[:,i]=x;
                currentcap[1,i]=(C/dt)*(V[3,i]-V[3,i-1]-V[4,i]+V[4,i-1]);
                 voltagecap[1,i]=V[3,i]-V[4,i];
                currentind[1,i]=(dt/L)*(V[1,i]-V[2,i]+V[1,i-1]-V[2,i-1])+currentind[1,i-1];
                if V[2,i]-V[3,i]>0 || V[6,i]>0
                   a=1;b=0;
                elseif V[6,i]<=0
                     a=0; b=1;
                else 
                  a=0; b=1;
                end
                if -V[3,i]>0 || V[7,i]>0
                   c=1; e=0;
                elseif V[7,i]<=0
                     c=0; e=1
               else
                    c=0; e=1;
                end
               if -V[4,i]>0 || V[8,i]>0
                  f=1; g=0;
               elseif V[8,i]<=0
                  f=0; g=1
               else
                  f=0; g=1;
               end
               if -V[4,i]-V[2,i]>0 || V[9,i]>0
                  h=1; n=0;
               elseif V[9,i]<=0
                  h=0; n=1
               else 
                  h=0; n=1;
               end  
            c2=c2+1;
end


voltagecap=voltagecap';
plot(t,voltagecap)

zdenek_hurak · July 30, 2021, 8:20pm

Before I encounter any for-loop-related problems, I encounter the problem that c is not defined. See the third line in your definition of Y.

Ronis_BR · July 30, 2021, 8:45pm

Yes, c is not defined (Julia is type sensitive, so c != C). Also, that for loop is not right. Where is the end ?

sajad_ghassemi · August 1, 2021, 3:52pm

Thanks for considering my problem.
Yeah you’re right. I defined but still loop doesn’t work

sajad_ghassemi · August 1, 2021, 3:54pm

Thanksfor the response,
at the end of the loop there is end. every if has one end and last end is related to loop

nilshg · August 1, 2021, 4:28pm

What is the actual error you are seeing? Make sure that any code you post is reproducible so that others can help you. In trying to run your code, I got c not defined, when defining c, I got the same error for a, then b, then e, and so on until n. After defining all variables, I got a mismatching dimensions error in the creation of the Y matrix:

Y = [dt/L -dt/L 0 0 1 0 0 0 0
             -dt/L dt/L+1/Roff -c/Roff -a/Roff 0 1 0 0 -1
              0 -c/Roff (1/R)+(C/dt)+1/Roff -1/R-C/dt 0 -1 -1 0 0
              0 -a/Roff -1/R-C/dt 1/R+C/dt+1/Roff  0 0 0 0 1 1
              1 0 0 0 0 0 0 0 0
              0 a -a 0 0 b 0 0 0
              0 0 -c 0 0 0 e 0 0
              0 0 0 f 0 0 0 g 0
              0 -h 0 0 0 0 0 0 n]

If you count carefully, you see that row 4 has 10 entries, while rows 1 to 3 only have 9 entries.

So in short there are loads of errors with your current code, but it’s not clear whether any of them are your actual problem or just come from the incomplete snippet you’ve posted.

sajad_ghassemi · August 1, 2021, 5:09pm

Sorry dear for my silly mistake
this is the new code:

using Plots
plotly()
Plots.PlotlyBackend()
currentcap=zeros(1,250001);
voltagecap=zeros(1,250001);
currentind=zeros(1,250001);
V= zeros(9,250004);
t_end = 250e-3;
dt = 1e-6;
t = 0:dt:t_end;
vp=20;
freq=50;
vs = vp*sin.(2*π*freq.* t')
L=1e-3;
R=1000;
C=1e-6;
Roff=10^10;
i=2;
c1=0;
c2=0;
a=0; b=1; f=0; g=1; c=0; e=1; h=0; n=1;
  for i=2:250001
         c1=c1+1;
         Y = [dt/L -dt/L 0 0 1 0 0 0 0
             -dt/L dt/L+1/Roff -c/Roff -a/Roff 0 1 0 0 -1
              0 -c/Roff (1/R)+(C/dt)+1/Roff -1/R-C/dt 0 -1 -1 0 0
              0 -a/Roff -1/R-C/dt 1/R+C/dt+1/Roff  0 0 0 1 1
              1 0 0 0 0 0 0 0 0
              0 a -a 0 0 b 0 0 0
              0 0 -c 0 0 0 e 0 0
              0 0 0 f 0 0 0 g 0
              0 -h 0 0 0 0 0 0 n];
              x1=V[5,i]-currentind[1,i-1]; x2=currentind[1,i-1]-V[6,i]+V[9,i]; x3=V[6,i]-currentcap[1,i-1]; x4=currentcap[1,i]-V[8,i]-V[9,i] ;
              I =[x1
                  x2
                  x3
                  x4
                  vp*sin.(2*π*freq.*t'[1,i])
                  0
                  0
                  0
                  0];
                x=inv(Y)*I;
                V[:,i]=x;
                currentcap[1,i]=(C/dt)*(V[3,i]-V[3,i-1]-V[4,i]+V[4,i-1]);
                 voltagecap[1,i]=V[3,i]-V[4,i];
                currentind[1,i]=(dt/L)*(V[1,i]-V[2,i]+V[1,i-1]-V[2,i-1])+currentind[1,i-1];
                if V[2,i]-V[3,i]>0 || V[6,i]>0
                   a=1;b=0;
                elseif V[6,i]<=0
                     a=0; b=1;
                else
                  a=0; b=1;
                end
                if -V[3,i]>0 || V[7,i]>0
                   c=1; e=0;
                elseif V[7,i]<=0
                     c=0; e=1
               else
                    c=0; e=1;
                end
               if V[4,i]>0 || V[8,i]>0
                  f=1; g=0;
               elseif V[8,i]<=0
                  f=0; g=1
               else
                  f=0; g=1;
               end
               if V[4,i]-V[2,i]>0 || V[9,i]>0
                  h=1; n=0;
               elseif V[9,i]<=0
                  h=0; n=1
               else
                  h=0; n=1;
               end
            c2=c2+1;
   end


voltagecap=voltagecap';
plot(t,voltagecap)

through checking c1 it shows the loop executed only 26430, while the it should be equal to 250001.

kristoffer.carlsson · August 1, 2021, 5:33pm

Have you tried to successively remove code until the problem disappears (or until the issue is obvious)?

sajad_ghassemi · August 1, 2021, 5:39pm

actually, I wrote the code in MATLAB and it works ( this is a simple code).
But Ido not know why whenever I wrote my code in Julia. I came accross many silly errors. I agree about my silly mistake like mismatch number of columns or not defining a,b,…n. but now everything is similar to code in MATLAB.

kristoffer.carlsson · August 1, 2021, 5:40pm

I meant regarding

If you take away a piece of your code, does this still happen? If so, you can keep taking away things until the answer might be obvious.

DNF · August 1, 2021, 5:47pm

Maybe you could post your Matlab code as well?

sajad_ghassemi · August 1, 2021, 5:51pm

yeah. nice suggestion. I will try this.
You know dear friend. I am simulating a Full-Bridge Rectifier through Modified Augmented Nodal Analysis (MANA) through programming in Julia. since I was familiar with MATLAB, I wrote in MATLAB first and then I wanted to wrote in Julia. My supervisor told me I should solve the equation through iteration. but still I have problem with understanding the concpet of MANA as well his algorithm in his paper.
he told me I should use another loop in first loop to solve the equation but I haven’t found yet
any way thanks for the suggestion

sajad_ghassemi · August 1, 2021, 5:56pm

nice way but I didn’t want to waste people’s time.

ToucheSir · August 1, 2021, 6:06pm

Not at all! If anything, it might be easier/faster to help translate the matlab code (assuming it’s self-contained and runs).

DNF · August 1, 2021, 6:11pm

I think it will be faster and easier to fix the Julia code if we can see the Matlab code.

sajad_ghassemi · August 1, 2021, 6:16pm

okay. it is very easy for me to show my code in MATLAB

clc
clear all

currentcap=zeros(1,250001);
voltagecap=zeros(1,250001);
currentind=zeros(1,250001);
V= zeros(9,250004);
t_end = 250e-3;                            %simulation time [sec.]
dt = 1e-6;                                 %simulation time-step [sec.]
t = 0:dt:t_end;
vp = 20;                                   %peak excitation voltage [V]
w = 2*pi*50;                               %angular frequency [rad/sec]
vs = vp*sin(w*t);                          % voltage source
L=1e-3;                                    % inductance [Henri]
R=10000;                                      % resistance [ohm] 
C=1e-6;                                    % Capacitor [Farad]
Roff=10^(10);
a=0; b=1; f=0; g=1; c=0; e=1; h=0; n=1;    % all the diodes are off 

  for i=2:250001 
      
          %build Y matrix
          Y = [dt/(L), -dt/(L), 0, 0, 1, 0, 0, 0, 0;
              -dt/(L), dt/(L)+1/Roff, -c/Roff, -a/Roff, 0, 1, 0, 0, -1;
               0, -c/Roff, (1/R)+((C)/dt)+1/Roff, -(1/R)-((C)/dt), 0, -1, -1, 0, 0;
               0, -a/Roff, -(1/R)-((C)/dt), (1/R)+((C)/dt)+1/Roff, 0, 0, 0, 1, 1;
               1, 0, 0, 0, 0, 0, 0, 0, 0;
               0, a,-a, 0, 0, b, 0, 0, 0;
               0, 0,-c, 0, 0, 0, e, 0, 0;
               0, 0, 0,-f, 0, 0, 0, g, 0;
               0,-h, 0, h, 0, 0, 0, 0, n];
           %using current kirchhoff law in every nodes and building current matrix    
           x1=V(5,i)-currentind(1,i-1); x2=currentind(1,i-1)-V(6,i)+V(9,i); x3=V(6,i)-currentcap(1,i)+V(7,i); x4= currentcap(1,i)-V(8,i)-V(9,i);     
           
           I =[x1;
               x2;
               x3;
               x4;
               vp*sin(w*t(1,i));
               0;
               0;
               0;
               0];
            
            % calculate voltage of every node as well as current passing through diodes
               x=inv(Y)*I;
               V(:,i)=x;
               
            % calculating the current and voltage of capacitor and inductor
              currentcap(1,i)=((C)/dt)*[V(3,i)-V(4,i)-V(3,i-1)+V(4,i-1)];
              voltagecap(1,i)=V(3,i)-V(4,i);
              currentind(1,i)=(dt/(L))*(V(1,i)-V(2,i)+V(1,i-1)-V(2,i-1))+currentind(1,i-1);
               
            % update diode-1 status  
               if V(2,i)-V(3,i)>0 || V(6,i)>0   % check the voltage across diode 1 or its current
                  a=1; b=0;                     % diode 1 become on
               elseif V(6,i)<=0                 % check the current passing diode 1
                  a=0; b=01;                    % diode 1 become off
               else
                  a=0; b=1;                     % diode 1 become off
               end
               
            % update diode-2 status 
               if 0-V(3,i)>0  || V(7,i)>0   % check the voltage across diode2 or its current
                  c=1; e=0;                 % diode2 become on
               elseif V(7,i)<=0             % check the current passing through diode2
                  c=0; e=1;                 % diode2 become off
                  else
                  c=0; e=1;                 % diode2 become off
               end
               
            % update diode-3 status 
               if V(4,i)>0  || V(8,i)>0     % check the voltage across diode3 or its current
                  f=1; g=0;                 % diode3 become on
               elseif V(8,i)<=0             % check the current passing diode3
                  f=0; g=1;                 % diode3 become off
               else
                  f=0; g=1;                 % diode3 become off
               end
               
             % update diode-4 status 
               if V(4,i)-V(2,i)>0 || V(9,i)>0   % check the voltage across diode4 or its current
                  h=1; n=0;                     % diode4 become on
               elseif V(9,i)<=0                 % check the current passing diode4
                  h=0; n=1;                     % diode4 become off
               else
                  h=0; n=1;                     % diode4 become off
               end               
  end
V=V';
t=t';
currentcap=currentcap';
voltagecap=voltagecap';
plot(t,voltagecap)

DNF · August 1, 2021, 9:03pm

You have not been sufficiently careful when translating your code. There are several discrepancies between your Matlab and Julia codes. At least I found these:

Matlab:

Y = [dt/(L), -dt/(L), 0, 0, 1, 0, 0, 0, 0;
              ...
               0, 0, 0,-f, 0, 0, 0, g, 0;
               0,-h, 0, h, 0, 0, 0, 0, n];

Julia, notice the errors in the last two rows:

Y = [dt/L -dt/L 0 0 1 0 0 0 0
              ...
              0 0 0 f 0 0 0 g 0
              0 -h 0 0 0 0 0 0 n];

Matlab code:

x3=V(6,i)-currentcap(1,i)+V(7,i);

Julia code:

x3=V[6,i]-currentcap[1,i-1];

I also noticed that you introduce c1 and c2 and increment them for each iteration, but you don’t use them for anything, and they are not present in your Matlab code.

There may be more, but I suggest that you do this part.

I’ll give one more piece of advice: Don’t write

x=inv(Y)*I;

This is bad numerically. Use

x = Y \ I;

instead. (This is the same in Matlab, use Y \ I)

sajad_ghassemi · August 4, 2021, 6:48pm

I really appreciate the time you look carefuly at my code dear friend
I changed it but loop still haven’t worked. I defined c1 and c2 to see wheter my code inters the loop or not? if yes how many times it iterates?

using Plots
plotly()
Plots.PlotlyBackend()
currentcap=zeros(1,250001);
voltagecap=zeros(1,250001);
currentind=zeros(1,250001);
V= zeros(9,250004);
t_end = 250e-3;
dt = 1e-6;
t = 0:dt:t_end;
vp=20;
freq=50;
vs = vp*sin.(2*π*freq.* t')
L=1e-3;
R=1000;
C=1e-6;
Roff=10^10;
i=2;
c1=0;
c2=0;
a=0; b=1; f=0; g=1; c=0; e=1; h=0; n=1;
  for i=2:250001
         c1=c1+1;
         Y = [dt/L -dt/L 0 0 1 0 0 0 0 1
             -dt/L dt/L+1/Roff -c/Roff -a/Roff 0 1 0 0 -1
              0 -c/Roff (1/R)+(C/dt)+1/Roff -1/R-C/dt 0 -1 -1 0 0
              0 -a/Roff -1/R-C/dt 1/R+C/dt+1/Roff  0 0 0 1 1
              1 0 0 0 0 0 0 0 0
              0 a -a 0 0 b 0 0 0
              0 0 -c 0 0 0 e 0 0
              0 0 0 f 0 0 0 g 0
              0 -h 0 h 0 0 0 0 n];
              x1=V[5,i]-currentind[1,i-1]; x2=currentind[1,i-1]-V[6,i]+V[9,i]; x3=V[6,i]-currentcap[1,i-1]+V[7,i]; x4=currentcap[1,i]-V[8,i]-V[9,i] ;
              I =[x1
                  x2
                  x3
                  x4
                  vp*sin.(2*π*freq.*t'[1,i])
                  0
                  0
                  0
                  0];
                x=I/Y;
                V[:,i]=x;
                currentcap[1,i]=(C/dt)*(V[3,i]-V[3,i-1]-V[4,i]+V[4,i-1]);
                 voltagecap[1,i]=V[3,i]-V[4,i];
                currentind[1,i]=(dt/L)*(V[1,i]-V[2,i]+V[1,i-1]-V[2,i-1])+currentind[1,i-1];
                if V[2,i]-V[3,i]>0 || V[6,i]>0
                   a=1;b=0;
                elseif V[6,i]<=0
                     a=0; b=1;
                else
                  a=0; b=1;
                end
                if -V[3,i]>0 || V[7,i]>0
                   c=1; e=0;
                elseif V[7,i]<=0
                     c=0; e=1
               else
                    c=0; e=1;
                end
               if V[4,i]>0 || V[8,i]>0
                  f=1; g=0;
               elseif V[8,i]<=0
                  f=0; g=1
               else
                  f=0; g=1;
               end
               if V[4,i]-V[2,i]>0 || V[9,i]>0
                  h=1; n=0;
               elseif V[9,i]<=0
                  h=0; n=1
               else
                  h=0; n=1;
               end
            c2=c2+1;
   end


voltagecap=voltagecap';
plot(t,voltagecap)

zdenek_hurak · August 4, 2021, 7:08pm

Just one thing upon visual inspection: in your previous code you solved the linear system by

x=inv(Y)*I;

Then you were advised by @DNF not to do it by explicitly computing the inverse matrix and instead solve the linear system in the “matlab style”

x = Y \ I;

Now in your new code you have

x=I/Y;

That is quite a difference from what you were advised.

By the way, I recommend avoiding the symbol I for your variable. You may want to use it for the identity matrix

julia> using LinearAlgebra

julia> 1.0I(2)
2×2 Diagonal{Float64, Vector{Float64}}:
 1.0   ⋅ 
  ⋅   1.0

sajad_ghassemi · August 4, 2021, 7:16pm

Thanks , nice advice