run takes too long!

3 次查看(过去 30 天)
arian hoseini
arian hoseini 2023-1-21
i have this code that it runs and answers are correct but it takes so long...why?is there any solution?
clc
clear;
% percent of the relay failure rate
SE=[0.99 0.9 0.8 0.5 0];
ME=[0 0 0 0 0];
eta=0.1;
for k=1:5
syms Tc
% permanent fault (lpf)
l1=2/8760;
% temporary faults (ltf)
l2=15/8760;
% relay failure (lp)
l3=0.01/8760;
% mal-trip failure rate (lext)
l4=0.00001/8760;
% reclosing switching rate (lrct)
l5=10800;
% main switching rate (lmct)
l6=30857.14;
% back-up switching rate (lbct)
l7=8640;
% relay instantaneous mal-trip (lrt-op)
l8=0.001/Tc;
% common cause failure rate of C and P (lcp)
l9=0.000001;
% manual switching rate to restore C (litc)
l10=0.5;
% manual switching rate to restore X (litx)
l11=0.5;
% failure rate of self-cheking system (lsc)
l12=0.002/8760;
% failure rato of monitoring circuit (lmn)
% l13=0.0005/8760;
% routine inseption rates (lrt)
l19=1/Tc;
% number of relay tested by self-cheking (msc)
m1=720;
% number of relay tested by inspected routine test (mrt)
m2=1;
% repair rates (mp)
m3=1;
% repaire rate of C (mc)
m4=1;
% portion of relay failure rate by self-cheking (lpsc)
l14=(1-eta)*l3*SE(k);
% portion of relay failure rate by monitoring (lpmn)
l15=0.0005/8760;
% portion of relay failure rate not detected (lprt)
l16=((1-eta)*l3*(1-SE(k)-ME(k)));
% portion of relay potential mal-trip failure rate (lpt-sc)
l17=eta*l3*SE(k);
% portion of relay potential mal-trip failure rate (lpt-rt)
l18=eta*l3*(1-SE(k)-ME(k));
A(1,1)=1-(l19+l1+l2+l12+l16+l14+l15+l18+l17+l4+l9);
A(1,2)=l1;A(1,4)=l2;A(1,6)=l12;A(1,7)=l19;A(1,8)=l16;A(1,9)=l14;A(1,10)=l15;
A(1,11)=l18;A(1,12)=l17;A(1,13)=l4;A(1,15)=l9;
A(2,2)=1-l6;
A(2,3)=l6;
A(3,3)=1-(l3+m4);
A(3,1)=m4;
A(3,17)=l3;
A(4,4)=1-l6;
A(4,5)=l6;
A(5,5)=1-l5;
A(5,1)=l5;
A(6,6)=1-(l1+l2+m1);
A(6,1)=m1;A(6,15)=(l1+l2);
A(7,7)=1-(m2+l1+l8);
A(7,1)=m2;A(7,13)=l8;A(7,15)=l1;
A(8,8)=1-(l19+l1);
A(8,10)=l19;A(8,15)=l1;
A(9,9)=1-(l12+l1+l2);
A(9,10)=l12;A(9,15)=l1+l2;
A(10,10)=1-(l1+l2+m3);
A(10,1)=m3;A(10,15)=l1+l2;
A(11,11)=1-(l19+l1+l3);
A(11,10)=l19;A(11,13)=l3;A(11,15)=l1;
A(12,12)=1-(l12+l3+l1+l2);
A(12,10)=l12;A(12,13)=l3;A(12,15)=l1+l2;
A(13,13)=1-l6;
A(13,14)=l6;
A(14,14)=1-(l10+l10);
A(14,1)=l10;A(14,10)=l10;
A(15,15)=1-l7;
A(15,16)=l7;
A(16,16)=1-l11;
A(16,17)=l11;
A(17,17)=1-(m3+m4);
A(17,3)=m3;A(17,10)=m4;
A(17,:)=1;
B=[0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1];
% format long
P=inv(A)*B;
% R=VPA(P);
% disp(P)
% disp('P1=')
% disp(P(1,1))
% disp('P2=')
% disp(P(2,1))
% disp('P3=')
% disp(P(3,1))
% disp('P4=')
% disp(P(4,1))
% disp('P5=')
% disp(P(5,1))
% disp('P6=')
% disp(P(6,1))
% disp('P7=')
% disp(P(7,1))
% disp('P8=')
% disp(P(8,1))
% disp('P9=')
% disp(P(9,1))
% disp('P10=')
% disp(P(10,1))
% disp('P11=')
% disp(P(11,1))
% disp('P12=')
% disp(P(12,1))
% disp('P13=')
% disp(P(13,1))
% disp('P14=')
% disp(P(14,1))
% disp('P15=')
% disp(P(15,1))
% disp('P16=')
% disp(P(16,1))
% disp('P17=')
% disp(P(17,1))
PI(k)=P(1);
PII(k)=P(2,1)+P(3,1)+P(4,1)+P(5,1);
PIII(k)=P(6,1)+P(7,1)+P(8,1)+P(9,1)+P(10,1)+P(11,1)+P(12,1);
PIV(k)=P(15,1)+P(16,1)+P(17,1);
PV(k)=P(13,1)+P(14,1);
% PI=double(PI);
end
for c=1:1:500
%%
% Pi1(c)=subs(PI(1),Tc,c);
% Pi1(c)=double(Pi1(c));
%
% Pi2(c)=subs(PI(2),Tc,c);
% Pi2(c)=double(Pi2(c));
%
% Pi3(c)=subs(PI(3),Tc,c);
% Pi3(c)=double(Pi3(c));
%
% Pi4(c)=subs(PI(4),Tc,c);
% Pi4(c)=double(Pi4(c));
%
% Pi5(c)=subs(PI(5),Tc,c);
% Pi5(c)=double(Pi5(c));
%%
Piii1(c)=subs(PIII(1),Tc,c);
Piii1(c)=double(Piii1(c));
Piii2(c)=subs(PIII(2),Tc,c);
Piii2(c)=double(Piii2(c));
Piii3(c)=subs(PIII(3),Tc,c);
Piii3(c)=double(Piii3(c));
Piii4(c)=subs(PIII(4),Tc,c);
Piii4(c)=double(Piii4(c));
Piii5(c)=subs(PIII(5),Tc,c);
Piii5(c)=double(Piii5(c));
%%
Piv1(c)=subs(PIV(1),Tc,c);
Piv1(c)=double(Piv1(c));
Piv2(c)=subs(PIV(2),Tc,c);
Piv2(c)=double(Piv2(c));
Piv3(c)=subs(PIV(3),Tc,c);
Piv3(c)=double(Piv3(c));
Piv4(c)=subs(PIV(4),Tc,c);
Piv4(c)=double(Piv4(c));
Piv5(c)=subs(PIV(5),Tc,c);
Piv5(c)=double(Piv5(c));
% Pv(c)=subs(PV,Tc,c);
% Pv(c)=double(Pv(c));
end
% figure(1)
% plot(Pi1(150:end));
% hold on
% plot(Pi2(150:end));
% hold on
% plot(Pi3(150:end));
% hold on
% plot(Pi4(150:end));
% hold on
% plot(Pi5(150:end));
figure(2)
% plot(Pii(150:end));
% hold on
% figure(3)
plot(Piii1(150:end));
hold on
plot(Piii2(150:end));
hold on
plot(Piii3(150:end));
hold on
plot(Piii4(150:end));
hold on
plot(Piii5(150:end));
figure(3)
plot(Piv1(150:end));
hold on
plot(Piv2(150:end));
hold on
plot(Piv3(150:end));
hold on
plot(Piv4(150:end));
hold on
plot(Piv5(150:end));
% figure(5)
% plot(Pv(150:end));
% hold on

