using Accord.Math;
|
using Accord.Math.Optimization;
|
using DevExpress.Data.Extensions;
|
using MathNet.Numerics;
|
using MathNet.Numerics.LinearAlgebra;
|
using MathNet.Numerics.LinearRegression;
|
using MathNet.Numerics.Statistics;
|
|
namespace IStation.Test
|
{
|
|
/// <summary>
|
/// 泵曲线数据融合校正器（用于合并样条曲线和实测数据并进行优化）
|
/// 增加了数据聚合和动态多项式拟合阶数的功能
|
/// </summary>
|
public class PumpCurveDataFusionCorrectorHelper2
|
{
|
# region 配置参数
|
// 异常值检测阈值（基于标准正态分布，2.5对应98.76%置信区间）
|
public const double ZScoreThreshold = 2;
|
// 过渡区域宽度（单位与输入数据相同，控制曲线衔接平滑度）
|
public const double TransitionWidth = 500;
|
// 多项式拟合阶数（三次平衡拟合能力与过拟合风险）
|
public const int PolynomialDegree = 3;
|
// 优化器最大迭代次数（Cobyla算法收敛阈值）
|
private const int OptimizationMaxIterations = 1000;
|
|
//动态多项式拟合阶数（用于动态调整）
|
public int DynamicPolyDegree = 2;
|
|
# endregion
|
|
/// <summary>
|
/// 主校正方法（入口函数）
|
/// </summary>
|
/// <param name="splineX">样条曲线X坐标</param>
|
/// <param name="splineY">样条曲线Y坐标</param>
|
/// <param name="measuredXAllFiltered">原始实测X坐标</param>
|
/// <param name="measuredYAllFiltered">原始实测Y坐标</param>
|
/// <returns>
|
/// mergedX: 融合后的X坐标
|
/// mergedY: 初步融合的Y坐标
|
/// optimizedX: 优化后的X坐标采样
|
/// optimizedY: 优化后的Y坐标结果
|
/// </returns>
|
public (double[] mergedX, double[] mergedY, double[] optimizedX, double[] optimizedY) Corrent(
|
double[] splineX, double[] splineY,
|
double[] measuredXAll, double[] measuredYAll)
|
{
|
|
if (measuredXAll.Length != measuredYAll.Length)
|
{
|
return default;
|
}
|
|
if (measuredXAll.Length>10000)
|
{
|
|
}
|
|
var minX = measuredXAll[0];
|
var maxX = measuredXAll[^1];
|
|
var fv = 0.15;
|
var threshold = minX + fv * (maxX - minX);
|
var startIndex = measuredXAll.FindIndex(x => x > threshold);
|
|
var measuredXAllFiltered = measuredXAll.Skip(startIndex).ToArray();
|
var measuredYAllFiltered = measuredYAll.Skip(startIndex).ToArray();
|
|
|
|
if (measuredXAllFiltered.Length < 5)
|
{
|
Console.WriteLine("异常值过滤后数据过小");
|
return default;
|
}
|
|
#region 步骤1：稳健回归与异常值过滤（最小二乘法）
|
// 执行简单线性回归获取基准线
|
(double intercept, double slope) = SimpleRegression.Fit(measuredXAllFiltered, measuredYAllFiltered);
|
|
// 计算预测值和残差
|
double[] predictedY = measuredXAllFiltered.Select(x => slope * x + intercept).ToArray();
|
double[] residuals = measuredYAllFiltered.Zip(predictedY, (m, p) => m - p).ToArray();
|
|
// 基于Z-Score的异常值过滤
|
var zScoreCalculator = new ZScore(residuals);
|
bool[] validMask = zScoreCalculator.Scores
|
.Select(score => Math.Abs(score) < ZScoreThreshold).ToArray();
|
|
// 过滤得到有效实测数据
|
var validData = measuredXAllFiltered.Zip(measuredYAllFiltered, validMask)
|
.Where(t => t.Third)
|
.Select(t => (t.First, t.Second)).ToList();
|
double[] measuredXValid = validData.Select(t => t.First).ToArray();
|
double[] measuredYValid = validData.Select(t => t.Second).ToArray();
|
#endregion
|
|
|
|
if (measuredXValid.Length < 5)
|
{
|
Console.WriteLine("异常值过滤后数据过小");
|
return default;
|
}
|
|
return ProcessData(splineX, splineY, measuredXValid, measuredYValid,
|
PolynomialDegree, TransitionWidth);
|
}
|
|
/// <summary>
|
/// 核心数据处理流程（包含数据融合、过渡处理和优化）
|
/// </summary>
|
private (double[] mergedX, double[] mergedY, double[] optimizedX, double[] optimizedY) ProcessData(
|
double[] splineX, double[] splineY,
|
double[] measuredXValid, double[] measuredYValid,
|
int polyDegree, double transitionWidth)
|
{
|
|
/* labels, areas = classify_points(filtered_first_x, filtered_first_y, method = 'dbscan',
|
params={ 'eps': 0.25, 'min_samples': 15})
|
# 动态多项式次数调整
|
Config.Dynamic_POLY_DEGREE = 2 if len(areas) > 1 else 1*/
|
int min_samples = (int)(measuredXValid.Length * 0.02);
|
min_samples = min_samples < 5 ? 5 : min_samples;
|
var pointClassifier = new PointClassifier();
|
var (labels, areas) = pointClassifier.ClassifyPoints(measuredXValid, measuredYValid, "dbscan", new Dictionary<string, object>
|
{
|
{ "eps", 0.3 },
|
{ "min_samples",(int)min_samples}
|
});
|
|
// 动态多项式次数调整
|
DynamicPolyDegree = areas.Count > 1 ? 2 : 1;
|
|
|
if (areas.Count<1)
|
{
|
throw new Exception("areas.Count<1");
|
}
|
|
var measured_x = new List<double>();
|
var measured_y = new List<double>();
|
|
foreach (var item in areas)
|
{
|
var x_min = item.XRange.Min;
|
var x_max = item.XRange.Max;
|
var y_min = item.YRange.Min;
|
var y_max = item.YRange.Max;
|
var mask = measuredXValid
|
.Select((x, i) => (x >= x_min && x <= x_max && measuredYValid[i] >= y_min && measuredYValid[i] <= y_max))
|
.ToArray();
|
for (int i = 0; i < mask.Length; i++)
|
{
|
var x = mask[i];
|
if (x)
|
{
|
measured_x.Add(measuredXValid[i]);
|
measured_y.Add(measuredYValid[i]);
|
}
|
}
|
}
|
|
#region
|
|
|
double minX = measuredXValid.Min(), maxX = measuredXValid.Max();
|
var (mergedX, mergedY) = MergeCurves(splineX, splineY, measured_x.ToArray(), measured_y.ToArray(), DynamicPolyDegree);
|
#endregion
|
|
#region 步骤5：带约束优化
|
|
// 数据标准化处理
|
double xMean = mergedX.Average();
|
double xStd = ArrayStatistics.PopulationStandardDeviation(mergedX);
|
double[] xNorm = mergedX.Select(x => (x - xMean) / xStd).ToArray();
|
|
// 构建三次多项式优化问题
|
var X = Matrix<double>.Build.DenseOfRowArrays(
|
xNorm.Select(x => new[] { Math.Pow(x, 3), x * x, x, 1.0 }).ToArray()
|
);
|
var initialParams = MultipleRegression.QR(X, Vector<double>.Build.DenseOfArray(mergedY));
|
|
|
|
|
// 设置单调性约束（导数非负）
|
var constraints = GenerateLinspace(xNorm.Min(), xNorm.Max(), 50)
|
.Select(x => new NonlinearConstraint(
|
numberOfVariables: 4,
|
function: p => -Derivative(p, x), // 导数>=0
|
shouldBe: ConstraintType.GreaterThanOrEqualTo,
|
value: 0
|
)).ToList();
|
|
// 执行COBYLA优化
|
var optimizer = new Cobyla(
|
function: BuildObjectiveFunction(xNorm, mergedY),
|
constraints: constraints.ToArray()
|
)
|
{ MaxIterations = OptimizationMaxIterations };
|
|
bool success = optimizer.Minimize(initialParams.ToArray());
|
double[] optimizedParams = optimizer.Solution;
|
#endregion
|
|
|
|
var a = optimizedParams.Select(x => Math.Round(x, 6)).ToList();
|
var json = JsonHelper.Object2FormatJson(a);
|
// 生成最终优化曲线
|
double[] optimizedX = GenerateLinspace(
|
Math.Min(splineX.Min(), minX),
|
Math.Max(splineX.Max(), maxX),
|
300
|
);
|
double[] optimizedY = optimizedX
|
.Select(x => CubicModel(optimizedParams, (x - xMean) / xStd))
|
.ToArray();
|
|
return (mergedX, mergedY, optimizedX, optimizedY);
|
}
|
|
# region 工具方法
|
|
/// <summary>
|
/// 三次多项式模型（降幂排列）
|
/// </summary>
|
private static double CubicModel(double[] param, double x) =>
|
param[0] * x * x * x + param[1] * x * x + param[2] * x + param[3];
|
|
/// <summary>
|
/// 多项式导数（用于单调性约束）
|
/// </summary>
|
private static double Derivative(double[] param, double x) =>
|
3 * param[0] * x * x + 2 * param[1] * x + param[2];
|
|
|
/// <summary>
|
/// 构建目标函数（均方误差）
|
/// </summary>
|
private NonlinearObjectiveFunction BuildObjectiveFunction(double[] xNorm, double[] y) =>
|
new NonlinearObjectiveFunction(4, p =>
|
xNorm.Select((x, i) => Math.Pow(CubicModel(p, x) - y[i], 2)).Average()
|
);
|
|
|
|
/// <summary>
|
/// 多项式求值（霍纳法则）
|
/// </summary>
|
private static double EvaluatePolynomial(double[] coeffs, double x) =>
|
coeffs.Aggregate((sum, c) => sum * x + c);
|
|
/// <summary>
|
/// 生成等间距采样点
|
/// </summary>
|
private static double[] GenerateLinspace(double start, double end, int numPoints) =>
|
Enumerable.Range(0, numPoints)
|
.Select(i => start + i * (end - start) / (numPoints - 1))
|
.ToArray();
|
# endregion
|
|
# region 过渡处理
|
|
|
/// <summary>
|
/// 生成合并的 x 值
|
/// </summary>
|
public static double[] GenerateMergedX(double[] splineX, double[] measuredXValid, int numPoints)
|
{
|
// 生成等间距的点
|
double[] linspace = GenerateLinspace(measuredXValid.Min(), measuredXValid.Max(), numPoints);
|
// 拼接并去重排序
|
var concat = splineX.Concat(linspace).ToArray();
|
var merged = concat.Distinct().OrderBy(x => x).ToArray();
|
return merged;
|
}
|
|
# endregion
|
|
# region 辅助类
|
|
/// <summary>
|
/// Z分数计算器（用于异常值检测）
|
/// </summary>
|
private class ZScore
|
{
|
public readonly double[] Scores;
|
|
public ZScore(double[] data)
|
{
|
double mean = data.Average();
|
double std = data.StandardDeviation();
|
Scores = data.Select(x => (x - mean) / std).ToArray();
|
}
|
}
|
|
#endregion
|
|
|
|
public static (double[] CorrectedX, double[] CorrectedY) MergeCurves(
|
double[] originalX, double[] originalY,
|
double[] measuredX, double[] measuredY,int degree)
|
{
|
// 实测数据不足处理
|
if (measuredX.Length < 3)
|
{
|
Console.WriteLine("实测数据不足，使用原始曲线");
|
return (originalX, originalY);
|
}
|
|
// 多项式拟合实测数据（系数反转保持与numpy一致）
|
double[] polyCoeff = Fit.Polynomial(measuredX, measuredY, degree).Reverse().ToArray();
|
|
// 创建修正曲线副本
|
double[] correctedX = originalX.Select(x => x).ToArray();
|
double[] correctedY = originalY.Select(x => x).ToArray();
|
|
// 计算实测数据范围
|
double measuredMin = measuredX.Min();
|
double measuredMax = measuredX.Max();
|
|
// 核心区域修正
|
for (int i = 0; i < correctedX.Length; i++)
|
{
|
if (correctedX[i] >= measuredMin && correctedX[i] <= measuredMax)
|
correctedY[i] = Polyval(polyCoeff, correctedX[i]);
|
}
|
|
// 过渡区域处理
|
double transition = TransitionWidth;
|
|
// 左过渡区域
|
double[] originalPolyCoeff = Fit.Polynomial(originalX, originalY, 3).Reverse().ToArray();
|
double xLeftStart = Math.Max(measuredMin - transition, originalX.Min());
|
double yLeftStart = Polyval(originalPolyCoeff, xLeftStart);
|
double yLeftEnd = Polyval(polyCoeff, measuredMin);
|
|
for (int i = 0; i < correctedX.Length; i++)
|
{
|
double x = correctedX[i];
|
if (x >= xLeftStart && x < measuredMin)
|
correctedY[i] = LinearInterp(x, xLeftStart, yLeftStart, measuredMin, yLeftEnd);
|
}
|
|
// 右过渡区域
|
double xRightEnd = Math.Min(measuredMax + transition, originalX.Max());
|
double yRightStart = Polyval(polyCoeff, measuredMax);
|
double yRightEnd = Polyval(originalPolyCoeff, xRightEnd);
|
|
for (int i = 0; i < correctedX.Length; i++)
|
{
|
double x = correctedX[i];
|
if (x > measuredMax && x <= xRightEnd)
|
correctedY[i] = LinearInterp(x, measuredMax, yRightStart, xRightEnd, yRightEnd);
|
}
|
|
return (correctedX, correctedY);
|
}
|
|
// 多项式求值（支持降序系数）
|
private static double Polyval(double[] coefficients, double x)
|
{
|
double result = 0;
|
int degree = coefficients.Length - 1;
|
for (int i = 0; i < coefficients.Length; i++)
|
result += coefficients[i] * Math.Pow(x, degree - i);
|
return result;
|
}
|
|
// 线性插值
|
private static double LinearInterp(double x, double x0, double y0, double x1, double y1)
|
{
|
return x1 != x0
|
? y0 + (y1 - y0) * (x - x0) / (x1 - x0)
|
: y0;
|
}
|
|
}
|
|
|
|
|
|
|
public class PointClassifier
|
{
|
public class DenseArea
|
{
|
public (double Min, double Max) XRange { get; set; }
|
public (double Min, double Max) YRange { get; set; }
|
}
|
|
public (int[] Labels, List<DenseArea> DenseAreas) ClassifyPoints(
|
double[] x,
|
double[] y,
|
string method = "dbscan",
|
Dictionary<string, object> parameters = null)
|
{
|
|
if (x.Length != y.Length)
|
throw new ArgumentException("x 和 y 的长度必须相同");
|
|
// 数据标准化
|
double[,] data = new double[x.Length, 2];
|
for (int i = 0; i < x.Length; i++)
|
{
|
data[i, 0] = x[i];
|
data[i, 1] = y[i];
|
}
|
|
// 计算均值和标准差
|
double[] mean = new double[2];
|
double[] std = new double[2];
|
for (int j = 0; j < 2; j++)
|
{
|
var jList = Enumerable.Range(0, data.GetLength(0)).Select(i => data[i, j]).ToArray();
|
mean[j] = jList.Mean();
|
std[j] = Math.Sqrt(jList.Variance());
|
}
|
|
// 标准化数据
|
double[,] scaledData = new double[x.Length, 2];
|
for (int i = 0; i < x.Length; i++)
|
{
|
scaledData[i, 0] = (data[i, 0] - mean[0]) / std[0];
|
scaledData[i, 1] = (data[i, 1] - mean[1]) / std[1];
|
}
|
|
int[] labels = new int[x.Length];
|
List<DenseArea> denseAreas = new List<DenseArea>();
|
|
// 方法选择
|
if (method == "dbscan")
|
{
|
int min_samples = (int)(x.Length * 0.02);
|
min_samples = min_samples < 5 ? 5 : min_samples;
|
// DBSCAN参数
|
double eps = parameters != null && parameters.ContainsKey("eps") ? (double)parameters["eps"] : 0.3;
|
int minSamples = parameters != null && parameters.ContainsKey("min_samples")
|
? (int)parameters["min_samples"]
|
: min_samples;
|
|
// 计算距离矩阵
|
double[,] distanceMatrix = DistanceMatrix(scaledData);
|
|
// 执行 DBSCAN
|
var dbscan = new DBSCAN(eps, minSamples);
|
labels = dbscan.Run(scaledData);
|
}
|
else if (method == "kde")
|
{
|
// 核密度估计
|
int bandwidth = parameters != null && parameters.ContainsKey("bandwidth")
|
? (int)parameters["bandwidth"]
|
: 1;
|
double percentile = parameters != null && parameters.ContainsKey("percentile")
|
? (double)parameters["percentile"]
|
: 70.0;
|
|
// 计算密度
|
double[] density = new double[x.Length];
|
for (int i = 0; i < x.Length; i++)
|
{
|
double sum = 0;
|
for (int j = 0; j < x.Length; j++)
|
{
|
double distance = Math.Sqrt(
|
Math.Pow(scaledData[i, 0] - scaledData[j, 0], 2) +
|
Math.Pow(scaledData[i, 1] - scaledData[j, 1], 2));
|
sum += Math.Exp(-distance * distance / (2 * bandwidth * bandwidth));
|
}
|
density[i] = sum / (x.Length * 2 * Math.PI * bandwidth * bandwidth);
|
}
|
|
// 计算阈值
|
double threshold = GetPercentile(density, percentile);
|
|
// 标记密集区域
|
labels = density.Select(d => d >= threshold ? 1 : 0).ToArray();
|
}
|
else if (method == "grid")
|
{
|
// 动态网格划分
|
int xBins = parameters != null && parameters.ContainsKey("x_bins")
|
? (int)parameters["x_bins"]
|
: 30;
|
int yBins = parameters != null && parameters.ContainsKey("y_bins")
|
? (int)parameters["y_bins"]
|
: 20;
|
double percentile = parameters != null && parameters.ContainsKey("percentile")
|
? (double)parameters["percentile"]
|
: 75.0;
|
|
// 创建网格
|
double[] xBinsArray = CreateBins(x, xBins);
|
double[] yBinsArray = CreateBins(y, yBins);
|
|
// 计算网格密度
|
int[,] gridDensity = new int[xBins, yBins];
|
for (int i = 0; i < x.Length; i++)
|
{
|
int xBin = Array.BinarySearch(xBinsArray, x[i]);
|
int yBin = Array.BinarySearch(yBinsArray, y[i]);
|
if (xBin >= 0 && xBin < xBins && yBin >= 0 && yBin < yBins)
|
gridDensity[xBin, yBin]++;
|
}
|
|
// 计算阈值
|
double threshold = GetPercentile(gridDensity.Cast<double>().ToArray(), percentile);
|
|
// 标记密集网格
|
for (int i = 0; i < x.Length; i++)
|
{
|
int xBin = Array.BinarySearch(xBinsArray, x[i]);
|
int yBin = Array.BinarySearch(yBinsArray, y[i]);
|
if (xBin >= 0 && xBin < xBins && yBin >= 0 && yBin < yBins)
|
labels[i] = gridDensity[xBin, yBin] > threshold ? 1 : 0;
|
else
|
labels[i] = 0;
|
}
|
}
|
else
|
{
|
throw new ArgumentException("不支持的检测方法");
|
}
|
|
// 提取密集区域
|
var uniqueLabels = labels.Distinct();
|
foreach (var label in uniqueLabels)
|
{
|
if (label == -1) continue; // 排除噪声点
|
|
bool[] mask = labels.Select(l => l == label).ToArray();
|
if (mask.Count(m => m) > x.Length * 0.05) // 过滤小簇
|
{
|
double[,] clusterPoints = new double[mask.Count(m => m), 2];
|
int idx = 0;
|
for (int i = 0; i < mask.Length; i++)
|
{
|
if (mask[i])
|
{
|
clusterPoints[idx, 0] = data[i, 0];
|
clusterPoints[idx, 1] = data[i, 1];
|
idx++;
|
}
|
}
|
|
var xList = Enumerable.Range(0, clusterPoints.GetLength(0)).Select(i => clusterPoints[i, 0]).ToArray();
|
var yList = Enumerable.Range(0, clusterPoints.GetLength(0)).Select(i => clusterPoints[i, 1]).ToArray();
|
|
|
double xMin = xList.Min();
|
double xMax = xList.Max();
|
double yMin = yList.Min();
|
double yMax = yList.Max();
|
|
denseAreas.Add(new DenseArea
|
{
|
XRange = (xMin, xMax),
|
YRange = (yMin, yMax)
|
});
|
}
|
}
|
|
return (labels, denseAreas);
|
}
|
|
private double[] CreateBins(double[] data, int bins)
|
{
|
double min = data.Min();
|
double max = data.Max();
|
double step = (max - min) / bins;
|
return Enumerable.Range(0, bins + 1).Select(i => min + i * step).ToArray();
|
}
|
|
private double GetPercentile(double[] data, double percentile)
|
{
|
Array.Sort(data);
|
int index = (int)Math.Ceiling(percentile / 100.0 * data.Length) - 1;
|
return data[index];
|
}
|
|
private double[,] DistanceMatrix(double[,] data)
|
{
|
int n = data.GetLength(0);
|
double[,] matrix = new double[n, n];
|
for (int i = 0; i < n; i++)
|
{
|
for (int j = i + 1; j < n; j++)
|
{
|
double dx = data[i, 0] - data[j, 0];
|
double dy = data[i, 1] - data[j, 1];
|
matrix[i, j] = matrix[j, i] = Math.Sqrt(dx * dx + dy * dy);
|
}
|
}
|
return matrix;
|
}
|
|
}
|
|
|
public class DBSCAN
|
{
|
private readonly double eps;
|
private readonly int minPts;
|
|
public DBSCAN(double eps, int minPts)
|
{
|
this.eps = eps;
|
this.minPts = minPts;
|
}
|
|
public int[] Run(double[,] points)
|
{
|
int n = points.GetLength(0);
|
int[] labels = new int[n];
|
Array.Fill(labels, -2); // 初始化为未处理状态
|
|
int clusterId = -1;
|
|
for (int i = 0; i < n; i++)
|
{
|
if (labels[i] != -2) continue;
|
|
var neighbors = GetNeighbors(points, i);
|
|
if (neighbors.Count < minPts)
|
{
|
labels[i] = -1; // 标记为噪声点
|
continue;
|
}
|
|
clusterId++;
|
ExpandCluster(points, i, neighbors, clusterId, labels);
|
}
|
|
return labels;
|
}
|
|
private void ExpandCluster(double[,] points, int pointIndex, List<int> neighbors, int clusterId, int[] labels)
|
{
|
labels[pointIndex] = clusterId;
|
|
for (int i = 0; i < neighbors.Count; i++)
|
{
|
int neighborIndex = neighbors[i];
|
if (labels[neighborIndex] == -1)
|
{
|
labels[neighborIndex] = clusterId;
|
}
|
else if (labels[neighborIndex] == -2)
|
{
|
labels[neighborIndex] = clusterId;
|
var newNeighbors = GetNeighbors(points, neighborIndex);
|
if (newNeighbors.Count >= minPts)
|
{
|
neighbors.AddRange(newNeighbors);
|
}
|
}
|
}
|
}
|
|
private List<int> GetNeighbors(double[,] points, int pointIndex)
|
{
|
List<int> neighbors = new List<int>();
|
int n = points.GetLength(0);
|
for (int i = 0; i < n; i++)
|
{
|
double distance = EuclideanDistance(points, pointIndex, i);
|
if (distance < eps)
|
{
|
neighbors.Add(i);
|
}
|
}
|
return neighbors;
|
}
|
|
private double EuclideanDistance(double[,] points, int pointIndex1, int pointIndex2)
|
{
|
double dx = points[pointIndex1, 0] - points[pointIndex2, 0];
|
double dy = points[pointIndex1, 1] - points[pointIndex2, 1];
|
return Math.Sqrt(dx * dx + dy * dy);
|
}
|
}
|
|
}
|
