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