YXS
/
nmWTAI-Platform
8
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

feat(nmNum): 增加 LM 固定对数时间窗口与分段误差诊断

Browse Source
feature/AutoFit-Optimize-20260914
lvjunjie 3 weeks ago
parent 4583e657cf
commit 90b6c5914a
2 changed files with 162 additions and 9 deletions
Show Stats Download Patch File Download Diff File

21
Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h
Unescape Escape View File

@ -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<double> residualVector;
QVector<AutoFitTimeWindowLM> timeWindows;
double verticalCommonBias;
double verticalLoss;
bool verticalReliable;
@ -140,13 +155,14 @@ private:
const AutoFitObjectiveBreakdownLM* objectiveBreakdown = nullptr);
QVector<double> buildTraceParameterVector(const QVector<double>& selectedParameters) const;
void emitRunSummary(bool success, StopReasonLM finalReason);
void emitTimeWindowDiagnostics(const AutoFitObjectiveBreakdownLM& breakdown);
bool validateParameters(const QVector<double>& parameters) const;
bool validateLogLogData(const QVector<QVector<double> >& logLogData) const;
bool validateInitialValues() const;
bool validateSolverResult(const QVector<QVector<double> >& result) const;
double calculateLogLogCurveError(const QVector<QVector<double> >& target,
const QVector<QVector<double> >& result) const;
const QVector<QVector<double> >& result);
private:
bool m_isRunning;
@ -172,6 +188,9 @@ private:
QVector<double> m_parameterUpper;
QVector<int> m_enabledParamIndices;
QVector<QVector<double> > m_targetLogLogData;
// 首次有效评价后固定，后续候选必须覆盖同一时间区间。
double m_comparisonTimeMin;
double m_comparisonTimeMax;
QString m_targetWellName;
int m_maxIterations;

150
Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp
Unescape Escape View File

@ -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<AutoFitTimeWindowLM> calculateAutoFitTimeWindows(
const QVector<double>& residualVector, double timeMin, double timeMax)
{
// 残差前后两半分别为压力和导数，已包含各占一半及采样点数的归一化。
const int pointCount = residualVector.size() / 2;
const double logMin = qLn(timeMin);
const double logSpan = qLn(timeMax) - logMin;
QVector<AutoFitTimeWindowLM> 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<double>(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<QVector<double> >
// 目标曲线由界面层从目标井 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<double> nmCalculationAutoFitLM::buildTraceParameterVector(const QVector<double>& 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<QVector<double>>
double nmCalculationAutoFitLM::calculateLogLogCurveError(
const QVector<QVector<double> >& target,
const QVector<QVector<double> >& result) const
const QVector<QVector<double> >& 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<double> commonX(numPoints);
QVector<double> 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(

Write Preview
Loading…
Cancel
Save