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

feat(nmNum): 新增 LM 高度与形状预调整的三阶段拟合

Browse Source 
- 优先根据压力与压力导数的上下偏差调整渗透率，对齐曲线高度
- 在固定对数时间网格计算双曲线斜率残差，优先改善形状并限制整体偏离
- 按形状改善与剩余预算切换阶段，第三阶段沿用仅按整体误差接受候选的 LM
- 联合缓存数值与形状灵敏度，复用阶段间模型并同步已接受参数及曲线
- 迭代候选求解失败后交由信赖域缩步，避免相同参数重复重试
- 补充三阶段拟合日志、跟踪诊断字段及中文翻译
feature/AutoFit-Optimize-20260914
lvjunjie 3 weeks ago
parent e027bf14d6
commit 57c3d613f6
4 changed files with 364 additions and 114 deletions
Show Stats Download Patch File Download Diff File

BIN
Bin/Config/Lang/cn/nmNum_cn.qm
View File

Binary file not shown.

72
Bin/Config/Lang/cn/nmNum_cn.ts
Unescape Escape View File

@ -977,8 +977,8 @@ Reason: %1</source>
<translation>连续 %1 次无有效改善，正在重建灵敏度模型进行确认</translation>
</message>
<message>
<source>Effective improvement threshold: max(%1, %2% of baseline error); %3 consecutive ineffective steps trigger convergence confirmation</source>
<translation>有效改善阈值：取 %1 与基准误差的 %2% 中较大值；连续 %3 次无有效改善后进行收敛确认</translation>
<source>Total-stage effective improvement threshold: max(%1, %2% of baseline error); %3 consecutive ineffective steps trigger convergence confirmation</source>
<translation>整体阶段有效改善阈值：取 %1 与基准误差的 %2% 中较大值；连续 %3 次无有效改善后进行收敛确认</translation>
</message>
<message>
<source>Sensitivity probe accepted: error reduced to %1</source>
@ -1116,6 +1116,74 @@ Reason: %1</source>
<source>Target well name is empty</source>
<translation>目标井名称为空</translation>
</message>
<message>
<source>LM stage 1: align curve height using permeability only (up to %1 evaluations).</source>
<translation>LM 阶段一：仅调整渗透率对齐曲线高度（最多求解 %1 次）。</translation>
</message>
<message>
<source>Permeability alignment: k=%1, height error=%2, shape error=%3, result=%4</source>
<translation>渗透率对齐：k=%1，上下误差=%2，形状误差=%3，结果=%4</translation>
</message>
<message>
<source>LM stage 2: optimize pressure and derivative shape; stop after 3 ineffective steps.</source>
<translation>LM 阶段二：优先调整压力和压力导数形状，连续 3 步无明显改善后切换。</translation>
</message>
<message>
<source>LM stage 3: original LM fitting; accept by total error only.</source>
<translation>LM 阶段三：按原有 LM 拟合，仅依据整体误差接受调整。</translation>
</message>
<message>
<source>Sensitivity probe accepted: total error=%1</source>
<translation>采用灵敏度试算点：整体误差=%1</translation>
</message>
<message>
<source>Shape stage ended: %1</source>
<translation>形状阶段结束：%1</translation>
</message>
<message>
<source>3 consecutive steps without effective shape improvement</source>
<translation>连续 3 步形状没有明显改善</translation>
</message>
<message>
<source>reserve remaining iterations and evaluations for total fitting</source>
<translation>为整体拟合保留剩余迭代和求解预算</translation>
</message>
<message>
<source>no valid shape sensitivity model</source>
<translation>没有有效的形状灵敏度模型</translation>
</message>
<message>
<source>no feasible shape descent step</source>
<translation>没有满足约束的形状下降步</translation>
</message>
<message>
<source>Permeability alignment ended: %1</source>
<translation>渗透率高度对齐结束：%1</translation>
</message>
<message>
<source>height evaluation budget reached</source>
<translation>达到高度调整求解预算</translation>
</message>
<message>
<source>stopped by user</source>
<translation>用户停止</translation>
</message>
<message>
<source>pressure and derivative height directions conflict</source>
<translation>压力与压力导数的上下调整方向冲突</translation>
</message>
<message>
<source>curve height is approximately aligned</source>
<translation>曲线高度已大致对齐</translation>
</message>
<message>
<source>2 consecutive steps without effective height improvement</source>
<translation>连续 2 步上下偏差没有明显改善</translation>
</message>
<message>
<source>permeability reached its bound</source>
<translation>渗透率已到达范围边界</translation>
</message>
</context>
<context>
<name>nmCalculationSolver</name>

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

