|
using Accord.Math.Optimization;
|
using MathNet.Numerics;
|
using MathNet.Numerics.Interpolation;
|
using MathNet.Numerics.LinearAlgebra;
|
using MathNet.Numerics.LinearRegression;
|
using MathNet.Numerics.Statistics;
|
|
namespace IStation.Test
|
{
|
|
/// <summary>
|
/// 泵曲线数据融合校正器（用于合并样条曲线和实测数据并进行优化）
|
/// </summary>
|
public class PumpCurveDataFusionCorrectorHelper
|
{
|
# region 配置参数
|
// 异常值检测阈值（基于标准正态分布，2.5对应98.76%置信区间）
|
public const double ZScoreThreshold = 2.5;
|
// 过渡区域宽度（单位与输入数据相同，控制曲线衔接平滑度）
|
public const double TransitionWidth = 500;
|
// 多项式拟合阶数（三次平衡拟合能力与过拟合风险）
|
public const int PolynomialDegree = 3;
|
// 优化器最大迭代次数（Cobyla算法收敛阈值）
|
private const int OptimizationMaxIterations = 1000;
|
# endregion
|
|
/// <summary>
|
/// 主校正方法（入口函数）
|
/// </summary>
|
/// <param name="splineX">样条曲线X坐标</param>
|
/// <param name="splineY">样条曲线Y坐标</param>
|
/// <param name="measuredXAll">原始实测X坐标</param>
|
/// <param name="measuredYAll">原始实测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)
|
{
|
# region 步骤1：稳健回归与异常值过滤（最小二乘法）
|
// 执行简单线性回归获取基准线
|
(double intercept, double slope) = SimpleRegression.Fit(measuredXAll, measuredYAll);
|
|
// 计算预测值和残差
|
double[] predictedY = measuredXAll.Select(x => slope * x + intercept).ToArray();
|
double[] residuals = measuredYAll.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 = measuredXAll.Zip(measuredYAll, 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)
|
{
|
# region 步骤2：多项式拟合与初步融合
|
// 对有效实测数据进行多项式拟合
|
var polyFunc = Polynomial.Fit(measuredXValid, measuredYValid, polyDegree);
|
double[] polyCoeffs = polyFunc.Coefficients.Reverse().ToArray(); // 注意系数顺序转换
|
|
// 生成融合后的X坐标序列（合并样条和实测范围）
|
double[] mergedX = GenerateMergedX(splineX, measuredXValid, 200);
|
|
// 创建样条插值器并生成初始融合曲线
|
var splineInterpolator = LinearSpline.InterpolateSorted(splineX, splineY);
|
//double[] mergedY = mergedX.Select(x =>
|
// x < splineX[0] ? splineY[0] : // 左外推
|
// x > splineX.Last() ? splineY.Last() : // 右外推
|
// splineInterpolator.Interpolate(x) // 插值区域
|
//).ToArray();
|
|
double[] mergedY = mergedX.Select(x =>
|
//x < splineX[0] ? splineY[0] : // 左外推
|
//x > splineX.Last() ? splineY.Last() : // 右外推
|
splineInterpolator.Interpolate(x) // 插值区域
|
).ToArray();
|
# endregion
|
|
# region 步骤3：核心区域修正
|
// 在实测数据范围内应用多项式拟合结果
|
double minX = measuredXValid.Min(), maxX = measuredXValid.Max();
|
bool[] coreMask = mergedX.Select(x => x >= minX && x <= maxX).ToArray();
|
for (int i = 0; i < mergedX.Length; i++)
|
{
|
if (coreMask[i])
|
mergedY[i] = EvaluatePolynomial(polyCoeffs, mergedX[i]);
|
}
|
# endregion
|
|
# region 步骤4：过渡区域处理
|
// 对左右过渡区域进行平滑处理
|
var splinePoly = Polynomial.Fit(splineX, splineY, polyDegree);
|
double[] splineCoeffs = splinePoly.Coefficients.Reverse().ToArray();
|
|
// 左过渡区处理（样条->多项式）
|
ApplyTransition(
|
(minX - transitionWidth, minX),
|
splineCoeffs,
|
polyCoeffs,
|
mergedX,
|
ref mergedY
|
);
|
|
// 右过渡区处理（多项式->样条）
|
ApplyTransition(
|
(maxX, maxX + transitionWidth),
|
polyCoeffs,
|
splineCoeffs,
|
mergedX,
|
ref mergedY
|
);
|
# 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
|
|
// 生成最终优化曲线
|
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>
|
/// 应用过渡混合（在两个多项式之间渐变）
|
/// </summary>
|
private static void ApplyTransition(
|
(double start, double end) range,
|
double[] baseCoeffs,
|
double[] targetCoeffs,
|
double[] xValues,
|
ref double[] yValues)
|
{
|
for (int i = 0; i < xValues.Length; i++)
|
{
|
double x = xValues[i];
|
if (x < range.start || x > range.end) continue;
|
|
double weight = (x - range.start) / (range.end - range.start);
|
double baseVal = EvaluatePolynomial(baseCoeffs, x);
|
double targetVal = EvaluatePolynomial(targetCoeffs, x);
|
yValues[i] = (1 - weight) * baseVal + weight * targetVal;
|
}
|
}
|
|
/// <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
|
}
|
}
|
