apacsys-vaporrecoverymonitor-document/程序设计文档/油气回收智能诊断系统/guzhangzhenduan_26.4.2.m

clear; clc; close all;

% V、P分别代表回气流速和回气压力前20个采样点信号数据
P = [0.2	0	0.6	2.4	14	14.5	13.6	3.9	13.4	9.5	7.2	6.9	14.1	2.1	9.5	8.2	12.6	11.9	19.2	4.2]; % 压力信号
V = [0.2	0	0.6	2.4	14	14.5	13.6	3.9	13.4	9.5	7.2	6.9	14.1	2.1	9.5	8.2	12.6	11.9	19.2	4.2];% 流速信号
%% 特征值
% 1. 平均值 (Mean)
mean_value_v = mean(V );
mean_value_p = mean(P);
% 2. 方差 (Variance)
variance_value_v = var(V );
variance_value_p = var(P);
% 3. 标准差 (Standard Deviation)
std_value_v = std(V );
std_value_p = std(P);
% 4. 最大值 (Maximum)
max_value_v = max(V (1:10));
max_value_p = max(P (1:10));
% 5. 最小值 (Minimum)
min_value_v = min(V );
min_value_p = min(P );
% 定义实时工况特征数据
% Y0 = [115.70,2832.562, 53.221, 226,126.9]; 
Y01 = [mean_value_p, variance_value_p,std_value_p,max_value_p,min_value_p]; %压力传感
Y02 = [mean_value_v, variance_value_v,std_value_v,max_value_v,min_value_v]; %流速传感


%%回气压力历史数据，回气流速也一样
P1 = [0.2	0	0.6	2.4	14	14.5	13.6	3.9	13.4	9.5	7.2	6.9	14.1	2.1	9.5	8.2	12.6	11.9	19.2	4.2;
    0.2	0	0.6	2.4	14	14.5	13.6	3.9	13.4	9.5	7.2	6.9	14.1	2.1	9.5	8.2	12.6	11.9	19.2	4.2;
    0.2	0	0.6	2.4	14	14.5	13.6	3.9	13.4	9.5	7.2	6.9	14.1	2.1	9.5	8.2	12.6	11.9	19.2	4.2;
    0.2	0	0.6	2.4	14	14.5	13.6	3.9	13.4	9.5	7.2	6.9	14.1	2.1	9.5	8.2	12.6	11.9	19.2	4.2;
    0.2	0	0.6	2.4	14	14.5	13.6	3.9	13.4	9.5	7.2	6.9	14.1	2.1	9.5	8.2	12.6	11.9	19.2	4.2;
    0.2	0	0.6	2.4	14	14.5	13.6	3.9	13.4	9.5	7.2	6.9	14.1	2.1	9.5	8.2	12.6	11.9	19.2	4.2;
    0.2	0	0.6	2.4	14	14.5	13.6	3.9	13.4	9.5	7.2	6.9	14.1	2.1	9.5	8.2	12.6	11.9	19.2	4.2;
    0.2	0	0.6	2.4	14	14.5	13.6	3.9	13.4	9.5	7.2	6.9	14.1	2.1	9.5	8.2	12.6	11.9	19.2	4.2;
    0.2	0	0.6	2.4	14	14.5	13.6	3.9	13.4	9.5	7.2	6.9	14.1	2.1	9.5	8.2	12.6	11.9	19.2	4.2;
    0.2	0	0.6	2.4	14	14.5	13.6	3.9	13.4	9.5	7.2	6.9	14.1	2.1	9.5	8.2	12.6	11.9	19.2	4.2]
% 获取矩阵的行数和列数
[n, m] = size(P1);
%  计算每一行的平均值
mean_mean_p = zeros(n, 1); 
mean_variance_p=zeros(n, 1); 
mean_std_p=zeros(n, 1); 
mean_max_p=zeros(n, 1); 
mean_min_p=zeros(n, 1); 



for i = 1:n
    % 对第 i 行的所有列求和，再除以元素个数
    mean_mean_p(i) = mean(P1(i, :));
    mean_variance_p(i) = var(P1(i, :));
    mean_std_p(i) = std(P1(i, :));
    mean_max_p(i) =  max(P1(i, :));
    mean_min_p(i) = min(P1(i, :));
end
% 3. 计算10行的总平均值
 total_mean_p = mean(mean_mean_p);
    total_variance_p = mean(mean_variance_p);
    total_std_p = mean(mean_std_p);
    total_max_p =  mean(mean_max_p);
    total_min_p = mean(mean_min_p);
%%
% 定义比较序列（不同故障数据，每行为一种故障）压力传感
Y1 = [ total_mean_p, total_variance_p, total_std_p,total_max_p,total_min_p;    % 正常1   此处的数值就需要历史加油过程中单笔交易中回气流速和回气压力前20个采样点信号数据，也是求5个特征值
    17.95, 248.8711,15.775,1,1;    % 轻微泄露2
     23.3900,338.2757, 18.3923,1,1;    % 一般泄露3
     12.5850, 130.3287, 11.4162,1,1];   % 严重泄露4
 % 定义比较序列（不同故障数据，每行为一种故障）流速传感