@ -25,8 +25,8 @@ struct AutoFitTimeWindowLM {
{}
};
// 双对数曲线误差分解。total 是候选接受和排序的统一依据；分层模式使用全部目标点，
// 时间窗口和残差用于 Fisher 选参，其余诊断量用于解释曲线失配。
// 双对数曲线误差分解。预调整使用纵向偏差，形状阶段使用双曲线斜率残差，
// 整体阶段只使用 total 接受候选；时间窗口和残差用于 Fisher 选参。
struct AutoFitObjectiveBreakdownLM {
bool valid;
double total;
@ -46,6 +46,10 @@ struct AutoFitObjectiveBreakdownLM {
double horizontalLoss;
bool horizontalReliable;
bool registrationAmbiguous;
// 固定 log-time 网格上的双曲线斜率残差，独立于数值残差的采样层级。
QVector<double> shapeResiduals;
double pressureVerticalBias;
double derivativeVerticalBias;
double shapeLoss;
double lateDerivativeSlopeBias;
double lateDerivativeTrendLoss;
@ -66,6 +70,8 @@ struct AutoFitObjectiveBreakdownLM {
, horizontalLoss(std::numeric_limits<double>::quiet_NaN())
, horizontalReliable(false)
, registrationAmbiguous(false)
, pressureVerticalBias(0.0)
, derivativeVerticalBias(0.0)
, shapeLoss(std::numeric_limits<double>::quiet_NaN())
, lateDerivativeSlopeBias(std::numeric_limits<double>::quiet_NaN())
, lateDerivativeTrendLoss(std::numeric_limits<double>::quiet_NaN())
@ -134,7 +140,7 @@ private:
AutoFitObjectiveBreakdownLM* breakdown,
QVector<QVector<double> >* curve,
int* elapsedMs);
double evaluateFitness(const QVector<double>& parameters);
double evaluateFitness(const QVector<double>& parameters, bool retrySolver = true);
void applyParametersToDataManager(const QVector<double>& parameters);
void updateReservoirParameters(const QVector<double>& parameters);

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

@ -276,14 +276,6 @@ static double fromTrustRegionCoordinate(double coordinate,
return qMax(lower, qMin(upper, value));
}
enum TrustRegionErrorComponent
{
TRUST_REGION_VERTICAL_COMPONENT = 0,
TRUST_REGION_HORIZONTAL_COMPONENT,
TRUST_REGION_SHAPE_COMPONENT,
TRUST_REGION_TOTAL_COMPONENT
};
// 一次真实求解的完整快照。除了参数和总误差，还保存内部坐标、诊断分量和
// 双对数曲线，因此拒绝候选后可以完整恢复上一个已接受工作点。
struct TrustRegionEvaluation
@ -320,9 +312,21 @@ static bool trustRegionResidualsValid(
return false;
}
}
if(breakdown.shapeResiduals.size() != 146 || !isFiniteNumber(breakdown.shapeLoss)) return false;
for(int i = 0; i < breakdown.shapeResiduals.size(); ++i) {
if(!isFiniteNumber(breakdown.shapeResiduals[i])) return false;
}
return true;
}
// 两类残差在同一次真实评价中获得，联合缓存使阶段切换不必重算灵敏度。
static QVector<double> trustRegionFullResidual(const AutoFitObjectiveBreakdownLM& objective)
{
QVector<double> residual = objective.residualVector;
residual += objective.shapeResiduals;
return residual;
}
// 计算向量二范数的平方，避免在只比较能量或计算正规方程时反复开方。
static double trustRegionSquaredNorm(const QVector<double>& values)
{
@ -348,50 +352,6 @@ static double trustRegionDotProduct(const QVector<double>& left,
return sum;
}
// trace 和运行日志使用稳定的英文标识，便于现有离线脚本继续按字段筛选。
static QString trustRegionComponentName(int component)
{
if(component == TRUST_REGION_VERTICAL_COMPONENT) {
return "vertical";
}
if(component == TRUST_REGION_HORIZONTAL_COMPONENT) {
return "horizontal";
}
if(component == TRUST_REGION_SHAPE_COMPONENT) {
return "shape";
}
return "total";
}
// 三类损失量纲一致，直接选择当前最大的可靠分量；都很小时退回总残差梯度。
static int trustRegionDominantComponent(
const AutoFitObjectiveBreakdownLM& breakdown,
double diagnosisThreshold)
{
int component = TRUST_REGION_TOTAL_COMPONENT;
double largestLoss = diagnosisThreshold;
if(breakdown.verticalReliable &&
isFiniteNumber(breakdown.verticalLoss) &&
breakdown.verticalLoss > largestLoss) {
component = TRUST_REGION_VERTICAL_COMPONENT;
largestLoss = breakdown.verticalLoss;
}
if(breakdown.horizontalReliable &&
!breakdown.registrationAmbiguous &&
isFiniteNumber(breakdown.horizontalLoss) &&
breakdown.horizontalLoss > largestLoss) {
component = TRUST_REGION_HORIZONTAL_COMPONENT;
largestLoss = breakdown.horizontalLoss;
}
if(isFiniteNumber(breakdown.shapeLoss) &&
breakdown.shapeLoss > largestLoss) {
component = TRUST_REGION_SHAPE_COMPONENT;
}
return component;
}
// 求解选中参数对应的阻尼正规方程。参数最多七维，使用带部分主元的
// 高斯消元处理该小矩阵，并在主元退化时明确返回失败。
static bool solveTrustRegionLinearSystem(
@ -1043,7 +1003,8 @@ void nmCalculationAutoFitLM::writeTraceHeader()
}
cols << "sampling_mode" << "sampling_stride" << "sampling_points"
<< "full_target_points" << "layer_objective";
<< "full_target_points" << "layer_objective"
<< "pressure_vertical_bias" << "derivative_vertical_bias";
QTextStream out(&m_traceFile);
out << cols.join(",") << "\n";
}
@ -1080,7 +1041,11 @@ void nmCalculationAutoFitLM::writeTraceMetaFile()
QTextStream out(&metaFile);
out << "{\n";
out << " \"schema_version\": 4,\n";
out << " \"schema_version\": 9,\n";
out << " \"strategy\": \"permeability_height_then_shape_then_original_lm\",\n";
out << " \"shape_metric\": \"pressure_and_derivative_log_slopes_81_points_lag_8\",\n";
out << " \"shape_stage_total_tolerance\": \"max(0.02, 25% of stage entry total)\",\n";
out << " \"total_stage_shape_constraint\": false,\n";
out << " \"trace_type\": \"finite_difference_lm_trust_region\",\n";
out << " \"run_id\": " << jsonEscape(m_traceRunId) << ",\n";
out << " \"created_at\": "
@ -1196,7 +1161,9 @@ void nmCalculationAutoFitLM::writeTraceRow(
<< QString::number(objectiveBreakdown ? objectiveBreakdown->samplingStride : m_samplingStride)
<< (objectiveBreakdown ? QString::number(objectiveBreakdown->residualVector.size() / 2) : QString())
<< (objectiveBreakdown ? QString::number(objectiveBreakdown->fullPointCount) : QString())
<< (objectiveBreakdown ? traceNumber(objectiveBreakdown->layerError) : QString());
<< (objectiveBreakdown ? traceNumber(objectiveBreakdown->layerError) : QString())
<< (objectiveBreakdown ? traceNumber(objectiveBreakdown->pressureVerticalBias) : QString())
<< (objectiveBreakdown ? traceNumber(objectiveBreakdown->derivativeVerticalBias) : QString());
QTextStream out(&m_traceFile);
out << cols.join(",") << "\n";
m_traceFile.flush();
@ -1809,7 +1776,7 @@ bool nmCalculationAutoFitLM::evaluateTrustRegionPoint(
// 真实评价次数和耗时，同时要求当前层残差、误差结构和结果曲线均有效。
QTime timer;
timer.start();
*fitness = evaluateFitness(parameters);
*fitness = evaluateFitness(parameters, false);
*elapsedMs = timer.elapsed();
*breakdown = m_lastObjectiveBreakdown;
*curve = m_lastEvaluatedLogLogData;
@ -1845,7 +1812,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
const double minimumCoordinateStep = 1.0e-5;
const double minimumTrustRadius = 2.0e-3;
const double maximumTrustRadius = 0.30;
const double diagnosisThreshold = 1.0e-5;
// 误差下降至少达到绝对 1e-5 且相对当前有效基准 0.2% 才算有效改善。
// 更小的下降仍保留为最佳解，但不能反复清除停滞状态、延长拟合时间。
const double effectiveRelativeImprovement = 2.0e-3;
@ -1869,7 +1835,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
StopReasonLM stopReason = LM_MAX_ITERATIONS;
// jacobian 的行对应当前采样层的残差，列对应用户勾选的参数。
// Fisher 直接复用残差 Jacobian，上下/左右/形状诊断仅保留用于结果说明。
// Fisher 按当前阶段取数值或形状残差行，完整 Jacobian 在阶段之间复用。
QVector<QVector<double> > jacobian;
QVector<bool> jacobianColumnValid(dimensions, false);
@ -1909,7 +1875,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
m_lastEvaluatedLogLogData = evaluation.curve;
};
// 只有真实总误差更小的工作点才能发布为全局最优；曲线和诊断快照必须
// 各阶段只发布通过当前目标和约束检查的工作点；曲线和诊断快照必须
// 与参数同步更新，防止界面显示或最终精英保护使用错配的数据。
auto publishAcceptedPoint = [&](const TrustRegionEvaluation& evaluation) {
m_globalBestPosition = evaluation.parameters;
@ -1976,9 +1942,101 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
.arg(current.fitness, 0, 'e', 4)
.arg(maximumEvaluations));
// 高度预调整只改变用户勾选的渗透率。ln(k) 的初始变化由有符号高度差
// 给出，真实求解后用割线估计修正；拒绝时缩步，不让其他参数补偿高度。
const int permeabilityColumn = m_enabledParamIndices.indexOf(0);
const int heightBudget = qMin(4, qMax(0, (maximumEvaluations - m_totalEvaluations - dimensions - 2) / 6));
const double heightShapeLimit = current.breakdown.shapeLoss + qMax(0.01, 0.20 * current.breakdown.shapeLoss);
double heightBaseline = current.breakdown.verticalLoss;
int heightStagnation = 0;
double heightScale = 1.0;
double heightSlope = -1.0;
if(permeabilityColumn >= 0 && current.fitness >= m_targetError) {
emit logMessageGenerated(tr("LM stage 1: align curve height using permeability only (up to %1 evaluations).")
.arg(heightBudget));
for(int trial = 0; trial < heightBudget && processPauseAndStop(); ++trial) {
const double bias = current.breakdown.verticalCommonBias;
if(!current.breakdown.verticalReliable || qAbs(bias) <= 0.01 || heightStagnation >= 2) break;
const double k = current.parameters[permeabilityColumn];
if(k <= 0.0 || m_parameterLower[0] <= 0.0 || m_parameterUpper[0] <= m_parameterLower[0]) break;
const double logRange = qLn(m_parameterUpper[0]) - qLn(m_parameterLower[0]);
const double change = qBound(-0.30 * logRange, -bias / heightSlope * heightScale, 0.30 * logRange);
TrustRegionEvaluation candidate;
candidate.parameters = current.parameters;
candidate.parameters[permeabilityColumn] = qBound(m_parameterLower[0], k * qExp(change), m_parameterUpper[0]);
candidate.coordinates = coordinatesFromParameters(candidate.parameters);
const double actualStep = qLn(candidate.parameters[permeabilityColumn] / k);
if(qAbs(actualStep) < 1.0e-5) break;
candidate.valid = evaluateTrustRegionPoint(candidate.parameters, &candidate.fitness,
&candidate.breakdown, &candidate.curve, &candidate.elapsedMs);
// 两条曲线不能靠一高一低互相抵消，也不能为对齐高度严重破坏形状。
const bool accepted = candidate.valid && candidate.breakdown.verticalReliable &&
candidate.breakdown.verticalLoss < current.breakdown.verticalLoss &&
candidate.breakdown.shapeLoss <= heightShapeLimit;
if(candidate.valid) {
const double slope = (candidate.breakdown.verticalCommonBias - bias) / actualStep;
if(slope < -0.05 && isFiniteNumber(slope)) heightSlope = slope;
}
writeTraceRow(-1, 0, "permeability_height", candidate.parameters, candidate.fitness,
candidate.valid, candidate.elapsedMs, accepted ? "accepted_height" : "rejected_height",
candidate.valid ? &candidate.breakdown : nullptr);
if(accepted) {
current = candidate;
publishAcceptedPoint(current);
} else {
heightScale *= 0.5;
}
restoreEvaluationState(current);
if(heightBaseline - current.breakdown.verticalLoss >= qMax(1.0e-4, 0.01 * heightBaseline)) {
heightBaseline = current.breakdown.verticalLoss;
heightStagnation = 0;
} else ++heightStagnation;
emit logMessageGenerated(tr("Permeability alignment: k=%1, height error=%2, shape error=%3, result=%4")
.arg(candidate.parameters[permeabilityColumn], 0, 'g', 6)
.arg(candidate.breakdown.verticalLoss, 0, 'e', 4).arg(candidate.breakdown.shapeLoss, 0, 'e', 4)
.arg(accepted ? tr("accepted") : tr("rejected")));
}
QString heightReason = tr("height evaluation budget reached");
if(m_shouldStop) heightReason = tr("stopped by user");
else if(!current.breakdown.verticalReliable) heightReason = tr("pressure and derivative height directions conflict");
else if(current.breakdown.verticalLoss <= 0.01) heightReason = tr("curve height is approximately aligned");
else if(heightStagnation >= 2) heightReason = tr("2 consecutive steps without effective height improvement");
else if(current.parameters[permeabilityColumn] <= m_parameterLower[0] ||
current.parameters[permeabilityColumn] >= m_parameterUpper[0]) heightReason = tr("permeability reached its bound");
emit logMessageGenerated(tr("Permeability alignment ended: %1").arg(heightReason));
}
if(m_shouldStop) return LM_USER_STOPPED;
// 形状阶段不设达标阈值；连续三步无有效改善或耗用约三分之一预算即切换。
// 允许整体误差适度回升，但上限固定在高度对齐后，禁止逐步放宽导致漂移。
bool shapeStage = m_maxIterations >= 3 && current.fitness >= m_targetError;
const double shapeValueLimit = current.fitness + qMax(0.02, 0.25 * current.fitness);
const int shapeEvaluationDeadline = m_totalEvaluations + qMax(0,
(maximumEvaluations - m_totalEvaluations - dimensions - 2) / 3);
const int shapeIterationLimit = qMax(1, m_maxIterations / 3);
auto stageError = [&](const TrustRegionEvaluation& point) -> double {
return shapeStage ? point.breakdown.shapeLoss : point.fitness;
};
auto acceptable = [&](const TrustRegionEvaluation& point, const TrustRegionEvaluation& base) -> bool {
if(!point.valid) return false;
return shapeStage
? point.breakdown.shapeLoss < base.breakdown.shapeLoss && point.fitness <= shapeValueLimit
// 预调整结束后恢复原 LM：有效候选只按整体误差下降接受。
: point.fitness < base.fitness;
};
auto acceptPoint = [&](const TrustRegionEvaluation& point) {
current = point;
publishAcceptedPoint(current);
restoreEvaluationState(current);
};
emit logMessageGenerated(shapeStage
? tr("LM stage 2: optimize pressure and derivative shape; stop after 3 ineffective steps.")
: tr("LM stage 3: original LM fitting; accept by total error only."));
// 有效改善始终相对“上一次有效改善后的误差”累计判断，避免一连串微小
// 下降每次都清零计数；累计达到门槛后才开始新的有效改善基准。
double effectiveImprovementBaseline = current.fitness;
double effectiveImprovementBaseline = stageError(current);
int fullDataRejections = 0;
bool samplingRefreshFailed = false;
auto promoteSampling = [&](bool complete) -> bool {
@ -2014,7 +2072,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
attemptedWindows.fill(false);
globalFallbackAttempted = false;
fullDataRejections = 0;
effectiveImprovementBaseline = current.fitness;
effectiveImprovementBaseline = stageError(current);
emit logMessageGenerated(tr("LM sampling refined: %1 / %2 target points; full-target error: %3")
.arg(current.breakdown.residualVector.size() / 2)
.arg(current.breakdown.fullPointCount).arg(current.fitness, 0, 'e', 4));
@ -2028,11 +2086,31 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
.arg(current.breakdown.residualVector.size() / 2).arg(current.breakdown.fullPointCount)
: tr("LM sampling: fixed 80 points (original mode)."));
auto enterTotalStage = [&](const QString& reason) {
// 保留当前曲线和完整 J，只重置阶段停滞状态、阻尼及信赖半径。
rebuildRequested = jacobian.isEmpty() || consecutiveSolverFailures > 0;
modelRebuiltAtMinimumRadius = false;
shapeStage = false;
effectiveImprovementBaseline = current.fitness;
consecutiveIneffectiveSteps = 0;
consecutiveRejectedSteps = 0;
consecutiveSolverFailures = 0;
stagnationConfirmationRequested = false;
attemptedWindows.fill(false);
globalFallbackAttempted = false;
trustRadius = 0.12;
damping = 0.01;
emit logMessageGenerated(tr("Shape stage ended: %1").arg(reason));
emit logMessageGenerated(tr("LM stage 3: original LM fitting; accept by total error only."));
writeTraceRow(m_currentIteration, -1, "stage_switch", current.parameters, current.fitness,
true, 0, "shape_to_total", &current.breakdown);
};
auto registerEffectiveImprovement = [&](double fitness) -> bool {
const double requiredImprovement = qMax(
effectiveAbsoluteImprovement,
shapeStage ? 1.0e-4 : effectiveAbsoluteImprovement,
qAbs(effectiveImprovementBaseline) *
effectiveRelativeImprovement);
(shapeStage ? 0.01 : effectiveRelativeImprovement));
const double improvement = effectiveImprovementBaseline - fitness;
if(improvement < requiredImprovement) {
return false;
@ -2051,6 +2129,10 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
// 完整一轮仍无改善时沿用原有重建确认，避免某个难处理窗口提前终止拟合。
auto recordIneffectiveStep = [&]() -> bool {
++consecutiveIneffectiveSteps;
if(shapeStage) {
if(consecutiveIneffectiveSteps >= maximumIneffectiveSteps) enterTotalStage(tr("3 consecutive steps without effective shape improvement"));
return false;
}
if(!globalFallbackAttempted) {
return false;
}
@ -2075,7 +2157,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
};
emit logMessageGenerated(
tr("Effective improvement threshold: max(%1, %2% of baseline error); "
tr("Total-stage effective improvement threshold: max(%1, %2% of baseline error); "
"%3 consecutive ineffective steps trigger convergence confirmation")
.arg(effectiveAbsoluteImprovement, 0, 'e', 2)
.arg(effectiveRelativeImprovement * 100.0, 0, 'f', 2)
@ -2090,7 +2172,8 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
// 求解失败时才补算反方向，因此初次建模通常每个参数只增加一次真实求解。
auto rebuildSensitivity = [&]() -> bool {
const TrustRegionEvaluation base = current;
const int residualCount = base.breakdown.residualVector.size();
const QVector<double> baseResidual = trustRegionFullResidual(base.breakdown);
const int residualCount = baseResidual.size();
if(residualCount <= 0) {
return false;
}
@ -2129,8 +2212,16 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
? preferredSign : -preferredSign;
double availableRoom = direction > 0.0
? positiveRoom : negativeRoom;
double deltaMagnitude = qMin(
finiteDifferenceStep, availableRoom);
const int parameterIndex = m_enabledParamIndices[column];
const double lower = m_parameterLower[parameterIndex];
const double upper = m_parameterUpper[parameterIndex];
double localStep = finiteDifferenceStep;
if(shapeStage && upper > lower && parameterIndex == 1) {
localStep = qMin(localStep, 0.02 * qMax(0.1, qAbs(base.parameters[column])) / (upper - lower));
} else if(shapeStage && useTrustRegionLogScale(parameterIndex, lower, upper)) {
localStep = qMin(localStep, qLn(1.05) / (qLn(upper) - qLn(lower)));
}
double deltaMagnitude = qMin(localStep, availableRoom);
if(deltaMagnitude < minimumCoordinateStep) {
continue;
}
@ -2158,7 +2249,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
probe.fitness,
probe.valid,
probe.elapsedMs,
decision,
(shapeStage ? "shape_" : "total_") + decision,
probe.valid ? &probe.breakdown : nullptr);
if(!probe.valid) {
@ -2169,25 +2260,23 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
double delta = probe.coordinates[column] -
base.coordinates[column];
if(qAbs(delta) < minimumCoordinateStep ||
probe.breakdown.residualVector.size() != residualCount) {
trustRegionFullResidual(probe.breakdown).size() != residualCount) {
restoreEvaluationState(base);
continue;
}
// 第 column 列是固定残差向量相对内部参数坐标的有限差分：
// J[:,column] = (r_probe-r_base)/delta。
const QVector<double> probeResidual = trustRegionFullResidual(probe.breakdown);
for(int row = 0; row < residualCount; ++row) {
jacobian[row][column] =
(probe.breakdown.residualVector[row] -
base.breakdown.residualVector[row]) / delta;
jacobian[row][column] = (probeResidual[row] - baseResidual[row]) / delta;
}
jacobianColumnValid[column] = true;
columnBuilt = true;
if(probe.fitness < base.fitness &&
(!bestProbe.valid ||
probe.fitness < bestProbe.fitness)) {
if(acceptable(probe, base) &&
(!bestProbe.valid || stageError(probe) < stageError(bestProbe))) {
bestProbe = probe;
bestProbeColumn = column;
bestProbeDelta = delta;
@ -2214,12 +2303,10 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
acceptedStep[bestProbeColumn] = bestProbeDelta;
updateTrustRegionJacobian(
&jacobian,
base.breakdown.residualVector,
bestProbe.breakdown.residualVector,
trustRegionFullResidual(base.breakdown),
trustRegionFullResidual(bestProbe.breakdown),
acceptedStep);
current = bestProbe;
publishAcceptedPoint(current);
restoreEvaluationState(current);
acceptPoint(bestProbe);
writeTraceRow(m_currentIteration,
bestProbeColumn,
"trust_region_sensitivity_accept",
@ -2227,10 +2314,10 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
current.fitness,
true,
0,
"accepted_cached_probe",
shapeStage ? "shape_accepted_cached_probe" : "total_accepted_cached_probe",
&current.breakdown);
emit logMessageGenerated(
tr("Sensitivity probe accepted: error reduced to %1")
tr("Sensitivity probe accepted: total error=%1")
.arg(current.fitness, 0, 'e', 4));
} else {
restoreEvaluationState(current);
@ -2264,7 +2351,11 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
if(!processPauseAndStop()) {
break;
}
if(m_layeredSampling && m_samplingStride > 1) {
if(shapeStage && (iteration >= shapeIterationLimit ||
m_totalEvaluations + (rebuildRequested ? dimensions : 0) >= shapeEvaluationDeadline)) {
enterTotalStage(tr("reserve remaining iterations and evaluations for total fitting"));
}
if(!shapeStage && m_layeredSampling && m_samplingStride > 1) {
const bool reserveFinalBudget = maximumEvaluations - m_totalEvaluations <= 2 * (dimensions + 1);
const int layerDeadline = qMax(1, m_maxIterations * (m_samplingStride == 4 ? 1 : 2) / 3);
if(reserveFinalBudget || iteration >= layerDeadline) {
@ -2276,6 +2367,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
const bool confirmingStagnation =
stagnationConfirmationRequested;
if(!rebuildSensitivity()) {
if(shapeStage && !m_shouldStop) { enterTotalStage(tr("no valid shape sensitivity model")); rebuildRequested = true; continue; }
if(promoteSampling(false)) continue;
stopReason = m_shouldStop
? LM_USER_STOPPED
@ -2292,7 +2384,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
break;
}
const bool rebuildEffective =
registerEffectiveImprovement(current.fitness);
registerEffectiveImprovement(stageError(current));
if(confirmingStagnation && !rebuildEffective) {
if(promoteSampling(false)) continue;
emit logMessageGenerated(
@ -2304,12 +2396,22 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
}
// 每轮从最新 J 和当前残差重算窗口 Fisher，包含有限差分与割线更新的变化。
const QVector<double> objectiveResidual = shapeStage
? current.breakdown.shapeResiduals : current.breakdown.residualVector;
const int rowOffset = shapeStage ? current.breakdown.residualVector.size() : 0;
const QVector<QVector<double> > objectiveJacobian = jacobian.mid(rowOffset, objectiveResidual.size());
QVector<double> objectiveCoordinates = current.breakdown.sampleCoordinates;
if(shapeStage) {
objectiveCoordinates.clear();
for(int i = 0; i < 73; ++i) objectiveCoordinates.append((i + 4.0) / 80.0);
}
const QVector<AutoFitTimeWindowLM> objectiveWindows = shapeStage
? calculateAutoFitTimeWindows(objectiveResidual, m_comparisonTimeMin, m_comparisonTimeMax,
objectiveCoordinates, autoFitLogTimeWeights(objectiveCoordinates))
: current.breakdown.timeWindows;
const QVector<TrustRegionFisher> information = buildTrustRegionFisher(
jacobian, current.breakdown.residualVector, jacobianColumnValid,
current.breakdown.sampleCoordinates);
objectiveJacobian, objectiveResidual, jacobianColumnValid, objectiveCoordinates);
const TrustRegionFisher& global = information[kAutoFitTimeWindowCount];
const int dominantComponent = trustRegionDominantComponent(
current.breakdown, diagnosisThreshold); // 仅用于现有诊断日志。
QVector<int> selectedColumns;
QVector<double> coordinateStep;
double predictedReduction = 0.0;
@ -2319,7 +2421,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
// 全部窗口都处理不动后，再用全局 Fisher 作一次补充选参。
while(selectedColumns.isEmpty()) {
selectedWindow = nextTrustRegionWindow(
current.breakdown.timeWindows, attemptedWindows);
objectiveWindows, attemptedWindows);
if(selectedWindow >= 0) {
attemptedWindows[selectedWindow] = true;
} else {
@ -2347,6 +2449,38 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
&step, &reduction)) {
continue;
}
// 仅形状预调整限制整体偏离；第三阶段直接使用原 LM 步长，
// 不预测形状上限，也不因形状变化缩短候选步长。
if(shapeStage) {
auto predictedGuard = [&](double scale) -> double {
double energy = 0.0;
for(int row = 0; row < current.breakdown.residualVector.size(); ++row) {
double value = current.breakdown.residualVector[row];
for(int column = 0; column < dimensions; ++column)
value += scale * jacobian[row][column] * step[column];
energy += value * value;
}
return qSqrt(energy);
};
const double currentGuard = predictedGuard(0.0);
const double predictedLimit = currentGuard + 0.9 * qMax(0.0, shapeValueLimit - currentGuard);
if(predictedGuard(1.0) > predictedLimit) {
double low = 0.0, high = 1.0;
for(int search = 0; search < 32; ++search) {
const double middle = 0.5 * (low + high);
if(predictedGuard(middle) <= predictedLimit) low = middle;
else high = middle;
}
const double stepScale = low;
for(int column = 0; column < dimensions; ++column) step[column] *= stepScale;
if(qSqrt(trustRegionSquaredNorm(step)) < minimumCoordinateStep) continue;
reduction = -trustRegionDotProduct(global.gradient, step);
for(int a = 0; a < dimensions; ++a)
for(int b = 0; b < dimensions; ++b)
reduction -= 0.5 * step[a] * global.matrix[a][b] * step[b];
if(!isFiniteNumber(reduction) || reduction <= 1.0e-14) continue;
}
}
const double tolerance = 1.0e-12 * qMax(predictedReduction, reduction);
if(selectedColumns.isEmpty() || reduction > predictedReduction + tolerance ||
(qAbs(reduction - predictedReduction) <= tolerance &&
@ -2363,6 +2497,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
// 只有窗口候选与全局回退均无方向，才收缩半径并进入原有重建/收敛处理。
if(selectedColumns.isEmpty()) {
if(shapeStage) { enterTotalStage(tr("no feasible shape descent step")); continue; }
if(trustRadius <= minimumTrustRadius * 1.01 &&
modelRebuiltAtMinimumRadius) {
if(promoteSampling(false)) continue;
@ -2390,7 +2525,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
for(int i = 0; i < selectedColumns.size(); ++i) {
selectedParameterIndices << QString::number(m_enabledParamIndices[selectedColumns[i]]);
}
const QString selectionName = (selectedWindow >= 0
const QString selectionName = QString(shapeStage ? "shape_" : "total_") + (selectedWindow >= 0
? QString("window_%1").arg(selectedWindow + 1) : QString("global")) +
"_params_" + selectedParameterIndices.join("_");
@ -2438,35 +2573,31 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
consecutiveSolverFailures = 0;
// 有效候选即使最终被拒绝，也提供了一条真实割线，可用于修正下一轮
// 局部模型；是否成为新工作点仍只由下面的 total 严格比较决定。
// 局部模型；是否成为新工作点由当前阶段的目标和约束共同决定。
const AutoFitObjectiveBreakdownLM oldBreakdown = current.breakdown;
updateTrustRegionJacobian(
&jacobian,
oldBreakdown.residualVector,
candidate.breakdown.residualVector,
trustRegionFullResidual(oldBreakdown),
trustRegionFullResidual(candidate.breakdown),
coordinateStep);
// reductionRatio 衡量局部线性模型的可信度：接近 1 表示预测准确；
// 值较小表示虽然可能下降，但模型低估了非线性，需要收紧下一步。
double actualReduction = m_layeredSampling
? 0.5 * (trustRegionSquaredNorm(current.breakdown.residualVector) -
trustRegionSquaredNorm(candidate.breakdown.residualVector))
: 0.5 * (current.fitness * current.fitness - candidate.fitness * candidate.fitness);
double actualReduction = 0.5 * (trustRegionSquaredNorm(objectiveResidual) -
trustRegionSquaredNorm(shapeStage ? candidate.breakdown.shapeResiduals : candidate.breakdown.residualVector));
double reductionRatio = actualReduction / predictedReduction;
bool accepted = candidate.fitness < current.fitness;
bool accepted = acceptable(candidate, current);
// 粗层认为下降而完整数据不认可时累计，连续两次就提前加密。
if(m_layeredSampling && m_samplingStride > 1 && !accepted && actualReduction > 0.0) {
if(!shapeStage && m_layeredSampling && m_samplingStride > 1 && !accepted && actualReduction > 0.0) {
++fullDataRejections;
} else {
fullDataRejections = 0;
}
QString componentName = trustRegionComponentName(dominantComponent);
QString componentName = shapeStage ? "shape" : "total";
if(accepted) {
// 真实总误差下降后才正式替换 current，并同步发布参数、曲线和诊断。
// 当前阶段接受候选后同步发布参数、曲线和诊断。
// 模型预测可靠时减小阻尼并可扩大半径，预测较差时保守收缩。
current = candidate;
publishAcceptedPoint(current);
restoreEvaluationState(current);
acceptPoint(candidate);
++acceptedSinceRebuild;
movementSinceRebuild += stepNorm;
consecutiveRejectedSteps = 0;
@ -2505,7 +2636,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
// 候选只要更优就继续作为 current 保存；是否足以解除停滞，则统一
// 相对上一次有效改善基准判断。拒绝和微小改善都会累计无效次数。
const bool effectiveImprovement =
registerEffectiveImprovement(current.fitness);
registerEffectiveImprovement(stageError(current));
if(!effectiveImprovement && recordIneffectiveStep()) {
stopReason = LM_LOCAL_OPTIMUM;
}
@ -2587,7 +2718,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
return LM_MAX_ITERATIONS;
}
double nmCalculationAutoFitLM::evaluateFitness(const QVector<double>& parameters)
double nmCalculationAutoFitLM::evaluateFitness(const QVector<double>& parameters, bool retrySolver)
{
// LM 候选评价函数，也是自动拟合最核心的闭环：
// 1. 校验候选参数是否在用户设置的上下界和基本物理范围内；
@ -2708,9 +2839,9 @@ double nmCalculationAutoFitLM::evaluateFitness(const QVector<double>& parameters
}
// 4. 运行求解器。真实求解器偶发失败时允许重试，避免一次 DLL 调用异常
// 直接让整个粒子评价失败。
// 直接让初始评价失败；迭代候选不重复相同参数，交给信赖域缩步。
QVector<QVector<double>> solverResult;
const int maxRetries = 2;
const int maxRetries = retrySolver ? 2 : 0;
bool solverSuccess = false;
for(int retry = 0; retry <= maxRetries; ++retry) {
@ -3406,6 +3537,49 @@ double nmCalculationAutoFitLM::calculateLogLogCurveError(
return invalidLoss;
}
auto populateStageMetrics = [&](AutoFitObjectiveBreakdownLM* objective) -> bool {
// 81 个固定对数时间点只用于插值评价，不增加求解点数或 DLL 调用。
// 斜率用跨度为全时域 10% 的差分，避免相邻点噪声；两条曲线等权。
const int count = 81;
const int lag = 8;
const double span = qLn(overlapMaxX) - qLn(overlapMinX);
QVector<double> residuals[2];
double biases[2] = {0.0, 0.0};
for(int component = 0; component < 2; ++component) {
for(int i = 0; i < count; ++i) {
const double time = i == 0 ? overlapMinX : (i == count - 1 ? overlapMaxX
: qExp(qLn(overlapMinX) + span * i / (count - 1)));
double targetValue = 0.0, resultValue = 0.0;
if(!interpolateLogValue(component == 0 ? targetPressure : targetDerivative,
time, &targetValue) ||
!interpolateLogValue(component == 0 ? resultPressure : resultDerivative,
time, &resultValue)) return false;
const double residual = resultValue - targetValue;
residuals[component].append(residual);
// 对数时间梯形权重等价于互补窗口加权求和，避免密集段主导高度。
biases[component] += residual * ((i == 0 || i == count - 1) ? 0.5 : 1.0)
/ (count - 1);
}
}
objective->pressureVerticalBias = biases[0];
objective->derivativeVerticalBias = biases[1];
objective->verticalCommonBias = 0.5 * (biases[0] + biases[1]);
objective->verticalLoss = qAbs(objective->verticalCommonBias);
objective->verticalReliable = !(biases[0] * biases[1] < 0.0 &&
qMin(qAbs(biases[0]), qAbs(biases[1])) > 0.01);
objective->shapeResiduals.clear();
const double scale = qSqrt(0.5 / (count - lag));
for(int component = 0; component < 2; ++component) {
for(int i = 0; i < count - lag; ++i) {
objective->shapeResiduals.append(scale *
(residuals[component][i + lag] - residuals[component][i]) /
(span * lag / (count - 1)));
}
}
objective->shapeLoss = qSqrt(trustRegionSquaredNorm(objective->shapeResiduals));
return isFiniteNumber(objective->shapeLoss);
};
if(m_layeredSampling) {
// 完整基准只取公共范围内的目标原始时间点；模拟点数不改变评价标准。
QVector<double> times;
@ -3469,6 +3643,7 @@ double nmCalculationAutoFitLM::calculateLogLogCurveError(
}
breakdown.layerError = qSqrt(trustRegionSquaredNorm(breakdown.residualVector));
breakdown.valid = isFiniteNumber(breakdown.total) && breakdown.total < 1.0e9;
if(!populateStageMetrics(&breakdown)) return invalidLoss;
if(!trustRegionResidualsValid(breakdown)) {
return invalidLoss;
}
@ -4104,7 +4279,7 @@ double nmCalculationAutoFitLM::calculateLogLogCurveError(
}
// LM 总目标等于固定残差向量的二范数；压力和导数各占一半能量。
// 上下、左右和形状分量不参与候选排序与接受。
// total 的定义不变；阶段形状和高度指标在下方以固定网格单独计算。
breakdown.total = qSqrt(
0.5 * breakdown.pressureLoss * breakdown.pressureLoss +
0.5 * breakdown.derivativeLoss * breakdown.derivativeLoss);
@ -4122,6 +4297,7 @@ double nmCalculationAutoFitLM::calculateLogLogCurveError(
writeTraceMetaFile();
}
}
if(!populateStageMetrics(&breakdown)) return invalidLoss;
m_lastObjectiveBreakdown = breakdown;
DEBUG_OUT(

Write Preview
Loading…
Cancel
Save