How do I Vectorize the for loop containing if else statements?

Hi!

I have a following code snippet containing for loop and if-else statements, I wish to change it into vectorized version for better computational speed.

for k=2:1:length(L)

if TREG(k) >= L(k)

if TREG(k) > L(k)

P_ch(k)=((TREG(k)-L(k))*(BT_eff));

E_ch(k)=P_ch(k).*1;

d1(k)=(SOC_bat_max-SOC_bat(k-1));

if E_ch(k) <= (d1(k))

SOC_bat(k)=SOC_bat(k-1)+E_ch(k);

if SOC_bat(k) > SOC_bat_max

SOC_bat(k)=SOC_bat_max;

E_dump(k)=E_ch(k)-(d1(k));

else

E_dump(k)=0;

end

else

SOC_bat(k)=SOC_bat_max;

E_dump(k)=E_ch(k)-(d1(k));

end

time1(k)=1;

else

SOC_bat(k)=SOC_bat(k-1);

end

elseif L(k) > TREG(k)

P_dh(k)=((L(k)-TREG(k))/(BT_eff));

E_dh(k)=P_dh(k).*1;

diff(k)=SOC_bat(k-1)-SOC_bat_min;

if (diff(k)) >= E_dh(k)

E_dh1(k)=E_dh(k);

SOC_bat(k)= (SOC_bat(k-1))-E_dh(k);

time1(k)=2;

elseif (diff(k)) < E_dh(k)

E_dh1(k)=E_dh(k);

if diff(k)~=0

SOC_bat(k)= E_dh(k)-(diff(k));

else

SOC_bat(k)=SOC_bat(k-1);

end

if E_dh(k)>0

LSBG(k)=E_dh(k);

RPOG(k)=(P_grid*1)-E_dh(k);

P_grid_out(k)=LSBG(k)*1;

end

SOC_bat(k)=SOC_bat(k);

end

if SOC_bat(k) < SOC_bat_min

SOC_bat(k)=SOC_bat_min;

elseif SOC_bat(k) > SOC_bat_max

SOC_bat(k)=SOC_bat_max;

end

end

end

I have tried the following version:

if (TREG > L) == 1

P_ch(TREG > L)=0;

P_ch(:)=((TREG-L)*(BT_eff));

E_ch=P_ch.*1;

d1(TREG >= L & TREG > L)=(SOC_bat_max-SOC_bat);

SOC_bat(TREG >= L & TREG > L & E_ch <= (d1))=SOC_bat+E_ch(k);

SOC_bat(TREG >= L & TREG > L & E_ch <= (d1) & SOC_bat > SOC_bat_max)=SOC_bat_max;

E_dump(TREG >= L & TREG > L & E_ch <= (d1) & SOC_bat > SOC_bat_max)=E_ch-(d1);

E_dump(TREG >= L & TREG > L & E_ch <= (d1) & SOC_bat < SOC_bat_max)=0;

SOC_bat(TREG >= L & TREG > L & E_ch >= (d1))=SOC_bat_max;

E_dump(TREG >= L & TREG > L & E_ch >= (d1))=E_ch-(d1);

SOC_bat(TREG < L)=SOC_bat;

%---------------------------------------------------------------------------

elseif (L > TREG) == 1

%---------------------------------------------------------------------------

P_dh(:)=((L-TREG)./(BT_eff));

P_dh(L > TREG)=0;

E_dh=P_dh.*1;

diff(L > TREG)=(SOC_bat-SOC_bat_min);

SOC_bat(diff>E_dh)=SOC_bat-E_dh(k);

SOC_bat(L > TREG & diff < E_dh & diff ~=0)= E_dh-diff;

SOC_bat(L > TREG & diff < E_dh & diff ==0)= SOC_bat;

LSBG(L > TREG & diff < E_dh & E_dh >0)=E_dh;

RPOG(L > TREG & diff < E_dh & E_dh >0)=(P_grid.*1)-E_dh;