Y2 = [ 36.5250, 365.1199, 19.1081,1,1;    % 正常1
    17.95, 248.8711,15.775,1,1;    % 轻微泄露2
     23.3900,338.2757, 18.3923,1,1;    % 一般泄露3
     12.5850, 130.3287, 11.4162,1,1];   % 严重泄露4

% 设置分辨系数
rho = 0.5;

% 执行灰色关联分析
[degree_1, order] = grey_relation_analysis1(Y01, Y1, rho);%压力传感-灰色关联度结果
[degree_2, order] = grey_relation_analysis1(Y02, Y2, rho);%流速传感-灰色关联度结果

%% 证据理论（Dempster-Shafer Theory）故障诊断
% 辨识框架：Theta = {F1, F2, F3}，分别表示三种故障类型
% 焦元：各子集的基本概率分配（BPA）

% 分别获取气路压力和流速传感器的BPA
bpa1=degree_1;%压力传感
bpa2=degree_2;%流速传感

% 打印原始BPA
disp('压力传感器 BPA:');
disp(bpa1);
disp('流速传感器 BPA:');
disp(bpa2);

% 执行证据组合
combined = dempster_combine(bpa1, bpa2);

% 显示融合结果
disp('融合后BPA:');
disp(combined);

% 故障决策（选择最高置信度）
[~, idx] = max([combined.F1, combined.F2, combined.F3, combined.F4]);
faults = {'F1', 'F2', 'F3','F4'};
fprintf('\n最终诊断结果：%s故障（置信度%.2f%%）\n', ...
    faults{idx}, 100*(combined.(faults{idx})));

%%

function [grey_relation_degree, order] = grey_relation_analysis1(Y0, Y, rho)
% 灰色关联分析用于故障诊断
% 输入：
%   Y0 - 参考序列（行或列向量）
%   Y  - 比较序列矩阵，每行代表一个比较序列
%   rho - 分辨系数（可选，默认0.5）
% 输出：
%   grey_relation_degree - 各比较序列的关联度
%   order - 关联度排序（从高到低）

% 参数检查与默认值设置
if nargin < 3
    rho = 0.5; % 默认分辨系数
end

% 确保Y0是行向量
Y0 = Y0(:)';

% 获取比较序列参数
[n, m] = size(Y);
m0 = length(Y0);

% 检查序列长度一致性
if m ~= m0
    error('参考序列和比较序列的长度必须一致');
end

% 数据预处理：初值化（每元素除以序列第一个元素）
% Y0_normalized = Y0 / Y0(1);
% Y_normalized = Y ./ Y(:,1); % 按行归一化
%%
% 数据预处理：均值化（每元素除以序列平均值）
Y0_mean = mean(Y0); % 计算参考序列均值
Y0_normalized = Y0 / Y0_mean; % 参考序列均值化

Y_mean = mean(Y, 2); % 计算各比较序列的均值（按行计算）
Y_normalized = Y ./ Y_mean; % 比较序列均值化（自动广播）

%%

% 计算绝对差值矩阵
diff = abs(Y0_normalized - Y_normalized);

% 计算全局最小差和最大差
delta_min = min(diff(:));
delta_max = max(diff(:));

% 计算关联系数矩阵
relations = (delta_min + rho * delta_max) ./ (diff + rho * delta_max);

% 计算关联度（按行求平均）
grey_relation_degree = mean(relations, 2);

% 按关联度降序排序
[~, order] = sort(grey_relation_degree, 'descend');
end
%% 步骤2：Dempster组合规则实现
function combined_bpa = dempster_combine(bpa1, bpa2)
    % 定义所有可能的焦元组合
    sets = {{'F1'}, {'F2'}, {'F3'},{'F4'}};
    
    % 计算冲突系数K
    K = 0;
    for i = 1:length(sets)
        for j = 1:length(sets)
            if isempty(intersect(sets{i}, sets{j}))%C = intersect(A,B)，其中A和B是要找出交集的两个数组   isempty是MATLAB中的一个函数，用于判断一个变量是否为空。
                K = K + bpa1(i)* bpa2(j);
            end
        end
    end
    
    % 计算合并后的BPA
    combined_bpa = struct();
    for k = 1:length(sets)
        sum_val = 0;
        for i = 1:length(sets)
            for j = 1:length(sets)
                if isequal(union(sets{i}, sets{j}), sets{k})%如果你想计算两个或多个向量的并集，你可以使用union函数。这个函数会合并这些向量，并去除重复的元素
                    sum_val = sum_val + bpa1(i) * bpa2(j);
                end
            end
        end
        combined_bpa.(strjoin(sets{k},'')) = sum_val / (1 - K);
    end
end

% 辅助函数：获取BPA值
function val = get_value(bpa, set)
    if length(set) == 1
        val = bpa.(set{1});
    else
        val = bpa.Theta;
    end
end