Why is my code running for forever?

I'm trying to implement the code for the following system shown in the graph. However, when my N gets larger than 10, it's becoming extremely large. And if I tried N=20 or 30, it just runs for whole night without any result. There's no bug so I really find it hard to modify. The idealized graph should be a sine wave with big amplitudes at first and then attenuate to zero. Could anybody help me check and to modify it?

Thanks a lot!

 % code
clear all;
close all;
syms s; %constants
K1 = 0.4;
K2 = 0.4;
B1 = 0.42;
B2 = 0.4;
D1 = 0.008;
D2 = 0.008;
R1 = 3.00;
R2 = 3.33;
Tg1 = 0.08;
Tg2 = 0.08;
Tt1 = 0.4;
Tt2 = 0.4;
H1 = 0.1667; %value for 2H
H2 = 0.1500;
Pl = 1;
G1 = 0.08;
G2 = 0.08;
T1 = 0.4;
T2 = 0.3;
T12 = 0.2;
T21 = 0.2;
t = 0;
fdev1 = 0;
fdev2 = 0;
suma = 0;
sumb = 0;
N=30;
val_t = zeros(N,1);
val_fdev1 = zeros(N,1);
val_fdev2 = zeros(N,1);
%transfer functions of each part
k1 = ilaplace(K1./s,t);
k2 = ilaplace(K2./s,t);
rot1 = ilaplace(1./(D1 + H1*s),t);
rot2 = ilaplace(1./(D2 + H2*s),t);
c = ilaplace(2*pi/s,t);
for i=1:1:N

 sum11 = suma*c+ fdev1*B1;
 sum12 = sum11*k1-fdev1*(1/R1);
 sum13 = sum12*Tg1-Pl-suma*c;
 fdev1 = sum13*rot1;

 sum21 = sumb*c+ fdev2*B2;
 sum22 = sum21*k2-fdev2*(1/R2);
 sum23 = sum22*Tg2-Pl-sumb*c;
 fdev2 = sum23*rot2;

 suma = T12*fdev1-T21*fdev2;
 sumb = T21*fdev2-T12*fdev1;

 Pl = 0.5*(1-cos(2*pi*t));

 val_t(i)=t;
 val_fdev1(i)=fdev1;
 val_fdev2(i)=fdev2;
 t = t+0.1;

 k1 = ilaplace(K1./s,t);
 k2 = ilaplace(K2./s,t);
 rot1 = ilaplace(1./(D1 + H1*s),t);
 rot2 = ilaplace(1./(D2 + H2*s),t); %rot = 1/(D + 2H*s)
 c = ilaplace(2*pi/s,t);
end
%plot the graph
x = linspace(min(val_t),max(val_t),500); %make the curve smooth
y1 = interp1(val_t,val_fdev1,x,'spline','extrap');
y2 = interp1(val_t,val_fdev2,x,'spline','extrap');
figure(1);
plot(val_t,val_fdev1,'ro')
hold on;
plot(val_t,val_fdev2,'bo')
hold on;
plot(x, y1,'-r','LineWidth', 2)
hold on;
plot(x, y2,'-b','LineWidth', 2)
hold off
grid
title('Two Interconnected Control Areas','fontweight','bold')
xlabel('Time(s)','fontweight','bold');
ylabel('Frequncy Deviation(Hz)','fontweight','bold')

Answer 1

[Stephan]

Hi,

your code is much too complicated:

1. You can calculate the values of all the Transfer functions in one step, by vectorizing the code. They all depend only on time and constants.
2. You can avoid the loop by getting rid of all the sums (sum11, sum12...) but putting all together in a equation for fdev1 and one for fdev2. Then in a seperate step you find out that there are solutions for both if you treat them as a System of equations. This i great, because now you can also calculate These values in a vectorized manner.

This is my version of the code:

syms s;

%constants
K1 = 0.4;
K2 = 0.4;
B1 = 0.42;
B2 = 0.4;
D1 = 0.008;
D2 = 0.008;
R1 = 3.00;
R2 = 3.33;
Tg1 = 0.08;
Tg2 = 0.08;
Tt1 = 0.4;
Tt2 = 0.4;
H1 = 0.1667; %value for 2H
H2 = 0.1500;
G1 = 0.08;
G2 = 0.08;
T1 = 0.4;
T2 = 0.3;
T12 = 0.2;
T21 = 0.2;
N=300;

% transfer functions of each part
t = linspace(0,(N/10),(N+1));
Pl = 0.5*(1-cos(2*pi*t));
k1 = double(ilaplace(K1./s,t));
k2 = double(ilaplace(K2./s,t));
rot1 = double(ilaplace(1./(D1 + H1*s),t));
rot2 = double(ilaplace(1./(D2 + H2*s),t));
c = double(ilaplace(2*pi/s,t));

