From 90b6c5914a6d46a304ad0f79120e60c95b950289 Mon Sep 17 00:00:00 2001 From: lvjunjie Date: Mon, 14 Sep 2026 14:30:26 +0800 Subject: [PATCH] =?UTF-8?q?feat(nmNum):=20=E5=A2=9E=E5=8A=A0=20LM=20?= =?UTF-8?q?=E5=9B=BA=E5=AE=9A=E5=AF=B9=E6=95=B0=E6=97=B6=E9=97=B4=E7=AA=97?= =?UTF-8?q?=E5=8F=A3=E4=B8=8E=E5=88=86=E6=AE=B5=E8=AF=AF=E5=B7=AE=E8=AF=8A?= =?UTF-8?q?=E6=96=AD?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../nmCalculation/nmCalculationAutoFitLM.h | 21 ++- .../nmCalculation/nmCalculationAutoFitLM.cpp | 150 +++++++++++++++++- 2 files changed, 162 insertions(+), 9 deletions(-) diff --git a/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h b/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h index f9fce51a..bcded5f3 100644 --- a/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h +++ b/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h @@ -11,6 +11,20 @@ #include "nmCalculation_global.h" +// 固定对数时间窗口的误差诊断。时间边界为不含重叠区的基础边界； +// rmsError 是局部加权均方根，energy 是对 total 平方的贡献。 +struct AutoFitTimeWindowLM { + double timeMin; + double timeMax; + double weightSum; + double rmsError; + double energy; + + AutoFitTimeWindowLM() + : timeMin(0.0), timeMax(0.0), weightSum(0.0), rmsError(0.0), energy(0.0) + {} +}; + // 双对数曲线误差分解。total 是 LM 候选接受和排序的唯一依据， // 其余诊断量用于有限差分灵敏度分析和信赖域选参。 struct AutoFitObjectiveBreakdownLM { @@ -19,6 +33,7 @@ struct AutoFitObjectiveBreakdownLM { double pressureLoss; double derivativeLoss; QVector residualVector; + QVector timeWindows; double verticalCommonBias; double verticalLoss; bool verticalReliable; @@ -140,13 +155,14 @@ private: const AutoFitObjectiveBreakdownLM* objectiveBreakdown = nullptr); QVector buildTraceParameterVector(const QVector& selectedParameters) const; void emitRunSummary(bool success, StopReasonLM finalReason); + void emitTimeWindowDiagnostics(const AutoFitObjectiveBreakdownLM& breakdown); bool validateParameters(const QVector& parameters) const; bool validateLogLogData(const QVector >& logLogData) const; bool validateInitialValues() const; bool validateSolverResult(const QVector >& result) const; double calculateLogLogCurveError(const QVector >& target, - const QVector >& result) const; + const QVector >& result); private: bool m_isRunning; @@ -172,6 +188,9 @@ private: QVector m_parameterUpper; QVector m_enabledParamIndices; QVector > m_targetLogLogData; + // 首次有效评价后固定，后续候选必须覆盖同一时间区间。 + double m_comparisonTimeMin; + double m_comparisonTimeMax; QString m_targetWellName; int m_maxIterations; diff --git a/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp b/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp index fee3e9bb..0cd0045f 100644 --- a/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp +++ b/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp @@ -27,6 +27,58 @@ #endif static const bool kAutoFitDiagnosticTraceEnabled = true; +static const int kAutoFitTimeWindowCount = 4; +static const double kAutoFitTimeWindowOverlapRatio = 0.20; + +static double autoFitTimeWindowWeight(double coordinate, int windowIndex) +{ + // 重叠区以基础边界为中心，总宽度为窗口宽度的 20%；相邻权重互补。 + const double width = 1.0 / kAutoFitTimeWindowCount; + const double overlap = width * kAutoFitTimeWindowOverlapRatio; + const double halfOverlap = overlap * 0.5; + const double pi = 3.14159265358979323846; + double weight = 1.0; + if(windowIndex > 0) { + double u = qBound(0.0, + (coordinate - windowIndex * width + halfOverlap) / overlap, + 1.0); + weight *= 0.5 * (1.0 - qCos(pi * u)); + } + if(windowIndex + 1 < kAutoFitTimeWindowCount) { + double u = qBound(0.0, + (coordinate - (windowIndex + 1) * width + halfOverlap) / overlap, + 1.0); + weight *= 0.5 * (1.0 + qCos(pi * u)); + } + return weight; +} + +static QVector calculateAutoFitTimeWindows( + const QVector& residualVector, double timeMin, double timeMax) +{ + // 残差前后两半分别为压力和导数，已包含各占一半及采样点数的归一化。 + const int pointCount = residualVector.size() / 2; + const double logMin = qLn(timeMin); + const double logSpan = qLn(timeMax) - logMin; + QVector windows(kAutoFitTimeWindowCount); + for(int k = 0; k < windows.size(); ++k) { + AutoFitTimeWindowLM& window = windows[k]; + window.timeMin = k == 0 ? timeMin + : qExp(logMin + logSpan * k / kAutoFitTimeWindowCount); + window.timeMax = k + 1 == windows.size() ? timeMax + : qExp(logMin + logSpan * (k + 1) / kAutoFitTimeWindowCount); + for(int i = 0; i < pointCount; ++i) { + const double weight = autoFitTimeWindowWeight( + static_cast(i) / (pointCount - 1), k); + window.weightSum += weight; + window.energy += weight * + (residualVector[i] * residualVector[i] + + residualVector[pointCount + i] * residualVector[pointCount + i]); + } + window.rmsError = qSqrt(window.energy * pointCount / window.weightSum); + } + return windows; +} static inline bool isFiniteNumber(double value) { @@ -466,6 +518,8 @@ nmCalculationAutoFitLM::nmCalculationAutoFitLM(QObject* parent) , m_isFinalizing(false) , m_currentIteration(0) , m_globalBestFitness(1e10) + , m_comparisonTimeMin(0.0) + , m_comparisonTimeMax(0.0) , m_maxIterations(100) , m_targetError(0.001) , m_totalEvaluations(0) @@ -549,6 +603,8 @@ void nmCalculationAutoFitLM::setTargetLogLogData(const QVector > // 目标曲线由界面层从目标井 history log-log 传入。 // 约定 targetData[0]=time，targetData[1]=pressure，targetData[2]=pressure derivative。 m_targetLogLogData = targetData; + m_comparisonTimeMin = 0.0; + m_comparisonTimeMax = 0.0; DEBUG_OUT(QString("Target LogLog data set: %1 arrays").arg(targetData.size())); if(targetData.size() >= 3) { @@ -641,6 +697,8 @@ void nmCalculationAutoFitLM::resetOptimizer() m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM(); m_userInitialLogLogData.clear(); m_userInitialObjectiveBreakdown = AutoFitObjectiveBreakdownLM(); + m_comparisonTimeMin = 0.0; + m_comparisonTimeMax = 0.0; m_currentIteration = 0; m_totalEvaluations = 0; m_successfulEvaluations = 0; @@ -753,6 +811,12 @@ void nmCalculationAutoFitLM::writeTraceHeader() << "late_trend_reliable" << "registration_ambiguous"; + for(int k = 0; k < kAutoFitTimeWindowCount; ++k) { + const QString prefix = QString("window_%1_").arg(k + 1); + cols << prefix + "time_min" << prefix + "time_max" + << prefix + "weight_sum" << prefix + "rms_error" << prefix + "energy"; + } + QTextStream out(&m_traceFile); out << cols.join(",") << "\n"; } @@ -789,7 +853,7 @@ void nmCalculationAutoFitLM::writeTraceMetaFile() QTextStream out(&metaFile); out << "{\n"; - out << " \"schema_version\": 2,\n"; + out << " \"schema_version\": 3,\n"; out << " \"trace_type\": \"finite_difference_lm_trust_region\",\n"; out << " \"run_id\": " << jsonEscape(m_traceRunId) << ",\n"; out << " \"created_at\": " @@ -805,6 +869,15 @@ void nmCalculationAutoFitLM::writeTraceMetaFile() out << " \"max_iterations\": " << m_maxIterations << ",\n"; out << " \"target_error\": " << jsonNumber(m_targetError) << "\n"; out << " },\n"; + // 区间尚未冻结时写 null；首次有效评价后重写元数据，保存实际窗口基准。 + out << " \"time_windows\": {\n"; + out << " \"count\": " << kAutoFitTimeWindowCount << ",\n"; + out << " \"overlap_ratio\": " << jsonNumber(kAutoFitTimeWindowOverlapRatio) << ",\n"; + out << " \"comparison_time_min\": " + << (m_comparisonTimeMin > 0.0 ? jsonNumber(m_comparisonTimeMin) : QString("null")) << ",\n"; + out << " \"comparison_time_max\": " + << (m_comparisonTimeMax > 0.0 ? jsonNumber(m_comparisonTimeMax) : QString("null")) << "\n"; + out << " },\n"; out << " \"parameters\": {\n"; out << " \"names\": " << jsonStringArray(parameterNames) << ",\n"; out << " \"enabled_indices\": " << jsonIntArray(m_enabledParamIndices) << ",\n"; @@ -875,6 +948,21 @@ void nmCalculationAutoFitLM::writeTraceRow( } } + // 无效候选也补齐窗口列，保证轨迹每行结构一致。 + for(int k = 0; k < kAutoFitTimeWindowCount; ++k) { + if(objectiveBreakdown && objectiveBreakdown->valid && + k < objectiveBreakdown->timeWindows.size()) { + const AutoFitTimeWindowLM& window = objectiveBreakdown->timeWindows[k]; + cols << traceNumber(window.timeMin) << traceNumber(window.timeMax) + << traceNumber(window.weightSum) << traceNumber(window.rmsError) + << traceNumber(window.energy); + } else { + for(int i = 0; i < 5; ++i) { + cols << QString(); + } + } + } + QTextStream out(&m_traceFile); out << cols.join(",") << "\n"; m_traceFile.flush(); @@ -893,6 +981,7 @@ void nmCalculationAutoFitLM::emitRunSummary(bool success, StopReasonLM finalReas .arg(m_totalEvaluations) .arg(m_successfulEvaluations) .arg(m_totalEvaluations - m_successfulEvaluations)); + emitTimeWindowDiagnostics(m_globalBestObjectiveBreakdown); if(!m_traceFilePath.isEmpty()) { emit logMessageGenerated(tr("Artifacts: trace=%1").arg(m_traceFilePath)); @@ -902,6 +991,26 @@ void nmCalculationAutoFitLM::emitRunSummary(bool success, StopReasonLM finalReas } } +void nmCalculationAutoFitLM::emitTimeWindowDiagnostics( + const AutoFitObjectiveBreakdownLM& breakdown) +{ + // 只汇报有效工作点；能量占比说明各时间段对全局误差的贡献。 + if(!breakdown.valid) { + return; + } + const double totalEnergy = breakdown.total * breakdown.total; + for(int k = 0; k < breakdown.timeWindows.size(); ++k) { + const AutoFitTimeWindowLM& window = breakdown.timeWindows[k]; + const double percentage = totalEnergy > 0.0 + ? 100.0 * window.energy / totalEnergy : 0.0; + emit logMessageGenerated( + tr("Time window %1 [%2, %3]: RMS=%4, energy=%5 (%6%)") + .arg(k + 1).arg(window.timeMin, 0, 'g', 6) + .arg(window.timeMax, 0, 'g', 6).arg(window.rmsError, 0, 'e', 4) + .arg(window.energy, 0, 'e', 4).arg(percentage, 0, 'f', 1)); + } +} + QVector nmCalculationAutoFitLM::buildTraceParameterVector(const QVector& selectedParameters) const { // 将 LM 内部使用的“启用参数向量”还原成完整 7 维参数向量。 @@ -1149,6 +1258,7 @@ bool nmCalculationAutoFitLM::startAutoFitting() m_globalBestLogLogData, 0, m_globalBestFitness); + emitTimeWindowDiagnostics(m_globalBestObjectiveBreakdown); } else { m_hasValidUserSolution = false; emit logMessageGenerated(tr("Initial solution evaluation failed")); @@ -1573,6 +1683,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() m_globalBestFitness = evaluation.fitness; m_globalBestObjectiveBreakdown = evaluation.breakdown; m_globalBestLogLogData = evaluation.curve; + emitTimeWindowDiagnostics(evaluation.breakdown); emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, m_currentIteration + 1, @@ -3071,10 +3182,10 @@ bool nmCalculationAutoFitLM::validateSolverResult(const QVector> double nmCalculationAutoFitLM::calculateLogLogCurveError( const QVector >& target, - const QVector >& result) const + const QVector >& result) { - // 主目标在目标与模拟曲线的公共时间范围内比较压力和导数残差；上下、左右 - // 和形状只负责诊断误差来源和选择参数，避免同一残差在 total 中被重复计算。 + // 首次有效评价后固定公共时间范围。窗口与上下、左右、形状只负责诊断， + // 主目标仍由完整压力和导数残差计算，避免窗口重叠造成重复计权。 // 整个计算过程均位于 log(time)-log(value) 坐标。 const double invalidLoss = 1.0e10; const double valueFloor = 1.0e-12; @@ -3200,13 +3311,21 @@ double nmCalculationAutoFitLM::calculateLogLogCurveError( return invalidLoss; } - // 与 PSO 保持一致：只在目标与模拟曲线的时间交集内比较，不再设置 - // 覆盖率门槛，也不对交集之外的首尾数据做外推。 - const double overlapMinX = qMax(targetMinX, resultMinX); - const double overlapMaxX = qMin(targetMaxX, resultMaxX); + // 仅首次有效评价使用交集建立基准；失败试算不能冻结区间，后续候选 + // 必须覆盖完整基准，不允许靠丢失首尾点缩小误差或改变窗口位置。 + const bool comparisonRangeFixed = m_comparisonTimeMin > 0.0; + const double overlapMinX = comparisonRangeFixed + ? m_comparisonTimeMin : qMax(targetMinX, resultMinX); + const double overlapMaxX = comparisonRangeFixed + ? m_comparisonTimeMax : qMin(targetMaxX, resultMaxX); if(overlapMinX >= overlapMaxX) { return invalidLoss; } + if(targetMinX > overlapMinX || targetMaxX < overlapMaxX || + resultMinX > overlapMinX || resultMaxX < overlapMaxX) { + DEBUG_OUT("Candidate does not cover the fixed LM comparison time range"); + return invalidLoss; + } QVector commonX(numPoints); QVector commonLogX(numPoints); @@ -3833,6 +3952,21 @@ double nmCalculationAutoFitLM::calculateLogLogCurveError( breakdown.valid = isFiniteNumber(breakdown.total) && breakdown.total >= 0.0; + if(breakdown.valid && breakdown.total < 1.0e9) { + breakdown.timeWindows = calculateAutoFitTimeWindows( + breakdown.residualVector, overlapMinX, overlapMaxX); + if(!comparisonRangeFixed) { + // 整个目标及诊断均有效后才提交基准，中点回退可重新建立区间。 + m_comparisonTimeMin = overlapMinX; + m_comparisonTimeMax = overlapMaxX; + writeTraceMetaFile(); + emit logMessageGenerated( + tr("Fixed LM comparison time range: [%1, %2]; %3 windows, %4% overlap") + .arg(overlapMinX, 0, 'g', 8).arg(overlapMaxX, 0, 'g', 8) + .arg(kAutoFitTimeWindowCount) + .arg(kAutoFitTimeWindowOverlapRatio * 100.0, 0, 'f', 0)); + } + } m_lastObjectiveBreakdown = breakdown; DEBUG_OUT(