diff --git a/pqs-advance/advance-boot/src/main/java/com/njcn/advance/responsibility/analysis/CanonicalCorrelationAnalysis.java b/pqs-advance/advance-boot/src/main/java/com/njcn/advance/responsibility/analysis/CanonicalCorrelationAnalysis.java index f0e1085a9..4d6a8bb98 100644 --- a/pqs-advance/advance-boot/src/main/java/com/njcn/advance/responsibility/analysis/CanonicalCorrelationAnalysis.java +++ b/pqs-advance/advance-boot/src/main/java/com/njcn/advance/responsibility/analysis/CanonicalCorrelationAnalysis.java @@ -9,107 +9,309 @@ import org.slf4j.LoggerFactory; /** * 典则相关分析类 * 实现典则相关系数的计算 - * + * * @author hongawen * @version 1.0 */ public class CanonicalCorrelationAnalysis { - + private static final Logger logger = LoggerFactory.getLogger(CanonicalCorrelationAnalysis.class); - + /** * 计算典则相关系数 * 对应C代码中的TransCancor函数 - * - * @param powerData 功率数据矩阵 [时间][节点] + * + * @param powerData 功率数据矩阵 [时间][节点] * @param harmonicData 谐波数据向量 - * @param windowSize 窗口大小 - * @param nodeCount 节点数量 + * @param windowSize 窗口大小 + * @param nodeCount 节点数量 * @return 典则相关系数 */ - public static float computeCanonicalCorrelation(float[][] powerData, float[] harmonicData, + public static float computeCanonicalCorrelation(float[][] powerData, float[] harmonicData, int windowSize, int nodeCount) { + logger.info("===== 开始典型相关分析 ====="); + logger.info("输入参数: windowSize={}, nodeCount={}", windowSize, nodeCount); + try { // 提取窗口数据 double[][] x = new double[windowSize][nodeCount]; double[] y = new double[windowSize]; - + + // ===== 数据质量统计（只统计，不影响计算） ===== + int nanCountPower = 0, infiniteCountPower = 0, zeroCountPower = 0; + int nanCountHarmonic = 0, infiniteCountHarmonic = 0, zeroCountHarmonic = 0; + double powerSum = 0, harmonicSum = 0; + double powerMin = Double.MAX_VALUE, powerMax = -Double.MAX_VALUE; + double harmonicMin = Double.MAX_VALUE, harmonicMax = -Double.MAX_VALUE; + for (int i = 0; i < windowSize; i++) { for (int j = 0; j < nodeCount; j++) { - x[i][j] = powerData[i][j]; + float val = powerData[i][j]; + x[i][j] = val; + + // 仅统计，不改变原逻辑 + if (Float.isNaN(val)) { + nanCountPower++; + } else if (Float.isInfinite(val)) { + infiniteCountPower++; + } else if (val == 0.0f) { + zeroCountPower++; + } else { + powerSum += val; + powerMin = Math.min(powerMin, val); + powerMax = Math.max(powerMax, val); + } + } + + float harmonicVal = harmonicData[i]; + y[i] = harmonicVal; + + // 仅统计，不改变原逻辑 + if (Float.isNaN(harmonicVal)) { + nanCountHarmonic++; + } else if (Float.isInfinite(harmonicVal)) { + infiniteCountHarmonic++; + } else if (harmonicVal == 0.0f) { + zeroCountHarmonic++; + } else { + harmonicSum += harmonicVal; + harmonicMin = Math.min(harmonicMin, harmonicVal); + harmonicMax = Math.max(harmonicMax, harmonicVal); } - y[i] = harmonicData[i]; } - + + // ===== 数据质量报告（只记录日志） ===== + int totalPowerCount = windowSize * nodeCount; + int totalHarmonicCount = windowSize; + + logger.info("功率数据质量分析:"); + logger.info(" 总数据点: {}", totalPowerCount); + logger.info(" NaN数量: {} ({:.2f}%)", nanCountPower, nanCountPower * 100.0 / totalPowerCount); + logger.info(" 无穷大数量: {} ({:.2f}%)", infiniteCountPower, infiniteCountPower * 100.0 / totalPowerCount); + logger.info(" 零值数量: {} ({:.2f}%)", zeroCountPower, zeroCountPower * 100.0 / totalPowerCount); + logger.info(" 有效数据范围: [{}, {}]", powerMin == Double.MAX_VALUE ? "N/A" : powerMin, + powerMax == -Double.MAX_VALUE ? "N/A" : powerMax); + + logger.info("谐波数据质量分析:"); + logger.info(" 总数据点: {}", totalHarmonicCount); + logger.info(" NaN数量: {} ({:.2f}%)", nanCountHarmonic, nanCountHarmonic * 100.0 / totalHarmonicCount); + logger.info(" 无穷大数量: {} ({:.2f}%)", infiniteCountHarmonic, infiniteCountHarmonic * 100.0 / totalHarmonicCount); + logger.info(" 零值数量: {} ({:.2f}%)", zeroCountHarmonic, zeroCountHarmonic * 100.0 / totalHarmonicCount); + logger.info(" 有效数据范围: [{}, {}]", harmonicMin == Double.MAX_VALUE ? "N/A" : harmonicMin, + harmonicMax == -Double.MAX_VALUE ? "N/A" : harmonicMax); + + // 只记录警告，不停止计算 + if (nanCountPower > 0 || infiniteCountPower > 0) { + logger.warn("功率数据包含异常值！NaN: {}, Infinite: {}", nanCountPower, infiniteCountPower); + } + if (nanCountHarmonic > 0 || infiniteCountHarmonic > 0) { + logger.warn("谐波数据包含异常值！NaN: {}, Infinite: {}", nanCountHarmonic, infiniteCountHarmonic); + } + // 计算协方差矩阵 SXX + logger.info("===== 开始协方差计算 ====="); double[][] sxxMatrix = MathUtils.covarianceMatrix(x, windowSize, nodeCount); - + + // ===== SXX矩阵诊断（只记录日志）===== + double sxxTrace = 0; + double sxxFrobeniusNorm = 0; + boolean sxxHasNaN = false, sxxHasInfinite = false; + + for (int i = 0; i < nodeCount; i++) { + sxxTrace += sxxMatrix[i][i]; + for (int j = 0; j < nodeCount; j++) { + double val = sxxMatrix[i][j]; + if (Double.isNaN(val)) { + sxxHasNaN = true; + } + if (Double.isInfinite(val)) { + sxxHasInfinite = true; + } + sxxFrobeniusNorm += val * val; + } + } + sxxFrobeniusNorm = Math.sqrt(sxxFrobeniusNorm); + + logger.info("SXX矩阵诊断:"); + logger.info(" 维度: {}x{}", nodeCount, nodeCount); + logger.info(" 迹(trace): {}", sxxTrace); + logger.info(" Frobenius范数: {}", sxxFrobeniusNorm); + logger.info(" 包含NaN: {}", sxxHasNaN); + logger.info(" 包含无穷大: {}", sxxHasInfinite); + logger.info(" 对角线元素: {}", java.util.Arrays.toString( + java.util.stream.IntStream.range(0, nodeCount) + .mapToDouble(i -> sxxMatrix[i][i]) + .toArray())); + // 计算协方差 SYY double syyMatrix = MathUtils.covariance(y, y, windowSize); + logger.info("SYY协方差: {}", syyMatrix); + if (Math.abs(syyMatrix) < HarmonicConstants.MIN_COVARIANCE) { + logger.warn("SYY过小 ({}), 调整为最小值: {}", syyMatrix, HarmonicConstants.MIN_COVARIANCE); syyMatrix = HarmonicConstants.MIN_COVARIANCE; } - + // 计算协方差向量 SXY double[] sxyVector = MathUtils.covarianceVector(x, y, windowSize, nodeCount); - + + // ===== SXY向量诊断（只记录日志）===== + double sxyNorm = 0; + boolean sxyHasNaN = false, sxyHasInfinite = false; + for (double val : sxyVector) { + if (Double.isNaN(val)) { + sxyHasNaN = true; + } + if (Double.isInfinite(val)) { + sxyHasInfinite = true; + } + sxyNorm += val * val; + } + sxyNorm = Math.sqrt(sxyNorm); + + logger.info("SXY向量诊断:"); + logger.info(" 长度: {}", sxyVector.length); + logger.info(" L2范数: {}", sxyNorm); + logger.info(" 包含NaN: {}", sxyHasNaN); + logger.info(" 包含无穷大: {}", sxyHasInfinite); + logger.info(" 向量值: {}", java.util.Arrays.toString(sxyVector)); + // 使用Apache Commons Math进行矩阵运算 + logger.info("===== 开始矩阵分解 ====="); RealMatrix sxx = new Array2DRowRealMatrix(sxxMatrix); RealVector sxy = new ArrayRealVector(sxyVector); - + // 计算 SXX^(-1) + logger.info("准备计算SXX逆矩阵..."); DecompositionSolver solver = new LUDecomposition(sxx).getSolver(); RealMatrix invSxx; - + + // 检查矩阵奇异性 + double sxxDet = new LUDecomposition(sxx).getDeterminant(); + logger.info("SXX矩阵行列式: {}", sxxDet); + + if (Math.abs(sxxDet) < 1e-15) { + logger.warn("SXX矩阵几乎奇异 (det={})", sxxDet); + } + if (!solver.isNonSingular()) { // 如果矩阵奇异，使用伪逆 logger.warn("SXX matrix is singular, using pseudo-inverse"); - SingularValueDecomposition svd = new SingularValueDecomposition(sxx); - invSxx = svd.getSolver().getInverse(); + try { + SingularValueDecomposition svd = new SingularValueDecomposition(sxx); + invSxx = svd.getSolver().getInverse(); + } catch (Exception svdException) { + logger.error("SVD pseudo-inverse failed, using regularized inverse", svdException); + // 添加正则化项 + RealMatrix identity = MatrixUtils.createRealIdentityMatrix(sxx.getRowDimension()); + RealMatrix regularized = sxx.add(identity.scalarMultiply(1e-10)); + invSxx = new LUDecomposition(regularized).getSolver().getInverse(); + } } else { invSxx = solver.getInverse(); } - + // 计算 U = SXX^(-1) * SXY * (1/SYY) * SXY' RealVector temp = invSxx.operate(sxy); double scale = 1.0 / syyMatrix; RealMatrix uMatrix = temp.outerProduct(sxy).scalarMultiply(scale); - - // 计算特征值 - EigenDecomposition eigenDecomposition = new EigenDecomposition(uMatrix); - double[] eigenvalues = eigenDecomposition.getRealEigenvalues(); - - // 找最大特征值 + + // 计算特征值 - 添加数值稳定性处理 double maxEigenvalue = 0.0; - for (double eigenvalue : eigenvalues) { - maxEigenvalue = Math.max(maxEigenvalue, Math.abs(eigenvalue)); + + try { + // 首先检查矩阵是否有效 + double[][] uMatrixData = uMatrix.getData(); + boolean hasNaN = false; + boolean hasInfinite = false; + + for (int i = 0; i < uMatrixData.length; i++) { + for (int j = 0; j < uMatrixData[i].length; j++) { + if (Double.isNaN(uMatrixData[i][j])) { + hasNaN = true; + } + if (Double.isInfinite(uMatrixData[i][j])) { + hasInfinite = true; + } + } + } + + if (hasNaN || hasInfinite) { + logger.warn("U matrix contains NaN or Infinite values, returning 0"); + return 0.0f; + } + + // 检查矩阵条件数 + SingularValueDecomposition svdCheck = new SingularValueDecomposition(uMatrix); + double conditionNumber = svdCheck.getConditionNumber(); + + if (conditionNumber > 1e12) { + logger.warn("U matrix is ill-conditioned (condition number: {}), using SVD approach", conditionNumber); + + // 使用SVD方法获取最大奇异值的平方作为最大特征值 + double[] singularValues = svdCheck.getSingularValues(); + if (singularValues.length > 0) { + maxEigenvalue = singularValues[0] * singularValues[0]; + } + } else { + // 正常的特征值分解 + EigenDecomposition eigenDecomposition = new EigenDecomposition(uMatrix); + double[] eigenvalues = eigenDecomposition.getRealEigenvalues(); + + // 找最大特征值 + for (double eigenvalue : eigenvalues) { + maxEigenvalue = Math.max(maxEigenvalue, Math.abs(eigenvalue)); + } + } + + } catch (Exception eigenException) { + logger.warn("EigenDecomposition failed, trying alternative approach: {}", eigenException.getMessage()); + + // 备用方案：使用SVD方法 + try { + SingularValueDecomposition svd = new SingularValueDecomposition(uMatrix); + double[] singularValues = svd.getSingularValues(); + if (singularValues.length > 0) { + maxEigenvalue = singularValues[0] * singularValues[0]; + } + } catch (Exception svdException) { + logger.error("Both EigenDecomposition and SVD failed, returning 0", svdException); + return 0.0f; + } } - + // 典则相关系数是最大特征值的平方根 - double canonicalCorr = Math.sqrt(maxEigenvalue); - + double canonicalCorr = Math.sqrt(Math.abs(maxEigenvalue)); + // 限制在[0,1]范围内 if (canonicalCorr > 1.0) { canonicalCorr = 1.0; } - + + logger.info("===== 典型相关分析计算完成 ====="); + logger.info("最大特征值: {}", maxEigenvalue); + logger.info("典型相关系数: {}", canonicalCorr); + logger.info("是否被截断到1.0: {}", canonicalCorr >= 1.0); + return (float) canonicalCorr; - + } catch (Exception e) { logger.error("Error computing canonical correlation", e); + logger.error("异常详情: {}", e.getMessage()); + logger.error("异常类型: {}", e.getClass().getSimpleName()); return 0.0f; } } - + /** * 滑动窗口计算典则相关系数序列 * 对应C代码中的SlideCanCor函数 - * - * @param powerData 功率数据矩阵 [时间][节点] + * + * @param powerData 功率数据矩阵 [时间][节点] * @param harmonicData 谐波数据向量 - * @param windowSize 窗口大小 - * @param nodeCount 节点数量 - * @param dataLength 数据总长度 + * @param windowSize 窗口大小 + * @param nodeCount 节点数量 + * @param dataLength 数据总长度 * @return 典则相关系数序列 */ public static float[] slidingCanonicalCorrelation(float[][] powerData, float[] harmonicData, @@ -118,62 +320,62 @@ public class CanonicalCorrelationAnalysis { if (slideLength <= 0) { throw new IllegalArgumentException("Data length must be greater than window size"); } - + float[] slideCanCor = new float[slideLength]; - + logger.info("Starting sliding canonical correlation analysis, slide length: {}", slideLength); - + for (int i = 0; i < slideLength; i++) { // 提取窗口数据 float[][] windowPower = new float[windowSize][nodeCount]; float[] windowHarmonic = new float[windowSize]; - + for (int j = 0; j < windowSize; j++) { System.arraycopy(powerData[i + j], 0, windowPower[j], 0, nodeCount); windowHarmonic[j] = harmonicData[i + j]; } - + // 计算当前窗口的典则相关系数 - slideCanCor[i] = computeCanonicalCorrelation(windowPower, windowHarmonic, - windowSize, nodeCount); - + slideCanCor[i] = computeCanonicalCorrelation(windowPower, windowHarmonic, + windowSize, nodeCount); + if (i % 10 == 0) { logger.debug("Processed window {}/{}", i, slideLength); } } - + logger.info("Sliding canonical correlation analysis completed"); - + return slideCanCor; } - + /** * 计算包含/不包含背景的动态相关系数 * 对应C代码中的SlideCor函数 - * - * @param powerData 功率数据（单个节点） + * + * @param powerData 功率数据（单个节点） * @param harmonicData 谐波数据 - * @param slideCanCor 滑动典则相关系数 - * @param windowSize 窗口大小 + * @param slideCanCor 滑动典则相关系数 + * @param windowSize 窗口大小 * @return 动态相关系数序列 */ public static float[] slidingCorrelation(float[] powerData, float[] harmonicData, - float[] slideCanCor, int windowSize) { + float[] slideCanCor, int windowSize) { int slideLength = slideCanCor.length; float[] slideCor = new float[slideLength]; - + for (int i = 0; i < slideLength; i++) { float[] tempPower = new float[windowSize]; float[] tempHarmonic = new float[windowSize]; - + for (int j = 0; j < windowSize; j++) { tempPower[j] = powerData[i + j]; tempHarmonic[j] = harmonicData[i + j] * slideCanCor[i]; } - + slideCor[i] = MathUtils.pearsonCorrelation(tempHarmonic, tempPower, windowSize); } - + return slideCor; } } \ No newline at end of file diff --git a/pqs-advance/advance-boot/src/main/java/com/njcn/advance/responsibility/calculator/ResponsibilityCalculator.java b/pqs-advance/advance-boot/src/main/java/com/njcn/advance/responsibility/calculator/ResponsibilityCalculator.java index 26e131786..ed4fe1533 100644 --- a/pqs-advance/advance-boot/src/main/java/com/njcn/advance/responsibility/calculator/ResponsibilityCalculator.java +++ b/pqs-advance/advance-boot/src/main/java/com/njcn/advance/responsibility/calculator/ResponsibilityCalculator.java @@ -1,5 +1,6 @@ package com.njcn.advance.responsibility.calculator; +import com.njcn.advance.responsibility.analysis.CanonicalCorrelationAnalysis; import org.slf4j.Logger; import org.slf4j.LoggerFactory; @@ -201,8 +202,18 @@ public class ResponsibilityCalculator { double[] HKSum = new double[columnCount]; // 对应C代码中的 VectorXd HKSum - 使用double精度 int exceedCount = 0; // 对应C代码中的 coutt - logger.debug("sumResponsibility参数: threshold={}, windowSize={}, columnCount={}, tl_num={}, slideLength={}", - threshold, windowSize, columnCount, tl_num, slideLength); + // ===== 添加详细调试日志 ===== + logger.info("======= 开始 sumResponsibility 计算 ======="); + logger.info("输入参数:"); + logger.info(" threshold(阈值): {}", threshold); + logger.info(" windowSize(窗口大小): {}", windowSize); + logger.info(" columnCount(列数): {}", columnCount); + logger.info(" tl_num(总数据长度): {}", tl_num); + logger.info(" slideLength(滑动长度): {}", slideLength); + logger.info(" responsibilityData维度: {}x{}", + responsibilityData != null ? responsibilityData.length : "null", + responsibilityData != null && responsibilityData.length > 0 ? responsibilityData[0].length : "null"); + logger.info(" harmonicData长度: {}", harmonicData != null ? harmonicData.length : "null"); // 数据验证 if (harmonicData == null) { @@ -218,7 +229,7 @@ public class ResponsibilityCalculator { // 关键验证：检查数组长度是否充足 // C代码中Udata长度是TL，循环遍历slg=TL-width个元素 - logger.debug("数据验证: slideLength={}, harmonicData.length={}, responsibilityData.length={}", + logger.info("数据验证: slideLength={}, harmonicData.length={}, responsibilityData.length={}", slideLength, harmonicData.length, responsibilityData.length); if (harmonicData.length < slideLength) { @@ -237,13 +248,52 @@ public class ResponsibilityCalculator { String.format("责任数据行数不足: 需要%d, 实际%d", slideLength, responsibilityData.length)); } + // ===== 添加数据分布统计 ===== + logger.info("谐波数据分析:"); + float harmonicMin = Float.MAX_VALUE, harmonicMax = Float.MIN_VALUE; + double harmonicSum = 0; + int preliminaryExceedCount = 0; + for (int i = 0; i < Math.min(slideLength, harmonicData.length); i++) { + float val = harmonicData[i]; + harmonicMin = Math.min(harmonicMin, val); + harmonicMax = Math.max(harmonicMax, val); + harmonicSum += val; + if (val > threshold) { + preliminaryExceedCount++; + } + } + double harmonicAvg = harmonicSum / Math.min(slideLength, harmonicData.length); + logger.info(" 谐波数据范围: [{}, {}]", harmonicMin, harmonicMax); + logger.info(" 谐波数据平均值: {}", harmonicAvg); + logger.info(" 设定阈值: {}", threshold); + logger.info(" 初步统计超限个数: {}/{} ({:.2f}%)", + preliminaryExceedCount, Math.min(slideLength, harmonicData.length), + preliminaryExceedCount * 100.0 / Math.min(slideLength, harmonicData.length)); + + // ===== 责任数据分析 ===== + logger.info("责任数据分析(检查前5行的归一化情况):"); + for (int i = 0; i < Math.min(5, responsibilityData.length); i++) { + float rowSum = 0; + for (int j = 0; j < responsibilityData[i].length; j++) { + rowSum += responsibilityData[i][j]; + } + logger.info(" 第{}行和: {} (期望值: 1.0, 偏差: {})", i, rowSum, Math.abs(rowSum - 1.0f)); + } + // 统计超限时段的责任 - 对应C代码行437-449 // 重要：C代码中有一个设计缺陷：coutt在每个j循环中被重置， // 但最后计算百分比时使用的是最后一次j循环的coutt值 // 为了严格保持一致，我们也要复现这个逻辑 + logger.info("===== 开始循环计算每列的累加值 ====="); + int[] exceedCountPerColumn = new int[columnCount]; // 记录每列的超限次数，用于调试 + for (int j = 0; j < columnCount; j++) { HKSum[j] = 0; exceedCount = 0; // 对应C代码行440: coutt = 0; + logger.info("开始计算第{}列 (共{}列)", j, columnCount); + + double columnSum = 0; // 用于调试 + int columnExceedCount = 0; // 用于调试 for (int i = 0; i < slideLength; i++) { // 添加越界检查 @@ -257,31 +307,64 @@ public class ResponsibilityCalculator { // 对应C代码行443-447 if (harmonicData[i] > threshold) { - HKSum[j] += responsibilityData[i][j]; // 对应C代码行445 + double currentResponsibility = responsibilityData[i][j]; + HKSum[j] += currentResponsibility; // 对应C代码行445 exceedCount++; // 对应C代码行446 + columnSum += currentResponsibility; + columnExceedCount++; + + // 只打印前几个超限情况的详细信息 + if (columnExceedCount <= 3) { + logger.info(" 时刻i={}: 谐波值={} > 阈值={}, 责任值={}, 累加到{}", + i, harmonicData[i], threshold, currentResponsibility, HKSum[j]); + } } } + + exceedCountPerColumn[j] = columnExceedCount; + logger.info("第{}列计算完成: 累加值={}, 超限次数={}", j, HKSum[j], columnExceedCount); } // 注意：这里exceedCount保留的是最后一列（j=columnCount-1）的超限次数 // 这与C代码的行为一致 - logger.debug("最终exceedCount={} (来自最后一列的计算)", exceedCount); + logger.info("===== 循环计算完成 ====="); + logger.info("最终exceedCount={} (来自最后一列的计算)", exceedCount); + logger.info("各列超限次数对比: {}", java.util.Arrays.toString(exceedCountPerColumn)); + logger.info("各列累加值: {}", java.util.Arrays.toString(HKSum)); // 计算平均责任百分比 - 对应C代码行453-459 + logger.info("===== 开始计算最终百分比 ====="); for (int i = 0; i < columnCount; i++) { sumData[i] = 0; // 对应C代码行454 } + double totalPercentage = 0; // 用于统计总和 for (int i = 0; i < columnCount; i++) { if (exceedCount > 0) { // 对应C代码行458: arrHKsum[i] = 100 * (HKSum(i)/coutt); // 使用double进行计算，然后转换为float - sumData[i] = (float)(100.0 * (HKSum[i] / (double)exceedCount)); + double percentage = 100.0 * (HKSum[i] / (double)exceedCount); + sumData[i] = (float)percentage; + totalPercentage += percentage; + + logger.info("节点{}: 累加值={}, 除以超限次数={}, 百分比={}%", + i, HKSum[i], exceedCount, percentage); + } else { + logger.warn("节点{}: 超限次数为0，百分比设为0", i); } } - logger.info("Exceeded count: {}, average responsibilities calculated", exceedCount); + logger.info("===== 计算结果汇总 ====="); + logger.info("使用的超限次数(分母): {}", exceedCount); + logger.info("各节点百分比: {}", java.util.Arrays.toString(sumData)); + logger.info("百分比总和: {}% (期望100%)", totalPercentage); + logger.info("偏差: {}%", Math.abs(totalPercentage - 100.0)); + if (Math.abs(totalPercentage - 100.0) > 1.0) { + logger.warn("!!!注意!!! 百分比总和偏离100%超过1%，可能存在问题"); + } + + logger.info("======= sumResponsibility 计算完成 ======="); return sumData; } @@ -330,7 +413,7 @@ public class ResponsibilityCalculator { } // 计算该节点的动态相关系数 - float[] nodeCorr = com.njcn.advance.responsibility.analysis.CanonicalCorrelationAnalysis + float[] nodeCorr = CanonicalCorrelationAnalysis .slidingCorrelation(nodePower, harmonicData, canonicalCorr, windowSize); // 存储结果