Scientific computing loops and iterations

Hello

I am trying to simulate a process with the help of an m.file. The process is containing multiple structures created previously, although i think the looping sequence I am using has a fault in the way i structure it.

The process is as follows

Initially there is the main if loop which if satisfied it will enclose the consequential simulations:

if Demand <= Ewave; % this is the main limitation to be satisfied
 count = length(Ewave) % Based on the length of the input file to create a counter for the loops
    for count = (1:count-1)
        if E_comp>P_comp
            m_air_flow = (P_comp_joule *E_comp_joule*n_isc*n_mech)/(Cp_air*T_amb*((Press_r^inverse)-1));
            T_out_comp = T_amb*(1+(1/n_isc)*((Press_r^inverse)-1));
        elseif E_comp<=P_comp
            m_air_flow = (P_comp_joule*E_comp_joule *n_isc*n_mech)/(Cp_air*T_amb*(Press_r^inverse));
            T_out_comp = T_amb*(1+(1/n_isc)*(Press_r^inverse));
        end
    end
    E_to_store = E_comp*Default.n_comp; 
    E_to_store_joule = E_to_store * 1000000;
    count = length(E_to_store_joule)
        for count = (1:count-1)
% This is the main portion of interests, the initial values are to calculated once for the 1st iteration
            M_cavern(count) = E_to_store_joule(count) / (Cp_air*T_amb*(Press_r^inverse)-1);
            P_cavern(count) = Default.P_cavern; 
            T_cavern(count) = Default.T_cavern;       
% Then the 1st iteration results are used, and from that point on the counter has to ensure that the next iteration will use the previous quantity to calculate the new quantity              
            M_cavern(count+1) = M_cavern(count)+(m_air_flow-m_air_flow*Default.dm);
            P_cavern(count+1) = P_cavern(count) - (Default.dP*Press_r);
            T_cavern(count+1) = T_cavern(count) - (Default.dt*T_cavern(count));
        end
    Essmin = Load.N_ex *ho; 
    Vssmin_MWh =Essmin/DOD; 
    Vssmin_joule = Vssmin_MWh*3600000;
    mass_flow_guar = (Vssmin_joule *P_cavern*DOD)/(ho*T_cavern*R_g); 
            if count==1
                Tt3(count) = T_cavern+(T_cavern.*n_pre); 
                Tt4(count) = (Cp_air/Cp_gas)*(lamda*Ma/((lamda*Ma)+1))*...
                    Tt3a+(Hu/(((lamda*Ma)+1)*Cp_gas));
                if Tt4(count)=<Ttmax;
                    Tt5(count) = Tt4(count)*(1-n_isT*(1-(1/(expansion_ratio^inverse))));
                end
            end
end

As you can see, my focus is to create a loop sequence which i can expand later on in various loops, the initial quantities are calculated (1st) then until the end of the (Nth) iteration the process should be 1st,2nd=(2nd-1*,,) .... Nth=(Nth-1* ...)

Can anyone help me, in suggesting or if i has faulty set-up the looping sequence.

thank you in advance

Answer 1

[Geoff Hayes]

George - I did notice a couple of pieces of your code which I'll just highlight here. In your for loops you do the following

 count = length(Ewave) % or something similar
 for count = (1:count-1)

While the above works, it does seem a little confusing to name your index variable to that of one that you are indexing/iterating over. Consider instead

 for u=1:count-1

or something similar.

In your first for loop, you do the following

 for count = (1:count-1)
        if E_comp>P_comp
            m_air_flow = (P_comp_joule *E_comp_joule*n_isc*n_mech)/(Cp_air*T_amb*...
                    ((Press_r^inverse)-1));
            T_out_comp = T_amb*(1+(1/n_isc)*((Press_r^inverse)-1));
        elseif E_comp<=P_comp
            m_air_flow = (P_comp_joule*E_comp_joule *n_isc*n_mech)/(Cp_air*T_amb*...
                     (Press_r^inverse));
            T_out_comp = T_amb*(1+(1/n_isc)*(Press_r^inverse));
        end
  end

Note how none of the iterations of the loop depends upon count, and so you are calculating the same values for m_air_flow and T_out_comp at each iteration. Is this loop necessary, or should you be making use of some other variable (array) in these equations that would depend upon count?

Note also how T_out_comp is not used in your block of code (at least not in the portion that you have included in your question). Should it be used somewhere?

