run takes too long!

4 ビュー (過去 30 日間)
arian hoseini
arian hoseini 2023 年 1 月 21 日
コメント済み: 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
  1 件のコメント
arian hoseini
arian hoseini 2023 年 1 月 21 日
thank u sir

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeTime Series Objects についてさらに検索

タグ

製品


リリース

R2016b

Community Treasure Hunt

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

Start Hunting!

Translated by