From 68c79378bbe04e6dbfca1112c168a48f8ac086fc Mon Sep 17 00:00:00 2001 From: lvjunjie Date: Sun, 20 Sep 2026 15:03:48 +0800 Subject: [PATCH] =?UTF-8?q?feat(nmNum):=20=E4=BC=98=E5=8C=96=20LM=20?= =?UTF-8?q?=E5=89=8D=E6=9C=9F=E5=BD=A2=E6=80=81=E5=BC=95=E5=AF=BC=E4=B8=8E?= =?UTF-8?q?=E5=88=86=E9=98=B6=E6=AE=B5=E5=8F=82=E6=95=B0=E8=B0=83=E6=95=B4?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 按初始曲线间距固定井储和表皮方向，小步起调，形状变差时同向扩步 - 前期采用双曲线斜率形状误差，放大第二阶段灵敏度试探步长 - 渗透率预调整仅看上下偏移，表皮搜索下限设为 0 - 补充间距与方向诊断日志，保留第三阶段原有逻辑 --- .../nmCalculation/nmCalculationAutoFitLM.h | 5 +- .../nmCalculation/nmCalculationAutoFitLM.cpp | 244 +++++++++++++++--- Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp | 11 +- 3 files changed, 217 insertions(+), 43 deletions(-) diff --git a/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h b/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h index 579abe33..5e6d028d 100644 --- a/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h +++ b/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h @@ -48,8 +48,9 @@ struct AutoFitObjectiveBreakdownLM { bool registrationAmbiguous; // 固定 log-time 网格上的双曲线斜率残差，独立于数值残差的采样层级。 QVector shapeResiduals; - // 前期数值残差为 81 点，平行程度残差为 73 个斜率区间，均含第一窗口权重。 - // 平行误差比较模拟与目标各自的压力—导数斜率差，不要求模拟自身斜率差为零。 + // 前期数值残差为 81 点，形状残差为 73 个斜率区间，均含第一窗口权重。 + // earlyParallelResiduals/Loss 沿用历史字段名，现为压力、导数分别匹配目标的形状误差。 + // earlyParallelBias 仍记录相对斜率偏差，仅供诊断。 QVector earlyValueResiduals; QVector earlyParallelResiduals; double earlyValueLoss; diff --git a/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp b/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp index 2af9406a..1aeccba9 100644 --- a/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp +++ b/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp @@ -328,7 +328,7 @@ static bool trustRegionResidualsValid( return true; } -// 数值、形状和前期平行程度残差共用一次求解，联合缓存供各子阶段复用。 +// 数值、整体形状和前期形状残差共用一次求解，联合缓存供各子阶段复用。 static QVector trustRegionFullResidual(const AutoFitObjectiveBreakdownLM& objective) { QVector residual = objective.residualVector; @@ -348,6 +348,46 @@ static double trustRegionSquaredNorm(const QVector& 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 weights(count, 0.0); + double weightSum = 0.0; + for(int i = 0; i < count; ++i) { + weights[i] = autoFitTimeWindowWeight(static_cast(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& residual) @@ -361,7 +401,7 @@ static double trustRegionEarlyShapeEnergy(const QVector& residual) return energy; } -// 比较模拟与目标各自的压力—导数斜率差；误差为零表示两组曲线的相对走势一致。 +// 前期复用整体形状的双曲线斜率残差，只按第一窗口重新加权并归一化。 // 单独平移任一曲线不改变此指标；前期数值误差仅作诊断，不参与井储、表皮验收。 static void populateEarlyWellboreMetrics(AutoFitObjectiveBreakdownLM* objective, const QVector targetLogs[2], const QVector 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& selected, const QVector& coordinates, double direction, + double damping, double trustRadius, double minimumStep, + QVector* 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& coordinates, const QVector& originalStep, - double maximumRadius, QVector* expandedStep, double* prediction) + double maximumRadius, QVector* 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), + ¤t.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", ¤t.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", ¤t.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")), ¤t.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 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 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; diff --git a/Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp b/Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp index 1f0e2f3c..9204f0e5 100644 --- a/Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp +++ b/Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp @@ -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)