apacsys-vaporrecoverymonitor-document/程序设计文档/油气回收智能诊断系统/实验进油数据/ai_modify_new2.txt

clc; clear;

%% ---- 压力矩阵（P_F1~P_F4）----
P_F1 = [ ...
    0 0 0 14 36 54 61 60 56 52 57 55 50 45 41 37 33 31 29 28;
    0 0 2 27 68 91 99 89 75 65 56 55 61 62 59 55 51 53 59 59;
    0 0 5 41 75 94 104 107 111 109 97 83 71 65 61 59 57 56 55 55];

P_F2 = [ ...
    0 0 0 51 86 129 114 95 73 69 66 81 99 97 75 60 53 47 44 46;
    0 0 7 96 149 138 130 139 134 129 126 125 123 121 120 115 100 70 60 52;
    0 0 19 115 150 141 136 140 135 128 127 126 123 121 120 119 118 118 117 117];

P_F3 = [ ...
    0 0 0 75 196 154 185 162 144 147 139 134 130 127 125 123 121 120 120 119;
    0 0 10 152 184 187 173 163 142 140 137 131 130 126 125 123 123 121 121 120;
    0 0 23 190 183 185 156 157 141 139 135 131 127 125 124 122 122 121 120 120];

P_F4 = [ ...
    0 0 0 101 226 226 183 170 156 143 141 141 137 135 131 129 126 124 123 122;
    1 1 12 176 226 211 200 168 149 142 139 134 132 129 127 124 122 121 120 119;
    0 0 21 226 226 194 164 152 148 146 139 136 131 129 126 125 123 122 120 120];

%% ---- 流速矩阵（V_F1~V_F4）----
V_F1 = [ ...
    0.3 0 1.4 17.4 28.3 35.3 37.2 37.4 35.4 35.3 36.2 35.6 34 32.6 30.8 29 27.7 26.3 25.4 24.8;
    0.3 0 4.0 24.5 39.2 45.9 47.8 46.0 41.7 39.0 36.1 36.1 37.8 37.7 37.3 35.4 34.8 35.0 37.1 37.3;
    0.3 0 7.6 29.8 41.4 47.1 48.6 50.5 50.2 50.7 47.5 44.3 40.8 39.5 37.7 37.3 36.3 36.1 35.7 35.9];

V_F2 = [ ...
    0.3 0 0.5 3.1 21.8 31.7 37.7 37.8 28.8 2.7 7.8 25.4 35.1 40.2 37 33.9 32.6 31.4 30.7 31.2;
    0.4 0 0.2 12.7 30.2 40.1 43.9 25.4 6.4 11.2 32.4 39.4 40.9 45.3 45.7 46.8 44.1 40.4 36.2 34.6;
    0.3 0 0.1 19.6 32.2 40.8 44.1 19.1 5 10.6 19.2 32.8 41.3 45.9 45.7 47 47.5 47.4 48.1 48.1];

V_F3 = [ ...
    0.3 0 0.3 0.9 7.9 27.5 17.7 18.8 10.1 5.9 6.3 7.4 15.2 25 23.9 24.6 24.3 30.1 46.1 47.7;
    0.3 0 0.2 1.7 23.1 25.3 17.7 13.8 9.3 3.2 9.1 9.1 12.7 15.5 24.3 24.6 33.8 45.3 46.5 47.4;
    0.4 0 0.6 3.6 27.9 23.3 16.7 11.4 10 4.4 4.5 11.9 12.9 20 23.8 32.3 44.3 46 47.3 47.2];

V_F4 = [ ...
    0.3 0 0.9 1.8 6.3 19.6 17.6 19 16.6 10.3 6.6 9.9 11.3 13.4 18.1 21.9 26.9 24.5 24.7 26.4;
    0.2 0 0.2 1.5 10.7 16.4 18.4 16.7 12.2 8.3 8.2 12.3 13.9 17.7 21.7 25.8 27 30.7 33.2 34.1;
    0.3 0 0.4 3.3 15.8 16.1 16.6 17.4 11.7 5.6 8.4 9.9 12.3 16.3 18.7 22.9 27.7 23.1 22 23.5];

