RBTS_F4序贯.zip_RBTS_可靠性_可靠性 RBTS_序贯_序贯RBTS可靠性

function [saifi,saidi,caidi,asui,aens]=pdfd( dataline,dataload) %UNTITLED Summary of this function goes here % Detailed explanation goes here % change the input mttf=1./dataline(:,4)*8760;%hours mttr=dataline(:,5); state=0; maxhour=3000*8760;%3000year i=1; num_line=size(dataline,1); component_up=[]; component_down=[]; % get the componet state while 1 component_up=[component_up,exprnd(mttf)]; component_down=[component_down,exprnd(mttr)]; least_sum_time=min(sum(component_up+component_down,2)); if least_sum_time>maxhour break; end i=i+1; end % get the intial system state(not consider the operator) component_history=zeros(size(component_up,1),size(component_up,2)*2); component_history(:,1:2:end)=component_up; component_history(:,2:2:end)=component_down; sys_state=combinestate(component_history); % change the sys_state (considering the operator) sys_state=change(sys_state); % get the load_point_operating_histories lpoh=get_lpoh(sys_state,dataline); %% get the three reliability indices of each load point % get the distribution of outage rate and unavailability in each load point p_outagerate=[]; p_unavailability=[]; lpoh_num=size(lpoh,2); first=1; for i=1:lpoh_num last=i; while 1 state_time=sum(lpoh(end,first:last)); if state_time>8760 temp_outagerate=sum(lpoh(1:(end-1),first:last)<1e-3,2); temp_unavailability=(lpoh(1:(end-1),first:last)<1e-3)*(lpoh(end,first:last))'; temp_unavailability=temp_unavailability-(lpoh(1:(end-1),last)<1e-3)*(state_time-8760); p_unavailability=[p_unavailability,temp_unavailability]; p_outagerate=[p_outagerate,temp_outagerate]; lpoh(end,last)=state_time-8760; first=last; end if lpoh(end,last)<8760 break; end end end % discrete the unavailibility interval=0.5; longest_unav=max(max(p_unavailability)); longest_unav=round(longest_unav)+1; for i=0:interval:longest_unav idx=find(abs(p_unavailability-i)<interval/2); p_unavailability(idx)=i; end %% get the distribution of reliability indicies of the system p_saifi=dataload(:,3)'*p_outagerate/sum(dataload(:,3)); p_saidi=dataload(:,3)'*p_unavailability/sum(dataload(:,3)); p_caidi=p_saidi./p_saifi; p_asui=p_saidi/8760; p_aens=dataload(:,2)'*p_unavailability/sum(dataload(:,3));%(MWh/customer/yr) %average value saifi=mean(p_saifi); saidi=mean(p_saidi); caidi=mean(p_caidi); asui=mean(p_asui); aens=mean(p_aens); % plot these indicies % figure(1); % histogram(p_saifi); figure(1); ksdensity(p_saifi); figure(2); ksdensity(p_saidi); figure(3); ksdensity(p_caidi); figure(4); ksdensity(p_asui); figure(5); ksdensity(p_aens); end function [lpoh]=get_lpoh(sys_state,dataline) lpoh=[]; sys_state_num=size(sys_state,2); for i=1:sys_state_num adjacency_matrix=sparse(dataline(:,1),dataline(:,2),sys_state(1:(end-1),i),54,54)+sparse(dataline(:,2),dataline(:,1),sys_state(1:(end-1),i),54,54); c=checkc(adjacency_matrix,1); elapesed_time=sys_state(end,i); temp=[]; for j=32:54 if c(j)==1 temp= [temp;1]; else temp=[temp;0]; end end temp_lpoh=[temp;elapesed_time]; lpoh=[lpoh,temp_lpoh]; end end function [ sys_state ] = combinestate( m ) %UNTITLED Summary of this function goes here initialstate=ones(size(m,1),1); sys_state=[]; % Detailed explanation goes here leasttime=min(sum(m,2)); minustime=0; while (leasttime-minustime)>1e-3% tolerance [value,where]=min(m(:,1)); temp=[initialstate;value]; sys_state=[sys_state,temp]; %minus the elapsed time m(:,1)=m(:,1)-value; %reorder the corresponding element m(where,1:end-1)=m(where,2:end); m(where,end)=0; initialstate(where)=~initialstate(where); minustime=minustime+value; end end function [sys_state]= change(sys_state)%%different topologe ,different content add_state=[]; total_state_num=size(sys_state,2); for i=1:total_state_num fault_line=find(sys_state(:,i)<1e-3); if any(fault_line==1)||any(fault_line==3)||any(fault_line==5)||any(fault_line==7)||any(fault_line==9)||any(fault_line==11)||any(fault_line==14)||any(fault_line==17)||any(fault_line==29) sys_state(1,i)=0;% open F1 else if any(fault_line==20)||any(fault_line==23)||any(fault_line==25)||any(fault_line==27) sys_state(18,i)=0;% open F2 end if any(fault_line==50)||any(fault_line==52) sys_state(48,i)=0;%open F4 end if any(fault_line==39)||any(fault_line==42)||any(fault_line==44)||any(fault_line==46) sys_state(37,i)=0;%open F3 end if any(fault_line==31)||any(fault_line==34)||any(fault_line==36) if sys_state(end,i)>1 sys_state(1,i)=0;%open F1 for 1 hour %need add one sys_state temp=sys_state(:,i); temp(1)=1; temp(29)=0; temp(end)=temp(end)-1; temp_addstate=[i;temp]; add_state=[add_state,temp_addstate]; else sys_state(1,i)=0;%open F1 end end end end %add sys_state if needed if ~isempty(add_state) add_num=size(add_state,2); for j=1:add_num k=add_state(1,j); sys_state(:,end+1)=ones(size(sys_state,1),1); sys_state(:,(k+2):end)=sys_state(:,(k+1):(end-1)); sys_state(:,k+1)=add_state(2:end,j); end end end function [ c] = checkc( am,k,tol ) %check the connectivity of the nondirection_graph. %check the connectivity between one point and other points. %Input: % am: % n x n Adjacency Matrix % k: % the number of the point % tol(optional): % tolerance %Output: % c: % n x 1 the connectivity between the point k and other points % (the element is 1 ,if connected;otherwise 0.) % Copyright (c) 2016 by Chongqing University Electrical Engineering College(CUEEC) % by Tang Yi % initializing parameters if nargin <3 tol=1e-6; if nargin <2 k=1; end end n=size(am,2); i=zeros(n,1); i(k)=1e6; % reconstruct the adjacency matrix change_am=am*-1+sparse(1:n,1:n,(sum(am,2))')+eye(n); % calculate the voltage v=change_am\i; % check the connectivity c=v>tol; % if the input is just one parameter ,then just check the whole graph if nargin==1 if sum(c)==n c=1; else c=0; end end end