# How to replace a "for loop" with vectorization?

219 views (last 30 days)
Ismaeel on 30 Apr 2018
Edited: Ismaeel on 8 Jun 2019
Hi everyone, I have built in Simulink a user-defined function block for solving a set of algebraic nonlinear equations (with sine and cosine). It works alright. However, I used for loop to construct the equations, which is very slow when the system size increases. It takes up to one hour to solve a system consisting of 136 equations with a reasonable accuracy and step size. I have read that vectorizing can enormously speed up the simulation. I was not able to do so, however. Any help is greatly appreciated. This is my code (part of it):
function outt = DAE_Full01(x)
n=68;m=16;
%======= Defining the variables=====================
S=1;
for i1=1:m
Id(i1)=x(S);S=S+1;
end
for i2=1:m
Iq(i2)=x(S);S=S+1;
end
for i3=1:n
V(i3)=x(S);S=S+1;
end
for i4=1:n
TH(i4)=x(S);S=S+1;
end
%===End of defining the variables===============
%=== Construct the equations ================
for i6=1:m
SE1(i6,:)=Edp(i6)-V(i6)*sin(D(i6)-TH(i6))-Rs(i6)*Id(i6)+Xqp(i6)*Iq(i6);
SE2(i6,:)=Eqp(i6)-V(i6)*cos(D(i6)-TH(i6))-Rs(i6)*Iq(i6)-Xdp(i6)*Id(i6);
end
.
.
.
[xx ]=solve...
Where Rs, Xdp, Xqp are constant values. Id, Iq, Edp, Eqp, and D are inputs to the block (variables) with size m. x0 is a vector of initial values already defined. How one can replace the above for loop with vectorization form?

Stephan on 3 May 2018
Edited: Stephan on 3 May 2018
Hi,
I finished ;-)
You could change the code as follows:
1.) To replace the nested loop this is needed:
[b, a] = ndgrid(1:m, 1:n);
2.) I made a function for this calculation because the Performance of a function against a script is usually better. I tested both of it and found the function faster. Call the function this way:
[PV1, PV2] = vectorized_run(a, b, m,n, V, TH, Iq, Id, IC, Yabs, Yang, D);
3.) Due to matlab wants defintions of a function behind all other code, define the function as follows:
function [PV1, PV2] = vectorized_run(a, b, m, n, V, TH, Iq, Id, IC, Yabs, Yang, D)
% The first step saves calculation time, because the first part of both nested
% loops is the same for calculating sum1 and sum2 - same for PV_common
sum_common = V(b).*V(a) .* Yabs(1:m,1:n);
sum1 = sum((sum_common .* cos(TH(b)-TH(a)-Yang(1:m,1:n))),2);
sum2 = sum((sum_common .* sin(TH(b)-TH(a)-Yang(1:m,1:n))),2);
PV_common = (Id(1:m) .* V(1:m));
PV1 = diag(PV_common .* (sin(D(1:m)-TH(1:m))) + Iq(1:m).*V(1:m) .* (cos(D(1:m)-TH(1:m))) - IC(1:m,5) - sum1(1:m));
PV2 = diag(PV_common .* (cos(D(1:m)-TH(1:m))) - Iq(1:m).*V(1:m) .* (sin(D(1:m)-TH(1:m))) - IC(1:m,6) - sum2(1:m));
end
Thats it ;-)
When i tested with different sizes of the problem the vectorized code took 33...40% less time to calculate PV1 and PV2. I guess that's a pretty good result.
Please test for your purposes and give me a Feedback to this.
Best regards
Stephan
Stephan on 6 May 2018
Edited: Stephan on 6 May 2018
Hi Ismael,
i have some pretty good news for you. i could solve the problems regarding the global variables and improve the performance of sum3 and sum4 so that we got this result:
To reach this all funtions are now in only .m-file. The overall execution time reaced <1s and the part for solving DAE is now clearly below 600ms - in the attached test it is 576ms.
You should still keep both variants, because my computer today obviously has a good day and the variant of previously also achieved very similar results. I can not say which of the two variants works better in the end for larger problems - please test. Therefore, I will send you 3 files by email:
1.) Main_vec.m & DAE_vec.m
These now additonaly include the improved code for sum3 and sum4 for optimal performance, but still work with global variables.
2.) solve_DAE_fast.m
This file contains all the improvements that have been incorporated into the other variant, but eliminates the global variables. If you want to know more about how this works, have a look at this link:
So my suggestions for the improvement of your code are implemented and now I really have nothing more to improve ... (or maybe even later? ;-)
Best regards and keep me up to date please
Stephan

Stephan on 1 May 2018
Edited: Stephan on 3 May 2018
Hi,
You can find a good start to this topic in the link below.
You assign the decision variables x(1) ... x(168) in the 4 loops by increasing S successively. Instead you could omit S and take advantage of the known quantities of m and n like this:
Id(1:m)=x(1:m); % x(1) ... x(16)
Iq(1:m)=x(m+1:2*m); % x(17) ... x(32)
V(1:n)=x(2*m+1:2*m+n); % x(33) ... x(100)
TH(1:n)=x(2*m+1+n:2*m+2*n); % x(101) ... x(168))
I hope that's what you wanted to achieve in your code.
The equations below in the code should also work without the for loop:
SE1(1:m,:)=Edp(1:m)-V(1:m)*sin(D(1:m)-TH(1:m))-Rs(1:m)*Id(1:m)+Xqp(1:m)*Iq(1:m);
SE2(1:m,:)=Eqp(1:m)-V(1:m)*cos(D(1:m)-TH(1:m))-Rs(1:m)*Iq(1:m)-Xdp(1:m)*Id(1:m);
Since I could not test the code, I would be very grateful for any feedback as to whether it works. I have found that with every help that I offer here in the forum, I am learning and would be happy if you contributed with your feedback.
I would also be interested to know if you have seen an improvement in the performance of your code and how large it is.
Best regards
Stephan
Stephan on 3 May 2018
Edited: Stephan on 3 May 2018
Hi,
i could find the issue... the sums sum1 and sum2 work fine now. i saved about 50% calculation time for the loops for the m=68 and n=16 dimensions. Excited to hear about the runtime behavior for larger problems.
This works for both sums, but still i have to look at the PV values.
The matrices and vectors i have made are only for comparision of the results with your formulas to eliminate possible errors that i could have made. If my calculation brings the same values as your calculation i made no errors. I guess that should work for this purpose.
I keep you updated ;-)
Best regards
Stephan

### Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by