P_grid_out(L > TREG & diff < E_dh & E_dh >0)=LSBG.*1;

SOC_bat(L > TREG & diff < E_dh)=SOC_bat;

SOC_bat(L > TREG & SOC_bat < SOC_bat_min)=SOC_bat_min;

SOC_bat(L > TREG & SOC_bat > SOC_bat_max)=SOC_bat_max;

%---------------------------------------------------------------------------

end

But I am stuck this error, I understand the error. However, I have not found any way to resolve it.

Error @ above code snippet(vectorized version) ( Line: P_dh(:)=((L-TREG)./(BT_eff)) ):

"Error using gpuArray/subsasgn

In an assignment A(:) = B, the number of elements in A and B must be the same."

Please help..

Thak You

Answer 1

[dpb]

In

if (TREG > L) == 1

the "== 1" is superfluous, remove it. The result of a logical expression is 0/1 and if compares to 1 for true so there's no need for the explicit value.

if (TREG > L)

is True IFF (if and only if) ALL elements of the logical array are true -- iow, only if every element of TREG is > every corresponding element of L. This probably isn't so in your data, so this branch will never get executed.

The solution is to use logical addressing/indexing with the vector

P_ch=zeros(size(L));                            % initialize entire array
isGreater = (TREG>L);                           % logical vector; true where elements satisfy
diffT_L=abs(TREG-L);                            % precompute the used differences
P_ch(isGreater)=diffT_L(isGreater)*BT_eff;      % assign the TREG>L elements
P_ch(~isGreater)=diffT_L(~isGreater)/BT_eff;    % assign the TREG<=L elements
...

NB: In your original you do not assign all elements of either P_ch nor the P_dh arrays; I've made the presumption they're the same except for the two cases.

Carry on in similar fashion for the rest of the code; assign the elements of the arrays as needed based on the logical addressing vector.

[Muhammad Abdullah Malik]

Thank you for your response. I am trying your approach and I am getting stuck as to what may be the correct vectorized alternative to SOC_bat(k-1) i.e. In loop version, I can simply set the index value to 'k-1' to get previous step value but in vectors how do I achieve that?.

Code Line: "diff(k)=SOC_bat(k-1)-SOC_bat_min;"

[dpb]

Well, some things simply can't be vectorized; others can be, but the indexing gets to be so convoluted it may not be worth it for code maintainability.

We don't know what SOC_bat and SOC_bat_min are -- are they static array and the minimum of that array? If so, you can write an indexing expression for them; the loop went from k=2:N so k-1 goes from 1:N-1...

diff=[0;SOC_bat(1:N-1)-SOC_bat_min];

If these are also modified during the loop, however, then you simply may not be able to fully vectorize the code..."stuff happens".

[Muhammad Abdullah Malik]

SOC_bat_min is the minimum value that is explicitly defined and yes it is from the array. Moreover, SOC_bat is not a static array as in case for loop version you can see that at each time step a value is added into the array so for each timestep it is modified. So if that is the case then, it won't be vectorized?

How about if i replace the for loop with parfor loop, will that may result in speed benefit or not?

[dpb]

Speed, maybe, but parfor is order independent so squential execution here is necessary.

I've not dug into the code enough to try to see if there is anyway to work around the looping construct or not, but such cases are those in which it may simply not be feasible otherwise without, at least, more machinations than otherwise.

Loops in and of themselves are not necessarily all that time-consuming in MATLAB, anyways...if you are sure to preallocate and store into predined arrays, combined with the JIT compiler/optimizer, performance may not be all that bad anyway.

The thing to do always before worrying much about performance is to code first for clarity and legibility, ensure correctness of the algorithm/implementation, then see if it is fast enough. If performance simply must be improved, then before anything else run the profiler to ensure you know what are actually the performance bottlenecks -- it makes no sense to spend a great deal of time cutting a portion of a code by even 3X if overall it's only 10% of the total.

[Muhammad Abdullah Malik]

Thanks, you were of great help.