Microservice/pqs
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 1 Wiki Activity

波形解析代码移植

Browse Source
This commit is contained in:
hongawen
2025-09-15 09:54:46 +08:00
parent 74420e7107
commit c1bf8d79ce
2 changed files with 352 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files

320
pqs-advance/advance-boot/src/main/java/com/njcn/advance/responsibility/analysis/CanonicalCorrelationAnalysis.java
View File

@@ -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;
}
}

99
pqs-advance/advance-boot/src/main/java/com/njcn/advance/responsibility/calculator/ResponsibilityCalculator.java
View File

@@ -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);
// 存储结果