In your second for loop

 for count = (1:count-1)
     % This is the main portion of interests, the initial values are to calculated 
     % once for the 1st iteration
     M_cavern(count) = E_to_store_joule(count) / (Cp_air*T_amb*(Press_r^inverse)-1);
     P_cavern(count) = Default.P_cavern; 
     T_cavern(count) = Default.T_cavern;       

     % Then the 1st iteration results are used, and from that point on the counter 
     % has to ensure that the next iteration will use the previous quantity to 
     % calculate the new quantity              
     M_cavern(count+1) = M_cavern(count)+(m_air_flow-m_air_flow*Default.dm);
     P_cavern(count+1) = P_cavern(count) - (Default.dP*Press_r);
     T_cavern(count+1) = T_cavern(count) - (Default.dt*T_cavern(count));
 end

On each iteration, you update the count and count+1 elements of the M_cavern, P_cavern and T_cavern arrays. On the subsequent iteration, you update the count+1 and count+2 elements of these same arrays and so overwrite the second block of calculations from the previous iteration i.e. {1,2} get updated when count==1, {2,3} get updated when count==2, {3,4} get updated when count==3, etc. Is this intentional, or do you want to skip by two at each iteration

 for u=1:2:count-1

With the above change to the indexing in this loop, you would update {1,2} when u==1, {3,4} when u==3, etc.

Note that the arrays M_cavern, P_cavern and T_cavern are not used outside of this for loop.

The last if block is

if count==1
    Tt3(count) = T_cavern+(T_cavern.*n_pre); 
    Tt4(count) = (Cp_air/Cp_gas)*(lamda*Ma/((lamda*Ma)+1))*...
                  Tt3a+(Hu/(((lamda*Ma)+1)*Cp_gas));
    if Tt4(count)=<Ttmax;
        Tt5(count) = Tt4(count)*(1-n_isT*(1-(1/(expansion_ratio^inverse))));
    end
end

What does count refer to - the length of which array? What are Tt3,*Tt4*, and Tt5?

[George]

Thank you very much for your answer and patience to look the code. The presented code is part of a greater portion of work, which due to its size (and i thought its not appropriate to past a small beast :) as a question). That is why some quantities are not defined here, they have been defined earlier.

I used the count=length(variable name)

because i do not always know the length/ numbers of iteration that may be required

As far as your second remark, with the skip iterations, not i think i have got it wrong, i would like to go every iteration at a time The Tt3,Tt4,Tt5 are quantities that have to be also calulcated

may i ask how you think i should change it?

[Geoff Hayes]

Not knowing enough about what your code is supposed to do, or how it fits into the larger program, I can only review it and point out what strikes me as different or unexpected:

• naming the index variable in the for loop to be the same as the variable that it is indexing over
• the perhaps unnecessary first for loop as nothing depends upon the index
• unused T_out_comp
• the overwriting of the data in the M_cavern, P_cavern and T_cavern arrays on subsequent iterations
• the last if block and what is count supposed to be

For each of the above points, determine whether the code needs to be reviewed further or whether it is written appropriately given how you see the software is to behave.

[George]

It basically tries to simulate a process, at this point I am try to simulate the compression and storage stage of the medium that is why the M_cavern,P_cavern and T_cavern have an initial condition solver and then must be looped over and produced a re-computed for every timestep I have in-putted.

• The first for loop i wrote is if the initial requirement is satisfied, becasue based on the >=< the T_out_comp uses a different solution for m_air_flow and T_out_comp
• I do not know how to stop the over-writing of iteration :(
• You are absolutely right the last if is a fault of mine, i should write for first then the enclosed if in the iteration scheme ( i have update the example code), this last portion is for the temperature criterion check if it is not exceeding a specified threshold

I hope is a bit clearer now

[Geoff Hayes]

Since the first for loop doesn't depend on count, just remove it and use simply

 if E_comp>P_comp
     m_air_flow = (P_comp_joule *E_comp_joule*n_isc*n_mech)/...
         (Cp_air*T_amb*((Press_r^inverse)-1));
     T_out_comp = T_amb*(1+(1/n_isc)*((Press_r^inverse)-1))
 else
     m_air_flow = (P_comp_joule*E_comp_joule *n_isc*n_mech)/...
         (Cp_air*T_amb*(Press_r^inverse));
     T_out_comp = T_amb*(1+(1/n_isc)*(Press_r^inverse))
 end

As for the overwriting of the data on subsequent iterations, just change the step size from the default of 1 to 2. Replace (in the second for loop)

 for count = (1:count-1)

with

 for k=1:2:(count-1)

That way, you will update array indices 1 and 2 when k==1, 3 and 4 when k==3, etc. because now our step size is 2. Try stepping through the code (in the debugger) and observe these changes. (Note that I've switched the index from count to k.)

[George]

Thank you very much for your patience and help, form your first comments i changed the count=1:(count-1) with u , so as not to contradict.

I will try it to run in section by section, hopefully it will run and i think i have understood how you do a loop in reverse order, I will post the proper code when it run so other may have visual help.