From b549ce7fef279be46cf6dbc5991382289682e8d9 Mon Sep 17 00:00:00 2001 From: Giulio Masetti Date: Tue, 11 Jul 2023 16:08:11 +0200 Subject: [PATCH] Initial commit --- ACload/CTMC_CR.m | 41 ++++++++++++++++++++++ ACload/CTMC_SR.m | 29 ++++++++++++++++ ACload/compare.m | 71 ++++++++++++++++++++++++++++++++++++++ ACload/compare_lambda.m | 76 +++++++++++++++++++++++++++++++++++++++++ ACload/evaluate_CR.m | 34 ++++++++++++++++++ ACload/evaluate_SR.m | 35 +++++++++++++++++++ DCload/CTMC_CR.m | 32 +++++++++++++++++ DCload/CTMC_SR.m | 20 +++++++++++ DCload/compare.m | 71 ++++++++++++++++++++++++++++++++++++++ DCload/evaluate_CR.m | 34 ++++++++++++++++++ DCload/evaluate_SR.m | 35 +++++++++++++++++++ 11 files changed, 478 insertions(+) create mode 100644 ACload/CTMC_CR.m create mode 100644 ACload/CTMC_SR.m create mode 100644 ACload/compare.m create mode 100644 ACload/compare_lambda.m create mode 100644 ACload/evaluate_CR.m create mode 100644 ACload/evaluate_SR.m create mode 100644 DCload/CTMC_CR.m create mode 100644 DCload/CTMC_SR.m create mode 100644 DCload/compare.m create mode 100644 DCload/evaluate_CR.m create mode 100644 DCload/evaluate_SR.m diff --git a/ACload/CTMC_CR.m b/ACload/CTMC_CR.m new file mode 100644 index 0000000..ac3f117 --- /dev/null +++ b/ACload/CTMC_CR.m @@ -0,0 +1,41 @@ +function [nstates, Qav, Qre, pi0, av_indx] = CTMC_CR(nm,n,lambdaB,lambdaR,muB,muR) +%CTMC_CR infinitesimal generator matrix and initial probability vector +% of the Component Redundancy (CR) Continuous Time Markov Chain (CTMC) +% nm : number of hot standby (n+m>=n+1) +% n : number of required rectifiers or batteries to address load +% lambdaB : battery failure rate +% lambdaR : rectifier failure rate +% muB : battery recovery rate +% muR : rectifier recovery rate + +aus = nm-n+2;% number of states only considering rectifiers +nstates = (nm-n+2)^2; + +eaus = zeros(1,nm-n+2); +eaus(end) = 1; +av_indx = find(kron(ones(1,nm-n+2), eaus)+kron(eaus, ones(1,nm-n+2))<1); + +% define rectifiers' and batteries' behavior +Rr_failure = diag(lambdaR*[nm:-1:n],1); +Rb_failure = diag(lambdaB*[nm:-1:n],1); + +Rre = kron(Rr_failure, eye(aus, aus)) + kron(eye(aus, aus), Rb_failure); +Rre(av_indx(2:end),1) = muB; + +Rav = Rre; +Rav(2:end,1) = muB; + +% % Plotting reliability model +% D= digraph(Rre); +% Plot = plot(D); +% highlight(Plot,[5,10,15,20,21,22,23,24,25],'NodeColor','r'); +% highlight(Plot,'Edges',[inedges(D,5)',inedges(D,10)',inedges(D,15)',inedges(D,20)',... +% inedges(D,21)',inedges(D,22)',inedges(D,23)',inedges(D,24)',inedges(D,25)'],'EdgeColor','r'); + +Qav = Rav - diag(Rav*ones(nstates,1)); +Qre = Rre - diag(Rre*ones(nstates,1)); + +pi0 = zeros(1,nstates); +pi0(1) = 1; + +end \ No newline at end of file diff --git a/ACload/CTMC_SR.m b/ACload/CTMC_SR.m new file mode 100644 index 0000000..3fa127e --- /dev/null +++ b/ACload/CTMC_SR.m @@ -0,0 +1,29 @@ +function [Q,pi0] = CTMC_SR(lambdaB,lambdaR,muB,muR) +%CTMC_SR infinitesimal generator matrix and initial probability vector +% of the System Redundancy (SR) Continuous Time Markov Chain (CTMC) +% lambdaB : battery failure rate +% lambdaR : rectifier failure rate +% muB : battery recovery rate +% muR : rectifier recovery rate + +lambda = lambdaB + lambdaR; + +R = zeros(3,3); +R(1,2) = 2*lambda; +R(2,1) = 2*muB; +R(2,3) = lambda; + +% % Plot reliability model +% D= digraph(R); +% Plot = plot(D); +% highlight(Plot,3,'NodeColor','r'); +% highlight(Plot,'Edges',inedges(D,3)','EdgeColor','r'); + +R(3,2) = muB; + +Q = R - diag(R*ones(3,1)); + +pi0 = zeros(1,3); +pi0(1) = 1; + +end \ No newline at end of file diff --git a/ACload/compare.m b/ACload/compare.m new file mode 100644 index 0000000..cc2c86d --- /dev/null +++ b/ACload/compare.m @@ -0,0 +1,71 @@ +%% parameters setting +% System Redundancy (SR) +lambdaB_SR = 11.76e-6;% hours^(-1) +lambdaR_SR = 16.6e-6; +muB_SR = 0.5;% hours^(-1) +muR_SR = muB_SR; + +% Component Redundancy (CR) +lambdaB_CR = 17*lambdaB_SR; +lambdaR_CR = 17*lambdaR_SR; +muB_CR = 0.5;% hours^(-1) +muR_CR = muR_SR; + + +ModulePower = 10;% kW +minLoad = 40;% kW +loadStep = 20;% kW + +%% start evaluation +numberOfEvaluations = 6; + +MTTF_perc_impr = zeros(numberOfEvaluations,1); +MTBF_perc_impr = zeros(numberOfEvaluations,1); +MTBF_SR = zeros(numberOfEvaluations,1); +MTBF_CR = zeros(numberOfEvaluations,1); +Energy_perc_impr = zeros(numberOfEvaluations,1); +for i = 1 : numberOfEvaluations + Load = minLoad + i*loadStep;% kW, typical values from 50 to 100 kW + + n = ceil(Load/ModulePower); + nm = n+3;% nm is n+m + + PotTransformerMonolite=0.04*2*Load; + PotTransformerModular=0.04*Load; + + Energy_perc_impr(i) = 100; % estimate based only on transformers' consumption + + [A_SR,MTTF_SR,MTBF_SR_value] = evaluate_SR(nm,n,lambdaB_SR,lambdaR_SR,muB_SR,muR_SR); + [A_CR,MTTF_CR,MTBF_CR_value] = evaluate_CR(nm,n,lambdaB_CR,lambdaR_CR,muB_CR,muR_CR); + + MTBF_SR(i) = MTBF_SR_value; + MTBF_CR(i) = MTBF_CR_value; + + A_perc_impr = 100*(A_CR-A_SR)/A_SR; + MTTF_perc_impr(i) = 100*(MTTF_CR-MTTF_SR)/MTTF_SR; + MTBF_perc_impr(i) = 100*(MTBF_CR_value-MTBF_SR_value)/MTBF_SR_value; + +end + +%% Plot +x = minLoad : loadStep : (minLoad+(numberOfEvaluations-1)*loadStep); +zeroline=zeros(1,length(x)); + +fs = 14; +f1=figure(1); +f1.Position = [50 50 650 650]; +title('Percentage of improvement','FontSize',fs-2); +hold on; +xlabel('DC load [kW]','FontWeight','bold','FontSize',fs); +plot(x,MTBF_perc_impr,'Color',[0,0,0],'LineWidth',2); +plot(x,zeroline,'Color',[0.4,0.4,0.4]); + +f2=figure(2); +f2.Position = [50 50 650 650]; +title('Mean Time Between two Failures (MTBF)','FontSize',fs-2); +hold on; +xlabel('DC load [kW]','FontWeight','bold','FontSize',fs); +ylabel('MTBF [hours]','FontWeight','bold','FontSize',fs); +plot(x,MTBF_SR,"--",'Color',[0.5,0,0.5],'LineWidth',2); +plot(x,MTBF_CR,"-.",'Color',[0,0.7,0],'LineWidth',2); +legend('SR','CR','FontSize',fs); \ No newline at end of file diff --git a/ACload/compare_lambda.m b/ACload/compare_lambda.m new file mode 100644 index 0000000..1bfd951 --- /dev/null +++ b/ACload/compare_lambda.m @@ -0,0 +1,76 @@ +%% parameters setting +% System Redundancy (SR) +lambdaB_SR = 11.76e-6;% hours^(-1) +lambdaR_SR = 16.6e-6; +muB_SR = 0.5;% hours^(-1) +muR_SR = muB_SR; + +% % Component Redundancy (CR) +% lambdaB_CR = 17*lambdaB_SR; +% lambdaR_CR = 17*lambdaR_SR; +muB_CR = 0.5;% hours^(-1) +muR_CR = muR_SR; + + +ModulePower = 10;% kW +Load = 40 + 6*20;% kW, typical values from 50 to 100 kW + +min = 1; +step = 1; +max = 20; + +%% start evaluation +MTTF_perc_impr = zeros(max,1); +MTBF_perc_impr = zeros(max,1); +MTBF_SR = zeros(max,1); +MTBF_CR = zeros(max,1); +Energy_perc_impr = zeros(max,1); +for i = min : step : max + + % Component Redundancy (CR) + lambdaB_CR = i*lambdaB_SR; + lambdaR_CR = i*lambdaR_SR; + + n = ceil(Load/ModulePower); + nm = n+3;% nm is n+m + + PotTransformerMonolite=0.04*2*Load; + PotTransformerModular=0.04*Load; + + Energy_perc_impr(i) = 100; % estimate based only on transformers' consumption + + [A_SR,MTTF_SR,MTBF_SR_value] = evaluate_SR(nm,n,lambdaB_SR,lambdaR_SR,muB_SR,muR_SR); + [A_CR,MTTF_CR,MTBF_CR_value] = evaluate_CR(nm,n,lambdaB_CR,lambdaR_CR,muB_CR,muR_CR); + + MTBF_SR(i) = MTBF_SR_value; + MTBF_CR(i) = MTBF_CR_value; + + A_perc_impr = 100*(A_CR-A_SR)/A_SR; + MTTF_perc_impr(i) = 100*(MTTF_CR-MTTF_SR)/MTTF_SR; + MTBF_perc_impr(i) = 100*(MTBF_CR_value-MTBF_SR_value)/MTBF_SR_value; + +end + +%% Plot +x = min : step : (min+(max-1)*step); +zeroline=zeros(1,length(x)); +positive=find(MTBF_perc_impr>0); +negative=find(MTBF_perc_impr<=0); + +figure(1); +semilogy(x(positive),MTBF_perc_impr(positive),'Color',[0,0,0],'LineWidth',2); +hold on; +% semilogy(x(negative),MTBF_perc_impr(negative),'Color',[0,0,0],'LineWidth',2); +semilogy(x(negative),10.^(-log10(-MTBF_perc_impr(negative))),'Color',[0,0,0],'LineWidth',2) +% semilogy(x,zeroline,'Color',[0.4,0.4,0.4]); +title(strcat('Percentage of improvement, load=', num2str(Load),' kW')); +xlabel('r'); + +figure(2); +semilogy(x,MTBF_SR,"--",'Color',[0.5,0,0.5],'LineWidth',2); +hold on; +semilogy(x,MTBF_CR,"-.",'Color',[0,0.7,0],'LineWidth',2); +title(strcat('Mean Time Between two Failures (MTBF), load=', num2str(Load),' kW')); +xlabel('r'); +ylabel('MTBF [hours]'); +legend('SR','CR'); \ No newline at end of file diff --git a/ACload/evaluate_CR.m b/ACload/evaluate_CR.m new file mode 100644 index 0000000..879085d --- /dev/null +++ b/ACload/evaluate_CR.m @@ -0,0 +1,34 @@ +function [A,MTTF,MTBF] = evaluate_CR(nm,n,lambdaB,lambdaR,muB,muR) + +%% model definition +% notice: here last state is, by definition, one of the system failed states +[nstates,Qav,Qre,pi0,av_indx] = CTMC_CR(nm,n,lambdaB,lambdaR,muB,muR); + +%% solution of the availability model +tildeQ = Qav; +tildeQ(:,end) = ones(nstates,1); + +elast = zeros(1,nstates); +elast(end) = 1; + +% steady-state probability vector +pi = (tildeQ')\(elast'); + +% availability = MTTF/MTBF +A = sum(pi(av_indx)); + +%% solution of the reliability model +hatQ = Qre(av_indx,av_indx); + +hatpi0 = pi0(av_indx); + +tau = -(hatQ')\(hatpi0'); + +MTTF = (tau')*ones(length(av_indx),1); + +%% evaluate MTBF +MTBF = MTTF/A; + +fprintf("Percentage of up states is %f\n", 100.0*length(av_indx)/nstates); + +end \ No newline at end of file diff --git a/ACload/evaluate_SR.m b/ACload/evaluate_SR.m new file mode 100644 index 0000000..e823fb6 --- /dev/null +++ b/ACload/evaluate_SR.m @@ -0,0 +1,35 @@ +function [A, MTTF, MTBF] = evaluate_SR(nm,n,lambdaB,lambdaR,muB,muR) + +%% model definition +% notice: here last state is the system failed state +[Q,pi0] = CTMC_SR(lambdaB,lambdaR,muB,muR); + +%% solution of the availability model +tildeQ = Q; +tildeQ(:,end) = ones(3,1); + +en = zeros(1,3); +en(end) = 1; + +% steady-state probability vector +pi = (tildeQ')\(en'); + +% availability = MTTF/MTBF +A = sum(pi(1:end-1)); + +%% solution of the reliability model +hatQ = Q(1:end-1,1:end-1); + +hatpi0 = pi0(1:end-1); + +tau = -(hatQ')\(hatpi0'); + +MTTF = (tau')*ones(2,1); + +%% evaluate MTBF +MTBF = MTTF/A; + +%% closed formula for Availability +% Asys = 1 - (1-muB/(muB+lambdaB))^2*(1-muR/(muR+lambdaR))^2; + +end \ No newline at end of file diff --git a/DCload/CTMC_CR.m b/DCload/CTMC_CR.m new file mode 100644 index 0000000..0dd7349 --- /dev/null +++ b/DCload/CTMC_CR.m @@ -0,0 +1,32 @@ +function [nstates, Qav, Qre, pi0] = CTMC_CR(nm,n,lambdaB,lambdaR,muB,muR) +%CTMC_CR infinitesimal generator matrix and initial probability vector +% of the Component Redundancy (CR) CMTC +% nm : number of hot standy (n+m>=n+1) +% n : number of required rectifiers or batteries to address load +% lambdaB : battery failure rate +% lambdaR : rectifier failure rate +% muB : battery recovery rate +% muR : rectifier recovery rate + +nstates = nm-n+2;% number of states only considering batteries + +% define batteries' behavior +Rre = diag(lambdaB*[nm:-1:n],1); +Rre(2:end-1,1) = muB; + +% % Plotting reliability model +% D= digraph(Rre); +% Plot = plot(D); +% highlight(Plot,nstates,'NodeColor','r'); +% highlight(Plot,'Edges',inedges(D,nstates)','EdgeColor','r'); + +Rav = Rre; +Rav(end,1) = muB; + +Qav = Rav - diag(Rav*ones(nstates,1)); +Qre = Rre - diag(Rre*ones(nstates,1)); + +pi0 = zeros(1,nstates); +pi0(1) = 1; + +end \ No newline at end of file diff --git a/DCload/CTMC_SR.m b/DCload/CTMC_SR.m new file mode 100644 index 0000000..495783e --- /dev/null +++ b/DCload/CTMC_SR.m @@ -0,0 +1,20 @@ +function [Q,pi0] = CTMC_SR(lambdaB,lambdaR,muB,muR) +%CTMC_SR infinitesimal generator matrix and initial probability vector +% of the System Redundancy (SR) CMTC +% lambdaB : battery failure rate +% lambdaR : rectifier failure rate +% muB : battery recovery rate +% muR : rectifier recovery rate + +R = zeros(3,3); +R(1,2) = 2*lambdaB; +R(2,1) = muB;% if there are two independent teams 2*muB; +R(2,3) = lambdaB; +R(3,2) = muB; + +Q = R - diag(R*ones(3,1)); + +pi0 = zeros(1,3); +pi0(1) = 1; + +end \ No newline at end of file diff --git a/DCload/compare.m b/DCload/compare.m new file mode 100644 index 0000000..c0a68e2 --- /dev/null +++ b/DCload/compare.m @@ -0,0 +1,71 @@ +%% parameters setting +% System Redundancy (SR) +lambdaB_SR = 11.76e-6;% hours^(-1) +lambdaR_SR = 16.6e-6; +muB_SR = 0.5;% hours^(-1) +muR_SR = muB_SR; + +% Component Redundancy (CR) +lambdaB_CR = 17*lambdaB_SR; +lambdaR_CR = 17*lambdaR_SR; +muB_CR = 0.5;% hours^(-1) +muR_CR = muR_SR; + +ModulePower = 10;% kW +minLoad = 40;% kW +loadStep = 20;% kW + +%% start evaluation +numberOfEvaluations = 6; + +MTTF_perc_impr = zeros(numberOfEvaluations,1); +MTBF_perc_impr = zeros(numberOfEvaluations,1); +MTBF_SR = zeros(numberOfEvaluations,1); +MTBF_CR = zeros(numberOfEvaluations,1); +Energy_perc_impr = zeros(numberOfEvaluations,1); +for i = 1 : numberOfEvaluations + Load = minLoad + i*loadStep;% kW + + n = ceil(Load/ModulePower); + nm = n + 3; + + PotTransformerMonolite=0.04*2*Load; + PotTransformerModular=0.04*Load; + + Energy_perc_impr(i) = 100; % estimate only considering transformers + + [A_SR,MTTF_SR,MTBF_SR_value] = evaluate_SR(nm,n,lambdaB_SR,lambdaR_SR,muB_SR,muR_SR); + [A_CR,MTTF_CR,MTBF_CR_value] = evaluate_CR(nm,n,lambdaB_CR,lambdaR_CR,muB_CR,muR_CR); + + MTBF_SR(i) = MTBF_SR_value; + MTBF_CR(i) = MTBF_CR_value; + + A_perc_impr = 100*(A_CR-A_SR)/A_SR; + MTTF_perc_impr(i) = 100*(MTTF_CR-MTTF_SR)/MTTF_SR; + MTBF_perc_impr(i) = 100*(MTBF_CR_value-MTBF_SR_value)/MTBF_SR_value; + +end + +%% Plot +x = minLoad : loadStep : (minLoad+(numberOfEvaluations-1)*loadStep); + +zeroline=zeros(1,length(x)); + +fs = 14; +f1=figure(1); +f1.Position = [50 50 650 650]; +title('Percentage of improvement','FontSize',fs-2); +hold on; +xlabel('DC load [kW]','FontWeight','bold','FontSize',fs); +plot(x,MTBF_perc_impr,'Color',[0,0,0],'LineWidth',2); +plot(x,zeroline,'Color',[0.4,0.4,0.4]); + +f2=figure(2); +f2.Position = [50 50 650 650]; +title('Mean Time Between two Failures (MTBF)','FontSize',fs); +hold on; +xlabel('DC load [kW]','FontWeight','bold','FontSize',fs); +ylabel('MTBF [hours]','FontWeight','bold','FontSize',fs); +plot(x,MTBF_SR,"--",'Color',[0.5,0,0.5],'LineWidth',2); +plot(x,MTBF_CR,"-.",'Color',[0,0.7,0],'LineWidth',2); +legend('SR','CR','FontSize',fs); \ No newline at end of file diff --git a/DCload/evaluate_CR.m b/DCload/evaluate_CR.m new file mode 100644 index 0000000..79c018a --- /dev/null +++ b/DCload/evaluate_CR.m @@ -0,0 +1,34 @@ +function [A,MTTF,MTBF] = evaluate_CR(nm,n,lambdaB,lambdaR,muB,muR) + +%% model definition +% notice: here last state is, by definition, one of the system failed states +[nstates,Qav,Qre,pi0] = CTMC_CR(nm,n,lambdaB,lambdaR,muB,muR); + +%% solution of the availability model +tildeQ = Qav; +tildeQ(:,end) = ones(nstates,1); + +elast = zeros(1,nstates); +elast(end) = 1; + +% here only the last is a system failure state + +% steady-state probability vector +pi = (tildeQ')\(elast'); + +% availability = MTTF/MTBF +A = sum(pi(1:end-1)); + +%% solution of the reliability model +hatQ = Qre(1:end-1,1:end-1); + +hatpi0 = pi0(1:end-1); + +tau = -(hatQ')\(hatpi0'); + +MTTF = (tau')*ones(nstates-1,1); + +%% evaluate MTBF +MTBF = MTTF/A; + +end \ No newline at end of file diff --git a/DCload/evaluate_SR.m b/DCload/evaluate_SR.m new file mode 100644 index 0000000..e823fb6 --- /dev/null +++ b/DCload/evaluate_SR.m @@ -0,0 +1,35 @@ +function [A, MTTF, MTBF] = evaluate_SR(nm,n,lambdaB,lambdaR,muB,muR) + +%% model definition +% notice: here last state is the system failed state +[Q,pi0] = CTMC_SR(lambdaB,lambdaR,muB,muR); + +%% solution of the availability model +tildeQ = Q; +tildeQ(:,end) = ones(3,1); + +en = zeros(1,3); +en(end) = 1; + +% steady-state probability vector +pi = (tildeQ')\(en'); + +% availability = MTTF/MTBF +A = sum(pi(1:end-1)); + +%% solution of the reliability model +hatQ = Q(1:end-1,1:end-1); + +hatpi0 = pi0(1:end-1); + +tau = -(hatQ')\(hatpi0'); + +MTTF = (tau')*ones(2,1); + +%% evaluate MTBF +MTBF = MTTF/A; + +%% closed formula for Availability +% Asys = 1 - (1-muB/(muB+lambdaB))^2*(1-muR/(muR+lambdaR))^2; + +end \ No newline at end of file