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 CenterFile Exchange에서 Time Series Objects에 대해 자세히 알아보기

태그

제품


릴리스

R2016b

Community Treasure Hunt

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

Start Hunting!

Translated by