采纳的回答

Matt J
Matt J 2023-1-21
编辑:Matt J 2023-1-21
i have this code that it runs and answers are correct but it takes so long...why?is there any solution?
It is unnecessary to use sym variables here. It is much faster to use numeric variables:
clc
clear;
% percent of the relay failure rate
tic;
for Tc=1:500
out(Tc)=doIt(Tc);
end
toc
Elapsed time is 0.111667 seconds.
Piii=[out.PIII];
Piv=[out.PIV];
figure(2)
plot(Piii(:,150:end)'); legend
figure(3)
plot(Piv(:,150:end)'); legend
function out=doIt(Tc)
SE=[0.99 0.9 0.8 0.5 0];
ME=[0 0 0 0 0];
eta=0.1;
for k=1:5
% permanent fault (lpf)
l1=2/8760;
% temporary faults (ltf)
l2=15/8760;
% relay failure (lp)
l3=0.01/8760;
% mal-trip failure rate (lext)
l4=0.00001/8760;
% reclosing switching rate (lrct)
l5=10800;
% main switching rate (lmct)
l6=30857.14;
% back-up switching rate (lbct)
l7=8640;
% relay instantaneous mal-trip (lrt-op)
l8=0.001/Tc;
% common cause failure rate of C and P (lcp)
l9=0.000001;
% manual switching rate to restore C (litc)
l10=0.5;
% manual switching rate to restore X (litx)
l11=0.5;
% failure rate of self-cheking system (lsc)
l12=0.002/8760;
% failure rato of monitoring circuit (lmn)
% l13=0.0005/8760;
% routine inseption rates (lrt)
l19=1/Tc;
% number of relay tested by self-cheking (msc)
m1=720;
% number of relay tested by inspected routine test (mrt)
m2=1;
% repair rates (mp)
m3=1;
% repaire rate of C (mc)
m4=1;
% portion of relay failure rate by self-cheking (lpsc)
l14=(1-eta)*l3*SE(k);
% portion of relay failure rate by monitoring (lpmn)
l15=0.0005/8760;
% portion of relay failure rate not detected (lprt)
l16=((1-eta)*l3*(1-SE(k)-ME(k)));
% portion of relay potential mal-trip failure rate (lpt-sc)
l17=eta*l3*SE(k);
% portion of relay potential mal-trip failure rate (lpt-rt)
l18=eta*l3*(1-SE(k)-ME(k));
A(1,1)=1-(l19+l1+l2+l12+l16+l14+l15+l18+l17+l4+l9);
A(1,2)=l1;A(1,4)=l2;A(1,6)=l12;A(1,7)=l19;A(1,8)=l16;A(1,9)=l14;A(1,10)=l15;
A(1,11)=l18;A(1,12)=l17;A(1,13)=l4;A(1,15)=l9;
A(2,2)=1-l6;
A(2,3)=l6;
A(3,3)=1-(l3+m4);
A(3,1)=m4;
A(3,17)=l3;
A(4,4)=1-l6;
A(4,5)=l6;
A(5,5)=1-l5;
A(5,1)=l5;
A(6,6)=1-(l1+l2+m1);
A(6,1)=m1;A(6,15)=(l1+l2);
A(7,7)=1-(m2+l1+l8);
A(7,1)=m2;A(7,13)=l8;A(7,15)=l1;
A(8,8)=1-(l19+l1);
A(8,10)=l19;A(8,15)=l1;
A(9,9)=1-(l12+l1+l2);
A(9,10)=l12;A(9,15)=l1+l2;
A(10,10)=1-(l1+l2+m3);
A(10,1)=m3;A(10,15)=l1+l2;
A(11,11)=1-(l19+l1+l3);
A(11,10)=l19;A(11,13)=l3;A(11,15)=l1;
A(12,12)=1-(l12+l3+l1+l2);
A(12,10)=l12;A(12,13)=l3;A(12,15)=l1+l2;
A(13,13)=1-l6;
A(13,14)=l6;
A(14,14)=1-(l10+l10);
A(14,1)=l10;A(14,10)=l10;
A(15,15)=1-l7;
A(15,16)=l7;
A(16,16)=1-l11;
A(16,17)=l11;
A(17,17)=1-(m3+m4);
A(17,3)=m3;A(17,10)=m4;
A(17,:)=1;
B=[0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1];
% format long
P=inv(A)*B;
out.PI(k,1)=P(1);
out.PII(k,1)=P(2,1)+P(3,1)+P(4,1)+P(5,1);
out.PIII(k,1)=P(6,1)+P(7,1)+P(8,1)+P(9,1)+P(10,1)+P(11,1)+P(12,1);
out.PIV(k,1)=P(15,1)+P(16,1)+P(17,1);
out.PV(k,1)=P(13,1)+P(14,1);
% PI=double(PI);
end
end
  2 个评论
arian hoseini
arian hoseini 2023-1-21
thank u adn another question if u see the plot data1,2,3,4,...are very close and i saw something in articles an papers that they zoom on that part but the zoom plot is on the main plot...can u help me with that.
Matt J
Matt J 2023-1-21
Here is one way you can bring out the differences in the plots, but as you can see they are very tiny. I can't imagine what significance the differences could have:
for Tc=1:1:500
out(Tc)=doIt(Tc);
end
Piii=[out.PIII]; Piii=Piii(2:5,:)-Piii(1,:);
Piv=[out.PIV]; Piv=Piv(2:5,:)-Piv(1,:);
figure(2)
semilogy(Piii(:,150:end)'); legend
figure(3)
semilogy(Piv(:,150:end)'); legend

请先登录,再进行评论。

更多回答(1 个)

the cyclist
the cyclist 2023-1-21
Vectorizing the for loop gives a nice speedup. There are almost certainly other optimizations possible.
clc
clear;
% percent of the relay failure rate
SE=[0.99 0.9 0.8 0.5 0];
ME=[0 0 0 0 0];
eta=0.1;
for k=1:5
syms Tc
% permanent fault (lpf)
l1=2/8760;
% temporary faults (ltf)
l2=15/8760;
% relay failure (lp)
l3=0.01/8760;
% mal-trip failure rate (lext)
l4=0.00001/8760;
% reclosing switching rate (lrct)
l5=10800;
% main switching rate (lmct)
l6=30857.14;
% back-up switching rate (lbct)
l7=8640;
% relay instantaneous mal-trip (lrt-op)
l8=0.001/Tc;
% common cause failure rate of C and P (lcp)
l9=0.000001;
% manual switching rate to restore C (litc)
l10=0.5;
% manual switching rate to restore X (litx)
l11=0.5;
% failure rate of self-cheking system (lsc)
l12=0.002/8760;
% failure rato of monitoring circuit (lmn)
% l13=0.0005/8760;
% routine inseption rates (lrt)
l19=1/Tc;
% number of relay tested by self-cheking (msc)
m1=720;
% number of relay tested by inspected routine test (mrt)
m2=1;
% repair rates (mp)
m3=1;
% repaire rate of C (mc)
m4=1;
% portion of relay failure rate by self-cheking (lpsc)
l14=(1-eta)*l3*SE(k);
% portion of relay failure rate by monitoring (lpmn)
l15=0.0005/8760;
% portion of relay failure rate not detected (lprt)
l16=((1-eta)*l3*(1-SE(k)-ME(k)));
% portion of relay potential mal-trip failure rate (lpt-sc)
l17=eta*l3*SE(k);
% portion of relay potential mal-trip failure rate (lpt-rt)
l18=eta*l3*(1-SE(k)-ME(k));
A(1,1)=1-(l19+l1+l2+l12+l16+l14+l15+l18+l17+l4+l9);
A(1,2)=l1;A(1,4)=l2;A(1,6)=l12;A(1,7)=l19;A(1,8)=l16;A(1,9)=l14;A(1,10)=l15;
A(1,11)=l18;A(1,12)=l17;A(1,13)=l4;A(1,15)=l9;
A(2,2)=1-l6;
A(2,3)=l6;
A(3,3)=1-(l3+m4);
A(3,1)=m4;
A(3,17)=l3;
A(4,4)=1-l6;
A(4,5)=l6;
A(5,5)=1-l5;
A(5,1)=l5;
A(6,6)=1-(l1+l2+m1);
A(6,1)=m1;A(6,15)=(l1+l2);
A(7,7)=1-(m2+l1+l8);
A(7,1)=m2;A(7,13)=l8;A(7,15)=l1;
A(8,8)=1-(l19+l1);
A(8,10)=l19;A(8,15)=l1;
A(9,9)=1-(l12+l1+l2);
A(9,10)=l12;A(9,15)=l1+l2;
A(10,10)=1-(l1+l2+m3);
A(10,1)=m3;A(10,15)=l1+l2;
A(11,11)=1-(l19+l1+l3);
A(11,10)=l19;A(11,13)=l3;A(11,15)=l1;
A(12,12)=1-(l12+l3+l1+l2);
A(12,10)=l12;A(12,13)=l3;A(12,15)=l1+l2;
A(13,13)=1-l6;
A(13,14)=l6;
A(14,14)=1-(l10+l10);
A(14,1)=l10;A(14,10)=l10;
A(15,15)=1-l7;
A(15,16)=l7;
A(16,16)=1-l11;
A(16,17)=l11;
A(17,17)=1-(m3+m4);
A(17,3)=m3;A(17,10)=m4;
A(17,:)=1;
B=[0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1];
% format long
P=inv(A)*B;
% R=VPA(P);
% disp(P)
% disp('P1=')
% disp(P(1,1))
% disp('P2=')
% disp(P(2,1))
% disp('P3=')
% disp(P(3,1))
% disp('P4=')
% disp(P(4,1))
% disp('P5=')
% disp(P(5,1))
% disp('P6=')
% disp(P(6,1))
% disp('P7=')
% disp(P(7,1))
% disp('P8=')
% disp(P(8,1))
% disp('P9=')
% disp(P(9,1))
% disp('P10=')
% disp(P(10,1))
% disp('P11=')
% disp(P(11,1))
% disp('P12=')
% disp(P(12,1))
% disp('P13=')
% disp(P(13,1))
% disp('P14=')
% disp(P(14,1))
% disp('P15=')
% disp(P(15,1))
% disp('P16=')
% disp(P(16,1))
% disp('P17=')
% disp(P(17,1))
PI(k)=P(1);
PII(k)=P(2,1)+P(3,1)+P(4,1)+P(5,1);
PIII(k)=P(6,1)+P(7,1)+P(8,1)+P(9,1)+P(10,1)+P(11,1)+P(12,1);
PIV(k)=P(15,1)+P(16,1)+P(17,1);
PV(k)=P(13,1)+P(14,1);
% PI=double(PI);
end
cvec = 1:500;
Piii1=double(subs(PIII(1),Tc,cvec));
Piii2=double(subs(PIII(2),Tc,cvec));
Piii3=double(subs(PIII(3),Tc,cvec));
Piii4=double(subs(PIII(4),Tc,cvec));
Piii5=double(subs(PIII(5),Tc,cvec));
Piv1=double(subs(PIV(1),Tc,cvec));
Piv2=double(subs(PIV(2),Tc,cvec));
Piv3=double(subs(PIV(3),Tc,cvec));
Piv4=double(subs(PIV(4),Tc,cvec));
Piv5=double(subs(PIV(5),Tc,cvec));
% figure(1)
% plot(Pi1(150:end));
% hold on
% plot(Pi2(150:end));
% hold on
% plot(Pi3(150:end));
% hold on
% plot(Pi4(150:end));
% hold on
% plot(Pi5(150:end));
figure(2)
% plot(Pii(150:end));
% hold on
% figure(3)
plot(Piii1(150:end));
hold on
plot(Piii2(150:end));
hold on
plot(Piii3(150:end));
hold on
plot(Piii4(150:end));
hold on
plot(Piii5(150:end));
figure(3)
plot(Piv1(150:end));
hold on
plot(Piv2(150:end));
hold on
plot(Piv3(150:end));
hold on
plot(Piv4(150:end));
hold on
plot(Piv5(150:end));
% figure(5)
% plot(Pv(150:end));
% hold on

类别

Help CenterFile Exchange 中查找有关 Continuous Waveforms 的更多信息

标签

产品


版本

R2016b

Community Treasure Hunt

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

Start Hunting!

Translated by