Markov模型SIL验算学习

clc;
clear;
format long;

% 定义初始值
lambda_sd = 1.93E-07;
lambda_su = 0;
lambda_dd = 1.0E-06;  % 检测到的危险失效故障率
lambda_du = 2.73E-08; % 未检测到的危险失效故障率
E = 0.95;                   % 功能测试覆盖率
MTTR = 24;             % 故障修复时间 (小时)
TI = 8760;                % 检验测试时间间隔（小时）
LT = 20;               % 系统运行时间
SD = 24;                 % 系统启动时间



lambda_s = lambda_sd + lambda_su;  
lambda_d = lambda_du + lambda_dd; 
DC = lambda_dd / lambda_d;  
SFF = (lambda_s+lambda_dd)/(lambda_s+lambda_d);

Mu_sd = 1 / SD;        
Mu_0 = 1 / MTTR;        


LT_hour = LT*365*24;   


S0(1,:) = [1 0 0 0 0];  

% 定义转移矩阵
p = [1 - (lambda_s + lambda_dd + E * lambda_du + (1 - E) * lambda_du), lambda_s, lambda_dd, E * lambda_du, (1 - E) * lambda_du;
     Mu_sd, 1 - Mu_sd, 0, 0, 0;
     Mu_0, 0, 1 - Mu_0, 0, 0;
     0, 0, 0, 1, 0;
     0, 0, 0, 0, 1];


PFD = zeros(LT_hour, 1);  % 预分配PFD数组

for i = 1:(LT_hour/TI)  
    S(1,:) = S0;  % 定义初始状态
    for j = 1:TI  % 计算TI小时内每个小时的系统状态
        S(j+1,:) = S(j,:) * p;
    end
    for j = 1:TI
        PFD((i-1)*TI + j) = S(j,3) + S(j,4) + S(j,5);  % 累加 PFD 值
    end
    S0 = [(1 - S(j,5)) 0 0 0 S(j,5)];
end

PFDavg = sum(PFD) / LT_hour;  % LT年内每小时PFD值的平均值
figure;
plot(1:LT_hour, PFD);
xlabel('年');
ylabel('PFD');
title('每年的PFD值');
set(gca, 'YScale', 'log');  % 设置Y轴为对数刻度
grid on;
% 添加PFDavg横线
yline(PFDavg, '--r', 'LineWidth', 2, 'Label', ['PFDavg = ', num2str(PFDavg, '%.2e')]);

disp(['SFF为: ', num2str(SFF)]);
disp(['诊断覆盖率为: ', num2str(DC)]);
disp(['PFD值: ', num2str(PFDavg)]);