I can't run this code, please fix it

if true
  % This code start from line 4
end
N = 33; % number of components
lambda = [0.3979; 0.8209; 0.7666; 0.1206; 0.8577; ...
0.2586; 0.7049; 0.2429; 0.3521; 0.8118; ...
0.1850; 0.5135; 0.3920; 0.7194; 0.6776; ...
0.2065; 0.3197; 0.3226; 0.8566; 0.3348; ...
0.8950; 0.8833; 0.6400; 0.0412; 0.0518; ...
0.9458; 0.2257; 0.7303; 0.2191; 0.0101; ...
0.7205; 0.1289; 0.1327];
MTTR = [0.3546; 0.6970; 0.8490; 0.8724; 0.0411; ...
0.2098; 0.7382; 0.7379; 0.1978; 0.4534; ...
0.2299; 0.0704; 0.3979; 0.8555; 0.6809; ...
0.2954; 0.8536; 0.7195; 0.3405; 0.0495; ...
0.0174; 0.7846; 0.2554; 0.6597; 0.8496; ...
0.2965; 0.2238; 0.0066; 0.0684; 0.4306; ...
0.7953; 0.7759; 0.9673];
power = [0 55 0 16 0 0 0 5 34 0 0 4 40 7 45 15 146 82 ...
0 67 9 405 45 83 0 0 4 0 32 0 28 2 0 450];
users = ceil(power/2)-1;
%%failure history
duration = 1000; % years
interuption = 0;
outage_time = 0;
maxT = 0;
for i = 1:N
downT(i),upT(i)] = failure_history2(lambda(i),MTTR(i),duration);
cur = max(upT{i}(1:end));
if (maxT < cur)
maxT = cur;
end
interuption = interruption + length(downT{i})-1;
outage_time = outage_time + sum(upT{i}-downT{i});
% figure;
line(1:length(downT{i}),downT{i});
line(1:length(downT{i}),upT{i},'color','r','Marker','.');
end
% INDEXES calculation
average_interuption = interuption/duration;
average_outage_time = outage_time/duration;
customers_num = sum(users);
total_power = sum(power);
SAIFI = average_interuption*customers_num/customers_num;
SAIDI = average_outage_time*customers_num/customers_num;
CAIDI = SAIDI/SAIFI;
ASAI = (customers_num*8760 - average_outage_time*customers_num)/(customers_num*8760);
EUE = average_outage_time*total_power; % kW*hour
disp('-------------------------------');
disp('Step 1');
disp(['SAIFI ' num2str(SAIFI)]);
disp(['SAIDI ' num2str(SAIDI)]);
disp(['CAIDI ' num2str(CAIDI)]);
disp(['ASAI ' num2str(ASAI)]);
disp(['EUE ' num2str(EUE)]);

[Ameer Hamza]

You also need to explain the problem you are facing or the error given by MATLAB. Also, we don't have the function failure_history2() so we cannot run your code successfully. An obvious mistake which can be spotted is that you are missing a square bracket

[downT(i),upT(i)] = failure_history2(lambda(i),MTTR(i),duration);
^ This is missing in your code.

[Fajri Mardiansyah]

This code i'm copied from other research. Can you fix it?

[Ameer Hamza]

We don't know what does the function failure_history2() do. The code is probably fine. To run the code, you will need the function file.

[Fajri Mardiansyah]

thank you sir, :)

[Enio Soares]

Fajri Mardiansyah, I'm with the same problem. Did you resolve this problem? Did you implament some code or find other that help you?

[Walter Roberson]

Is this related to https://ww2.mathworks.cn/matlabcentral/answers/uploaded_files/233458/JPEE_2017080914092727(1).pdf ? The authors of that code do not appear to have posted the definition of failure_history2

[Haddijatou Touray]

Does anyone knwo where to fine the file 'failure_history2()' ?

[Walter Roberson]

Ask the authors https://www.researchgate.net/publication/319014373_Matlab_Code_to_Assess_the_Reliability_of_the_Smart_Power_Distribution_System_Using_Monte_Carlo_Simulation

[Joao Cardoso]

does anyone know how other similar bars?

[Image Analyst]

@Joao Cardoso, not sure what that means. How "other similar bars" do what? What bars? Do you have the same code as @Fajri Mardiansyah??? How/why? Where did you get it and why are you using it?

[Joao Cardoso]

Hello the author uses lines 834 and 832. How could I simulate in other lines?

[Walter Roberson]

@Joao Cardoso

I do not see any reference to 832 or 834 in the code? And there posted code is less than 100 lines long.

I do not think I understand what you are asking?

[DGM]

FWIW, this is from here:

https://www.scirp.org/journal/PaperInformation.aspx?PaperID=78287

I edited it to fix the parts where the website formatting broke things. It still won't run without the missing function anyway.

"lines" isn't code line numbers, but electrical distribution lines; similarly, "bars" is bus bars, or in the context of variable naming, switches.

