Код:
% "Моделювання надійності з"єднання"
clear all; %очистити робочий простір
tic; %час початку обчислень
%% Вихідні дані
n_v=11; %кількість вузлів мережі
p_i=ones(n_v,1); %масив імовірності відмови кожного вузла
p_i=p_i*5e-3; %вважаємо, що всі імовірності відмови однакові
q_i=1-p_i; %масив імовірності роботи кожного вузла
%% Аналітичний метод розрахунку
Q(1)=1-q_i(3,1)*q_i(6,1)*q_i(9,1); %шлях 1: вузли 3-6-9
Q(2)=1-q_i(5,1)*q_i(6,1)*q_i(9,1); %шлях 2: вузли 5-6-9
Q(3)=1-q_i(5,1)*q_i(8,1)*q_i(9,1); %шлях 3: вузли 5-8-9
Q(4)=1-q_i(5,1)*q_i(8,1)*q_i(10,1); %шлях 4: вузли 5-8-10
Q_a=prod(Q); %для паралельних шляхів - добуток імовірностей відмови (всі шляхи незалежні)
Q_eg=Q(1)*Q(4); %для дійсно незалежних шляхів (без врахування ребер 5-6 та 8-9)
%% Статистичне моделювання
X=zeros(n_v,1); %масив станів вузлів
N=round(100/Q_eg); %кількість дослідів для однієї оцінки імовірності відмови
No=1300; %і кількість цих оцінок
sh=10; %кількість діапазонів гістограми
Q_sm=[]; %масив для накопичення оцінок імовірностей відмови
for k=1:No
Q_s=0;
for count = 1 : N %цикл статистичного моделювання
r_n=rand(n_v,1); %генерація випадкових величин з рівномірним розподілом
X=X*0; %обнулення вектора Х
for (i=1:n_v) %Цикл генерації стану вузлів мережі
if( r_n(i,1) <= p_i(i,1)) %Якщо значення генератора менше рівне імовірності відмови
X(i,1)=1; %тоді даний вузол відмовив !
end
end
if (X(3,1)+X(6,1)+X(9,1) == 0) %Перевірка шляхів 1-4
continue;
end
if (X(5,1)+X(6,1)+X(9,1) == 0) %якщо хоча б один шлях не вийшов з ладу
continue;
end
if (X(5,1)+X(8,1)+X(9,1) == 0) %то система працює
continue;
end
if (X(5,1)+X(8,1)+X(10,1) == 0)
continue;
end
Q_s=Q_s+1; %інакше - реєструємо відмову
end
Q_sm=[Q_sm; Q_s/N]; %запис чергової оцінки відмови в масив накопичення
end
%% Обчислення статистичних параметрів (No оцінок Q) - статистичне моделювання
My_s=mean(Q_sm); % математичне очікування
Sy_s=std (Q_sm); % середньоквадратичне відхилення
Dy_s=Sy_s*Sy_s; % дисперсія
v_s=Sy_s/My_s; % коефіцієнт варіації (точність оцінки)
sig3_s=Sy_s*3; % Імовірність незв"язності вузлів і та j
[Ry_s, ix_s]=hist(Q_sm, sh); %Гістограма і середини інтервалів (статистичне моделювання)
ts=toc;
tic;
%% Метод зважування
koef=0.5/p_i(1,1); %підбираємо ваговий коефіцієнт так, щоб імовірність відмови вузла була 0.5
p_i_strix=koef*p_i; %і знаходимо нову імовірність відмови кожного з вузлів, перемноживши всі імовірності на ваговий коефіцієнт
q_i_strix=1-p_i_strix; %нові імовірності роботи вузлів
Q_tmp=(1-q_i_strix(1,1)^3)^2; %оціночна імовірність відмови для обчислення кількості дослідів
N_z=round(100/Q_tmp); %обчислення кількості дослідів для методу зважування
Q_sz=[]; %масив для накопичення оцінок (метод зважування)
% Розрахунок вагового коефіцієнта
PX=p_i(3,1)*p_i(6,1)*p_i(9,1); %розрахунок проведемо для вектора Х, що рівний
q_i(3,1)=1;q_i(6,1)=1;q_i(9,1)=1; % Х = { 0 0 1 0 0 1 0 0 1 0 0 }
tmp=prod(q_i); %і добутків імовірностей
PX=PX*tmp; % q1*q2* p3 *q4*q5* p6 *q7*q8* p9 *q10*q11
PsX=p_i_strix(3,1)*p_i_strix(6,1)*p_i_strix(9,1);
q_i_strix(3,1)=1;q_i_strix(6,1)=1;q_i_strix(9,1)=1;
tmp=prod(q_i_strix);
PsX=PsX*tmp;
w_koef=PX/PsX; %Ваговий коефіцієнт
for k=1:No
Q_z=0;
for count = 1 : N_z %цикл статистичного моделювання
r_n=rand(n_v,1); %генерація випадкових величин з рівномірним розподілом
X=X*0; %обнулення вектора Х
for (i=1:n_v) %Цикл генерації стану вузлів мережі
if( r_n(i,1) <= p_i_strix(i,1)) %Якщо значення генератора менше рівне імовірності відмови
X(i,1)=1; %тоді даний вузол відмовив !
end
end
if (X(3,1)+X(6,1)+X(9,1) == 0) %Перевірка шляхів 1-4
continue;
end
if (X(5,1)+X(6,1)+X(9,1) == 0) %якщо хоча б один шлях не вийшов з ладу
continue;
end
if (X(5,1)+X(8,1)+X(9,1) == 0) %то система працює
continue;
end
if (X(5,1)+X(8,1)+X(10,1) == 0)
continue;
end
Q_z=Q_z+1; %інакше - реєструємо відмову
end
Q_sz=[Q_sz; Q_z*w_koef]; %Імовірність відмови шляхом статистичного моделювання
end
Q_sz=Q_sz/N_z;
%% Обчислення статистичних параметрів (No оцінок Q) - зважування
My_z=mean(Q_sz); % математичне очікування
Sy_z=std (Q_sz); % середньоквадратичне відхилення
Dy_z=Sy_z*Sy_z; % дисперсія
v_z=Sy_z/My_z; % коефіцієнт варіації (точність оцінки)
sig3_z=Sy_z*3; % Імовірність незв"язності вузлів і та j
[Ry_z, ix_z]=hist(Q_sz, sh); %Гістограма і середини інтервалів (статистичне моделювання)
%% Побудова гістограми
f=figure;
plot(Ry_s, '*r');
hold on;
plot(Ry_s, 'r');
plot(Ry_z, '^b');
plot(Ry_z, 'b');
xlabel('Номер інтервалу');
ylabel('Частота попадання в інтервал');
title('Дослідження імовірності відмови \bfQ');
tz=toc; %час закінчення обчислень
%% Форматування результатів
Rez(1,1)=N; Rez(1,2)=N_z;
Rez(2,1)=N*No; Rez(2,2)=N_z*No;
Rez(3,1)=ts; Rez(3,2)=tz;
Rez(4,1)=My_s; Rez(4,2)=My_z;
Rez(5,1)=Sy_s; Rez(5,2)=Sy_z;
Rez(6,1)=Dy_s; Rez(6,2)=Dy_z;
Rez(7,1)=v_s; Rez(7,2)=v_z;
Rez(8,1)=3*Sy_s; Rez(8,2)=3*Sy_z;
%% Кінець