diff --git a/Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h b/Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h index e65d1ee..b3e201a 100644 --- a/Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h +++ b/Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h @@ -20,64 +20,66 @@ class nmDataWellBase; class QTimer; class QProcess; -// 双对数曲线误差分解，供误差诊断和后续参数调整读取。 -// -// 所有 pressure/derivative/shape 数值均是在 log(value) 空间计算的无量纲误差。 -// verticalBias* 保留正负号：正值表示模拟曲线整体高于目标，负值表示整体低于目标。 -// horizontalPhysicalShift 是 log(time) 方向的等效平移量，正值表示模拟曲线相对目标偏右。 -// total 仍是 PSO 当前使用的 fitness（用于粒子比较和收敛判断）；其余字段只描述误差 -// 来源，不会在本次改动中直接修改粒子参数，避免损失诊断和参数更新策略互相耦合。 +// 双对数曲线误差分解。该结构同时保存用于候选排序的主目标，以及用于判断 +// 曲线上下、左右和形状偏差的诊断量。total 是唯一的接受和排序依据，诊断量 +// 只参与信赖域选参，不能再次叠加到 total，否则会重复计算同一批曲线残差。 struct AutoFitObjectiveBreakdown { - // valid/total 是本次评价是否有效及其最终 fitness（用于排序和收敛判断）。 + // valid 表示本次曲线评价完整有效；无效评价统一保留 total=1e10。 + // pressureLoss 和 derivativeLoss 均在 log(value) 空间按固定网格计算。 bool valid; double total; - // pressureLoss/derivativeLoss 是压力和导数两条曲线的整体误差。 double pressureLoss; double derivativeLoss; - // vertical* 描述整体上下偏移；保留 bias 的符号以判断偏高或偏低。 - double verticalBiasPressure; - double verticalBiasDerivative; + // 固定目标网格上的 Huber 等效残差。非代理搜索使用它建立完整 Jacobian， + // 向量平方和与 total 的平方一致。 + QVector residualVector; + + // 上下偏差使用压力和导数残差共享的 Huber 稳健中心。 + // verticalCommonBias 为正表示模拟曲线整体偏高，为负表示整体偏低； + // verticalReliable=false 时仍保留数值，但不能据此确定参数调整方向。 + double verticalCommonBias; double verticalLoss; - // horizontal* 描述等效的对数时间偏移；physicalShift 为正表示模拟曲线相对目标向右 - // （时间延迟），为负表示向左。 - double horizontalShift; + bool verticalReliable; + + // 水平偏差在 log(time) 坐标中计算。physicalShift 为正表示模拟曲线相对 + // 目标偏右，即相同曲线特征在模拟结果中出现得更晚。 double horizontalPhysicalShift; double horizontalLoss; - // shapeLoss 是去除整体上下和左右偏移后剩余的曲线形状差异。 + bool horizontalReliable; + // true 表示当前曲线无法可靠区分上下和左右误差；此时禁止使用两类有符号 + // 诊断量选参，但去除公共中心后的 shapeLoss 仍可用于局部选参。 + bool registrationAmbiguous; + + // 去除稳健公共中心和可信左右偏差后剩余的整体形状误差；verticalReliable + // 只控制能否把公共中心解释为上下参数方向，不改变 shape 的中心化公式。 double shapeLoss; - // 早、中、晚分段误差用于定位误差主要出现在哪个时间阶段。 - double pressureEarlyLoss; - double pressureMiddleLoss; - double pressureLateLoss; - double derivativeEarlyLoss; - double derivativeMiddleLoss; - double derivativeLateLoss; - // coverage 是 50 点目标网格上的 min(有效点比例、连续 log-time 跨度比例)。 - // coveragePenalty 是归一化覆盖缺口的平方惩罚，并以 0.1 权重加入 total。 + + // 兼容现有 trace 列。当前非代理搜索不再单独识别或调度晚期分量。 + double lateDerivativeSlopeBias; + double lateDerivativeTrendLoss; + bool lateDerivativeTrendReliable; + + // 模拟曲线对目标固定网格的有效覆盖率，取覆盖点比例与连续 log-time + // 跨度比例中的较小值。低于损失函数门槛时本次评价直接无效。 double coverage; - double coveragePenalty; - // 无效评价使用 1e10 作为统一的“差解”标记；其他字段用 NaN 表示尚未得到诊断值。 AutoFitObjectiveBreakdown() : valid(false) , total(1.0e10) , pressureLoss(std::numeric_limits::quiet_NaN()) , derivativeLoss(std::numeric_limits::quiet_NaN()) - , verticalBiasPressure(std::numeric_limits::quiet_NaN()) - , verticalBiasDerivative(std::numeric_limits::quiet_NaN()) + , verticalCommonBias(std::numeric_limits::quiet_NaN()) , verticalLoss(std::numeric_limits::quiet_NaN()) - , horizontalShift(std::numeric_limits::quiet_NaN()) + , verticalReliable(false) , horizontalPhysicalShift(std::numeric_limits::quiet_NaN()) , horizontalLoss(std::numeric_limits::quiet_NaN()) + , horizontalReliable(false) + , registrationAmbiguous(false) , shapeLoss(std::numeric_limits::quiet_NaN()) - , pressureEarlyLoss(std::numeric_limits::quiet_NaN()) - , pressureMiddleLoss(std::numeric_limits::quiet_NaN()) - , pressureLateLoss(std::numeric_limits::quiet_NaN()) - , derivativeEarlyLoss(std::numeric_limits::quiet_NaN()) - , derivativeMiddleLoss(std::numeric_limits::quiet_NaN()) - , derivativeLateLoss(std::numeric_limits::quiet_NaN()) + , lateDerivativeSlopeBias(std::numeric_limits::quiet_NaN()) + , lateDerivativeTrendLoss(std::numeric_limits::quiet_NaN()) + , lateDerivativeTrendReliable(false) , coverage(std::numeric_limits::quiet_NaN()) - , coveragePenalty(std::numeric_limits::quiet_NaN()) {} }; @@ -95,6 +97,8 @@ struct AutoFitParticle { QVector velocity; // 速度 QVector bestPosition; // 真实求解器确认的个体最优位置 QVector guideBestPosition; // 仅用于速度更新的引导位置；不会参与真实 gbest/最终结果 + AutoFitObjectiveBreakdown currentObjectiveBreakdown; // 当前真实评价对应的误差分解 + AutoFitObjectiveBreakdown bestObjectiveBreakdown; // pbest 对应的误差分解 double fitness; // 当前适应度 double bestFitness; // 真实求解器确认的个体最优适应度 double guideBestObjective; // guideBestPosition 对应的真实或代理目标值 @@ -200,20 +204,26 @@ private: void loadOptimizationConfig(); void loadParameterBounds(); - // ===== PSO核心算法 ===== + // ===== 自动拟合核心算法 ===== // - // 主流程： - // 1. extractUserInitialValues(): 从当前项目数据中取用户已有初始解； - // 2. initializeSwarm(): 根据初始解和上下界生成粒子群； - // 3. updateParticle(): 对单个粒子跑真实求解器并计算误差； - // 4. updateGlobalBest(): 只用真实求解器误差更新全局最优； - // 5. updateVelocityAndPosition(): 按 PSO 公式推进下一代粒子。 + // 代理开启时保留原 PSO 筛选流程；代理关闭时使用真实求解器驱动的 + // 诊断灵敏度信赖域搜索，不依赖 pbest/gbest 速度公式。 void extractUserInitialValues(); void initializeSwarm(); void updateVelocityAndPosition(); double evaluateFitness(const QVector& parameters); void updateGlobalBest(); void updateParticle(int particleIndex); + // 非代理拟合入口：建立有限差分灵敏度，按诊断分量选择参数，再用有界 + // LM/信赖域产生候选；所有候选最终都由真实求解器总误差决定是否接受。 + StopReasonPSO runTrustRegionFitting(); + // 对一个信赖域候选执行完整真实评价，并一次性返回误差、诊断量、曲线和耗时。 + // 返回 false 表示求解失败、损失无效或用户已请求停止。 + bool evaluateTrustRegionPoint(const QVector& parameters, + double* fitness, + AutoFitObjectiveBreakdown* breakdown, + QVector >* curve, + int* elapsedMs); // ===== 参数应用方法 ===== // @@ -284,7 +294,8 @@ private: double surrogateObjective, const QString& screeningDecision, const QVector& pbestPosition, - double pbestObjective); + double pbestObjective, + const AutoFitObjectiveBreakdown* objectiveBreakdown = nullptr); void writeIterationTraceRows(); QVector buildTraceParameterVector(const QVector& selectedParameters) const; void resetRunSummary(); @@ -340,19 +351,21 @@ private: bool m_isRunning; // 当前是否有一次自动拟合正在运行。 bool m_shouldStop; // 用户停止标志；主循环和求解器等待循环会定期检查它。 bool m_isPaused; // 预留暂停标志；主循环中有暂停等待逻辑。 - int m_currentIteration; // 当前 PSO 迭代序号，从 0 开始。 + int m_currentIteration; // 当前自动拟合迭代序号，从 0 开始。 QString m_lastError; // 最近一次失败原因，供 UI 展示或日志排查。 - // ===== PSO数据 ===== + // ===== 优化状态数据 ===== QVector m_initialValues; // 当前模型中提取的用户初始值，顺序与 m_enabledParamIndices 一致。 QVector m_swarm; // 粒子群，每个粒子只保存启用参数维度。 - QVector m_globalBestPosition; // 全局最优参数，仍是启用参数向量。 + QVector m_globalBestPosition; // 真实求解器确认的当前最优参数。 double m_globalBestFitness; // 全局最优真实误差，越小越好。 double m_previousBestFitness; // 上一轮全局最优误差，用于自适应参数更新。 + AutoFitObjectiveBreakdown m_globalBestObjectiveBreakdown; // 真实 gbest 对应的误差分解。 QVector > m_lastEvaluatedLogLogData; // 最近一次真实求解得到的 result log-log 曲线。 QVector > m_globalBestLogLogData; // 当前全局最优对应的 result log-log 曲线。 mutable AutoFitObjectiveBreakdown m_lastObjectiveBreakdown; // 最近一次损失评价的误差分解。 QVector > m_userInitialLogLogData; // 用户初始解对应的 result log-log 曲线，用于精英保护。 + AutoFitObjectiveBreakdown m_userInitialObjectiveBreakdown; // 用户初始解对应的误差分解。 // ===== 优化配置 ===== // @@ -377,7 +390,7 @@ private: double m_socialParam; // 群体学习因子，控制粒子靠近全局 gbest 的程度。 // ===== 统计信息 ===== - int m_totalEvaluations; // 已调用真实求解器评价的粒子总数。 + int m_totalEvaluations; // 真实求解器评价总次数，包含粒子评价和方向试算。 int m_successfulEvaluations; // 真实求解器成功且误差有效的评价次数。 QVector m_convergenceHistory; // 每代全局最优误差历史，用于收敛判断。 @@ -393,7 +406,7 @@ private: // ===== 精英保护 ===== QVector m_userInitialSolution; // 用户初始解参数，若最终改进不足会恢复它。 double m_userInitialFitness; // 用户初始解真实误差。 - double m_improvementThreshold; // 最终结果相对初始解至少需要达到的改进阈值。 + double m_improvementThreshold; // 仅用于日志区分显著改进和微小改进。 bool m_hasValidUserSolution; // 初始解是否成功跑过真实求解器。 int m_consecutiveFailedIterations; // 连续失败迭代次数 @@ -422,8 +435,8 @@ private: // 这些字段只描述代理筛选和运行复盘，不参与 PSO 数学更新。 bool m_traceEnabled; // 是否写出 trace CSV/meta 文件。 QString m_traceRunId; // 本次运行 ID，作为 trace/candidate/score 文件名的一部分。 - QString m_traceFilePath; // pso_baseline_trace_.csv 完整路径。 - QString m_traceMetaFilePath; // pso_baseline_trace_.meta.json 完整路径。 + QString m_traceFilePath; // 本次自动拟合 trace CSV 的完整路径。 + QString m_traceMetaFilePath; // 与 trace 匹配的 meta JSON 完整路径。 QFile m_traceFile; // trace CSV 文件句柄。 bool m_surrogateScreeningEnabled; // 用户配置中的 PSO acceleration 开关。 unsigned int m_psoRandomSeed; // PSO 随机种子，也用于可复现 random audit。 diff --git a/Include/nmNum/nmSubWxs/nmWxAutomaticFitting.h b/Include/nmNum/nmSubWxs/nmWxAutomaticFitting.h index 30a9119..5b3848c 100644 --- a/Include/nmNum/nmSubWxs/nmWxAutomaticFitting.h +++ b/Include/nmNum/nmSubWxs/nmWxAutomaticFitting.h @@ -57,7 +57,7 @@ private: void renumberVisibleParameterRows(QTableWidget* table); void updateParameterVisibility(QTableWidget* table, NM_SOLVER_MODEL_TYPE eType); void initializeSuggestedParameterRanges(); - void updateRangeForParameter(int parameterIndex, double centerValue, bool afterFit); + void updateRangeForParameter(int parameterIndex, double centerValue); void setParameterRange(int parameterIndex, double minValue, double maxValue); bool getPhysicalParameterRange(int parameterIndex, double& minValue, double& maxValue); void normalizeSavedParameterRanges(); diff --git a/Src/nmNum/nmCalculation/nmCalculationAutoFitPSO.cpp b/Src/nmNum/nmCalculation/nmCalculationAutoFitPSO.cpp index ccc05ea..a4eb476 100644 --- a/Src/nmNum/nmCalculation/nmCalculationAutoFitPSO.cpp +++ b/Src/nmNum/nmCalculation/nmCalculationAutoFitPSO.cpp @@ -104,6 +104,17 @@ static inline bool isFiniteNumber(double value) #endif } +// 两个 Huber RMS 的差不能直接解释为被消除的独立误差。RMS 的平方才对应 +// 稳健能量，因此先计算 reduced^2-full^2，再开方恢复原量纲。这里用于分别 +// 提取“消除公共上下偏差”和“消除水平位移”实际减少的误差贡献。 +static double nestedRmsContribution(double reducedModelLoss, + double fullModelLoss) +{ + return qSqrt(qMax(0.0, + reducedModelLoss * reducedModelLoss - + fullModelLoss * fullModelLoss)); +} + static inline bool isInClosedRange(double value, double lower, double upper) { return isFiniteNumber(value) && value >= lower && value <= upper; @@ -813,10 +824,12 @@ void nmCalculationAutoFitPSO::resetOptimizer() m_globalBestPosition.clear(); m_globalBestFitness = 1e10; m_previousBestFitness = 1e10; + m_globalBestObjectiveBreakdown = AutoFitObjectiveBreakdown(); m_lastEvaluatedLogLogData.clear(); m_globalBestLogLogData.clear(); m_lastObjectiveBreakdown = AutoFitObjectiveBreakdown(); m_userInitialLogLogData.clear(); + m_userInitialObjectiveBreakdown = AutoFitObjectiveBreakdown(); m_currentIteration = 0; m_totalEvaluations = 0; m_successfulEvaluations = 0; @@ -844,9 +857,8 @@ void nmCalculationAutoFitPSO::setPSOTargetWellName(const QString& wellName) void nmCalculationAutoFitPSO::initializeTraceFile() { - // 创建本次 PSO 的可复盘文件： - // - pso_baseline_trace_.csv：逐代逐粒子的参数、真实误差、代理误差和筛选决策； - // - pso_baseline_trace_.meta.json：目标曲线、流量制度、参数上下界、PSO/代理配置。 + // 创建本次自动拟合的可复盘文件。代理 PSO 与非代理信赖域使用不同前缀， + // 防止代理回放脚本把信赖域记录误当成最新 PSO 粒子记录。 // // 代理模型评分脚本也会读取 meta.json，因此 trace meta 不是单纯日志，而是C++ 与 Python 代理模型之间的运行上下文契约。 if(!m_traceEnabled) { @@ -865,8 +877,13 @@ void nmCalculationAutoFitPSO::initializeTraceFile() return; } - m_traceFilePath = traceDir.absoluteFilePath(QString("pso_baseline_trace_%1.csv").arg(m_traceRunId)); - m_traceMetaFilePath = traceDir.absoluteFilePath(QString("pso_baseline_trace_%1.meta.json").arg(m_traceRunId)); + QString tracePrefix = isSurrogateScreeningEnabled() + ? "pso_baseline_trace" + : "trust_region_trace"; + m_traceFilePath = traceDir.absoluteFilePath( + QString("%1_%2.csv").arg(tracePrefix).arg(m_traceRunId)); + m_traceMetaFilePath = traceDir.absoluteFilePath( + QString("%1_%2.meta.json").arg(tracePrefix).arg(m_traceRunId)); m_traceFile.setFileName(m_traceFilePath); if(!m_traceFile.open(QIODevice::WriteOnly | QIODevice::Text)) { @@ -878,11 +895,12 @@ void nmCalculationAutoFitPSO::initializeTraceFile() writeTraceHeader(); writeTraceMetaFile(); - DEBUG_OUT(QString("PSO baseline trace initialized: %1").arg(m_traceFilePath)); - emit logMessageGenerated(tr("PSO baseline trace: %1").arg(m_traceFilePath)); + DEBUG_OUT(QString("Automatic fitting trace initialized: %1").arg(m_traceFilePath)); + emit logMessageGenerated(tr("Automatic fitting trace: %1").arg(m_traceFilePath)); if(!m_traceMetaFilePath.isEmpty()) { - emit logMessageGenerated(tr("PSO baseline trace meta: %1").arg(m_traceMetaFilePath)); + emit logMessageGenerated( + tr("Automatic fitting trace meta: %1").arg(m_traceMetaFilePath)); } emit logMessageGenerated(tr("PSO surrogate screening: %1, model=%2, keep=%3, audit=%4, warmup=%5, min_solver=%6") @@ -997,7 +1015,8 @@ void nmCalculationAutoFitPSO::writeTraceHeader() // - solver_objective 是真实求解器误差； // - surrogate_objective 是 Python 代理评分； // - screening_decision 说明该粒子为什么跑/不跑真实求解器； - // - pbest/gbest 字段用于离线复盘 PSO 更新是否只依赖真实误差。 + // - pbest/gbest 字段用于离线复盘 PSO 更新是否只依赖真实误差； + // - 末尾诊断字段记录同一次真实评价的分量误差，便于核对引导方向和接受结果。 if(!m_traceFile.isOpen()) { return; } @@ -1035,7 +1054,26 @@ void nmCalculationAutoFitPSO::writeTraceHeader() << "gbest_h" << "gbest_Ct" << "gbest_Cf" - << "enabled_param_indices"; + << "enabled_param_indices" + << "pressure_loss" + << "derivative_loss"; + if(!isSurrogateScreeningEnabled()) { + cols << "vertical_common_bias"; + } + cols << "vertical_loss"; + if(!isSurrogateScreeningEnabled()) { + cols << "vertical_reliable" + << "horizontal_physical_shift"; + } + cols << "horizontal_loss"; + if(!isSurrogateScreeningEnabled()) { + cols << "horizontal_reliable"; + } + cols << "shape_loss" + << "late_trend_loss" + << "late_slope_bias" + << "late_trend_reliable" + << "registration_ambiguous"; QTextStream out(&m_traceFile); out << cols.join(",") << "\n"; @@ -1077,8 +1115,14 @@ void nmCalculationAutoFitPSO::writeTraceMetaFile() QTextStream out(&metaFile); out << "{\n"; - out << " \"schema_version\": 1,\n"; - out << " \"trace_type\": \"pso_baseline_replay_meta\",\n"; + // 非代理 v5 增加带符号诊断列；代理 PSO 保留原 v3 字段和目标，避免改变 + // 已有模型的训练和回放契约。 + out << " \"schema_version\": " + << (isSurrogateScreeningEnabled() ? 3 : 5) << ",\n"; + out << " \"trace_type\": " + << jsonEscape(isSurrogateScreeningEnabled() + ? "pso_baseline_replay_meta" + : "diagnostic_trust_region_meta") << ",\n"; out << " \"run_id\": " << jsonEscape(m_traceRunId) << ",\n"; out << " \"created_at\": " << jsonEscape(QDateTime::currentDateTime().toString(Qt::ISODate)) << ",\n"; out << " \"trace_csv\": " << jsonEscape(QFileInfo(m_traceFilePath).fileName()) << ",\n"; @@ -1180,7 +1224,8 @@ void nmCalculationAutoFitPSO::writeTraceRow(int generation, double surrogateObjective, const QString& screeningDecision, const QVector& pbestPosition, - double pbestObjective) + double pbestObjective, + const AutoFitObjectiveBreakdown* objectiveBreakdown) { // 写一行 trace。generation=-1/particleIndex=-1 表示用户初始解； // 普通粒子行的 phase 为 particle_solver、particle_verified_cache 或 particle_not_evaluated。 @@ -1234,6 +1279,35 @@ void nmCalculationAutoFitPSO::writeTraceRow(int generation, << traceParamAt(gbestParams, 6) << csvEscape(enabledIndices.join(";")); + if(objectiveBreakdown && objectiveBreakdown->valid) { + cols << traceNumber(objectiveBreakdown->pressureLoss) + << traceNumber(objectiveBreakdown->derivativeLoss); + if(!isSurrogateScreeningEnabled()) { + cols << traceNumber(objectiveBreakdown->verticalCommonBias); + } + cols << traceNumber(objectiveBreakdown->verticalLoss); + if(!isSurrogateScreeningEnabled()) { + cols << QString::number(objectiveBreakdown->verticalReliable ? 1 : 0) + << traceNumber(objectiveBreakdown->horizontalPhysicalShift); + } + cols << traceNumber(objectiveBreakdown->horizontalLoss); + if(!isSurrogateScreeningEnabled()) { + cols << QString::number(objectiveBreakdown->horizontalReliable ? 1 : 0); + } + cols << traceNumber(objectiveBreakdown->shapeLoss) + << traceNumber(objectiveBreakdown->lateDerivativeTrendLoss) + << traceNumber(objectiveBreakdown->lateDerivativeSlopeBias) + << QString::number(objectiveBreakdown->lateDerivativeTrendReliable ? 1 : 0) + << QString::number(objectiveBreakdown->registrationAmbiguous ? 1 : 0); + } else { + // 未运行真实求解器或评价无效时保持列数一致，诊断字段写空值。 + int diagnosticColumnCount = + isSurrogateScreeningEnabled() ? 9 : 13; + for(int i = 0; i < diagnosticColumnCount; ++i) { + cols << QString(); + } + } + QTextStream out(&m_traceFile); out << cols.join(",") << "\n"; m_traceFile.flush(); @@ -1265,9 +1339,11 @@ void nmCalculationAutoFitPSO::writeIterationTraceRows() particle.lastEvaluationSuccess, particle.lastEvaluationElapsedMs, particle.surrogateObjective, - particle.screeningDecision, - particle.bestPosition, - particle.bestFitness); + particle.screeningDecision, + particle.bestPosition, + particle.bestFitness, + particle.evaluatedThisIteration + ? &particle.currentObjectiveBreakdown : nullptr); } } @@ -2855,21 +2931,19 @@ void nmCalculationAutoFitPSO::loadParameterBounds() .arg(m_enabledParamIndices.size())); } -// ==================== PSO算法核心方法 ==================== +// ==================== 自动拟合核心方法 ==================== bool nmCalculationAutoFitPSO::startAutoFitting() { - // 自动拟合的总入口。可以把这个函数当成 PSO 的“运行剧本”： - // 读取配置 -> 校验输入 -> 评价用户初始解 -> 初始化粒子群 -> - // 按代循环评价粒子 -> 更新全局最优 -> 判断停止 -> 保存结果。 + // 自动拟合总入口：代理开启时保留原 PSO 筛选流程；代理关闭时改走 + // 诊断灵敏度信赖域搜索。两条路径共用初始解评价、真实求解器和结果写回。 + StopReasonPSO finalReason = PSO_CONTINUE_OPTIMIZATION; + bool useParticleSwarm = false; + if(m_isRunning) { m_lastError = "Auto fitting is already running"; return false; } - // 发送初始化日志 - //emit logMessageGenerated(tr("=== PSO Automatic Fitting Started ===")); - emit logMessageGenerated(tr("Algorithm: Particle Swarm Optimization")); - try { // 从 DataManager 读取界面保存的自动拟合配置。 // 本类不直接依赖 UI 控件，便于后续从脚本或其他入口复用。 @@ -2878,6 +2952,11 @@ bool nmCalculationAutoFitPSO::startAutoFitting() return false; } + useParticleSwarm = isSurrogateScreeningEnabled(); + emit logMessageGenerated(useParticleSwarm + ? tr("Algorithm: Particle Swarm Optimization") + : tr("Algorithm: Diagnostic Trust-Region Search")); + if(m_simulationMode) { // 调试/演示用快速路径，不调用真实求解器。正式工况通常不走这里。 DEBUG_OUT("=== SIMULATION MODE ACTIVATED ==="); @@ -3007,6 +3086,8 @@ bool nmCalculationAutoFitPSO::startAutoFitting() m_globalBestPosition = m_userInitialSolution; m_userInitialLogLogData = m_lastEvaluatedLogLogData; m_globalBestLogLogData = m_userInitialLogLogData; + m_userInitialObjectiveBreakdown = m_lastObjectiveBreakdown; + m_globalBestObjectiveBreakdown = m_userInitialObjectiveBreakdown; emit logMessageGenerated(tr("Initial solution evaluation successful")); emit logMessageGenerated(tr("Initial Error: %1").arg(m_userInitialFitness, 0, 'e', 4)); @@ -3025,9 +3106,11 @@ bool nmCalculationAutoFitPSO::startAutoFitting() m_userInitialFitness < 1e9, initialEvalElapsedMs, std::numeric_limits::quiet_NaN(), - "initial_solution", - m_hasValidUserSolution ? m_userInitialSolution : QVector(), - m_userInitialFitness); + "initial_solution", + m_hasValidUserSolution ? m_userInitialSolution : QVector(), + m_userInitialFitness, + m_hasValidUserSolution + ? &m_userInitialObjectiveBreakdown : nullptr); } catch(...) { m_hasValidUserSolution = false; emit logMessageGenerated(tr("Exception during initial solution evaluation")); @@ -3037,12 +3120,13 @@ bool nmCalculationAutoFitPSO::startAutoFitting() m_initialValues = savedInitialValues; } - // 初始化粒子群。粒子维度等于用户勾选的参数数量，而不是固定 11 维。 - if(kUseFixedPsoSeed) { - emit logMessageGenerated(tr("PSO random seed: %1 ").arg(m_psoRandomSeed)); - } else { - emit logMessageGenerated(tr("PSO random seed: %1 ").arg(m_psoRandomSeed)); - } + if(useParticleSwarm) { + // 初始化粒子群。粒子维度等于用户勾选的参数数量，而不是固定 11 维。 + if(kUseFixedPsoSeed) { + emit logMessageGenerated(tr("PSO random seed: %1 ").arg(m_psoRandomSeed)); + } else { + emit logMessageGenerated(tr("PSO random seed: %1 ").arg(m_psoRandomSeed)); + } initializeSwarm(); emit logMessageGenerated(tr("Swarm initialized: %1 particles, %2 dimensions").arg(m_swarmSize).arg(getEnabledParameterCount())); @@ -3191,9 +3275,7 @@ bool nmCalculationAutoFitPSO::startAutoFitting() .arg(currentSuccessRate * 100, 0, 'f', 1).arg(m_currentIteration + 1)); } - // 更新全局最优。updateGlobalBest() 只读取粒子的真实 bestFitness， - // 不使用代理模型的 surrogateObjective。 - //double previousGlobalBest = m_globalBestFitness; + // 只从真实求解器确认的粒子 bestFitness 更新全局最优，不使用代理分数。 updateGlobalBest(); // 记录本代所有粒子的真实/代理误差和筛选决策，用于复盘和排障。 writeIterationTraceRows(); @@ -3282,11 +3364,40 @@ bool nmCalculationAutoFitPSO::startAutoFitting() } } + finalReason = analyzeOptimizationStatus(); + } else { + // 非代理路径不初始化粒子，也不使用 pbest/gbest 速度更新。 + finalReason = runTrustRegionFitting(); + } + // 最终结果验证和保护 validateAndProtectFinalResult(); + if(!useParticleSwarm && !m_globalBestPosition.isEmpty() && + m_globalBestObjectiveBreakdown.valid) { + // 精英保护可能恢复用户初始解，最终行必须在保护之后写入，确保 trace + // 中最后记录的就是实际回写 DataManager 的参数，而非最后一次接受候选。 + writeTraceRow(m_currentIteration, -1, + "trust_region_final", + m_globalBestPosition, + m_globalBestFitness, + m_globalBestFitness < 1.0e9, + -1, + std::numeric_limits::quiet_NaN(), + "final_result", + m_globalBestPosition, + m_globalBestFitness, + &m_globalBestObjectiveBreakdown); + } + + if(m_globalBestFitness < m_targetError) { + finalReason = PSO_TARGET_ACHIEVED; + } + } catch(const std::exception& e) { - m_lastError = QString(tr("Critical exception in PSO main loop: %1")).arg(e.what()); + m_lastError = useParticleSwarm + ? QString(tr("Critical exception in PSO main loop: %1")).arg(e.what()) + : QString(tr("Critical exception in automatic fitting: %1")).arg(e.what()); emit logMessageGenerated(tr("CRITICAL ERROR: %1").arg(e.what())); closeTraceFile(); cleanupTemporaryDirectory(); @@ -3294,7 +3405,9 @@ bool nmCalculationAutoFitPSO::startAutoFitting() emit fittingFinished(false, m_lastError); return false; } catch(...) { - m_lastError = QString(tr("Unknown critical exception in PSO main loop")); + m_lastError = useParticleSwarm + ? QString(tr("Unknown critical exception in PSO main loop")) + : QString(tr("Unknown critical exception in automatic fitting")); emit logMessageGenerated(tr("CRITICAL ERROR: Unknown exception in PSO main loop")); closeTraceFile(); cleanupTemporaryDirectory(); @@ -3373,43 +3486,68 @@ bool nmCalculationAutoFitPSO::startAutoFitting() // 判断系统确定最终结果 bool success; QString message; - StopReasonPSO finalReason = analyzeOptimizationStatus(); if(finalReason == PSO_TARGET_ACHIEVED) { success = true; message = QString(tr("Target achieved. Best error: %1, Iterations: %2")) .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1); - emit logMessageGenerated(tr("=== PSO OPTIMIZATION SUCCESSFUL ===")); + emit logMessageGenerated(useParticleSwarm + ? tr("=== PSO OPTIMIZATION SUCCESSFUL ===") + : tr("=== AUTOMATIC FITTING SUCCESSFUL ===")); } else if(finalReason == PSO_TRUE_CONVERGENCE) { success = true; - message = QString(tr("PSO optimization converged to stable solution. Best error: %1, Iterations: %2")) - .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1); - emit logMessageGenerated(tr("=== PSO OPTIMIZATION CONVERGED ===")); + message = useParticleSwarm + ? QString(tr("PSO optimization converged to stable solution. Best error: %1, Iterations: %2")) + .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1) + : QString(tr("Automatic fitting converged to a stable solution. Best error: %1, Iterations: %2")) + .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1); + emit logMessageGenerated(useParticleSwarm + ? tr("=== PSO OPTIMIZATION CONVERGED ===") + : tr("=== AUTOMATIC FITTING CONVERGED ===")); } else if(finalReason == PSO_LOCAL_OPTIMUM) { success = true; - message = QString(tr("PSO optimization trapped in local optimum. Best error: %1, Iterations: %2")) - .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1); - emit logMessageGenerated(tr("=== PSO OPTIMIZATION - LOCAL OPTIMUM ===")); + message = useParticleSwarm + ? QString(tr("PSO optimization trapped in local optimum. Best error: %1, Iterations: %2")) + .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1) + : QString(tr("Automatic fitting reached a local optimum. Best error: %1, Iterations: %2")) + .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1); + emit logMessageGenerated(useParticleSwarm + ? tr("=== PSO OPTIMIZATION - LOCAL OPTIMUM ===") + : tr("=== AUTOMATIC FITTING - LOCAL OPTIMUM ===")); } else if(finalReason == PSO_MAX_ITERATIONS) { success = true; message = QString(tr("Max iterations reached. Best error: %1, Iterations: %2")) .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1); - emit logMessageGenerated(tr("=== PSO OPTIMIZATION - MAX ITERATIONS ===")); + emit logMessageGenerated(useParticleSwarm + ? tr("=== PSO OPTIMIZATION - MAX ITERATIONS ===") + : tr("=== AUTOMATIC FITTING - MAX ITERATIONS ===")); } else if(finalReason == PSO_USER_STOPPED) { success = true; message = QString(tr("Best error: %1, Iterations: %2")) .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1); - emit logMessageGenerated(tr("=== PSO OPTIMIZATION STOPPED BY USER ===")); + emit logMessageGenerated(useParticleSwarm + ? tr("=== PSO OPTIMIZATION STOPPED BY USER ===") + : tr("=== AUTOMATIC FITTING STOPPED BY USER ===")); } else if(finalReason == PSO_CONSECUTIVE_FAILURES) { success = false; - message = QString(tr("PSO optimization failed due to consecutive failures. Best error: %1, Iterations: %2")) - .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1); - emit logMessageGenerated(tr("=== PSO OPTIMIZATION FAILED ===")); + message = useParticleSwarm + ? QString(tr("PSO optimization failed due to consecutive failures. Best error: %1, Iterations: %2")) + .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1) + : QString(tr("Automatic fitting failed due to consecutive failures. Best error: %1, Iterations: %2")) + .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1); + emit logMessageGenerated(useParticleSwarm + ? tr("=== PSO OPTIMIZATION FAILED ===") + : tr("=== AUTOMATIC FITTING FAILED ===")); } else { success = false; - message = QString(tr("PSO optimization ended unexpectedly. Best error: %1, Iterations: %2")) - .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1); - emit logMessageGenerated(tr("=== PSO OPTIMIZATION - UNKNOWN END ===")); + message = useParticleSwarm + ? QString(tr("PSO optimization ended unexpectedly. Best error: %1, Iterations: %2")) + .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1) + : QString(tr("Automatic fitting ended unexpectedly. Best error: %1, Iterations: %2")) + .arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1); + emit logMessageGenerated(useParticleSwarm + ? tr("=== PSO OPTIMIZATION - UNKNOWN END ===") + : tr("=== AUTOMATIC FITTING - UNKNOWN END ===")); } if(!finalFullSolverSucceeded) { @@ -3544,6 +3682,8 @@ void nmCalculationAutoFitPSO::initializeSwarm() particle.velocity.resize(dimensions); particle.bestPosition.resize(dimensions); particle.guideBestPosition.resize(dimensions); + particle.currentObjectiveBreakdown = AutoFitObjectiveBreakdown(); + particle.bestObjectiveBreakdown = AutoFitObjectiveBreakdown(); particle.bestFitness = 1e10; particle.guideBestObjective = 1e10; particle.guideBestFromSurrogate = false; @@ -3609,7 +3749,1201 @@ void nmCalculationAutoFitPSO::initializeSwarm() particle.bestPosition = particle.position; particle.guideBestPosition = particle.position; + + // 第一个粒子可能直接复用用户初始解，因此同步保存初始解的误差分解， + // 后续误差引导只能使用真实求解器确认过的 breakdown。 + if(i == 0 && m_hasValidUserSolution) { + particle.currentObjectiveBreakdown = m_userInitialObjectiveBreakdown; + particle.bestObjectiveBreakdown = m_userInitialObjectiveBreakdown; + } + } +} + +// 信赖域搜索统一在 [0, 1] 内部坐标工作。正值参数使用对数坐标，使内部相同步长 +// 表示近似相同的相对变化，避免 k、C、Ct、Cf 等跨数量级参数被线性尺度支配； +// skin 可为负数、Swi 的物理意义是线性比例，因此二者保持有界线性坐标。 +static bool useTrustRegionLogScale(int parameterIndex, double lower, double upper) +{ + return parameterIndex != 1 && parameterIndex != 7 && + lower > 0.0 && upper > lower; +} + +static double toTrustRegionCoordinate(double value, + int parameterIndex, + double lower, + double upper) +{ + // 所有进入优化器的物理值先投影到用户上下界，再转换成无量纲坐标。 + // 这样有限差分步长、信赖半径和参数间相关性可以在统一尺度上比较。 + value = qMax(lower, qMin(upper, value)); + + if(useTrustRegionLogScale(parameterIndex, lower, upper)) { + return (qLn(value) - qLn(lower)) / (qLn(upper) - qLn(lower)); + } + + return upper > lower ? (value - lower) / (upper - lower) : 0.0; +} + +static double fromTrustRegionCoordinate(double coordinate, + int parameterIndex, + double lower, + double upper) +{ + // 候选内部坐标先限制在 [0,1]，再执行上述映射的逆变换，保证写回 + // DataManager 的参数始终位于用户设置的物理范围内。 + coordinate = qMax(0.0, qMin(1.0, coordinate)); + + if(useTrustRegionLogScale(parameterIndex, lower, upper)) { + return qExp(qLn(lower) + coordinate * (qLn(upper) - qLn(lower))); } + + return lower + coordinate * (upper - lower); +} + +enum TrustRegionErrorComponent +{ + TRUST_REGION_VERTICAL_COMPONENT = 0, + TRUST_REGION_HORIZONTAL_COMPONENT, + TRUST_REGION_SHAPE_COMPONENT, + TRUST_REGION_TOTAL_COMPONENT +}; + +// 一次真实求解的完整快照。除了参数和总误差，还保存内部坐标、诊断分量和 +// 双对数曲线，因此拒绝候选后可以完整恢复上一个已接受工作点。 +struct TrustRegionEvaluation +{ + QVector parameters; + QVector coordinates; + AutoFitObjectiveBreakdown breakdown; + QVector > curve; + double fitness; + int elapsedMs; + bool valid; + + TrustRegionEvaluation() + : fitness(1.0e10) + , elapsedMs(-1) + , valid(false) + {} +}; + +// LM 只使用固定长度、全部有限的稳健残差。代理路径不会进入本搜索器。 +static bool trustRegionResidualsValid( + const AutoFitObjectiveBreakdown& breakdown) +{ + // 损失函数固定使用 50 个压力点和 50 个导数点。严格校验长度，避免 + // Jacobian 沿用旧维度后访问另一候选的短残差向量。 + if(!breakdown.valid || breakdown.residualVector.size() != 100) { + return false; + } + + for(int i = 0; i < breakdown.residualVector.size(); ++i) { + if(!isFiniteNumber(breakdown.residualVector[i])) { + return false; + } + } + return true; +} + +// 计算向量二范数的平方，避免在只比较能量或计算正规方程时反复开方。 +static double trustRegionSquaredNorm(const QVector& values) +{ + double sum = 0.0; + for(int i = 0; i < values.size(); ++i) { + sum += values[i] * values[i]; + } + return sum; +} + +// 计算同维向量内积；维度不一致表示局部模型无效，返回零让调用方放弃修正。 +static double trustRegionDotProduct(const QVector& left, + const QVector& right) +{ + if(left.size() != right.size()) { + return 0.0; + } + + double sum = 0.0; + for(int i = 0; i < left.size(); ++i) { + sum += left[i] * right[i]; + } + 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 AutoFitObjectiveBreakdown& 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; +} + +// 求解选中参数对应的阻尼正规方程。上下和左右诊断量保留方向；形状没有 +// 天然正负，因此使用 shapeLoss 对参数的局部导数。参数最多八维，使用带 +// 部分主元的高斯消元即可处理该小矩阵，并在主元退化时明确返回失败。 +static bool solveTrustRegionLinearSystem( + QVector > matrix, + QVector rightHandSide, + QVector* solution) +{ + if(!solution || matrix.isEmpty() || + matrix.size() != rightHandSide.size()) { + return false; + } + + const int size = matrix.size(); + for(int i = 0; i < size; ++i) { + if(matrix[i].size() != size) { + return false; + } + } + + for(int column = 0; column < size; ++column) { + int pivotRow = column; + double pivotMagnitude = qAbs(matrix[column][column]); + for(int row = column + 1; row < size; ++row) { + double magnitude = qAbs(matrix[row][column]); + if(magnitude > pivotMagnitude) { + pivotMagnitude = magnitude; + pivotRow = row; + } + } + if(pivotMagnitude <= 1.0e-14) { + return false; + } + + if(pivotRow != column) { + qSwap(matrix[pivotRow], matrix[column]); + qSwap(rightHandSide[pivotRow], rightHandSide[column]); + } + + for(int row = column + 1; row < size; ++row) { + double factor = matrix[row][column] / + matrix[column][column]; + matrix[row][column] = 0.0; + for(int nextColumn = column + 1; + nextColumn < size; ++nextColumn) { + matrix[row][nextColumn] -= + factor * matrix[column][nextColumn]; + } + rightHandSide[row] -= factor * rightHandSide[column]; + } + } + + solution->fill(0.0, size); + for(int row = size - 1; row >= 0; --row) { + double value = rightHandSide[row]; + for(int column = row + 1; column < size; ++column) { + value -= matrix[row][column] * (*solution)[column]; + } + double pivot = matrix[row][row]; + if(qAbs(pivot) <= 1.0e-14) { + return false; + } + (*solution)[row] = value / pivot; + if(!isFiniteNumber((*solution)[row])) { + return false; + } + } + return true; +} + +// 计算两个 Jacobian 列向量的绝对余弦相似度。接近 1 表示两个参数在当前 +// 工作点对曲线的影响几乎相同，联合调整容易产生不可辨识方向。 +static double trustRegionJacobianColumnCorrelation( + const QVector >& jacobian, + int leftColumn, + int rightColumn) +{ + double product = 0.0; + double leftNorm = 0.0; + double rightNorm = 0.0; + for(int row = 0; row < jacobian.size(); ++row) { + if(leftColumn >= jacobian[row].size() || + rightColumn >= jacobian[row].size()) { + return 1.0; + } + double left = jacobian[row][leftColumn]; + double right = jacobian[row][rightColumn]; + product += left * right; + leftNorm += left * left; + rightNorm += right * right; + } + + if(leftNorm <= 1.0e-20 || rightNorm <= 1.0e-20) { + return 0.0; + } + return qAbs(product) / qSqrt(leftNorm * rightNorm); +} + +// 每次接受一个真实候选后，使用满足最新割线条件的秩一修正更新完整残差 +// Jacobian。这样模型吸收了刚得到的真实变化，又不必立即逐参数重新试算。 +static void updateTrustRegionJacobian( + QVector >* jacobian, + const QVector& oldResidual, + const QVector& newResidual, + const QVector& coordinateStep) +{ + if(!jacobian || jacobian->size() != oldResidual.size() || + oldResidual.size() != newResidual.size()) { + return; + } + + double denominator = trustRegionSquaredNorm(coordinateStep); + if(denominator <= 1.0e-12) { + return; + } + + for(int row = 0; row < jacobian->size(); ++row) { + if((*jacobian)[row].size() != coordinateStep.size()) { + return; + } + + double predictedChange = 0.0; + for(int column = 0; column < coordinateStep.size(); ++column) { + predictedChange += + (*jacobian)[row][column] * coordinateStep[column]; + } + double correction = + (newResidual[row] - oldResidual[row] - predictedChange) / + denominator; + for(int column = 0; column < coordinateStep.size(); ++column) { + (*jacobian)[row][column] += + correction * coordinateStep[column]; + } + } +} + +// 对上下偏差、左右偏差和形状损失的梯度执行同样的割线秩一修正，使诊断 +// 选参模型与完整残差 Jacobian 保持在同一个已接受工作点。 +static void updateTrustRegionScalarGradient( + QVector* gradient, + double oldValue, + double newValue, + const QVector& coordinateStep) +{ + if(!gradient || gradient->size() != coordinateStep.size() || + !isFiniteNumber(oldValue) || !isFiniteNumber(newValue)) { + return; + } + + double denominator = trustRegionSquaredNorm(coordinateStep); + if(denominator <= 1.0e-12) { + return; + } + + double predictedChange = trustRegionDotProduct( + *gradient, coordinateStep); + double correction = + (newValue - oldValue - predictedChange) / denominator; + for(int i = 0; i < gradient->size(); ++i) { + (*gradient)[i] += correction * coordinateStep[i]; + } +} + +bool nmCalculationAutoFitPSO::evaluateTrustRegionPoint( + const QVector& parameters, + double* fitness, + AutoFitObjectiveBreakdown* breakdown, + QVector >* curve, + int* elapsedMs) +{ + if(!fitness || !breakdown || !curve || !elapsedMs || m_shouldStop) { + return false; + } + + // evaluateFitness() 会写入 DataManager 并调用真实求解器。这里统一统计 + // 真实评价次数和耗时，同时严格要求固定残差、诊断结构和结果曲线均有效。 + QTime timer; + timer.start(); + *fitness = evaluateFitness(parameters); + *elapsedMs = timer.elapsed(); + *breakdown = m_lastObjectiveBreakdown; + *curve = m_lastEvaluatedLogLogData; + ++m_totalEvaluations; + + bool valid = isFiniteNumber(*fitness) && *fitness < 1.0e9 && + breakdown->valid && + trustRegionResidualsValid(*breakdown) && + !curve->isEmpty(); + if(valid) { + ++m_successfulEvaluations; + } + + return valid; +} + +StopReasonPSO nmCalculationAutoFitPSO::runTrustRegionFitting() +{ + const int dimensions = getEnabledParameterCount(); + if(dimensions <= 0 || m_enabledParamIndices.size() != dimensions) { + m_lastError = tr("No valid parameters are available for trust-region fitting"); + return PSO_OPTIMIZATION_FAILED; + } + + // 真实求解次数比“外层迭代次数”更能反映耗时。预算至少允许完成一次全参数 + // 灵敏度和两次候选评价，同时避免连续重建 Jacobian 导致运行时间失控。 + const int maximumEvaluations = qMax( + m_totalEvaluations + dimensions + 2, + qMax(20, m_maxIterations * 3)); + // 下列步长均位于归一化内部坐标：0.04 表示参数范围的 4%，信赖半径 + // 限制一次联合移动的二范数，相关性门槛用于排除响应近乎共线的参数。 + const double sensitivityStep = 0.04; + const double minimumCoordinateStep = 1.0e-5; + const double minimumTrustRadius = 2.0e-3; + const double maximumTrustRadius = 0.30; + const double columnCorrelationLimit = 0.995; + const double diagnosisThreshold = 1.0e-5; + + // damping 是 LM 阻尼；拒绝或预测失准时增大，真实下降与预测一致时减小。 + // 两组累计量控制 Jacobian 重建，避免长期使用已偏离当前工作点的局部模型。 + double trustRadius = 0.12; + double damping = 1.0e-2; + int consecutiveRejectedSteps = 0; + int consecutiveSolverFailures = 0; + int acceptedSinceRebuild = 0; + double movementSinceRebuild = 0.0; + bool rebuildRequested = true; + bool modelRebuiltAtMinimumRadius = false; + StopReasonPSO stopReason = PSO_MAX_ITERATIONS; + + // jacobian 的行对应固定 100 维稳健残差，列对应用户勾选的参数。 + // 三个 gradient 单独描述诊断分量对参数的局部变化，只用于本轮选参。 + QVector > jacobian; + QVector verticalGradient(dimensions, 0.0); + QVector horizontalGradient(dimensions, 0.0); + QVector shapeGradient(dimensions, 0.0); + QVector jacobianColumnValid(dimensions, false); + + // 参数向量的顺序始终与 m_enabledParamIndices 一致，不能按完整参数索引 + // 直接访问；下面两个转换函数集中维护这层映射关系。 + auto coordinatesFromParameters = [&](const QVector& parameters) + -> QVector { + QVector coordinates(dimensions, 0.0); + for(int i = 0; i < dimensions; ++i) { + int parameterIndex = m_enabledParamIndices[i]; + coordinates[i] = toTrustRegionCoordinate( + parameters[i], parameterIndex, + m_parameterLower[parameterIndex], + m_parameterUpper[parameterIndex]); + } + return coordinates; + }; + + auto parametersFromCoordinates = [&](const QVector& coordinates) + -> QVector { + QVector parameters(dimensions, 0.0); + for(int i = 0; i < dimensions; ++i) { + int parameterIndex = m_enabledParamIndices[i]; + parameters[i] = fromTrustRegionCoordinate( + coordinates[i], parameterIndex, + m_parameterLower[parameterIndex], + m_parameterUpper[parameterIndex]); + } + return parameters; + }; + + auto restoreEvaluationState = [&](const TrustRegionEvaluation& evaluation) { + // evaluateFitness() 会把试算参数写入 DataManager。无论候选是否接受， + // 下一次计算前都恢复到唯一的已接受工作点，防止失败试算污染后续求解。 + applyParametersToDataManager(evaluation.parameters); + m_lastObjectiveBreakdown = evaluation.breakdown; + m_lastEvaluatedLogLogData = evaluation.curve; + }; + + // 只有真实总误差更小的工作点才能发布为全局最优；曲线和诊断快照必须 + // 与参数同步更新，防止界面显示或最终精英保护使用错配的数据。 + auto publishAcceptedPoint = [&](const TrustRegionEvaluation& evaluation) { + m_previousBestFitness = m_globalBestFitness; + m_globalBestPosition = evaluation.parameters; + m_globalBestFitness = evaluation.fitness; + m_globalBestObjectiveBreakdown = evaluation.breakdown; + m_globalBestLogLogData = evaluation.curve; + emit bestCurveUpdated(m_targetLogLogData, + m_globalBestLogLogData, + m_currentIteration + 1, + m_globalBestFitness); + }; + + auto processPauseAndStop = [&]() -> bool { + while(m_isPaused && !m_shouldStop) { + QApplication::processEvents(); + msleep(100); + } + QApplication::processEvents(); + return !m_shouldStop; + }; + + // current 始终代表唯一已接受工作点。优先复用启动阶段已经真实验证的 + // 用户初始解，避免在信赖域入口重复调用一次昂贵求解器。 + TrustRegionEvaluation current; + if(m_hasValidUserSolution && + m_globalBestPosition.size() == dimensions && + trustRegionResidualsValid(m_globalBestObjectiveBreakdown) && + !m_globalBestLogLogData.isEmpty()) { + current.parameters = m_globalBestPosition; + current.coordinates = coordinatesFromParameters(current.parameters); + current.breakdown = m_globalBestObjectiveBreakdown; + current.curve = m_globalBestLogLogData; + current.fitness = m_globalBestFitness; + current.elapsedMs = 0; + current.valid = true; + } else { + // 用户初始解无效时只做一次确定性的范围中点回退；所有正值参数在对数 + // 坐标取中点，避免线性中点过分偏向跨数量级范围的上界。 + current.coordinates.fill(0.5, dimensions); + current.parameters = parametersFromCoordinates(current.coordinates); + current.valid = evaluateTrustRegionPoint( + current.parameters, + ¤t.fitness, + ¤t.breakdown, + ¤t.curve, + ¤t.elapsedMs); + writeTraceRow(-1, -1, + "trust_region_midpoint", + current.parameters, + current.fitness, + current.valid, + current.elapsedMs, + std::numeric_limits::quiet_NaN(), + current.valid ? "midpoint_valid" : "midpoint_invalid", + QVector(), + 1.0e10, + current.valid ? ¤t.breakdown : nullptr); + if(!current.valid) { + m_lastError = tr("The initial solution and parameter-range midpoint are both invalid"); + return m_shouldStop + ? PSO_USER_STOPPED + : PSO_OPTIMIZATION_FAILED; + } + publishAcceptedPoint(current); + } + + restoreEvaluationState(current); + m_convergenceHistory.append(current.fitness); + emit logMessageGenerated( + tr("Trust-region initial error: %1; evaluation budget: %2") + .arg(current.fitness, 0, 'e', 4) + .arg(maximumEvaluations)); + + if(current.fitness < m_targetError) { + return PSO_TARGET_ACHIEVED; + } + + // 在同一个真实工作点逐参数做单边差分。首选可用空间更大的方向；只有该方向 + // 求解失败时才补算反方向，因此初次建模通常每个参数只增加一次真实求解。 + auto rebuildSensitivity = [&]() -> bool { + const TrustRegionEvaluation base = current; + const int residualCount = base.breakdown.residualVector.size(); + if(residualCount <= 0) { + return false; + } + + jacobian = QVector >( + residualCount, QVector(dimensions, 0.0)); + verticalGradient.fill(0.0, dimensions); + horizontalGradient.fill(0.0, dimensions); + shapeGradient.fill(0.0, dimensions); + jacobianColumnValid.fill(false, dimensions); + + TrustRegionEvaluation bestProbe; + int bestProbeColumn = -1; + double bestProbeDelta = 0.0; + // 差分步长不超过参数范围的 4%，信赖域收缩后同步减小，但保留 0.5% + // 下限，避免步长太小使求解器数值噪声淹没真实灵敏度。 + const double finiteDifferenceStep = qMin( + sensitivityStep, + qMax(5.0e-3, trustRadius * 0.5)); + + for(int column = 0; + column < dimensions && + m_totalEvaluations < maximumEvaluations && + processPauseAndStop(); + ++column) { + // 单边差分优先选择离边界空间更大的方向；首方向求解无效时才反向 + // 补算，因此正常情况下每个参数只消耗一次真实求解。 + double positiveRoom = 1.0 - base.coordinates[column]; + double negativeRoom = base.coordinates[column]; + double preferredSign = positiveRoom >= negativeRoom ? 1.0 : -1.0; + bool columnBuilt = false; + + for(int directionAttempt = 0; + directionAttempt < 2 && + !columnBuilt && + m_totalEvaluations < maximumEvaluations; + ++directionAttempt) { + double direction = directionAttempt == 0 + ? preferredSign : -preferredSign; + double availableRoom = direction > 0.0 + ? positiveRoom : negativeRoom; + double deltaMagnitude = qMin( + finiteDifferenceStep, availableRoom); + if(deltaMagnitude < minimumCoordinateStep) { + continue; + } + + TrustRegionEvaluation probe; + probe.coordinates = base.coordinates; + probe.coordinates[column] += direction * deltaMagnitude; + probe.parameters = parametersFromCoordinates(probe.coordinates); + probe.valid = evaluateTrustRegionPoint( + probe.parameters, + &probe.fitness, + &probe.breakdown, + &probe.curve, + &probe.elapsedMs); + + QString decision = probe.valid + ? "sensitivity_valid" + : (directionAttempt == 0 + ? "sensitivity_retry_opposite" + : "sensitivity_invalid"); + writeTraceRow(m_currentIteration, + column, + "trust_region_sensitivity", + probe.parameters, + probe.fitness, + probe.valid, + probe.elapsedMs, + std::numeric_limits::quiet_NaN(), + decision, + base.parameters, + base.fitness, + probe.valid ? &probe.breakdown : nullptr); + + if(!probe.valid) { + restoreEvaluationState(base); + continue; + } + + double delta = probe.coordinates[column] - + base.coordinates[column]; + if(qAbs(delta) < minimumCoordinateStep || + probe.breakdown.residualVector.size() != residualCount) { + restoreEvaluationState(base); + continue; + } + + // 第 column 列是固定残差向量相对内部参数坐标的有限差分： + // J[:,column] = (r_probe-r_base)/delta。 + for(int row = 0; row < residualCount; ++row) { + jacobian[row][column] = + (probe.breakdown.residualVector[row] - + base.breakdown.residualVector[row]) / delta; + } + + // 有符号诊断量只有在基点和试算点都可靠时才能计算方向梯度； + // shapeLoss 无方向可靠性标志，始终记录其局部变化率。 + if(base.breakdown.verticalReliable && + probe.breakdown.verticalReliable && + !base.breakdown.registrationAmbiguous && + !probe.breakdown.registrationAmbiguous) { + verticalGradient[column] = + (probe.breakdown.verticalCommonBias - + base.breakdown.verticalCommonBias) / delta; + } + if(base.breakdown.horizontalReliable && + probe.breakdown.horizontalReliable && + !base.breakdown.registrationAmbiguous && + !probe.breakdown.registrationAmbiguous) { + horizontalGradient[column] = + (probe.breakdown.horizontalPhysicalShift - + base.breakdown.horizontalPhysicalShift) / delta; + } + shapeGradient[column] = + (probe.breakdown.shapeLoss - + base.breakdown.shapeLoss) / delta; + jacobianColumnValid[column] = true; + columnBuilt = true; + + if(probe.fitness < base.fitness && + (!bestProbe.valid || + probe.fitness < bestProbe.fitness)) { + bestProbe = probe; + bestProbeColumn = column; + bestProbeDelta = delta; + } + restoreEvaluationState(base); + } + } + + int validColumnCount = 0; + for(int i = 0; i < jacobianColumnValid.size(); ++i) { + if(jacobianColumnValid[i]) { + ++validColumnCount; + } + } + if(validColumnCount == 0 || m_shouldStop) { + restoreEvaluationState(base); + return false; + } + + // 灵敏度试算本身若找到更优真实解也应保留。所有列先基于同一个 base + // 建完，再用该已知割线把 Jacobian 平移到新工作点，避免边算边移动基点。 + if(bestProbe.valid && bestProbeColumn >= 0) { + QVector acceptedStep(dimensions, 0.0); + acceptedStep[bestProbeColumn] = bestProbeDelta; + updateTrustRegionJacobian( + &jacobian, + base.breakdown.residualVector, + bestProbe.breakdown.residualVector, + acceptedStep); + if(base.breakdown.verticalReliable && + bestProbe.breakdown.verticalReliable) { + updateTrustRegionScalarGradient( + &verticalGradient, + base.breakdown.verticalCommonBias, + bestProbe.breakdown.verticalCommonBias, + acceptedStep); + } + if(base.breakdown.horizontalReliable && + bestProbe.breakdown.horizontalReliable) { + updateTrustRegionScalarGradient( + &horizontalGradient, + base.breakdown.horizontalPhysicalShift, + bestProbe.breakdown.horizontalPhysicalShift, + acceptedStep); + } + updateTrustRegionScalarGradient( + &shapeGradient, + base.breakdown.shapeLoss, + bestProbe.breakdown.shapeLoss, + acceptedStep); + + current = bestProbe; + publishAcceptedPoint(current); + restoreEvaluationState(current); + m_convergenceHistory.append(current.fitness); + writeTraceRow(m_currentIteration, + bestProbeColumn, + "trust_region_sensitivity_accept", + current.parameters, + current.fitness, + true, + 0, + std::numeric_limits::quiet_NaN(), + "accepted_cached_probe", + current.parameters, + current.fitness, + ¤t.breakdown); + emit logMessageGenerated( + tr("Sensitivity probe accepted: error reduced to %1") + .arg(current.fitness, 0, 'e', 4)); + } else { + restoreEvaluationState(current); + } + + acceptedSinceRebuild = 0; + movementSinceRebuild = 0.0; + consecutiveRejectedSteps = 0; + rebuildRequested = false; + // 若重建过程中接受了试算点，当前模型已通过割线平移而不是在新点完整 + // 重算；再遇到最小半径停滞时仍允许做一次真正的新点重建。 + modelRebuiltAtMinimumRadius = + trustRadius <= minimumTrustRadius * 1.01 && + !bestProbe.valid; + emit logMessageGenerated( + tr("Sensitivity model rebuilt: %1/%2 parameter columns valid") + .arg(validColumnCount) + .arg(dimensions)); + return true; + }; + + int completedIterations = 0; + for(int iteration = 0; + iteration < m_maxIterations && + m_totalEvaluations < maximumEvaluations && + !m_shouldStop; + ++iteration) { + m_currentIteration = iteration; + completedIterations = iteration + 1; + + if(!processPauseAndStop()) { + break; + } + if(rebuildRequested) { + if(!rebuildSensitivity()) { + stopReason = m_shouldStop + ? PSO_USER_STOPPED + : PSO_LOCAL_OPTIMUM; + break; + } + if(current.fitness < m_targetError) { + stopReason = PSO_TARGET_ACHIEVED; + break; + } + if(m_totalEvaluations >= maximumEvaluations) { + stopReason = PSO_MAX_ITERATIONS; + break; + } + } + + // 先确定当前最突出的可靠诊断误差，用其梯度回答“哪些参数最能改善 + // 当前问题”；实际 LM 方向仍由完整残差梯度和 Jacobian 共同计算。 + int dominantComponent = trustRegionDominantComponent( + current.breakdown, diagnosisThreshold); + const QVector* componentGradient = nullptr; + if(dominantComponent == TRUST_REGION_VERTICAL_COMPONENT) { + componentGradient = &verticalGradient; + } else if(dominantComponent == TRUST_REGION_HORIZONTAL_COMPONENT) { + componentGradient = &horizontalGradient; + } else if(dominantComponent == TRUST_REGION_SHAPE_COMPONENT) { + componentGradient = &shapeGradient; + } + + // 主目标采用 0.5*||r||^2，其对参数的梯度为 J^T*r。这里不再叠加 + // vertical/horizontal/shape，保证诊断分量不会改变真实接受目标。 + QVector totalGradient(dimensions, 0.0); + for(int column = 0; column < dimensions; ++column) { + if(!jacobianColumnValid[column]) { + continue; + } + for(int row = 0; row < jacobian.size(); ++row) { + totalGradient[column] += + jacobian[row][column] * + current.breakdown.residualVector[row]; + } + } + + // 每轮最多联合调整三个灵敏参数。按当前诊断梯度绝对值由大到小选取， + // 并剔除 Jacobian 响应过度共线的列，降低弱可辨识参数互相补偿的风险。 + QVector selectedColumns; + QVector alreadyConsidered(dimensions, false); + for(int selection = 0; selection < qMin(3, dimensions); ++selection) { + int bestColumn = -1; + double bestScore = 0.0; + for(int column = 0; column < dimensions; ++column) { + if(alreadyConsidered[column] || + !jacobianColumnValid[column]) { + continue; + } + + double score = componentGradient + ? qAbs((*componentGradient)[column]) + : qAbs(totalGradient[column]); + if(!isFiniteNumber(score) || score <= bestScore) { + continue; + } + + bool excessivelyCorrelated = false; + for(int selectedIndex = 0; + selectedIndex < selectedColumns.size(); + ++selectedIndex) { + if(trustRegionJacobianColumnCorrelation( + jacobian, + column, + selectedColumns[selectedIndex]) > + columnCorrelationLimit) { + excessivelyCorrelated = true; + break; + } + } + if(!excessivelyCorrelated) { + bestColumn = column; + bestScore = score; + } + } + if(bestColumn < 0 || bestScore <= 1.0e-12) { + break; + } + selectedColumns.append(bestColumn); + alreadyConsidered[bestColumn] = true; + } + + // 诊断梯度接近零时，说明该分量在当前局部无法可靠选参，退回完整残差 + // 梯度，但接受标准仍然只有真实 total，诊断值不会重复计入目标函数。 + if(selectedColumns.isEmpty() && componentGradient) { + dominantComponent = TRUST_REGION_TOTAL_COMPONENT; + componentGradient = nullptr; + alreadyConsidered.fill(false, dimensions); + for(int selection = 0; + selection < qMin(3, dimensions); + ++selection) { + int bestColumn = -1; + double bestScore = 0.0; + for(int column = 0; column < dimensions; ++column) { + if(alreadyConsidered[column] || + !jacobianColumnValid[column]) { + continue; + } + double score = qAbs(totalGradient[column]); + if(score <= bestScore) { + continue; + } + bool excessivelyCorrelated = false; + for(int selectedIndex = 0; + selectedIndex < selectedColumns.size(); + ++selectedIndex) { + if(trustRegionJacobianColumnCorrelation( + jacobian, + column, + selectedColumns[selectedIndex]) > + columnCorrelationLimit) { + excessivelyCorrelated = true; + break; + } + } + if(!excessivelyCorrelated) { + bestColumn = column; + bestScore = score; + } + } + if(bestColumn < 0 || bestScore <= 1.0e-12) { + break; + } + selectedColumns.append(bestColumn); + alreadyConsidered[bestColumn] = true; + } + } + + // 当前局部没有可用方向时先缩小半径并重建灵敏度；只有已经在最小 + // 半径完整重建后仍无方向，才把它判定为局部最优。 + if(selectedColumns.isEmpty()) { + if(trustRadius <= minimumTrustRadius * 1.01 && + modelRebuiltAtMinimumRadius) { + stopReason = PSO_LOCAL_OPTIMUM; + break; + } + trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5); + damping = qMin(1.0e8, damping * 4.0); + rebuildRequested = true; + continue; + } + + // 在选中参数子空间构造 LM 正规方程： + // (J^T*J + damping*diag(J^T*J))*step = -J^T*r。 + // 对角缩放使不同参数列的灵敏度量级差异不会直接改变阻尼强弱。 + const int selectedCount = selectedColumns.size(); + QVector > normalMatrix( + selectedCount, QVector(selectedCount, 0.0)); + QVector rightHandSide(selectedCount, 0.0); + for(int left = 0; left < selectedCount; ++left) { + int leftColumn = selectedColumns[left]; + rightHandSide[left] = -totalGradient[leftColumn]; + for(int right = 0; right < selectedCount; ++right) { + int rightColumn = selectedColumns[right]; + for(int row = 0; row < jacobian.size(); ++row) { + normalMatrix[left][right] += + jacobian[row][leftColumn] * + jacobian[row][rightColumn]; + } + } + double diagonalScale = qMax( + 1.0e-10, normalMatrix[left][left]); + normalMatrix[left][left] += damping * diagonalScale; + } + + QVector selectedStep; + bool solved = solveTrustRegionLinearSystem( + normalMatrix, rightHandSide, &selectedStep); + QVector coordinateStep(dimensions, 0.0); + if(solved) { + for(int i = 0; i < selectedCount; ++i) { + coordinateStep[selectedColumns[i]] = selectedStep[i]; + } + } + + double stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep)); + if(!solved || !isFiniteNumber(stepNorm) || + stepNorm < minimumCoordinateStep) { + // 正规方程退化时使用投影最速下降方向，仍只移动本轮已选择的参数。 + coordinateStep.fill(0.0, dimensions); + double gradientNormSquared = 0.0; + for(int i = 0; i < selectedCount; ++i) { + int column = selectedColumns[i]; + double stepDirection = -totalGradient[column]; + if((current.coordinates[column] <= minimumCoordinateStep && + stepDirection < 0.0) || + (current.coordinates[column] >= + 1.0 - minimumCoordinateStep && + stepDirection > 0.0)) { + stepDirection = 0.0; + } + coordinateStep[column] = stepDirection; + gradientNormSquared += stepDirection * stepDirection; + } + double gradientNorm = qSqrt(gradientNormSquared); + if(gradientNorm > minimumCoordinateStep) { + double scale = trustRadius / gradientNorm; + for(int i = 0; i < selectedCount; ++i) { + int column = selectedColumns[i]; + coordinateStep[column] *= scale; + } + } + stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep)); + } + + // LM 解只给出局部模型建议方向；若超出当前信赖半径，保持方向不变并 + // 等比例截短，避免一次试算离开 Jacobian 有效的局部区域。 + if(stepNorm > trustRadius && stepNorm > 0.0) { + double scale = trustRadius / stepNorm; + for(int i = 0; i < coordinateStep.size(); ++i) { + coordinateStep[i] *= scale; + } + } + + // 将 LM 步长投影到用户给定的参数范围，实际用于预测下降的也是投影后步长。 + QVector candidateCoordinates = current.coordinates; + for(int i = 0; i < dimensions; ++i) { + candidateCoordinates[i] = qBound( + 0.0, + current.coordinates[i] + coordinateStep[i], + 1.0); + coordinateStep[i] = candidateCoordinates[i] - + current.coordinates[i]; + } + stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep)); + + // 用线性模型 r_new ~= r_current + J*step 预测稳健残差，再用平方能量 + // 的下降量与真实候选下降量比较，作为调整阻尼和半径的依据。 + QVector predictedResidual = + current.breakdown.residualVector; + for(int row = 0; row < jacobian.size(); ++row) { + for(int column = 0; column < dimensions; ++column) { + predictedResidual[row] += + jacobian[row][column] * coordinateStep[column]; + } + } + double predictedReduction = 0.5 * + (trustRegionSquaredNorm(current.breakdown.residualVector) - + trustRegionSquaredNorm(predictedResidual)); + + // 无实际移动或模型预测不下降时没有必要调用昂贵求解器。将它按一次 + // 拒绝处理，并在连续发生后重建灵敏度，防止继续沿失效模型试算。 + if(stepNorm < minimumCoordinateStep || + !isFiniteNumber(predictedReduction) || + predictedReduction <= 1.0e-14) { + trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5); + damping = qMin(1.0e8, damping * 4.0); + ++consecutiveRejectedSteps; + if(consecutiveRejectedSteps >= 2) { + if(trustRadius <= minimumTrustRadius * 1.01 && + modelRebuiltAtMinimumRadius) { + stopReason = PSO_LOCAL_OPTIMUM; + break; + } + rebuildRequested = true; + } + continue; + } + + TrustRegionEvaluation candidate; + candidate.coordinates = candidateCoordinates; + candidate.parameters = parametersFromCoordinates(candidate.coordinates); + candidate.valid = evaluateTrustRegionPoint( + candidate.parameters, + &candidate.fitness, + &candidate.breakdown, + &candidate.curve, + &candidate.elapsedMs); + + if(!candidate.valid) { + // 求解失败的候选不能改变 current。先完整恢复上一个已接受参数和 + // 对应误差快照，再缩小信赖域；连续失败达到上限才终止整个拟合。 + ++consecutiveSolverFailures; + ++consecutiveRejectedSteps; + trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5); + damping = qMin(1.0e8, damping * 4.0); + writeTraceRow(m_currentIteration, + -1, + "trust_region_candidate", + candidate.parameters, + candidate.fitness, + false, + candidate.elapsedMs, + std::numeric_limits::quiet_NaN(), + "solver_invalid", + current.parameters, + current.fitness, + nullptr); + restoreEvaluationState(current); + if(consecutiveRejectedSteps >= 2) { + rebuildRequested = true; + } + if(consecutiveSolverFailures >= m_maxConsecutiveFailures) { + stopReason = PSO_CONSECUTIVE_FAILURES; + break; + } + continue; + } + + consecutiveSolverFailures = 0; + // 有效候选即使最终被拒绝，也提供了一条真实割线，可用于修正下一轮 + // 局部模型；是否成为新工作点仍只由下面的 total 严格比较决定。 + const AutoFitObjectiveBreakdown oldBreakdown = current.breakdown; + updateTrustRegionJacobian( + &jacobian, + oldBreakdown.residualVector, + candidate.breakdown.residualVector, + coordinateStep); + if(oldBreakdown.verticalReliable && + candidate.breakdown.verticalReliable && + !oldBreakdown.registrationAmbiguous && + !candidate.breakdown.registrationAmbiguous) { + updateTrustRegionScalarGradient( + &verticalGradient, + oldBreakdown.verticalCommonBias, + candidate.breakdown.verticalCommonBias, + coordinateStep); + } + if(oldBreakdown.horizontalReliable && + candidate.breakdown.horizontalReliable && + !oldBreakdown.registrationAmbiguous && + !candidate.breakdown.registrationAmbiguous) { + updateTrustRegionScalarGradient( + &horizontalGradient, + oldBreakdown.horizontalPhysicalShift, + candidate.breakdown.horizontalPhysicalShift, + coordinateStep); + } + updateTrustRegionScalarGradient( + &shapeGradient, + oldBreakdown.shapeLoss, + candidate.breakdown.shapeLoss, + coordinateStep); + + // reductionRatio 衡量局部线性模型的可信度：接近 1 表示预测准确； + // 值较小表示虽然可能下降，但模型低估了非线性，需要收紧下一步。 + double actualReduction = 0.5 * + (current.fitness * current.fitness - + candidate.fitness * candidate.fitness); + double reductionRatio = actualReduction / predictedReduction; + bool accepted = candidate.fitness < current.fitness; + QString componentName = trustRegionComponentName(dominantComponent); + + if(accepted) { + // 真实总误差下降后才正式替换 current，并同步发布参数、曲线和诊断。 + // 模型预测可靠时减小阻尼并可扩大半径，预测较差时保守收缩。 + current = candidate; + publishAcceptedPoint(current); + restoreEvaluationState(current); + ++acceptedSinceRebuild; + movementSinceRebuild += stepNorm; + consecutiveRejectedSteps = 0; + m_convergenceHistory.append(current.fitness); + + if(reductionRatio > 0.75) { + damping = qMax(1.0e-8, damping * 0.5); + if(stepNorm >= trustRadius * 0.8) { + trustRadius = qMin( + maximumTrustRadius, trustRadius * 1.6); + } + } else if(reductionRatio > 0.25) { + damping = qMax(1.0e-8, damping * 0.8); + } else { + damping = qMin(1.0e8, damping * 2.0); + trustRadius = qMax( + minimumTrustRadius, trustRadius * 0.75); + } + + if(acceptedSinceRebuild >= 6 || + movementSinceRebuild >= 0.30) { + rebuildRequested = true; + } + modelRebuiltAtMinimumRadius = false; + } else { + // 拒绝时 candidate 只保留在 trace 中，DataManager 和内存状态都恢复 + // 到 current。连续拒绝说明割线模型可能失真，因此请求重新试算灵敏度。 + ++consecutiveRejectedSteps; + damping = qMin(1.0e8, damping * 4.0); + trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5); + restoreEvaluationState(current); + if(consecutiveRejectedSteps >= 2) { + rebuildRequested = true; + } + } + + writeTraceRow(m_currentIteration, + -1, + "trust_region_candidate", + candidate.parameters, + candidate.fitness, + true, + candidate.elapsedMs, + std::numeric_limits::quiet_NaN(), + accepted + ? QString("accepted_%1").arg(componentName) + : QString("rejected_%1").arg(componentName), + current.parameters, + current.fitness, + &candidate.breakdown); + + emit logMessageGenerated( + tr("Iteration %1: focus=%2, parameters=%3, error=%4, result=%5") + .arg(iteration + 1) + .arg(componentName) + .arg(selectedColumns.size()) + .arg(candidate.fitness, 0, 'e', 4) + .arg(accepted ? tr("accepted") : tr("rejected"))); + emit progressUpdated(iteration + 1, m_globalBestFitness); + + if(current.fitness < m_targetError) { + stopReason = PSO_TARGET_ACHIEVED; + break; + } + if(trustRadius <= minimumTrustRadius * 1.01 && + consecutiveRejectedSteps >= 2) { + if(modelRebuiltAtMinimumRadius) { + stopReason = PSO_LOCAL_OPTIMUM; + break; + } + rebuildRequested = true; + } + } + + if(completedIterations > 0) { + m_currentIteration = completedIterations - 1; + } + restoreEvaluationState(current); + + if(m_shouldStop) { + return PSO_USER_STOPPED; + } + if(current.fitness < m_targetError) { + return PSO_TARGET_ACHIEVED; + } + if(stopReason == PSO_CONSECUTIVE_FAILURES || + stopReason == PSO_LOCAL_OPTIMUM || + stopReason == PSO_OPTIMIZATION_FAILED) { + return stopReason; + } + return PSO_MAX_ITERATIONS; } void nmCalculationAutoFitPSO::updateParticle(int particleIndex) @@ -3647,6 +4981,7 @@ void nmCalculationAutoFitPSO::updateParticle(int particleIndex) // 直接复用真实误差和曲线，避免一次重复 DLL 调用且不改变 PSO 数学状态。 particle.fitness = m_userInitialFitness; particle.currentLogLogData = m_userInitialLogLogData; + particle.currentObjectiveBreakdown = m_userInitialObjectiveBreakdown; particle.lastEvaluationElapsedMs = 0; particle.evaluatedThisIteration = true; particle.lastEvaluationSuccess = true; @@ -3657,6 +4992,7 @@ void nmCalculationAutoFitPSO::updateParticle(int particleIndex) evalTimer.start(); particle.fitness = evaluateFitness(particle.position); particle.currentLogLogData = m_lastEvaluatedLogLogData; + particle.currentObjectiveBreakdown = m_lastObjectiveBreakdown; particle.lastEvaluationElapsedMs = evalTimer.elapsed(); particle.evaluatedThisIteration = true; particle.lastEvaluationSuccess = (particle.fitness < 1e9); @@ -3698,6 +5034,7 @@ void nmCalculationAutoFitPSO::updateParticle(int particleIndex) particle.bestFitness = particle.fitness; particle.bestPosition = particle.position; particle.bestLogLogData = particle.currentLogLogData; + particle.bestObjectiveBreakdown = particle.currentObjectiveBreakdown; if(!preserveSurrogateGuide) { particle.guideBestPosition = particle.position; @@ -3746,6 +5083,7 @@ void nmCalculationAutoFitPSO::updateGlobalBest() m_globalBestFitness = particle.bestFitness; m_globalBestPosition = particle.bestPosition; m_globalBestLogLogData = particle.bestLogLogData; + m_globalBestObjectiveBreakdown = particle.bestObjectiveBreakdown; globalBestUpdated = true; emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, m_currentIteration + 1, m_globalBestFitness); } @@ -3765,6 +5103,7 @@ void nmCalculationAutoFitPSO::updateGlobalBest() m_globalBestFitness = m_userInitialFitness; m_globalBestPosition = m_userInitialSolution; m_globalBestLogLogData = m_userInitialLogLogData; + m_globalBestObjectiveBreakdown = m_userInitialObjectiveBreakdown; emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, m_currentIteration + 1, m_globalBestFitness); } } @@ -3841,6 +5180,7 @@ void nmCalculationAutoFitPSO::updateVelocityAndPosition() int paramIndex = m_enabledParamIndices[j]; double range = m_parameterUpper[paramIndex] - m_parameterLower[paramIndex]; double maxVel = range * VELOCITY_LIMIT_FACTOR; + particle.velocity[j] = qMax(-maxVel, qMin(maxVel, particle.velocity[j])); // 更新位置 @@ -4650,10 +5990,8 @@ void nmCalculationAutoFitPSO::saveOptimizationResult() void nmCalculationAutoFitPSO::validateAndProtectFinalResult() { - // 最终精英保护。 - // PSO 是随机启发式算法，某些工况下可能没有找到比用户初始模型更好的结果。 - // 这里用初始解误差与最终全局最优误差做比较，如果改进不足，就恢复初始解， - // 避免“自动拟合”把已有模型调坏。 + // 最终精英保护只阻止无效结果或真正变差的结果。任何真实误差下降都应保留， + // 不能再用固定百分比门槛把已经找到的更优解恢复成初始值。 if(!m_hasValidUserSolution) { emit logMessageGenerated(tr("No initial solution for elite protection")); return; @@ -4668,28 +6006,36 @@ void nmCalculationAutoFitPSO::validateAndProtectFinalResult() emit logMessageGenerated(tr("Comparing results: Initial=%1, Final=%2") .arg(initialFitness, 0, 'e', 4).arg(finalFitness, 0, 'e', 4)); - // 计算改进程度 - double improvement = initialFitness - finalFitness; - double relativeImprovement = improvement / qMax(1e-10, qAbs(initialFitness)); - - emit logMessageGenerated(tr("Improvement: %1 (%2%)") - .arg(improvement, 0, 'e', 4).arg(relativeImprovement * 100, 0, 'f', 2)); + bool finalValid = isFiniteNumber(finalFitness) && + finalFitness < 1.0e9 && + m_globalBestPosition.size() == + m_userInitialSolution.size() && + !m_globalBestLogLogData.isEmpty() && + m_globalBestObjectiveBreakdown.valid; + if(finalValid) { + double improvement = initialFitness - finalFitness; + double relativeImprovement = + improvement / qMax(1.0e-10, qAbs(initialFitness)); + emit logMessageGenerated(tr("Improvement: %1 (%2%)") + .arg(improvement, 0, 'e', 4) + .arg(relativeImprovement * 100, 0, 'f', 2)); + } - if(relativeImprovement < m_improvementThreshold) { - emit logMessageGenerated(tr("Elite protection triggered: insufficient improvement")); - emit logMessageGenerated(tr("Threshold: %1%, Actual: %2%") - .arg(m_improvementThreshold * 100, 0, 'f', 2) - .arg(relativeImprovement * 100, 0, 'f', 4)); + if(!finalValid || finalFitness > initialFitness) { + emit logMessageGenerated( + tr("Elite protection triggered: final result is invalid or worse than initial")); emit logMessageGenerated(tr("Restoring initial solution as final result")); m_globalBestFitness = initialFitness; m_globalBestPosition = m_userInitialSolution; m_globalBestLogLogData = m_userInitialLogLogData; + m_globalBestObjectiveBreakdown = m_userInitialObjectiveBreakdown; emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, m_currentIteration + 1, m_globalBestFitness); emit logMessageGenerated(tr("Initial solution restored successfully")); } else { - emit logMessageGenerated(tr("Final result validated - significant improvement achieved")); + emit logMessageGenerated( + tr("Final result validated - solution is not worse than initial")); } } @@ -4997,17 +6343,16 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError( const QVector >& target, const QVector >& result) const { - // 这里只负责“曲线比较和误差诊断”，不根据诊断结果直接修改任何拟合参数。 - // 调用方可以读取 m_lastObjectiveBreakdown 做诊断或展示；本函数本身不修改参数。 - // 在统一的对数时间网格上计算压力和导数残差，并拆分为上下、左右、形状误差。 + // 主目标只比较固定网格上的压力和导数残差；上下、左右和形状只负责诊断 + // 误差来源和选择参数，避免同一残差在 total 中被重复计算。整个计算过程均 + // 位于 log(time)-log(value) 坐标，因此得到的是相对尺度偏差而非原始压力量纲。 const double invalidLoss = 1.0e10; - const double valueFloor = 1.0e-12; // 导数接近零时的对数下限，避免 log(0)。 - const double huberDelta = qLn(1.2); // 约对应 20% 的相对偏差拐点。 - const double minimumCoverage = 0.95; // 点数比例和连续跨度比例都至少接近 95%。 - const int numPoints = 50; // 固定网格使不同候选的损失具有可比性。 - // 目标函数的主排序项为 0.5*pressureLoss + 0.5*derivativeLoss； - // coveragePenalty 只在接近覆盖边界时提供连续惩罚，上下、左右和形状分量 - // 会写入 m_lastObjectiveBreakdown，供后续按误差类型选择参数。 + const double valueFloor = 1.0e-12; + // Huber 转折点 qLn(1.2) 对应约 20% 的倍率偏差；小偏差保持平方惩罚， + // 更大的局部尖峰改为近似线性惩罚，避免单点支配整条曲线。 + const double huberDelta = qLn(1.2); + const double minimumCoverage = 0.95; + const int numPoints = 50; m_lastObjectiveBreakdown = AutoFitObjectiveBreakdown(); if(!validateLogLogData(target) || !validateLogLogData(result)) { @@ -5015,9 +6360,8 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError( } try { - // 清洗曲线并拆成压力、导数两条曲线。导数可以为负，所以统一使用绝对值 - // 进入双对数空间；时间和压力必须为正，否则无法进行对数插值。这里的清洗 - // 只丢弃无法比较的采样点，不改变原始曲线或求解器输出。 + // 无法比较的采样行先跳过；有限但非正的导数无法进入双对数空间， + // 当前数据又没有逐点有效掩码，因此遇到这种导数时判本次评价无效。 auto prepareCurve = [valueFloor](const QVector >& data, QVector* pressure, QVector* derivative) -> bool { @@ -5028,37 +6372,41 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError( } for(int i = 0; i < data[0].size(); ++i) { - if(!isFiniteNumber(data[0][i]) || !isFiniteNumber(data[1][i]) || - !isFiniteNumber(data[2][i]) || data[0][i] <= 0.0 || + if(!isFiniteNumber(data[0][i]) || + !isFiniteNumber(data[1][i]) || + !isFiniteNumber(data[2][i]) || + data[0][i] <= 0.0 || data[1][i] <= 0.0) { continue; } + if(data[2][i] <= 0.0) { + return false; + } pressure->append(QPointF(data[0][i], data[1][i])); - derivative->append(QPointF(data[0][i], - qMax(qAbs(data[2][i]), valueFloor))); + derivative->append( + QPointF(data[0][i], qMax(data[2][i], valueFloor))); } - // 插值要求时间严格递增。重复时间点保留排序后的最后一个值， - // 避免重复横坐标导致对数插值分母为零。压力和导数分别去重， - // 这样即使某条曲线存在重复时间点，也不会污染另一条曲线的插值。 + // 求解器输出可能不是严格升序，且同一时刻可能出现重复记录。 + // 插值前统一排序并让后出现的记录覆盖同时间旧值，保证横坐标严格递增。 auto sortAndUnique = [](QVector* curve) { - std::stable_sort(curve->begin(), curve->end(), - [](const QPointF& left, const QPointF& right) { - return left.x() < right.x(); - }); + std::stable_sort( + curve->begin(), curve->end(), + [](const QPointF& left, const QPointF& right) { + return left.x() < right.x(); + }); QVector unique; unique.reserve(curve->size()); - for(int i = 0; i < curve->size(); ++i) { - if(unique.isEmpty() || curve->at(i).x() > unique.last().x()) { + if(unique.isEmpty() || + curve->at(i).x() > unique.last().x()) { unique.append(curve->at(i)); } else { unique[unique.size() - 1] = curve->at(i); } } - *curve = unique; }; @@ -5071,44 +6419,77 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError( QVector targetDerivative; QVector resultPressure; QVector resultDerivative; - if(!prepareCurve(target, &targetPressure, &targetDerivative) || !prepareCurve(result, &resultPressure, &resultDerivative)) { return invalidLoss; } - // 在 log(time)-log(value) 空间做线性插值，而不是在线性坐标直接插值。 - // 这样可以保持双对数曲线的时间尺度和数量级特征；返回值是 - // log(abs(value))，后续残差因此可以直接解释为相对幅值误差。 - auto interpolateLogValue = [valueFloor](const QVector& curve, - double x, - double* value) -> bool { - if(!value || curve.size() < 2 || x < curve.first().x() || - x > curve.last().x() || x <= 0.0) { + // 在双对数坐标中插值。二分定位用于后面的多次水平配准试算。 + auto interpolateLogValue = [valueFloor]( + const QVector& curve, + double x, + double* value) -> bool { + if(!value || curve.size() < 2 || x <= 0.0 || + x < curve.first().x() || x > curve.last().x()) { return false; } - int right = 1; - - while(right < curve.size() && curve[right].x() < x) { - ++right; + int low = 0; + int high = curve.size() - 1; + while(low < high) { + int middle = low + (high - low) / 2; + if(curve[middle].x() < x) { + low = middle + 1; + } else { + high = middle; + } } - right = qMin(right, curve.size() - 1); - int left = qMax(0, right - 1); + int right = qBound(1, low, curve.size() - 1); + int left = right - 1; double leftLogX = qLn(curve[left].x()); double rightLogX = qLn(curve[right].x()); double denominator = rightLogX - leftLogX; - double leftLogY = qLn(qMax(qAbs(curve[left].y()), valueFloor)); - double rightLogY = qLn(qMax(qAbs(curve[right].y()), valueFloor)); + double leftLogY = + qLn(qMax(qAbs(curve[left].y()), valueFloor)); + double rightLogY = + qLn(qMax(qAbs(curve[right].y()), valueFloor)); if(qAbs(denominator) <= 1.0e-12) { *value = leftLogY; } else { double ratio = (qLn(x) - leftLogX) / denominator; - *value = leftLogY + ratio * (rightLogY - leftLogY); + *value = leftLogY + + ratio * (rightLogY - leftLogY); } + return isFiniteNumber(*value); + }; + // 覆盖率通过后若只缺少首尾少量点，用模拟曲线自身的端点斜率作短距离 + // 双对数外推。该外推只用于主损失的固定网格，不参与水平配准搜索。 + auto extrapolateEndpointLogValue = [valueFloor]( + const QVector& curve, + double x, + double* value) -> bool { + if(!value || curve.size() < 2 || x <= 0.0) { + return false; + } + int left = x < curve.first().x() + ? 0 + : curve.size() - 2; + int right = left + 1; + double leftLogX = qLn(curve[left].x()); + double rightLogX = qLn(curve[right].x()); + double denominator = rightLogX - leftLogX; + if(qAbs(denominator) <= 1.0e-12) { + return false; + } + double leftLogY = + qLn(qMax(qAbs(curve[left].y()), valueFloor)); + double rightLogY = + qLn(qMax(qAbs(curve[right].y()), valueFloor)); + double ratio = (qLn(x) - leftLogX) / denominator; + *value = leftLogY + ratio * (rightLogY - leftLogY); return isFiniteNumber(*value); }; @@ -5116,14 +6497,11 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError( const double targetMaxX = targetPressure.last().x(); const double resultMinX = resultPressure.first().x(); const double resultMaxX = resultPressure.last().x(); - if(targetMinX <= 0.0 || targetMaxX <= targetMinX || resultMinX <= 0.0 || resultMaxX <= resultMinX) { return invalidLoss; } - // 目标曲线的完整时间范围作为统一比较区间。若候选曲线覆盖不足， - // 后面的 coverage 检查会拒绝它，防止候选通过缩短时间范围来降低误差。 QVector commonX(numPoints); QVector commonLogX(numPoints); QVector targetLogPressure(numPoints); @@ -5131,34 +6509,34 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError( const double targetLogMinX = qLn(targetMinX); const double targetLogMaxX = qLn(targetMaxX); + // 固定使用目标曲线的完整 log-time 网格，候选之间不会因采样点不同而失去可比性。 for(int i = 0; i < numPoints; ++i) { double logX = targetLogMinX + static_cast(i) * - (targetLogMaxX - targetLogMinX) / (numPoints - 1); + (targetLogMaxX - targetLogMinX) / + (numPoints - 1); commonLogX[i] = logX; commonX[i] = qExp(logX); - if(!interpolateLogValue(targetPressure, commonX[i], - &targetLogPressure[i]) || - !interpolateLogValue(targetDerivative, commonX[i], - &targetLogDerivative[i])) { + if(!interpolateLogValue( + targetPressure, commonX[i], + &targetLogPressure[i]) || + !interpolateLogValue( + targetDerivative, commonX[i], + &targetLogDerivative[i])) { return invalidLoss; } } - // 残差采用“模拟减目标”，因此正值表示模拟曲线在对数幅值上高于目标， - // 负值表示模拟曲线偏低。NaN 表示该网格点不在模拟曲线支持范围内， - // 后续统计会自动跳过，但覆盖率检查仍会限制候选不能靠缺失数据降低损失。 - QVector pressureResidual(numPoints, - std::numeric_limits::quiet_NaN()); - QVector derivativeResidual(numPoints, - std::numeric_limits::quiet_NaN()); - QVector pressureSlope(numPoints, 0.0); - QVector derivativeSlope(numPoints, 0.0); + QVector pressureResidual( + numPoints, std::numeric_limits::quiet_NaN()); + QVector derivativeResidual( + numPoints, std::numeric_limits::quiet_NaN()); int firstSupported = -1; int lastSupported = -1; int supportedCount = 0; + // 残差定义为“模拟减目标”：正值表示模拟曲线偏高，负值表示偏低。 for(int i = 0; i < numPoints; ++i) { if(commonX[i] < resultMinX || commonX[i] > resultMaxX) { continue; @@ -5166,87 +6544,95 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError( double resultLogPressure = 0.0; double resultLogDerivative = 0.0; - - if(!interpolateLogValue(resultPressure, commonX[i], - &resultLogPressure) || - !interpolateLogValue(resultDerivative, commonX[i], - &resultLogDerivative)) { - continue; + if(!interpolateLogValue( + resultPressure, commonX[i], + &resultLogPressure) || + !interpolateLogValue( + resultDerivative, commonX[i], + &resultLogDerivative)) { + // 已位于结果时间范围内却无法插值说明数据存在内部断点，不能补线。 + return invalidLoss; } - pressureResidual[i] = resultLogPressure - targetLogPressure[i]; - derivativeResidual[i] = resultLogDerivative - targetLogDerivative[i]; + pressureResidual[i] = + resultLogPressure - targetLogPressure[i]; + derivativeResidual[i] = + resultLogDerivative - targetLogDerivative[i]; ++supportedCount; - if(firstSupported < 0) { firstSupported = i; } - lastSupported = i; } - // 同时使用点覆盖率和连续时间跨度覆盖率，避免只覆盖少数离散点也被判定为完整。 + // 覆盖率同时约束“有效点数量”和“连续时间跨度”。取两者较小值可避免 + // 点数很多但只集中在局部时段的候选被误认为覆盖充分。 AutoFitObjectiveBreakdown breakdown; breakdown.coverage = supportedCount > 0 - ? static_cast(supportedCount) / numPoints + ? static_cast(supportedCount) / + numPoints : 0.0; - if(firstSupported >= 0 && lastSupported >= firstSupported) { - double span = qMax(1.0e-12, targetLogMaxX - targetLogMinX); - double coveredSpan = commonLogX[lastSupported] - - commonLogX[firstSupported]; - breakdown.coverage = qMin(breakdown.coverage, - qMax(0.0, coveredSpan / span)); + double targetSpan = + qMax(1.0e-12, targetLogMaxX - targetLogMinX); + double coveredSpan = + commonLogX[lastSupported] - commonLogX[firstSupported]; + breakdown.coverage = qMin( + breakdown.coverage, + qMax(0.0, coveredSpan / targetSpan)); } - // 覆盖率越接近 1，惩罚越小；覆盖不足 minimumCoverage 时直接返回无效损失。 - double coverageGap = qMax(0.0, 1.0 - breakdown.coverage); - breakdown.coveragePenalty = - qPow(coverageGap / (1.0 - minimumCoverage), 2.0); - + // 覆盖率只判断候选是否有效，不再加入固定惩罚，避免所有误差被整体抬高。 if(breakdown.coverage < minimumCoverage) { breakdown.total = invalidLoss; m_lastObjectiveBreakdown = breakdown; return invalidLoss; } - // 目标曲线斜率用于把“残差随时间的系统性变化”解释为左右平移。 - // 斜率用对数坐标计算，与前面的插值空间保持一致；平坦区斜率接近零， - // 不会凭空制造水平偏移量。 + // 通过门槛后最多只缺少首尾少量目标点。按模拟曲线端点趋势补齐后， + // 每个候选仍在固定 50 点上计算 Huber 均值，不能靠少算难拟合端点获益。 for(int i = 0; i < numPoints; ++i) { - int left = i == 0 ? 0 : i - 1; - int right = i == numPoints - 1 ? numPoints - 1 : i + 1; - double denominator = commonLogX[right] - commonLogX[left]; - - if(qAbs(denominator) > 1.0e-12) { - pressureSlope[i] = - (targetLogPressure[right] - targetLogPressure[left]) / - denominator; - derivativeSlope[i] = - (targetLogDerivative[right] - targetLogDerivative[left]) / - denominator; - } - } - - // Huber RMS 在小残差区域保持平方损失，在异常点区域转为线性增长， - // 避免少量求解器异常点完全主导候选排序。这里没有除以目标值， - // 因为残差已经是 log(value) 差值，本身就是相对误差的表达；返回值是 - // Huber rho 均值的平方根，保持与 RMS 类似的尺度。 - auto huberRms = [huberDelta](const QVector& values, - int begin, - int end) -> double { + if(isFiniteNumber(pressureResidual[i]) && + isFiniteNumber(derivativeResidual[i])) { + continue; + } + double resultLogPressure = 0.0; + double resultLogDerivative = 0.0; + if(!extrapolateEndpointLogValue( + resultPressure, commonX[i], + &resultLogPressure) || + !extrapolateEndpointLogValue( + resultDerivative, commonX[i], + &resultLogDerivative)) { + return invalidLoss; + } + pressureResidual[i] = + resultLogPressure - targetLogPressure[i]; + derivativeResidual[i] = + resultLogDerivative - targetLogDerivative[i]; + } + + // Huber RMS 在小残差处保持平方损失，在异常点处转为线性增长。 + auto huberRmsAround = [huberDelta]( + const QVector& values, + int begin, + int end, + double center) -> double { double sum = 0.0; int count = 0; + int validBegin = qMax(0, begin); + int validEnd = + qMin(end, static_cast(values.size())); - for(int i = qMax(0, begin); - i < qMin(end, static_cast(values.size())); ++i) { + for(int i = validBegin; i < validEnd; ++i) { if(!isFiniteNumber(values[i])) { continue; } - double absoluteValue = qAbs(values[i]); + double difference = values[i] - center; + double absoluteValue = qAbs(difference); double rho = absoluteValue <= huberDelta - ? values[i] * values[i] + ? difference * difference : 2.0 * huberDelta * absoluteValue - huberDelta * huberDelta; sum += rho; @@ -5258,41 +6644,61 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError( : std::numeric_limits::quiet_NaN(); }; - // 用 Huber 加权迭代估计残差中心，作为整体上下偏移。相比普通平均值， - // 它对局部尖峰更稳健，同时保留偏高/偏低的方向信息。迭代只用于诊断， - // 不会把残差“校正”后再写回求解器结果。 - auto huberCenter = [huberDelta](const QVector& values) -> double { + auto huberRms = [&huberRmsAround]( + const QVector& values, + int begin, + int end) -> double { + return huberRmsAround(values, begin, end, 0.0); + }; + + // 将 Huber 能量转换成带符号的等效残差，使向量平方和与稳健损失一致。 + // 非代理 LM 直接对该固定长度向量建立 Jacobian。 + auto huberEquivalentResidual = [huberDelta](double value) -> double { + double absoluteValue = qAbs(value); + double rho = absoluteValue <= huberDelta + ? value * value + : 2.0 * huberDelta * absoluteValue - + huberDelta * huberDelta; + double magnitude = qSqrt(qMax(0.0, rho)); + return value < 0.0 ? -magnitude : magnitude; + }; + + // Huber 加权中心保留上下偏差的符号，并降低局部尖峰的影响。 + auto huberCenterRange = [huberDelta]( + const QVector& values, + int begin, + int end) -> double { + int validBegin = qMax(0, begin); + int validEnd = + qMin(end, static_cast(values.size())); double center = 0.0; int count = 0; - for(int i = 0; i < values.size(); ++i) { + for(int i = validBegin; i < validEnd; ++i) { if(isFiniteNumber(values[i])) { center += values[i]; ++count; } } - if(count == 0) { return std::numeric_limits::quiet_NaN(); } - center /= count; - // 固定最多 8 次迭代，控制每个候选的计算开销并保持结果稳定。 for(int iteration = 0; iteration < 8; ++iteration) { double weightedSum = 0.0; double weightTotal = 0.0; - - for(int i = 0; i < values.size(); ++i) { + for(int i = validBegin; i < validEnd; ++i) { if(!isFiniteNumber(values[i])) { continue; } double distance = qAbs(values[i] - center); - // 距离接近零时直接取权重 1，避免除零并保持中心点不被放大。 - double weight = distance <= huberDelta || distance < 1.0e-12 - ? 1.0 - : huberDelta / distance; + double weight = + distance <= huberDelta || + distance < 1.0e-12 + ? 1.0 + : huberDelta / distance; weightedSum += weight * values[i]; weightTotal += weight; } @@ -5300,133 +6706,549 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError( if(weightTotal <= 1.0e-12) { break; } - double nextCenter = weightedSum / weightTotal; if(qAbs(nextCenter - center) <= 1.0e-12) { center = nextCenter; break; } - center = nextCenter; } - return center; }; - // 整体压力/导数误差用于排序；verticalLoss 主要用于解释整体上下偏移。 - breakdown.pressureLoss = huberRms(pressureResidual, 0, numPoints); - breakdown.derivativeLoss = huberRms(derivativeResidual, 0, numPoints); - breakdown.verticalBiasPressure = huberCenter(pressureResidual); - breakdown.verticalBiasDerivative = huberCenter(derivativeResidual); - breakdown.verticalLoss = - 0.5 * (qAbs(breakdown.verticalBiasPressure) + - qAbs(breakdown.verticalBiasDerivative)); - - // 去除上下中心后，把残差投影到目标曲线斜率上估计左右偏移。 - // residual ~= -physicalShift * targetSlope，因此 physicalShift 取回归系数的相反数。 - // 这是局部一阶近似，用于判断方向和大小，不等同于再次优化时间轴。 - double horizontalNumerator = 0.0; - double horizontalDenominator = 0.0; + // 压力和导数合并后只求一个公共中心，表示两条曲线共同的上下位移。 + // 分别去中心会把压力与导数之间真实的相对形状差异一并消除。 + auto huberCommonCenterRange = [&huberCenterRange]( + const QVector& pressureValues, + const QVector& derivativeValues, + int begin, + int end) -> double { + QVector combined; + int validBegin = qMax(0, begin); + int validEnd = qMin( + end, + qMin(static_cast(pressureValues.size()), + static_cast(derivativeValues.size()))); + combined.reserve(2 * qMax(0, validEnd - validBegin)); + + for(int i = validBegin; i < validEnd; ++i) { + if(isFiniteNumber(pressureValues[i])) { + combined.append(pressureValues[i]); + } + if(isFiniteNumber(derivativeValues[i])) { + combined.append(derivativeValues[i]); + } + } + return huberCenterRange(combined, 0, combined.size()); + }; + + // 两个通道按能量等权合并，返回值与单通道 Huber RMS 保持同一量纲。 + auto jointHuberRmsAround = [&huberRmsAround]( + const QVector& pressureValues, + const QVector& derivativeValues, + int begin, + int end, + double center) -> double { + double pressureLoss = huberRmsAround( + pressureValues, begin, end, center); + double derivativeLoss = huberRmsAround( + derivativeValues, begin, end, center); + if(!isFiniteNumber(pressureLoss) || + !isFiniteNumber(derivativeLoss)) { + return std::numeric_limits::quiet_NaN(); + } + return qSqrt(0.5 * + (pressureLoss * pressureLoss + + derivativeLoss * derivativeLoss)); + }; + + // 主目标始终使用未做上下或左右校正的完整曲线误差。 + breakdown.pressureLoss = + huberRms(pressureResidual, 0, numPoints); + breakdown.derivativeLoss = + huberRms(derivativeResidual, 0, numPoints); + // 压力和导数各占一半权重。缩放后 residualVector 的二范数就是 + // sqrt(0.5 * pressureLoss^2 + 0.5 * derivativeLoss^2)。 + const double residualScale = qSqrt(0.5 / numPoints); + breakdown.residualVector.reserve(2 * numPoints); for(int i = 0; i < numPoints; ++i) { - if(!isFiniteNumber(pressureResidual[i]) || - !isFiniteNumber(derivativeResidual[i])) { + breakdown.residualVector.append( + residualScale * + huberEquivalentResidual(pressureResidual[i])); + } + for(int i = 0; i < numPoints; ++i) { + breakdown.residualVector.append( + residualScale * + huberEquivalentResidual(derivativeResidual[i])); + } + + const double logGridStep = + (targetLogMaxX - targetLogMinX) / + (numPoints - 1); + const double resultLogMinX = qLn(resultMinX); + const double resultLogMaxX = qLn(resultMaxX); + // 左右配准只在所有候选位移都共同覆盖的固定区间比较，至少保留 80% + // 目标点；每个 log-time 网格间隔再细分为 8 份，提高位移分辨率。 + const int minimumRegistrationPoints = + (numPoints * 4) / 5; + const int shiftSubdivisions = 8; + int maximumShiftIntervals = 4; + int registrationBegin = 0; + int registrationEnd = numPoints; + double maximumPhysicalShift = + maximumShiftIntervals * logGridStep; + + // 所有位移候选使用同一组目标点。结果范围不足时逐步缩小最大位移， + // 但用于配准的固定公共区间不得少于目标网格的 80%。 + auto updateRegistrationRange = [&](double maximumShift) { + registrationBegin = 0; + registrationEnd = numPoints; + while(registrationBegin < registrationEnd && + commonLogX[registrationBegin] - maximumShift < + resultLogMinX - 1.0e-12) { + ++registrationBegin; + } + while(registrationEnd > registrationBegin && + commonLogX[registrationEnd - 1] + maximumShift > + resultLogMaxX + 1.0e-12) { + --registrationEnd; + } + }; + + updateRegistrationRange(maximumPhysicalShift); + while(maximumShiftIntervals > 0 && + registrationEnd - registrationBegin < + minimumRegistrationPoints) { + --maximumShiftIntervals; + maximumPhysicalShift = + maximumShiftIntervals * logGridStep; + updateRegistrationRange(maximumPhysicalShift); + } + if(registrationEnd - registrationBegin < + minimumRegistrationPoints) { + return invalidLoss; + } + + // physicalShift 为正表示模拟曲线偏右；对齐时在目标时刻右侧读取模拟值。 + auto buildShiftResidual = [&]( + double physicalShift, + int compareBegin, + int compareEnd, + QVector* shiftedPressure, + QVector* shiftedDerivative) -> bool { + if(!shiftedPressure || !shiftedDerivative) { + return false; + } + + shiftedPressure->fill( + std::numeric_limits::quiet_NaN(), + numPoints); + shiftedDerivative->fill( + std::numeric_limits::quiet_NaN(), + numPoints); + + int validBegin = qMax(0, compareBegin); + int validEnd = qMin(numPoints, compareEnd); + for(int i = validBegin; i < validEnd; ++i) { + double shiftedLogX = + commonLogX[i] + physicalShift; + if(shiftedLogX < resultLogMinX - 1.0e-12 || + shiftedLogX > + resultLogMaxX + 1.0e-12) { + return false; + } + + double shiftedX = qExp( + qBound(resultLogMinX, + shiftedLogX, + resultLogMaxX)); + double resultLogPressure = 0.0; + double resultLogDerivative = 0.0; + if(!interpolateLogValue( + resultPressure, shiftedX, + &resultLogPressure) || + !interpolateLogValue( + resultDerivative, shiftedX, + &resultLogDerivative)) { + return false; + } + + (*shiftedPressure)[i] = + resultLogPressure - targetLogPressure[i]; + (*shiftedDerivative)[i] = + resultLogDerivative - targetLogDerivative[i]; + } + return true; + }; + + // 损失相同时优先选择绝对位移更小的候选，避免平坦 profile 在数值噪声 + // 下无故偏向搜索边界。 + auto isBetterProfileValue = []( + double loss, + double shift, + double bestLoss, + double bestShift) -> bool { + const double tolerance = 1.0e-12; + return loss < bestLoss - tolerance || + (qAbs(loss - bestLoss) <= tolerance && + qAbs(shift) < qAbs(bestShift)); + }; + + const int halfShiftStepCount = + maximumShiftIntervals * shiftSubdivisions; + const double physicalShiftStep = + logGridStep / shiftSubdivisions; + double zeroShiftCenteredLoss = + std::numeric_limits::quiet_NaN(); + double bestCenteredLoss = + std::numeric_limits::infinity(); + double bestPhysicalShift = 0.0; + int bestShiftStep = 0; + double bestPressureLoss = + std::numeric_limits::infinity(); + double bestPressureShift = 0.0; + double bestDerivativeLoss = + std::numeric_limits::infinity(); + double bestDerivativeShift = 0.0; + QVector profileLosses( + 2 * halfShiftStepCount + 1, + std::numeric_limits::quiet_NaN()); + QVector profileCommonBiases( + 2 * halfShiftStepCount + 1, + std::numeric_limits::quiet_NaN()); + QVector shiftedPressureResidual; + QVector shiftedDerivativeResidual; + + // 位移和公共上下偏移联合求解，避免“先扣上下还是先扣左右”的顺序依赖。 + for(int shiftStep = -halfShiftStepCount; + shiftStep <= halfShiftStepCount; + ++shiftStep) { + double physicalShift = + shiftStep * physicalShiftStep; + if(!buildShiftResidual( + physicalShift, + registrationBegin, + registrationEnd, + &shiftedPressureResidual, + &shiftedDerivativeResidual)) { continue; } - double pressureCentered = - pressureResidual[i] - breakdown.verticalBiasPressure; - double derivativeCentered = - derivativeResidual[i] - breakdown.verticalBiasDerivative; - horizontalNumerator += pressureSlope[i] * pressureCentered + - derivativeSlope[i] * derivativeCentered; - horizontalDenominator += pressureSlope[i] * pressureSlope[i] + - derivativeSlope[i] * derivativeSlope[i]; + double commonBias = huberCommonCenterRange( + shiftedPressureResidual, + shiftedDerivativeResidual, + registrationBegin, + registrationEnd); + double centeredLoss = jointHuberRmsAround( + shiftedPressureResidual, + shiftedDerivativeResidual, + registrationBegin, + registrationEnd, + commonBias); + double pressureBias = huberCenterRange( + shiftedPressureResidual, + registrationBegin, + registrationEnd); + double derivativeBias = huberCenterRange( + shiftedDerivativeResidual, + registrationBegin, + registrationEnd); + double pressureLoss = huberRmsAround( + shiftedPressureResidual, + registrationBegin, + registrationEnd, + pressureBias); + double derivativeLoss = huberRmsAround( + shiftedDerivativeResidual, + registrationBegin, + registrationEnd, + derivativeBias); + + if(!isFiniteNumber(centeredLoss) || + !isFiniteNumber(pressureLoss) || + !isFiniteNumber(derivativeLoss)) { + continue; + } + int profileIndex = shiftStep + halfShiftStepCount; + profileLosses[profileIndex] = centeredLoss; + profileCommonBiases[profileIndex] = commonBias; + if(shiftStep == 0) { + zeroShiftCenteredLoss = centeredLoss; + } + if(isBetterProfileValue( + centeredLoss, + physicalShift, + bestCenteredLoss, + bestPhysicalShift)) { + bestCenteredLoss = centeredLoss; + bestPhysicalShift = physicalShift; + bestShiftStep = shiftStep; + } + if(isBetterProfileValue( + pressureLoss, + physicalShift, + bestPressureLoss, + bestPressureShift)) { + bestPressureLoss = pressureLoss; + bestPressureShift = physicalShift; + } + if(isBetterProfileValue( + derivativeLoss, + physicalShift, + bestDerivativeLoss, + bestDerivativeShift)) { + bestDerivativeLoss = derivativeLoss; + bestDerivativeShift = physicalShift; + } } - breakdown.horizontalShift = horizontalDenominator > 1.0e-12 - ? horizontalNumerator / - horizontalDenominator - : 0.0; - breakdown.horizontalPhysicalShift = -breakdown.horizontalShift; - breakdown.horizontalLoss = qAbs(breakdown.horizontalShift); + if(!isFiniteNumber(zeroShiftCenteredLoss) || + !isFiniteNumber(bestCenteredLoss) || + !isFiniteNumber(bestPressureLoss) || + !isFiniteNumber(bestDerivativeLoss)) { + return invalidLoss; + } - // 从原始残差中扣除“整体上下 + 等效左右”两部分，剩余项才作为形状误差。 - // 因此 shapeLoss 较大而 vertical/horizontal 较小时，说明主要是曲率、拐点 - // 或导数变化趋势不一致，而不是简单的整体平移。 - QVector shapePressure(numPoints, - std::numeric_limits::quiet_NaN()); - QVector shapeDerivative(numPoints, - std::numeric_limits::quiet_NaN()); + // horizontalGain 是“允许水平位移”相对“固定零位移”减少的稳健能量。 + // 只有改善足够明显且最优点不是边界，才把位移解释为可靠左右偏差。 + double horizontalGain = nestedRmsContribution( + zeroShiftCenteredLoss, bestCenteredLoss); + double horizontalSignalThreshold = + qMax(1.0e-5, zeroShiftCenteredLoss * 0.02); + int bestProfileIndex = bestShiftStep + halfShiftStepCount; + double nearbyProfileLoss = + std::numeric_limits::infinity(); + int leftProfileIndex = + bestProfileIndex - shiftSubdivisions; + int rightProfileIndex = + bestProfileIndex + shiftSubdivisions; + if(leftProfileIndex >= 0 && + leftProfileIndex < profileLosses.size() && + isFiniteNumber(profileLosses[leftProfileIndex])) { + nearbyProfileLoss = qMin( + nearbyProfileLoss, + profileLosses[leftProfileIndex]); + } + if(rightProfileIndex >= 0 && + rightProfileIndex < profileLosses.size() && + isFiniteNumber(profileLosses[rightProfileIndex])) { + nearbyProfileLoss = qMin( + nearbyProfileLoss, + profileLosses[rightProfileIndex]); + } + double profileContrast = isFiniteNumber(nearbyProfileLoss) + ? nestedRmsContribution( + nearbyProfileLoss, + bestCenteredLoss) + : 0.0; + bool flatRegistrationProfile = + profileContrast <= horizontalSignalThreshold; + + // 平台曲线的 profile 也可能很平，但公共 bias 在各个位移下保持稳定， + // 此时仍能可靠判断上下。只有近优位移会明显改变 bias 才说明上下/左右不可辨识。 + double minimumNearOptimalBias = + std::numeric_limits::infinity(); + double maximumNearOptimalBias = + -std::numeric_limits::infinity(); + for(int i = 0; i < profileLosses.size(); ++i) { + if(isFiniteNumber(profileLosses[i]) && + isFiniteNumber(profileCommonBiases[i]) && + profileLosses[i] <= + bestCenteredLoss + horizontalSignalThreshold) { + minimumNearOptimalBias = qMin( + minimumNearOptimalBias, + profileCommonBiases[i]); + maximumNearOptimalBias = qMax( + maximumNearOptimalBias, + profileCommonBiases[i]); + } + } + double nearOptimalBiasSpread = + isFiniteNumber(minimumNearOptimalBias) && + isFiniteNumber(maximumNearOptimalBias) + ? maximumNearOptimalBias - minimumNearOptimalBias + : std::numeric_limits::infinity(); + bool commonBiasStable = nearOptimalBiasSpread <= + qMax(1.0e-4, huberDelta * 0.05); + bool horizontalAtBoundary = + halfShiftStepCount > 0 && + qAbs(bestShiftStep) == halfShiftStepCount; + // 压力和导数通道分别求出的最佳位移若方向相反或相差过大，说明一个 + // 单一水平平移无法解释两条曲线，此时标记配准歧义并禁用左右引导。 + bool pressureShiftDetected = + qAbs(bestPressureShift) >= + 0.5 * physicalShiftStep; + bool derivativeShiftDetected = + qAbs(bestDerivativeShift) >= + 0.5 * physicalShiftStep; + bool channelShiftConflict = + pressureShiftDetected && + derivativeShiftDetected && + (bestPressureShift * bestDerivativeShift < 0.0 || + qAbs(bestPressureShift - bestDerivativeShift) > + 2.0 * logGridStep); + + breakdown.horizontalLoss = horizontalGain; + breakdown.horizontalReliable = + maximumShiftIntervals > 0 && + !horizontalAtBoundary && + !channelShiftConflict && + !flatRegistrationProfile && + horizontalGain > horizontalSignalThreshold && + qAbs(bestPhysicalShift) >= + 0.5 * physicalShiftStep; + breakdown.registrationAmbiguous = + channelShiftConflict || + (qAbs(bestPhysicalShift) >= + 0.5 * physicalShiftStep && + !breakdown.horizontalReliable) || + (flatRegistrationProfile && !commonBiasStable); + breakdown.horizontalPhysicalShift = + breakdown.horizontalReliable + ? bestPhysicalShift + : 0.0; + + // 可信水平位移确定后，在该位移实际覆盖的最大区间重新计算上下和形状。 + int diagnosticBegin = 0; + int diagnosticEnd = numPoints; + while(diagnosticBegin < diagnosticEnd && + commonLogX[diagnosticBegin] + + breakdown.horizontalPhysicalShift < + resultLogMinX - 1.0e-12) { + ++diagnosticBegin; + } + while(diagnosticEnd > diagnosticBegin && + commonLogX[diagnosticEnd - 1] + + breakdown.horizontalPhysicalShift > + resultLogMaxX + 1.0e-12) { + --diagnosticEnd; + } + if(diagnosticEnd - diagnosticBegin < + minimumRegistrationPoints || + !buildShiftResidual( + breakdown.horizontalPhysicalShift, + diagnosticBegin, + diagnosticEnd, + &shiftedPressureResidual, + &shiftedDerivativeResidual)) { + return invalidLoss; + } + + double commonBias = huberCommonCenterRange( + shiftedPressureResidual, + shiftedDerivativeResidual, + diagnosticBegin, + diagnosticEnd); + double rawAlignedLoss = jointHuberRmsAround( + shiftedPressureResidual, + shiftedDerivativeResidual, + diagnosticBegin, + diagnosticEnd, + 0.0); + double centeredAlignedLoss = jointHuberRmsAround( + shiftedPressureResidual, + shiftedDerivativeResidual, + diagnosticBegin, + diagnosticEnd, + commonBias); + if(!isFiniteNumber(commonBias) || + !isFiniteNumber(rawAlignedLoss) || + !isFiniteNumber(centeredAlignedLoss)) { + return invalidLoss; + } + + // 原始对齐误差减去公共中心后的能量差定义为上下误差贡献。只有它相对 + // 当前对齐误差足够明显，且配准无歧义时，公共 bias 才可用于有符号选参。 + breakdown.verticalCommonBias = commonBias; + breakdown.verticalLoss = nestedRmsContribution( + rawAlignedLoss, centeredAlignedLoss); + breakdown.verticalReliable = + !breakdown.registrationAmbiguous && + breakdown.verticalLoss > + qMax(1.0e-5, rawAlignedLoss * 0.02); + + QVector shapePressure( + numPoints, std::numeric_limits::quiet_NaN()); + QVector shapeDerivative( + numPoints, std::numeric_limits::quiet_NaN()); + for(int i = diagnosticBegin; i < diagnosticEnd; ++i) { + if(isFiniteNumber(shiftedPressureResidual[i])) { + shapePressure[i] = + shiftedPressureResidual[i] - commonBias; + } + if(isFiniteNumber(shiftedDerivativeResidual[i])) { + shapeDerivative[i] = + shiftedDerivativeResidual[i] - commonBias; + } + } + + // shapeLoss 是去除可信左右位移和稳健公共中心后的剩余误差。 + double shapePressureLoss = huberRms( + shapePressure, diagnosticBegin, diagnosticEnd); + double shapeDerivativeLoss = huberRms( + shapeDerivative, diagnosticBegin, diagnosticEnd); + if(!isFiniteNumber(shapePressureLoss) || + !isFiniteNumber(shapeDerivativeLoss)) { + return invalidLoss; + } + breakdown.shapeLoss = qSqrt( + 0.5 * + (shapePressureLoss * shapePressureLoss + + shapeDerivativeLoss * shapeDerivativeLoss)); + + // 现阶段不识别或单独调度晚期流动段；保留字段只为了维持现有 trace 列。 + breakdown.lateDerivativeSlopeBias = 0.0; + breakdown.lateDerivativeTrendLoss = 0.0; + breakdown.lateDerivativeTrendReliable = false; - for(int i = 0; i < numPoints; ++i) { - if(isFiniteNumber(pressureResidual[i])) { - shapePressure[i] = pressureResidual[i] - - breakdown.verticalBiasPressure - - breakdown.horizontalShift * pressureSlope[i]; - } - - if(isFiniteNumber(derivativeResidual[i])) { - shapeDerivative[i] = derivativeResidual[i] - - breakdown.verticalBiasDerivative - - breakdown.horizontalShift * - derivativeSlope[i]; - } - } - - breakdown.shapeLoss = - 0.5 * (huberRms(shapePressure, 0, numPoints) + - huberRms(shapeDerivative, 0, numPoints)); - - // 将对数时间网格分成早、中、晚三段，用于定位误差集中出现的阶段。 - // 网格本身按 log(time) 均匀分布，所以三段对应的是时间数量级，而非原始 - // 线性时间长度，适合双对数试井曲线的早期/中期/晚期判读。 - const int segment1 = numPoints / 3; - const int segment2 = (2 * numPoints) / 3; - breakdown.pressureEarlyLoss = huberRms(pressureResidual, 0, segment1); - breakdown.pressureMiddleLoss = - huberRms(pressureResidual, segment1, segment2); - breakdown.pressureLateLoss = - huberRms(pressureResidual, segment2, numPoints); - breakdown.derivativeEarlyLoss = - huberRms(derivativeResidual, 0, segment1); - breakdown.derivativeMiddleLoss = - huberRms(derivativeResidual, segment1, segment2); - breakdown.derivativeLateLoss = - huberRms(derivativeResidual, segment2, numPoints); - - // 函数入口已将 m_lastObjectiveBreakdown 重置为无效状态，因此这里直接返回 - // invalidLoss 时不会把上一候选的误差分解误报给调用方。 - // 导数、压力或形状无法形成有效统计时，整个候选都视为无效，避免 NaN - // 进入粒子排序。 if(!isFiniteNumber(breakdown.pressureLoss) || !isFiniteNumber(breakdown.derivativeLoss) || + !isFiniteNumber(breakdown.verticalCommonBias) || + !isFiniteNumber(breakdown.verticalLoss) || + !isFiniteNumber(breakdown.horizontalPhysicalShift) || + !isFiniteNumber(breakdown.horizontalLoss) || !isFiniteNumber(breakdown.shapeLoss) || - !isFiniteNumber(breakdown.verticalLoss)) { + breakdown.residualVector.size() != 2 * numPoints) { return invalidLoss; } - // 当前总损失作为 fitness 用于粒子比较、收敛/停止判断；上下、左右、形状和 - // 分段分量先作为诊断输出，不在本次改动中直接参与参数更新。 - breakdown.total = 0.5 * breakdown.pressureLoss + - 0.5 * breakdown.derivativeLoss + - 0.1 * breakdown.coveragePenalty; - breakdown.valid = isFiniteNumber(breakdown.total) && - breakdown.total >= 0.0; + if(isSurrogateScreeningEnabled()) { + // 代理模型仍按原目标比较，不能用未参与训练的新损失改变候选排序。 + breakdown.total = + 0.5 * breakdown.pressureLoss + + 0.5 * breakdown.derivativeLoss; + } else { + // 非代理总目标等于固定 100 维 Huber 等效残差的二范数；压力和 + // 导数各占一半能量。上下、左右和形状分量不参与候选排序与接受。 + breakdown.total = qSqrt( + 0.5 * breakdown.pressureLoss * breakdown.pressureLoss + + 0.5 * breakdown.derivativeLoss * breakdown.derivativeLoss); + } + breakdown.valid = + isFiniteNumber(breakdown.total) && + breakdown.total >= 0.0; m_lastObjectiveBreakdown = breakdown; - DEBUG_OUT(QString("LogLog objective: pressure=%1, derivative=%2, vertical=%3, horizontal=%4, shape=%5, coverage=%6, total=%7") - .arg(breakdown.pressureLoss, 0, 'e', 4) - .arg(breakdown.derivativeLoss, 0, 'e', 4) - .arg(breakdown.verticalLoss, 0, 'e', 4) - .arg(breakdown.horizontalLoss, 0, 'e', 4) - .arg(breakdown.shapeLoss, 0, 'e', 4) - .arg(breakdown.coverage, 0, 'f', 4) - .arg(breakdown.total, 0, 'e', 4)); - - return breakdown.valid ? qMin(1.0e9, breakdown.total) : invalidLoss; + DEBUG_OUT( + QString("LogLog objective: pressure=%1, derivative=%2, vertical=%3, horizontal=%4, shape=%5, ambiguous=%6, shift=%7, coverage=%8, total=%9") + .arg(breakdown.pressureLoss, 0, 'e', 4) + .arg(breakdown.derivativeLoss, 0, 'e', 4) + .arg(breakdown.verticalLoss, 0, 'e', 4) + .arg(breakdown.horizontalLoss, 0, 'e', 4) + .arg(breakdown.shapeLoss, 0, 'e', 4) + .arg(breakdown.registrationAmbiguous) + .arg(breakdown.horizontalPhysicalShift, 0, 'e', 4) + .arg(breakdown.coverage, 0, 'f', 4) + .arg(breakdown.total, 0, 'e', 4)); + + return breakdown.valid + ? qMin(1.0e9, breakdown.total) + : invalidLoss; } catch(const std::exception& e) { - DEBUG_OUT(QString("Exception in LogLog error calculation: %1").arg(e.what())); + DEBUG_OUT( + QString("Exception in LogLog error calculation: %1") + .arg(e.what())); return invalidLoss; } catch(...) { DEBUG_OUT("Unknown exception in LogLog error calculation"); diff --git a/Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp b/Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp index f0f243e..2a527ee 100644 --- a/Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp +++ b/Src/nmNum/nmSubWxs/nmWxAutomaticFitting.cpp @@ -35,6 +35,14 @@ bool nmAutoFitUiIsFinite(double value) #endif } +// 物理边界可能由多个浮点配置量计算得到，界面十进制文本再转回 double 后会有 +// 末位舍入差。这里只放宽约 8 个机器精度，不放宽实际物理范围。 +bool nmAutoFitUiNearlyEqual(double left, double right) +{ + const double scale = qMax(std::fabs(left), std::fabs(right)); + return std::fabs(left - right) <= DBL_EPSILON * 8.0 * scale; +} + // 从系统参数表读取物理边界，读取失败时保留调用方提供的兜底边界。 bool nmAutoFitReadPhysicalRange(const char* parameterName, double fallbackMin, double fallbackMax, double& minValue, double& maxValue) @@ -337,9 +345,10 @@ void nmWxAutomaticFitting::setParameterRange(int parameterIndex, m_updatingParameterRanges = wasUpdatingRanges; } -// 根据初值生成首次或拟合后的建议范围，并始终限制在物理边界内。 +// 根据当前初值生成建议搜索范围，并始终截断在系统物理边界内。skin 使用 +// 加减固定宽度，其余正值参数使用倍率范围；该规则在首次加载和拟合完成后复用。 void nmWxAutomaticFitting::updateRangeForParameter(int parameterIndex, - double centerValue, bool afterFit) + double centerValue) { if(!m_parameterTable || parameterIndex < 0 || parameterIndex >= 8 || !nmAutoFitUiIsFinite(centerValue)) { @@ -371,12 +380,12 @@ void nmWxAutomaticFitting::updateRangeForParameter(int parameterIndex, double newMin = physicalMin; double newMax = physicalMax; if(parameterIndex == 1) { - const double skinHalfRange = afterFit ? 1.0 : 10.0; + const double skinHalfRange = 10.0; newMin = qMax(physicalMin, reference - skinHalfRange); newMax = qMin(physicalMax, reference + skinHalfRange); } else if(reference > 0.0 && !(parameterIndex == 7 && centerValue <= 0.0)) { - const double lowerFactor = afterFit ? 0.5 : 0.1; - const double upperFactor = afterFit ? 2.0 : 10.0; + const double lowerFactor = 0.1; + const double upperFactor = 10.0; newMin = qMax(physicalMin, reference * lowerFactor); newMax = qMin(physicalMax, reference * upperFactor); } else if(parameterIndex == 7) { @@ -409,7 +418,7 @@ void nmWxAutomaticFitting::initializeSuggestedParameterRanges() bool initialOk = false; const double initialValue = initialItem->text().toDouble(&initialOk); if(initialOk && nmAutoFitUiIsFinite(initialValue)) { - updateRangeForParameter(parameterIndex, initialValue, false); + updateRangeForParameter(parameterIndex, initialValue); } else { // 数据对象没有提供该初值时使用完整物理区间，不回退到旧的默认范围。 double physicalMin = 0.0; @@ -462,7 +471,7 @@ void nmWxAutomaticFitting::normalizeSavedParameterRanges() bool initialOk = false; const double initialValue = initialItem->text().toDouble(&initialOk); if(initialOk && nmAutoFitUiIsFinite(initialValue)) { - updateRangeForParameter(parameterIndex, initialValue, false); + updateRangeForParameter(parameterIndex, initialValue); } else if(physicalMax >= physicalMin) { setParameterRange(parameterIndex, physicalMin, physicalMax); } @@ -529,9 +538,12 @@ bool nmWxAutomaticFitting::validateParameterTable(QString& errorMessage, int par return false; } - if(minValue < physicalMin || minValue > physicalMax - || initialValue < physicalMin || initialValue > physicalMax - || maxValue < physicalMin || maxValue > physicalMax) { + if((minValue < physicalMin && !nmAutoFitUiNearlyEqual(minValue, physicalMin)) + || (minValue > physicalMax && !nmAutoFitUiNearlyEqual(minValue, physicalMax)) + || (initialValue < physicalMin && !nmAutoFitUiNearlyEqual(initialValue, physicalMin)) + || (initialValue > physicalMax && !nmAutoFitUiNearlyEqual(initialValue, physicalMax)) + || (maxValue < physicalMin && !nmAutoFitUiNearlyEqual(maxValue, physicalMin)) + || (maxValue > physicalMax && !nmAutoFitUiNearlyEqual(maxValue, physicalMax))) { errorMessage = tr("The values of %1 exceed the physical range [%2, %3].") .arg(tr(parameterNames[currentParameterIndex])) .arg(QString::number(physicalMin, 'g', 10)) @@ -1195,8 +1207,8 @@ void nmWxAutomaticFitting::onWellSelected(int index) m_parameterTable->item(2, 3)->setText(QString::number(wellboreStorageValue)); if(m_autoParameterRanges) { - updateRangeForParameter(1, skinValue, false); - updateRangeForParameter(2, wellboreStorageValue, false); + updateRangeForParameter(1, skinValue); + updateRangeForParameter(2, wellboreStorageValue); } } } @@ -1531,9 +1543,9 @@ void nmWxAutomaticFitting::updateBestParametersToTable() // 更新初始值 m_parameterTable->item(paramIndex, 3)->setText(QString::number(bestValue, 'g', 4)); - // 只有自动范围模式才根据拟合结果收窄下一轮搜索区间；手工范围由用户保留。 + // 自动范围模式下，以拟合结果为中心复用首次建范围的规则；手工范围由用户保留。 if(m_autoParameterRanges) { - updateRangeForParameter(paramIndex, bestValue, true); + updateRangeForParameter(paramIndex, bestValue); } }