As far as how to rearrange this to simulate at other nodes, I don't know. I'd have to wrap my head around the network model to even know which lines/elements are which, and I'd have to T&D coursework and study this paper for a while.

clear;

N = 33; % number of components

lambda = [0.3979; 0.8209; 0.7666; 0.1206; 0.8577; ...
    0.2586; 0.7049; 0.2429; 0.3521; 0.8118; ...
    0.1850; 0.5135; 0.3920; 0.7194; 0.6776; ...
    0.2065; 0.3197; 0.3226; 0.8566; 0.3348; ...
    0.8950; 0.8833; 0.6400; 0.0412; 0.0518; ...
    0.9458; 0.2257; 0.7303; 0.2191; 0.0101; ...
    0.7205; 0.1289; 0.1327];

MTTR = [0.3546; 0.6970; 0.8490; 0.8724; 0.0411; ...
    0.2098; 0.7382; 0.7379; 0.1978; 0.4534; ...
    0.2299; 0.0704; 0.3979; 0.8555; 0.6809; ...
    0.2954; 0.8536; 0.7195; 0.3405; 0.0495; ...
    0.0174; 0.7846; 0.2554; 0.6597; 0.8496; ...
    0.2965; 0.2238; 0.0066; 0.0684; 0.4306; ...
    0.7953; 0.7759; 0.9673];

power = [0 55 0 16 0 0 0 5 34 0 0 4 40 7 45 15 146 82 0 67 9 405 45 83 0 0 4 0 32 0 28 2 0 450];
users = ceil(power/2)-1;

%% failure history

duration = 1000; % years
interuption = 0;
outage_time = 0;
maxT = 0;

for i = 1:N
    
    [downT{i},upT{i}] = failure_history2(lambda(i),MTTR(i),duration);
    cur = max(upT{i}(1:end));
    
    if (maxT < cur)
        maxT = cur;
    end
    
    interuption = interruption + length(downT{i})-1;
    outage_time = outage_time + sum(upT{i}-downT{i});
    
    % figure;
    line(1:length(downT{i}),downT{i});
    line(1:length(downT{i}),upT{i},'color','r','Marker','.');
end

% INDEXES calculation
average_interuption = interuption/duration;
average_outage_time = outage_time/duration;
customers_num = sum(users);
total_power = sum(power);

SAIFI = average_interuption*customers_num/customers_num;
SAIDI = average_outage_time*customers_num/customers_num;
CAIDI = SAIDI/SAIFI;
ASAI = (customers_num*8760 - average_outage_time*customers_num)/ (customers_num*8760);
EUE = average_outage_time*total_power; % kW*hour

disp('-------------------------------');
disp('Step 1');
disp(['SAIFI ' num2str(SAIFI)]);
disp(['SAIDI ' num2str(SAIDI)]);
disp(['CAIDI ' num2str(CAIDI)]);
disp(['ASAI ' num2str(ASAI)]);
disp(['EUE ' num2str(EUE)]);

%% Step2 switches

switch_834_lines = [20 21 22 31 32];
switch_832_lines = [17 18 19 30 23 24 25 switch_834_lines];
switch_834_elements = [19 20 21 30 31];
switch_832_elements = [switch_834_elements 18 22 23 24 25 29 32];

interaptionXcustomer_number = 0;
interaptionXcustomer_duration = 0;
unservedEnergy = 0;
cur_t = 0;

while cur_t < maxT
    
    % find next failure time
    min_t = (duration + 1)*24*365;
    fail_line_number = 0;
    
    for i = 1:N
        
        clear indx;
        indx = find(downT{i} > cur_t);
        
        if (~isempty(indx))
            
            indx = indx(1);
            
            if (downT{i}(indx) < min_t)
                min_t = downT{i}(indx);
                fail_line_number = i;
                n = indx;
                cur_outage_time = upT{i}(indx) - downT{i}(indx);
            end
        end
    end
    
    cur_t = upT{fail_line_number}(n);
    in_switch_834 = ~isempty(find(switch_834_lines == fail_line_number));
    customer_834 = sum(users(switch_834_elements));
    power_834 = sum(power(switch_834_elements));
    in_switch_832 = ~isempty(find(switch_832_elements == fail_line_number));
    customer_832 = sum(users(switch_832_elements));
    power_832 = sum(power(switch_832_elements));
    
    % interaption calculation
    if (~in_switch_832)
        interaptionXcustomer_number = interaptionXcustomer_number + customers_num;
        interaptionXcustomer_duration = interaptionXcustomer_duration + customers_num*cur_outage_time;
        unservedEnergy = unservedEnergy + total_power*cur_outage_time;
    else
        if (in_switch_834)
            interaptionXcustomer_number = interaptionXcustomer_number + customer_834;
            interaptionXcustomer_duration = interaptionXcustomer_duration + customer_834*cur_outage_time;
            unservedEnergy = unservedEnergy + power_834*cur_outage_time;
        else
            interaptionXcustomer_number = interaptionXcustomer_number + customer_832;
            interaptionXcustomer_duration = interaptionXcustomer_duration + customer_832*cur_outage_time;
            unservedEnergy = unservedEnergy + power_832*cur_outage_time;
        end
    end
