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

feat(nmNum): 优化 LM 前期形态引导与分阶段参数调整

Browse Source 
- 按初始曲线间距固定井储和表皮方向，小步起调，形状变差时同向扩步
- 前期采用双曲线斜率形状误差，放大第二阶段灵敏度试探步长
- 渗透率预调整仅看上下偏移，表皮搜索下限设为 0
- 补充间距与方向诊断日志，保留第三阶段原有逻辑
feature/AutoFit-Optimize-20260914
lvjunjie 2 weeks ago
parent 7e5255203e
commit 68c79378bb
3 changed files with 217 additions and 43 deletions
Show Stats Download Patch File Download Diff File

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

@ -48,8 +48,9 @@ struct AutoFitObjectiveBreakdownLM {
bool registrationAmbiguous;
// 固定 log-time 网格上的双曲线斜率残差，独立于数值残差的采样层级。
QVector<double> shapeResiduals;
// 前期数值残差为 81 点，平行程度残差为 73 个斜率区间，均含第一窗口权重。
// 平行误差比较模拟与目标各自的压力—导数斜率差，不要求模拟自身斜率差为零。
// 前期数值残差为 81 点，形状残差为 73 个斜率区间，均含第一窗口权重。
// earlyParallelResiduals/Loss 沿用历史字段名，现为压力、导数分别匹配目标的形状误差。
// earlyParallelBias 仍记录相对斜率偏差，仅供诊断。
QVector<double> earlyValueResiduals;
QVector<double> earlyParallelResiduals;
double earlyValueLoss;

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

@ -328,7 +328,7 @@ static bool trustRegionResidualsValid(
return true;
}
// 数值、形状和前期平行程度残差共用一次求解，联合缓存供各子阶段复用。
// 数值、整体形状和前期形状残差共用一次求解，联合缓存供各子阶段复用。
static QVector<double> trustRegionFullResidual(const AutoFitObjectiveBreakdownLM& objective)
{
QVector<double> residual = objective.residualVector;
@ -348,6 +348,46 @@ static double trustRegionSquaredNorm(const QVector<double>& values)
return sum;
}
// 数值残差已含 sqrt(0.5 * 窗口权重)，相减并换底后得到 log10 间距的加权 RMS。
// 不减去平均间距，保留“两条曲线整体离得太远”的信息；共同纵移会自然抵消。
static double trustRegionEarlyGapLoss(const AutoFitObjectiveBreakdownLM& objective)
{
const int count = objective.earlyValueResiduals.size() / 2;
double energy = 0.0;
for(int i = 0; i < count; ++i) {
const double difference = objective.earlyValueResiduals[i] - objective.earlyValueResiduals[count + i];
energy += 2.0 * difference * difference;
}
return qSqrt(energy) / qLn(10.0);
}
// 前期带符号的平均间距差：正值表示模拟比目标分得更开，负值表示更靠近。
// 复用现有加权数值残差，恢复第一窗口的归一化权重，不用 RMS 推断偏差符号。
static double trustRegionEarlyGapBias(const AutoFitObjectiveBreakdownLM& objective)
{
const int count = objective.earlyValueResiduals.size() / 2;
QVector<double> weights(count, 0.0);
double weightSum = 0.0;
for(int i = 0; i < count; ++i) {
weights[i] = autoFitTimeWindowWeight(static_cast<double>(i) / (count - 1), 0) *
((i == 0 || i == count - 1) ? 0.5 : 1.0);
weightSum += weights[i];
}
double bias = 0.0;
for(int i = 0; i < count; ++i) {
bias += qSqrt(2.0 * weights[i] / weightSum) *
(objective.earlyValueResiduals[i] - objective.earlyValueResiduals[count + i]);
}
return bias / qLn(10.0);
}
// 按用户指定的形态规则定向：初始模拟间距大则增大井储/表皮，小则减小。
static double trustRegionEarlyGapDirection(double initialGapBias)
{
if(!isFiniteNumber(initialGapBias) || qAbs(initialGapBias) <= 1.0e-12) return 0.0;
return initialGapBias > 0.0 ? 1.0 : -1.0;
}
// 第一窗口的形状能量与局部 Fisher 使用相同的中心时间和重叠权重，
// 压力、导数同时参与；不重新插值或调用求解器。
static double trustRegionEarlyShapeEnergy(const QVector<double>& residual)
@ -361,7 +401,7 @@ static double trustRegionEarlyShapeEnergy(const QVector<double>& residual)
return energy;
}
// 比较模拟与目标各自的压力—导数斜率差；误差为零表示两组曲线的相对走势一致。
// 前期复用整体形状的双曲线斜率残差，只按第一窗口重新加权并归一化。
// 单独平移任一曲线不改变此指标；前期数值误差仅作诊断，不参与井储、表皮验收。
static void populateEarlyWellboreMetrics(AutoFitObjectiveBreakdownLM* objective,
const QVector<double> targetLogs[2], const QVector<double> resultLogs[2], double logTimeSpan)
@ -404,10 +444,11 @@ static void populateEarlyWellboreMetrics(AutoFitObjectiveBreakdownLM* objective,
const double relativeSlopeError = resultSlopeDifference - targetSlopeDifference;
const double weight = parallelWeights[i] / parallelWeightSum;
objective->earlyParallelBias += weight * relativeSlopeError;
// 两半各占一半能量以兼容窗口 Fisher；保留相对目标的斜率差残差符号。
const double residual = qSqrt(0.5 * weight) * relativeSlopeError;
objective->earlyParallelResiduals[i] = residual;
objective->earlyParallelResiduals[slopeCount + i] = residual;
// 整体残差已含 sqrt(0.5 / slopeCount)，换成第一窗口权重后仍保持双曲线等权。
// 压力、导数分别匹配目标；相对斜率偏差只保留为诊断，不参与前期验收。
const double scale = qSqrt(slopeCount * weight);
objective->earlyParallelResiduals[i] = scale * objective->shapeResiduals[i];
objective->earlyParallelResiduals[slopeCount + i] = scale * objective->shapeResiduals[slopeCount + i];
}
objective->earlyParallelLoss = qSqrt(trustRegionSquaredNorm(objective->earlyParallelResiduals));
}
@ -719,10 +760,31 @@ static bool buildTrustRegionFisherStep(
return projectAndPredict();
}
// 放大已被真实结果验证可靠的联合方向，仍受最大半径和参数边界限制。
// 前期每次只调一个参数：间距锁定符号，形状 LM 提供幅度和后续验收依据。
static bool buildEarlyGapGuidedStep(const TrustRegionFisher& shapeInformation,
const QVector<int>& selected, const QVector<double>& coordinates, double direction,
double damping, double trustRadius, double minimumStep,
QVector<double>* step, double* predictedReduction)
{
if(selected.size() != 1 || direction == 0.0) return false;
const int column = selected[0];
TrustRegionFisher guided = shapeInformation;
guided.gradient[column] = -direction * qAbs(shapeInformation.gradient[column]);
if(!buildTrustRegionFisherStep(guided, selected, coordinates, damping,
trustRadius, minimumStep, step, predictedReduction)) return false;
// 改方向后的预测必须用原形状模型重算，不能把上坡伪装成预测下降。
const double delta = (*step)[column];
*predictedReduction = -shapeInformation.gradient[column] * delta -
0.5 * shapeInformation.matrix[column][column] * delta * delta;
return isFiniteNumber(*predictedReduction);
}
// 放大当前方向，仍受传入的最大半径和参数边界限制。
// 普通形状步要求预测下降；前期初始探路可跨过预测上坡区，最终仍按真实形状验收。
static bool buildExpandedTrustRegionStep(const TrustRegionFisher& information,
const QVector<double>& coordinates, const QVector<double>& originalStep,
double maximumRadius, QVector<double>* expandedStep, double* prediction)
double maximumRadius, QVector<double>* expandedStep, double* prediction,
bool requirePredictedDescent = true)
{
const double norm = qSqrt(trustRegionSquaredNorm(originalStep));
if(norm <= 1.0e-12) return false;
@ -739,7 +801,7 @@ static bool buildExpandedTrustRegionStep(const TrustRegionFisher& information,
for(int i = 0; i < expandedStep->size(); ++i)
for(int j = 0; j < expandedStep->size(); ++j)
*prediction -= 0.5 * (*expandedStep)[i] * information.matrix[i][j] * (*expandedStep)[j];
return isFiniteNumber(*prediction) && *prediction > 1.0e-14;
return isFiniteNumber(*prediction) && (!requirePredictedDescent || *prediction > 1.0e-14);
}
// 每次得到有效真实候选后，使用满足最新割线条件的秩一修正更新完整残差
@ -1113,7 +1175,7 @@ void nmCalculationAutoFitLM::writeTraceHeader()
cols << "sampling_mode" << "sampling_stride" << "sampling_points"
<< "full_target_points" << "layer_objective"
<< "pressure_vertical_bias" << "derivative_vertical_bias" << "first_window_shape_loss"
<< "early_value_loss" << "early_parallel_loss" << "early_parallel_bias";
<< "early_value_loss" << "early_parallel_loss" << "early_parallel_bias" << "early_gap_loss" << "early_gap_bias";
QTextStream out(&m_traceFile);
out << cols.join(",") << "\n";
}
@ -1150,12 +1212,18 @@ void nmCalculationAutoFitLM::writeTraceMetaFile()
QTextStream out(&metaFile);
out << "{\n";
out << " \"schema_version\": 21,\n";
out << " \"schema_version\": 29,\n";
out << " \"strategy\": \"permeability_height_then_shape_then_joint_lm\",\n";
out << " \"height_acceptance\": \"reliable_vertical_loss_decrease; no_shape_loss_constraint\",\n";
out << " \"shape_priority\": \"storage_then_skin_then_shape_without_wellbore_then_optional_wellbore_recheck\",\n";
out << " \"shape_metric\": \"pressure_and_derivative_log_slopes_81_points_lag_8\",\n";
out << " \"early_parallel_metric\": \"pressure_derivative_log_slope_difference_vs_target_81_points_lag_8\",\n";
out << " \"early_wellbore_acceptance\": \"relative_slope_matching_decrease_only\",\n";
out << " \"early_parallel_metric\": \"pressure_and_derivative_log_slopes_first_window_81_points_lag_8\",\n";
out << " \"early_gap_metric\": \"weighted_rms_log10_pressure_derivative_ratio_error_first_window_81_points\",\n";
out << " \"early_gap_bias_metric\": \"weighted_mean_signed_log10_gap_simulation_minus_target_first_window\",\n";
out << " \"early_direction_policy\": \"initial_gap_bias_positive_increase_negative_decrease; direction_shared_and_locked_for_storage_and_skin; recompute_only_at_wellbore_recheck_entry\",\n";
out << " \"early_initial_step_policy\": \"after_gap_direction_selection; storage_ratio_cap_1.05; skin_local_scale_fraction_0.02; no_cached_probe_acceptance_before_first_shape_improvement\",\n";
out << " \"early_initial_expansion_policy\": \"before_first_accepted_shape_improvement; if_shape_worsens_double_coordinate_step_from_unchanged_base_up_to_parameter_bound; retain_base_on_failure\",\n";
out << " \"early_wellbore_acceptance\": \"locked_gap_direction_and_early_shape_decrease\",\n";
out << " \"early_wellbore_total_constraint\": false,\n";
out << " \"shape_stage_total_tolerance\": \"max(0.02, 25% of stage entry total)\",\n";
out << " \"shape_wellbore_parameters_frozen\": true,\n";
@ -1171,6 +1239,7 @@ void nmCalculationAutoFitLM::writeTraceMetaFile()
out << " \"total_stage_shape_constraint\": false,\n";
out << " \"total_parameter_selection\": \"all_valid_free_columns\",\n";
out << " \"skin_difference_policy\": \"local_scale_in_all_stages\",\n";
out << " \"stage2_difference_policy\": \"coordinate_step_cap_0.10_with_full_trust_radius; log_parameter_ratio_cap_1.50; skin_local_scale_fraction_0.10; total_stage_unchanged\",\n";
out << " \"difference_failure_policy\": \"halve_before_opposite_direction\",\n";
out << " \"trace_type\": \"finite_difference_lm_trust_region\",\n";
out << " \"run_id\": " << jsonEscape(m_traceRunId) << ",\n";
@ -1294,7 +1363,9 @@ void nmCalculationAutoFitLM::writeTraceRow(
? traceNumber(qSqrt(trustRegionEarlyShapeEnergy(objectiveBreakdown->shapeResiduals))) : QString())
<< (objectiveBreakdown && objectiveBreakdown->valid ? traceNumber(objectiveBreakdown->earlyValueLoss) : QString())
<< (objectiveBreakdown && objectiveBreakdown->valid ? traceNumber(objectiveBreakdown->earlyParallelLoss) : QString())
<< (objectiveBreakdown && objectiveBreakdown->valid ? traceNumber(objectiveBreakdown->earlyParallelBias) : QString());
<< (objectiveBreakdown && objectiveBreakdown->valid ? traceNumber(objectiveBreakdown->earlyParallelBias) : QString())
<< (objectiveBreakdown && objectiveBreakdown->valid ? traceNumber(trustRegionEarlyGapLoss(*objectiveBreakdown)) : QString())
<< (objectiveBreakdown && objectiveBreakdown->valid ? traceNumber(trustRegionEarlyGapBias(*objectiveBreakdown)) : QString());
QTextStream out(&m_traceFile);
out << cols.join(",") << "\n";
m_traceFile.flush();
@ -2079,7 +2150,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
// 给出，真实求解后用割线估计修正；拒绝时缩步，不让其他参数补偿高度。
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;
@ -2102,10 +2172,9 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
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;
candidate.breakdown.verticalLoss < current.breakdown.verticalLoss;
if(candidate.valid) {
const double slope = (candidate.breakdown.verticalCommonBias - bias) / actualStep;
if(slope < -0.05 && isFiniteNumber(slope)) heightSlope = slope;
@ -2166,6 +2235,23 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
// 先井储、再表皮；回检沿用相同顺序和前期指标，不回到形状阶段反复循环。
bool earlyShapeStage = shapeStage && wellboreParameterCount > 0;
int earlyShapeParameter = storageColumn >= 0 ? 2 : 1;
double initialEarlyGapBias = trustRegionEarlyGapBias(current.breakdown);
double earlyGapDirection = trustRegionEarlyGapDirection(initialEarlyGapBias);
bool earlyShapeHasImproved = false;
auto announceEarlyGapDirection = [&]() {
// 两个参数共用本轮入口的方向；日志保留初始偏差，避免误读为探针响应方向。
const QString directionName = earlyGapDirection > 0.0 ? "increase" :
(earlyGapDirection < 0.0 ? "decrease" : "matched");
writeTraceRow(m_currentIteration, m_enabledParamIndices.indexOf(earlyShapeParameter),
"early_gap_direction", current.parameters, current.fitness, true, 0,
QString("%1_initial_gap_bias_%2").arg(directionName).arg(initialEarlyGapBias, 0, 'g', 12),
&current.breakdown);
emit logMessageGenerated(tr("Early gap guidance: %1, initial signed gap=%2, direction=%3; fixed for storage and skin.")
.arg(earlyShapeParameter == 2 ? tr("wellbore storage") : tr("skin"))
.arg(initialEarlyGapBias, 0, 'g', 6)
.arg(earlyGapDirection > 0.0 ? tr("increase") :
(earlyGapDirection < 0.0 ? tr("decrease") : tr("matched"))));
};
auto parameterAllowedInStage = [&](int column) -> bool {
const int parameterIndex = m_enabledParamIndices[column];
// 候选和缓存灵敏度探针共用此限制，防止探针绕过形状阶段的冻结。
@ -2180,8 +2266,10 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
auto acceptable = [&](const TrustRegionEvaluation& point, const TrustRegionEvaluation& base) -> bool {
if(!point.valid) return false;
if(earlyShapeStage) {
// 初调保持原规则；回检还需保护整轮入口的总误差和形状，不能逐步放宽。
return point.breakdown.earlyParallelLoss < base.breakdown.earlyParallelLoss &&
// 缓存探针也只能沿间距确定的方向接受；停止仍看原来的前期形状改善。
const int column = m_enabledParamIndices.indexOf(earlyShapeParameter);
return earlyGapDirection * (point.coordinates[column] - base.coordinates[column]) > 0.0 &&
point.breakdown.earlyParallelLoss < base.breakdown.earlyParallelLoss &&
(!wellboreRecheck || (point.fitness <= recheckTotalLimit &&
point.breakdown.shapeLoss <= recheckShapeLimit));
}
@ -2191,6 +2279,8 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
: point.fitness < base.fitness;
};
auto acceptPoint = [&](const TrustRegionEvaluation& point) {
// 缓存探针和正式候选共用此标记；首次接受形状改善后，不再启用初始扩步探路。
if(earlyShapeStage) earlyShapeHasImproved = true;
current = point;
publishAcceptedPoint(current);
restoreEvaluationState(current);
@ -2200,8 +2290,9 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
: tr("LM stage 3: original LM fitting; accept by total error only."));
if(earlyShapeStage) {
emit logMessageGenerated(earlyShapeParameter == 2
? tr("Early adjustment: adjust wellbore storage, then skin, to match the early pressure-derivative slope difference of the target.")
: tr("Early adjustment: adjust skin to match the early pressure-derivative slope difference of the target."));
? tr("Early adjustment: storage then skin; pressure-derivative gap selects direction, early shape error controls acceptance and stopping.")
: tr("Early adjustment: skin; pressure-derivative gap selects direction, early shape error controls acceptance and stopping."));
announceEarlyGapDirection();
}
// 有效改善始终相对“上一次有效改善后的误差”累计判断，避免一连串微小
@ -2314,10 +2405,13 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
wellboreRecheck = true;
earlyShapeStage = true;
earlyShapeParameter = storageColumn >= 0 ? 2 : 1;
initialEarlyGapBias = trustRegionEarlyGapBias(current.breakdown);
earlyGapDirection = trustRegionEarlyGapDirection(initialEarlyGapBias);
earlyShapeHasImproved = false;
recheckTotalLimit = current.fitness + qMax(1.0e-4, 0.05 * current.fitness);
recheckShapeLimit = current.breakdown.shapeLoss + qMax(1.0e-4, 0.05 * current.breakdown.shapeLoss);
effectiveImprovementBaseline = current.breakdown.earlyParallelLoss;
// 其他参数已改变，回检入口重测两列，不能沿用初调时的响应方向。
// 其他参数已改变，回检入口重测形状灵敏度，方向由本轮初始间距重新确定。
jacobian.clear();
reuseSensitivityForNextStage();
rebuildRequested = true;
@ -2327,16 +2421,22 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
.arg(recheckTotalLimit, 0, 'g', 6).arg(recheckShapeLimit, 0, 'g', 6));
writeTraceRow(m_currentIteration, -1, "stage_switch", current.parameters, current.fitness,
true, 0, "shape_to_wellbore_recheck", &current.breakdown);
announceEarlyGapDirection();
};
auto finishEarlyShapeStage = [&](const QString& reason) {
if(earlyShapeParameter == 2 && skinColumn >= 0) {
earlyShapeParameter = 1;
earlyShapeHasImproved = false;
effectiveImprovementBaseline = current.breakdown.earlyParallelLoss;
reuseSensitivityForNextStage();
// 井储已变，重测表皮的形状灵敏度和步幅，方向仍保持本轮入口的判断。
rebuildRequested = true;
rebuildReason = QT_TR_NOOP("refresh shape sensitivity at skin adjustment entry");
emit logMessageGenerated(tr("Wellbore storage adjustment ended: %1; now adjusting skin.").arg(reason));
writeTraceRow(m_currentIteration, -1, "stage_switch", current.parameters, current.fitness,
true, 0, wellboreRecheck ? "recheck_storage_to_skin" : "storage_to_skin", &current.breakdown);
announceEarlyGapDirection();
return;
}
if(wellboreRecheck) {
@ -2481,11 +2581,11 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
TrustRegionEvaluation bestProbe;
int bestProbeColumn = -1;
double bestProbeDelta = 0.0;
// 差分步长不超过参数范围的 4%，信赖域收缩后同步减小，但保留 0.5%
// 下限，避免步长太小使求解器数值噪声淹没真实灵敏度。
// 第二阶段允许内部坐标范围的 10% 及完整信赖半径，避免较大比例试探被旧上限截小。
// 第三阶段仍为 4% 及半个信赖半径；均保留 0.5% 下限以减少数值噪声影响。
const double finiteDifferenceStep = qMin(
sensitivityStep,
qMax(5.0e-3, trustRadius * 0.5));
shapeStage ? 0.10 : sensitivityStep,
qMax(5.0e-3, trustRadius * (shapeStage ? 1.0 : 0.5)));
for(int column = 0;
column < dimensions &&
@ -2497,12 +2597,14 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
const double lower = m_parameterLower[parameterIndex];
const double upper = m_parameterUpper[parameterIndex];
if(upper <= lower) continue;
// 表皮在所有阶段均按当前物理尺度扰动，避免整体阶段退回范围的 4%。
// 第二阶段扩大形状试探范围；第三阶段仍沿用原来的局部差分尺度。
// 表皮包含零和负值，按物理尺度扰动；正值参数按对数比例限制幅度。
double localStep = finiteDifferenceStep;
if(parameterIndex == 1) {
localStep = qMin(localStep, 0.02 * qMax(0.1, qAbs(base.parameters[column])) / (upper - lower));
const double skinStepFraction = shapeStage ? 0.10 : 0.02;
localStep = qMin(localStep, skinStepFraction * 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)));
localStep = qMin(localStep, qLn(1.50) / (qLn(upper) - qLn(lower)));
}
double positiveRoom = 1.0 - base.coordinates[column];
double negativeRoom = base.coordinates[column];
@ -2580,7 +2682,9 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
columnBuilt = true;
// 其他参数仍建立灵敏度供后续复用，但不能绕过当前子阶段的选参限制。
if(parameterAllowedInStage(column) && acceptable(probe, base) &&
// 前期首次正式调整必须从小步开始，测形状灵敏度的大探针不能提前成为工作点。
if((!earlyShapeStage || earlyShapeHasImproved) &&
parameterAllowedInStage(column) && acceptable(probe, base) &&
(!bestProbe.valid || stageError(probe) < stageError(bestProbe))) {
bestProbe = probe;
bestProbeColumn = column;
@ -2628,7 +2732,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
(shapeStage ? "shape_accepted_cached_probe" : "total_accepted_cached_probe")),
&current.breakdown);
emit logMessageGenerated(
(earlyShapeStage ? tr("Sensitivity probe accepted: early relative-slope matching error=%1")
(earlyShapeStage ? tr("Sensitivity probe accepted: early shape error=%1")
: tr("Sensitivity probe accepted: total error=%1"))
.arg(earlyShapeStage ? current.breakdown.earlyParallelLoss : current.fitness, 0, 'e', 4));
} else {
@ -2664,7 +2768,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
}
if(wellboreRecheck && current.breakdown.earlyParallelLoss <=
earlyWellboreBaseline + qMax(1.0e-4, 0.10 * earlyWellboreBaseline))
enterTotalStage(tr("early relative-slope error restored within tolerance"));
enterTotalStage(tr("early shape error restored within tolerance"));
if(shapeStage && !earlyShapeStage && !jointShapeStarted) {
jointShapeStarted = true;
shapeRoundBaseline = current.breakdown.shapeLoss;
@ -2831,8 +2935,25 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
}
QVector<double> step;
double reduction = 0.0;
// 井储、表皮均由相对斜率残差的实际灵敏度确定方向，不预设参数增减。
const bool stepBuilt = buildTrustRegionFisherStep(stepInformation, columns, current.coordinates,
double earlyStepRadius = trustRadius;
if(earlyShapeStage && !earlyShapeHasImproved && columns.size() == 1) {
// 确定方向后先小步：井储最多 1.05 倍，表皮为局部尺度的 2%。
// 首次形状变差由后面的扩步循环处理，首次接受后恢复正常 LM 幅度。
const int column = columns[0];
const int parameterIndex = m_enabledParamIndices[column];
const double lower = m_parameterLower[parameterIndex];
const double upper = m_parameterUpper[parameterIndex];
const double smallStep = useTrustRegionLogScale(parameterIndex, lower, upper)
? qLn(1.05) / (qLn(upper) - qLn(lower))
: (parameterIndex == 1 ? 0.02 * qMax(0.1, qAbs(current.parameters[column]))
: 0.05 * qMax(1.0e-8, qAbs(current.parameters[column]))) / (upper - lower);
earlyStepRadius = qMin(trustRadius, smallStep);
}
// 井储、表皮沿间距确定的方向调整；其他阶段保持原来的 LM 求步。
const bool stepBuilt = earlyShapeStage
? buildEarlyGapGuidedStep(stepInformation, columns, current.coordinates, earlyGapDirection,
damping, earlyStepRadius, minimumCoordinateStep, &step, &reduction)
: buildTrustRegionFisherStep(stepInformation, columns, current.coordinates,
damping, trustRadius, minimumCoordinateStep, &step, &reduction);
if(!stepBuilt) continue;
// 仅联合形状阶段预测总误差上限；井储、表皮步长不受总误差裁剪。
@ -2930,7 +3051,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
candidate.coordinates = candidateCoordinates;
candidate.parameters = parametersFromCoordinates(candidate.coordinates);
if(earlyShapeStage && earlyShapeParameter == 2) {
emit logMessageGenerated(tr("Wellbore storage trial: relative-slope matching error=%1, C=%2 -> %3")
emit logMessageGenerated(tr("Wellbore storage trial: early shape error=%1, C=%2 -> %3")
.arg(current.breakdown.earlyParallelLoss, 0, 'g', 5)
.arg(current.parameters[storageColumn], 0, 'g', 6)
.arg(candidate.parameters[storageColumn], 0, 'g', 6));
@ -2942,6 +3063,49 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
&candidate.curve,
&candidate.elapsedMs);
// 初始形状变差时先沿锁定方向跨大步探路，不立即缩步或累计无效次数。
// 始终以未移动的 current 为基点；每次翻倍并裁到边界，差解不发布、不接受。
bool earlyInitialExpanded = false;
bool earlyInitialExpansionAtBound = false;
while(earlyShapeStage && !earlyShapeHasImproved && candidate.valid &&
candidate.breakdown.earlyParallelLoss > current.breakdown.earlyParallelLoss &&
processPauseAndStop()) {
QVector<double> expandedStep;
double expandedPrediction = 0.0;
if(!buildExpandedTrustRegionStep(stepInformation, current.coordinates, coordinateStep,
1.0, &expandedStep, &expandedPrediction, false)) {
earlyInitialExpansionAtBound = true;
break;
}
writeTraceRow(m_currentIteration, selectedColumns[0], "early_initial_expansion",
candidate.parameters, candidate.fitness, true, candidate.elapsedMs,
"shape_worse_expand_same_direction", &candidate.breakdown);
restoreEvaluationState(current);
TrustRegionEvaluation expanded;
expanded.coordinates = current.coordinates;
for(int column = 0; column < dimensions; ++column)
expanded.coordinates[column] += expandedStep[column];
expanded.parameters = parametersFromCoordinates(expanded.coordinates);
const int column = selectedColumns[0];
emit logMessageGenerated(tr("Early shape initially worsened: %1, trial=%2 -> %3; expanding in the same direction.")
.arg(earlyShapeParameter == 2 ? tr("wellbore storage") : tr("skin"))
.arg(candidate.parameters[column], 0, 'g', 6)
.arg(expanded.parameters[column], 0, 'g', 6));
expanded.valid = evaluateTrustRegionPoint(expanded.parameters, &expanded.fitness,
&expanded.breakdown, &expanded.curve, &expanded.elapsedMs);
candidate = expanded;
coordinateStep = expandedStep;
predictedReduction = expandedPrediction;
stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep));
earlyInitialExpanded = true;
}
if(earlyInitialExpanded) selectionName += "_initial_expanded";
if(m_shouldStop) {
restoreEvaluationState(current);
stopReason = LM_USER_STOPPED;
break;
}
if(!candidate.valid) {
// 求解失败的候选不能改变 current。先完整恢复上一个已接受参数和
// 对应误差快照，再缩小信赖域；连续失败达到上限才终止整个拟合。
@ -3108,9 +3272,9 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
} else if(componentName == "horizontal") {
componentDisplayName = tr("horizontal deviation");
} else if(componentName == "storage_parallel") {
componentDisplayName = tr("wellbore storage relative-slope matching");
componentDisplayName = tr("wellbore storage early shape");
} else if(componentName == "skin_parallel") {
componentDisplayName = tr("skin relative-slope matching");
componentDisplayName = tr("skin early shape");
} else if(componentName == "shape") {
componentDisplayName = tr("shape deviation");
} else if(componentName == "total") {
@ -3127,6 +3291,10 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
emit progressUpdated(shapeStage ? -1 : totalIterations, m_globalBestFitness);
// 先记录当前试调结果，再发布阶段切换，避免日志显示为未试调就结束。
if(earlyInitialExpansionAtBound && !accepted) {
finishEarlyShapeStage(tr("initial same-direction expansion reached the parameter bound without early shape improvement"));
continue;
}
const bool effectiveImprovement = registerEffectiveImprovement(stageError(current));
if(!effectiveImprovement && recordIneffectiveStep()) stopReason = LM_LOCAL_OPTIMUM;

11
Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp
Unescape Escape View File

@ -157,6 +157,11 @@ bool nmWxAutomaticFitting::getPhysicalParameterRange(int parameterIndex,
return false;
}
// 自动拟合的表皮搜索下界固定为 0，自动范围和手工输入校验使用同一边界。
if(parameterIndex == 1) {
minValue = 0.0;
}
if(parameterIndex == 4) {
double soi = reservoirData.getSoi().getValue().toDouble();
double sgi = reservoirData.getSgi().getValue().toDouble();
@ -233,8 +238,8 @@ void nmWxAutomaticFitting::setParameterRange(int parameterIndex,
m_updatingParameterRanges = wasUpdatingRanges;
}
// 根据当前初值生成建议搜索范围，并始终截断在系统物理边界内。skin 使用
// 加减固定宽度，其余正值参数使用倍率范围；该规则在首次加载和拟合完成后复用。
// 根据当前初值生成建议搜索范围，并始终截断在拟合边界内。skin 下界为 0，
// 上界按初值加固定宽度，其余正值参数使用倍率范围；首次加载和拟合完成后复用。
void nmWxAutomaticFitting::updateRangeForParameter(int parameterIndex,
double centerValue)
{
@ -271,7 +276,7 @@ void nmWxAutomaticFitting::updateRangeForParameter(int parameterIndex,
double newMax = physicalMax;
if(parameterIndex == 1) {
const double skinHalfRange = 10.0;
newMin = qMax(physicalMin, reference - skinHalfRange);
newMin = physicalMin;
newMax = qMin(physicalMax, reference + skinHalfRange);
} else if(reference > 0.0
&& !(parameterIndex == 4 && centerValue <= 0.0)

Write Preview
Loading…
Cancel
Save