P_lib = {P_F1, P_F2, P_F3, P_F4};
V_lib = {V_F1, V_F2, V_F3, V_F4};
faultNames = {'F1','F2','F3','F4'};

rho = 0.5;
max_window = 10;

%% =========================
% 1) 模板区
% =========================
% 最初算法模板（若你有固定Y1/Y2，直接替换这两个变量）
Y1_old = build_template_from_matrix_mean(P_lib, max_window);
Y2_old = build_template_from_matrix_mean(V_lib, max_window);

% 改进算法模板
Y1_new = build_template_from_matrix_mean(P_lib, max_window);
Y2_new = build_template_from_matrix_mean(V_lib, max_window);

%% =========================
% 2) 模式1：保留原来自定义 P、V
% =========================
% 你原来自定义的 P、V 放这里：
P_custom = [118 118 117 118 117 99 65 51 41 35 30 28 26 26 25 24 24 23 24 24];
V_custom = [46.6 48.3 48.4 49.1 48.6 46.7 40.2 35.3 31.6 29.5 27.1 26.1 25.1 25.4 24.4 24.9 24.9 23.8 24.1 24.5];

fprintf('\n==================== 模式1：自定义 P,V ====================\n');
res_custom = run_dual_alg(P_custom, V_custom, Y1_old, Y2_old, Y1_new, Y2_new, rho, max_window, faultNames);
print_compare_result('自定义样本', res_custom, faultNames);

%% =========================
% 3) 模式2：批量遍历库中每一行
% =========================
fprintf('\n==================== 模式2：批量样本 ====================\n');

allRows = {};
idxRow = 0;

for f = 1:4
    P_mat = P_lib{f};
    V_mat = V_lib{f};
    nRow = size(P_mat,1);

    for r = 1:nRow
        idxRow = idxRow + 1;
        tag = sprintf('P_%s第%d行 + V_%s第%d行', faultNames{f}, r, faultNames{f}, r);

        res = run_dual_alg(P_mat(r,:), V_mat(r,:), Y1_old, Y2_old, Y1_new, Y2_new, rho, max_window, faultNames);
        print_compare_result(tag, res, faultNames);

        allRows{idxRow,1} = tag;
        allRows{idxRow,2} = faultNames{res.old_idx};
        allRows{idxRow,3} = res.old_conf;
        allRows{idxRow,4} = faultNames{res.new_idx};
        allRows{idxRow,5} = res.new_conf;
        allRows{idxRow,6} = strcmp(faultNames{res.old_idx}, faultNames{res.new_idx});
        allRows{idxRow,7} = abs(res.new_conf - res.old_conf);
    end
end

%% =========================
% 4) 汇总对比表
% =========================
fprintf('\n==================== 汇总对比表 ====================\n');
fprintf('%-28s %-6s %-10s %-6s %-10s %-8s %-10s\n', ...
    '样本', 'Old', 'OldConf', 'New', 'NewConf', '一致?', '置信度差');
for i = 1:size(allRows,1)
    fprintf('%-28s %-6s %-10.2f %-6s %-10.2f %-8d %-10.2f\n', ...
        allRows{i,1}, allRows{i,2}, allRows{i,3}*100, ...
        allRows{i,4}, allRows{i,5}*100, allRows{i,6}, allRows{i,7}*100);
end

fprintf('\n完成：已同时保留"自定义样本"与"批量样本"两种输入，并输出两种算法对比结果。\n');

%% ========================= 函数区 =========================

function res = run_dual_alg(P, V, Y1_old, Y2_old, Y1_new, Y2_new, rho, max_window, faultNames)
Y01 = calc_feature_row(P, max_window);
Y02 = calc_feature_row(V, max_window);

% ------- 最初算法：不归一化BPA + 原始DS -------
deg1_old = grey_relation_analysis1(Y01, Y1_old, rho);
deg2_old = grey_relation_analysis1(Y02, Y2_old, rho);

comb_old = dempster_combine_original(deg1_old, deg2_old);
vec_old = [comb_old.F1, comb_old.F2, comb_old.F3, comb_old.F4];
[old_conf, old_idx] = max(vec_old);