% Speeding up claculations a bit by using function handles
fdev1_fun = @(Pl,R1,R2,T12,T21,Tg1,Tg2,c,k1,k2,rot1,rot2)(Pl.*R1.*(R2.*-5.0-Tg2.*5.0+R2.*k2.*2.0+R2.*T21.*c.*5.0+R2.*T21.*c.*rot1.*5.0).*-5.0e1)./(R1.*R2.*-2.5e2-R1.*Tg2.*2.5e2-R2.*Tg1.*2.5e2-Tg1.*Tg2.*2.5e2+R1.*R2.*k1.*1.05e2+R1.*R2.*k2.*1.0e2+R1.*Tg2.*k1.*1.05e2+R2.*Tg1.*k2.*1.0e2-R1.*R2.*k1.*k2.*4.2e1+R1.*R2.*T12.*T21.*c.^2.*2.5e2+R1.*R2.*T12.*T21.*c.^2.*rot1.*2.5e2+R1.*R2.*T12.*T21.*c.^2.*rot2.*2.5e2+R1.*R2.*T12.*T21.*c.^2.*rot1.*rot2.*2.5e2);
fdev2_fun = @(Pl,R1,R2,T12,T21,Tg1,Tg2,c,k1,k2,rot1,rot2)(Pl.*R2.*(R1.*-5.0e1-Tg1.*5.0e1+R1.*k1.*2.1e1+R1.*T12.*c.*5.0e1+R1.*T12.*c.*rot2.*5.0e1).*-5.0)./(R1.*R2.*-2.5e2-R1.*Tg2.*2.5e2-R2.*Tg1.*2.5e2-Tg1.*Tg2.*2.5e2+R1.*R2.*k1.*1.05e2+R1.*R2.*k2.*1.0e2+R1.*Tg2.*k1.*1.05e2+R2.*Tg1.*k2.*1.0e2-R1.*R2.*k1.*k2.*4.2e1+R1.*R2.*T12.*T21.*c.^2.*2.5e2+R1.*R2.*T12.*T21.*c.^2.*rot1.*2.5e2+R1.*R2.*T12.*T21.*c.^2.*rot2.*2.5e2+R1.*R2.*T12.*T21.*c.^2.*rot1.*rot2.*2.5e2);
% calculation fdev1 and fdev2
fdev1 = fdev1_fun(Pl,R1,R2,T12,T21,Tg1,Tg2,c,k1,k2,rot1,rot2);
fdev2 = fdev2_fun(Pl,R1,R2,T12,T21,Tg1,Tg2,c,k1,k2,rot1,rot2);

% Plot
figure(1);
plot(t, fdev1,'r')
hold on;
plot(t, fdev2,'b')
hold off
grid
title('Two Interconnected Control Areas','fontweight','bold')
xlabel('Time(s)','fontweight','bold');
ylabel('Frequncy Deviation(Hz)','fontweight','bold')

Using this code with N=300(!) took about 21 seconds on Matlab Online:

Elapsed time is 21.070690 seconds.

But the result is only partially what you expect - it are sin-waves but amplitudes getting bigger with time - so maybe there is an error in the equations(?):

.

You should check this.

-------------------------

How to get to the functions for fdev1 and fdev2?

Do this in a seperate .m-file:

Replace all sum11,sum12,sum13 by their RHS in the equation for fdev - then you get:

fdev1 = T12*fdev1-T21*fdev2*c+ fdev1*B1.*k1-fdev1*(1/R1)*Tg1-Pl-T12*fdev1-T21*fdev2*c.*rot1;

If you do the same for fdev2 you can create a System of equations with 2 unknowns. in fact all other values are known from step 1 or are constants. So do this way:

syms fdev1 fdev2 T12 T21 c k1 k2 R1 R2 Tg1 Tg2 Pl rot1 rot2

eqn1 = fdev1 == T12*fdev1-T21*fdev2*c+ fdev1*B1.*k1-fdev1*(1/R1)*Tg1-Pl-T12*fdev1-T21*fdev2*c.*rot1;
eqn2 = fdev2 == T21*fdev2-T12*fdev1.*c+ fdev2*B2.*k2-fdev2*(1/R2)*Tg2-Pl-T21*fdev2-T12*fdev1*c.*rot2;

[fdev1, fdev2] = solve([eqn1, eqn2], [fdev1, fdev2]);

fdev1_fun = matlabFunction(fdev1)
fdev2_fun = matlabFunction(fdev2)

The both function handles you got can be copied in your script and you are finished.

-------------------------

Can this result be improved? Yes!

Look at your transfer functions - if you convert them back to time domain they are either constants or only depending on time. Why do you calculate this expensive computations for every time step? No need - just get the vaues and do the same like above - vectorize it and let Matlab do what it can realy fast: Calculations on matrices and vectors. Your trnasfer functions in time domain are:

>> rot1

rot1 =

(10000*exp(-(80*t)/1667))/1667

>> rot2

rot2 =

(20*exp(-(4*t)/75))/3

>> c

c =

2*pi

>> Pl

Pl =

1/2 - cos(2*pi*t)/2

>> k1

k1 =

2/5

>> k2

k2 =

2/5

If we change the corresponding code to:

% results of the becktranfered transfer functions of each part
t = linspace(0,(N/10),(N+1));
Pl = 0.5.*(1-cos(2*pi.*t));
k1 = 2/5;
k2 = 2/5;
rot1_fun = @(t)(10000.*exp(-(80.*t)./1667))./1667;
rot1 = rot1_fun(t);
rot2_fun = @(t)(20.*exp(-(4.*t)./75))./3;
rot2 = rot2_fun(t);
c = 2*pi;

we get an execution time for N=300 on Matlab Online:

Elapsed time is 1.711962 seconds.

which is pretty nice.

Best regards

Stephan