end

% INDEXES calculation
SAIFI = interaptionXcustomer_number/customers_num/duration;
SAIDI = interaptionXcustomer_duration/customers_num/duration;
CAIDI = SAIDI/SAIFI;
ASAI = (customers_num*8760 - interaptionXcustomer_duration/duration)/ (customers_num*8760);
EUE = unservedEnergy/duration; % kW*hour

disp('-------------------------------');
disp('Step 2');
disp(['SAIFI ' num2str(SAIFI)]);
disp(['SAIDI ' num2str(SAIDI)]);
disp(['CAIDI ' num2str(CAIDI)]);
disp(['ASAI ' num2str(ASAI)]);
disp(['EUE ' num2str(EUE)]);

%% Step 3 switches from Step 2 and DGs

switch_834_lines = [20 21 22 31 32];
switch_832_lines = [17 18 19 30 23 24 25 switch_834_lines];
switch_834_elements = [19 20 21 30 31];
switch_832_elements = [switch_834_elements 18 22 23 24 25 29 32];

interaptionXcustomer_number = 0;
interaptionXcustomer_duration = 0;
unservedEnergy = 0;
cur_t = 0;

while cur_t < maxT
    
    % find next failure time
    min_t = (duration + 1)*24*365;
    fail_line_number = 0;
    
    for i = 1:N
        
        clear indx;
        indx = find(downT{i} > cur_t);
        
        if (~isempty(indx))
            
            indx = indx(1);
            
            if (downT{i}(indx) < min_t)
                min_t = downT{i}(indx);
                fail_line_number = i;
                n = indx;
                cur_outage_time = upT{i}(indx) - downT{i}(indx);
            end
            
        end
    end
    
    cur_t = upT{fail_line_number}(n);
    in_switch_834 = ~isempty(find(switch_834_lines == fail_line_number));
    customer_834 = sum(users(switch_834_elements));
    power_834 = sum(power(switch_834_elements));
    in_switch_832 = ~isempty(find(switch_832_elements == fail_line_number));
    customer_832 = sum(users(switch_832_elements));
    power_832 = sum(power(switch_832_elements));
    customer_800 = customers_num-customer_832;
    power_800 = total_power-power_832;
    
    % interaption calculation
    if (~in_switch_832)
        interaptionXcustomer_number = interaptionXcustomer_number + customer_800;
        interaptionXcustomer_duration = interaptionXcustomer_duration + customer_800*cur_outage_time;
        unservedEnergy = unservedEnergy + power_832*cur_outage_time;
    else
        if (in_switch_834)
            interaptionXcustomer_number = interaptionXcustomer_number + customer_834;
            interaptionXcustomer_duration = interaptionXcustomer_duration + customer_834*cur_outage_time;
            unservedEnergy = unservedEnergy + power_834*cur_outage_time;
        else
            interaptionXcustomer_number = interaptionXcustomer_number + customer_832;
            interaptionXcustomer_duration = interaptionXcustomer_duration + customer_832*cur_outage_time;
            unservedEnergy = unservedEnergy + power_832*cur_outage_time;
        end
    end
end

% INDEXES calculation
SAIFI = interaptionXcustomer_number/customers_num/duration;
SAIDI = interaptionXcustomer_duration/customers_num/duration;
CAIDI = SAIDI/SAIFI;
ASAI = (customers_num*8760 - interaptionXcustomer_duration/duration)/ (customers_num*8760);
EUE = unservedEnergy/duration; % kW*hour

disp('-------------------------------');
disp('Step 3');
disp(['SAIFI ' num2str(SAIFI)]);
disp(['SAIDI ' num2str(SAIDI)]);
disp(['CAIDI ' num2str(CAIDI)]);
disp(['ASAI ' num2str(ASAI)]);
disp(['EUE ' num2str(EUE)]);

[Hung Nguyen]

@DGM Did you resolve this problem and have function file ?

[Walter Roberson]

The code that was published in the paper is only a subset of the actual code. Unfortunately the actual code does not appear to be available. You could try writing to the authors to see if they are willing to make it available.

[DGM]

I do not have any more information other than what I posted. As Walter says, it's incomplete. All I did was clean up the formatting.

[myrto pieridou]

hi.. do you find the file?

[Walter Roberson]

No, failure_history2 has not been located.

[DGM]

For sake of the future: I have found nothing beyond what I posted. This is not a problem I would pursue out of personal curiosity, so it will remain that way unless someone else responds with substantial contributions of their own. At this point, the most likely route to those answers is to contact the authors of the original paper. If anyone does find the missing code and get it working, you're free to post an answer and get the appropriate credit for it.