% ------- 改进算法：归一化BPA + 修正DS(intersect) -------
deg1_new = grey_relation_analysis1(Y01, Y1_new, rho);
deg2_new = grey_relation_analysis1(Y02, Y2_new, rho);

bpa1 = deg1_new(:) / sum(deg1_new);
bpa2 = deg2_new(:) / sum(deg2_new);

comb_new = dempster_combine_fixed(bpa1, bpa2);
vec_new = [comb_new.F1, comb_new.F2, comb_new.F3, comb_new.F4];
[new_conf, new_idx] = max(vec_new);

res.old_vec = vec_old; res.old_conf = old_conf; res.old_idx = old_idx;
res.new_vec = vec_new; res.new_conf = new_conf; res.new_idx = new_idx;
res.same = strcmp(faultNames{old_idx}, faultNames{new_idx});
res.delta_conf = abs(new_conf - old_conf);
end

function print_compare_result(tag, res, faultNames)
fprintf('\n样本: %s\n', tag);
fprintf('  [最初算法] BPA=[%.4f %.4f %.4f %.4f], 诊断=%s, 置信度=%.2f%%\n', ...
    res.old_vec(1), res.old_vec(2), res.old_vec(3), res.old_vec(4), ...
    faultNames{res.old_idx}, res.old_conf*100);
fprintf('  [改进算法] BPA=[%.4f %.4f %.4f %.4f], 诊断=%s, 置信度=%.2f%%\n', ...
    res.new_vec(1), res.new_vec(2), res.new_vec(3), res.new_vec(4), ...
    faultNames{res.new_idx}, res.new_conf*100);
fprintf('  [对比] 是否一致=%d, 置信度差=%.2f%%\n', res.same, res.delta_conf*100);
end

function feat = calc_feature_row(x, max_window)
x = x(:)';
max_window = min(max_window, numel(x));
feat = [mean(x), var(x), std(x), max(x(1:max_window)), min(x)];
end

function Y = build_template_from_matrix_mean(libCell, max_window)
Y = zeros(4,5);
for k = 1:4
    X = libCell{k};
    n = size(X,1);
    tmp = zeros(n,5);
    for i = 1:n
        tmp(i,:) = calc_feature_row(X(i,:), max_window);
    end
    Y(k,:) = mean(tmp,1);
end
end

function [grey_relation_degree, order] = grey_relation_analysis1(Y0, Y, rho)
if nargin < 3, rho = 0.5; end
Y0 = Y0(:)';
[n,m] = size(Y);
if m ~= length(Y0), error('参考序列和比较序列长度必须一致'); end
Y0n = Y0 / mean(Y0);
Yn  = Y ./ mean(Y,2);
diff = abs(Yn - Y0n);
dmin = min(diff(:)); dmax = max(diff(:));
rel = (dmin + rho*dmax) ./ (diff + rho*dmax);
grey_relation_degree = mean(rel,2);
[~, order] = sort(grey_relation_degree,'descend');
end

function combined_bpa = dempster_combine_original(m1, m2)
K = 0;
for i=1:4
    for j=1:4
        if i~=j
            K = K + m1(i)*m2(j);
        end
    end
end
combined_bpa = struct();
combined_bpa.F1 = (m1(1)*m2(1)) / (1-K);
combined_bpa.F2 = (m1(2)*m2(2)) / (1-K);
combined_bpa.F3 = (m1(3)*m2(3)) / (1-K);
combined_bpa.F4 = (m1(4)*m2(4)) / (1-K);
end

function combined_bpa = dempster_combine_fixed(m1, m2)
sets = {{'F1'},{'F2'},{'F3'},{'F4'}};
K = 0;
for i=1:4
    for j=1:4
        if isempty(intersect(sets{i}, sets{j}))
            K = K + m1(i)*m2(j);
        end
    end
end
vals = zeros(1,4);
for k=1:4
    s = 0;
    for i=1:4
        for j=1:4
            interSet = intersect(sets{i}, sets{j});
            if isequal(interSet, sets{k})
                s = s + m1(i)*m2(j);
            end
        end
    end
    vals(k) = s/(1-K);
end
combined_bpa = struct('F1',vals(1),'F2',vals(2),'F3',vals(3),'F4',vals(4));
end