优化自动拟合，引入灵敏度分析与信赖域求解

Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h

@@ -20,64 +20,66 @@ class nmDataWellBase;
 class QTimer;
 class QProcess;
 
-// 双对数曲线误差分解，供误差诊断和后续参数调整读取。
+// 双对数曲线误差分解。该结构同时保存用于候选排序的主目标，以及用于判断
-//
+// 曲线上下、左右和形状偏差的诊断量。total 是唯一的接受和排序依据，诊断量
-// 所有 pressure/derivative/shape 数值均是在 log(value) 空间计算的无量纲误差。
+// 只参与信赖域选参，不能再次叠加到 total，否则会重复计算同一批曲线残差。
-// verticalBias* 保留正负号：正值表示模拟曲线整体高于目标，负值表示整体低于目标。
-// horizontalPhysicalShift 是 log(time) 方向的等效平移量，正值表示模拟曲线相对目标偏右。
-// total 仍是 PSO 当前使用的 fitness（用于粒子比较和收敛判断）；其余字段只描述误差
-// 来源，不会在本次改动中直接修改粒子参数，避免损失诊断和参数更新策略互相耦合。
 struct AutoFitObjectiveBreakdown {
-	// valid/total 是本次评价是否有效及其最终 fitness（用于排序和收敛判断）。
+	// valid 表示本次曲线评价完整有效；无效评价统一保留 total=1e10。
+	// pressureLoss 和 derivativeLoss 均在 log(value) 空间按固定网格计算。
 	bool valid;
 	double total;
-	// pressureLoss/derivativeLoss 是压力和导数两条曲线的整体误差。
 	double pressureLoss;
 	double derivativeLoss;
-	// vertical* 描述整体上下偏移；保留 bias 的符号以判断偏高或偏低。
+	// 固定目标网格上的 Huber 等效残差。非代理搜索使用它建立完整 Jacobian，
-	double verticalBiasPressure;
+	// 向量平方和与 total 的平方一致。
-	double verticalBiasDerivative;
+	QVector<double> residualVector;
+
+	// 上下偏差使用压力和导数残差共享的 Huber 稳健中心。
+	// verticalCommonBias 为正表示模拟曲线整体偏高，为负表示整体偏低；
+	// verticalReliable=false 时仍保留数值，但不能据此确定参数调整方向。
+	double verticalCommonBias;
 	double verticalLoss;
-	// horizontal* 描述等效的对数时间偏移；physicalShift 为正表示模拟曲线相对目标向右
+	bool verticalReliable;
-	// （时间延迟），为负表示向左。
+
-	double horizontalShift;
+	// 水平偏差在 log(time) 坐标中计算。physicalShift 为正表示模拟曲线相对
+	// 目标偏右，即相同曲线特征在模拟结果中出现得更晚。
 	double horizontalPhysicalShift;
 	double horizontalLoss;
-	// shapeLoss 是去除整体上下和左右偏移后剩余的曲线形状差异。
+	bool horizontalReliable;
+	// true 表示当前曲线无法可靠区分上下和左右误差；此时禁止使用两类有符号
+	// 诊断量选参，但去除公共中心后的 shapeLoss 仍可用于局部选参。
+	bool registrationAmbiguous;
+
+	// 去除稳健公共中心和可信左右偏差后剩余的整体形状误差；verticalReliable
+	// 只控制能否把公共中心解释为上下参数方向，不改变 shape 的中心化公式。
 	double shapeLoss;
-	// 早、中、晚分段误差用于定位误差主要出现在哪个时间阶段。
+
-	double pressureEarlyLoss;
+	// 兼容现有 trace 列。当前非代理搜索不再单独识别或调度晚期分量。
-	double pressureMiddleLoss;
+	double lateDerivativeSlopeBias;
-	double pressureLateLoss;
+	double lateDerivativeTrendLoss;
-	double derivativeEarlyLoss;
+	bool lateDerivativeTrendReliable;
-	double derivativeMiddleLoss;
+
-	double derivativeLateLoss;
+	// 模拟曲线对目标固定网格的有效覆盖率，取覆盖点比例与连续 log-time
-	// coverage 是 50 点目标网格上的 min(有效点比例、连续 log-time 跨度比例)。
+	// 跨度比例中的较小值。低于损失函数门槛时本次评价直接无效。
-	// coveragePenalty 是归一化覆盖缺口的平方惩罚，并以 0.1 权重加入 total。
 	double coverage;
-	double coveragePenalty;
 
-	// 无效评价使用 1e10 作为统一的“差解”标记；其他字段用 NaN 表示尚未得到诊断值。
 	AutoFitObjectiveBreakdown()
 		: valid(false)
 		, total(1.0e10)
 		, pressureLoss(std::numeric_limits<double>::quiet_NaN())
 		, derivativeLoss(std::numeric_limits<double>::quiet_NaN())
-		, verticalBiasPressure(std::numeric_limits<double>::quiet_NaN())
+		, verticalCommonBias(std::numeric_limits<double>::quiet_NaN())
-		, verticalBiasDerivative(std::numeric_limits<double>::quiet_NaN())
 		, verticalLoss(std::numeric_limits<double>::quiet_NaN())
-		, horizontalShift(std::numeric_limits<double>::quiet_NaN())
+		, verticalReliable(false)
 		, horizontalPhysicalShift(std::numeric_limits<double>::quiet_NaN())
 		, horizontalLoss(std::numeric_limits<double>::quiet_NaN())
+		, horizontalReliable(false)
+		, registrationAmbiguous(false)
 		, shapeLoss(std::numeric_limits<double>::quiet_NaN())
-		, pressureEarlyLoss(std::numeric_limits<double>::quiet_NaN())
+		, lateDerivativeSlopeBias(std::numeric_limits<double>::quiet_NaN())
-		, pressureMiddleLoss(std::numeric_limits<double>::quiet_NaN())
+		, lateDerivativeTrendLoss(std::numeric_limits<double>::quiet_NaN())
-		, pressureLateLoss(std::numeric_limits<double>::quiet_NaN())
+		, lateDerivativeTrendReliable(false)
-		, derivativeEarlyLoss(std::numeric_limits<double>::quiet_NaN())
-		, derivativeMiddleLoss(std::numeric_limits<double>::quiet_NaN())
-		, derivativeLateLoss(std::numeric_limits<double>::quiet_NaN())
 		, coverage(std::numeric_limits<double>::quiet_NaN())
-		, coveragePenalty(std::numeric_limits<double>::quiet_NaN())
 	{}
 };
 
@@ -95,6 +97,8 @@ struct AutoFitParticle {
 	QVector<double> velocity;      // 速度
 	QVector<double> bestPosition;  // 真实求解器确认的个体最优位置
 	QVector<double> guideBestPosition; // 仅用于速度更新的引导位置；不会参与真实 gbest/最终结果
+	AutoFitObjectiveBreakdown currentObjectiveBreakdown; // 当前真实评价对应的误差分解
+	AutoFitObjectiveBreakdown bestObjectiveBreakdown;    // pbest 对应的误差分解
 	double fitness;                // 当前适应度
 	double bestFitness;            // 真实求解器确认的个体最优适应度
 	double guideBestObjective;     // guideBestPosition 对应的真实或代理目标值
@@ -200,20 +204,26 @@ private:
 	void loadOptimizationConfig();
 	void loadParameterBounds();
 
-	// ===== PSO核心算法 =====
+	// ===== 自动拟合核心算法 =====
 	//
-	// 主流程：
+	// 代理开启时保留原 PSO 筛选流程；代理关闭时使用真实求解器驱动的
-	// 1. extractUserInitialValues(): 从当前项目数据中取用户已有初始解；
+	// 诊断灵敏度信赖域搜索，不依赖 pbest/gbest 速度公式。
-	// 2. initializeSwarm(): 根据初始解和上下界生成粒子群；
-	// 3. updateParticle(): 对单个粒子跑真实求解器并计算误差；
-	// 4. updateGlobalBest(): 只用真实求解器误差更新全局最优；
-	// 5. updateVelocityAndPosition(): 按 PSO 公式推进下一代粒子。
 	void extractUserInitialValues();
 	void initializeSwarm();
 	void updateVelocityAndPosition();
 	double evaluateFitness(const QVector<double>& parameters);
 	void updateGlobalBest();
 	void updateParticle(int particleIndex);
+	// 非代理拟合入口：建立有限差分灵敏度，按诊断分量选择参数，再用有界
+	// LM/信赖域产生候选；所有候选最终都由真实求解器总误差决定是否接受。
+	StopReasonPSO runTrustRegionFitting();
+	// 对一个信赖域候选执行完整真实评价，并一次性返回误差、诊断量、曲线和耗时。
+	// 返回 false 表示求解失败、损失无效或用户已请求停止。
+	bool evaluateTrustRegionPoint(const QVector<double>& parameters,
+			double* fitness,
+			AutoFitObjectiveBreakdown* breakdown,
+			QVector<QVector<double> >* curve,
+			int* elapsedMs);
 
 	// ===== 参数应用方法 =====
 	//
@@ -284,7 +294,8 @@ private:
 		double surrogateObjective,
 		const QString& screeningDecision,
 		const QVector<double>& pbestPosition,
-		double pbestObjective);
+		double pbestObjective,
+		const AutoFitObjectiveBreakdown* objectiveBreakdown = nullptr);
 	void writeIterationTraceRows();
 	QVector<double> buildTraceParameterVector(const QVector<double>& 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<double> m_initialValues;              // 当前模型中提取的用户初始值，顺序与 m_enabledParamIndices 一致。
 	QVector<AutoFitParticle> m_swarm;             // 粒子群，每个粒子只保存启用参数维度。
-	QVector<double> m_globalBestPosition;         // 全局最优参数，仍是启用参数向量。
+	QVector<double> m_globalBestPosition;         // 真实求解器确认的当前最优参数。
 	double m_globalBestFitness;                   // 全局最优真实误差，越小越好。
 	double m_previousBestFitness;                 // 上一轮全局最优误差，用于自适应参数更新。
+	AutoFitObjectiveBreakdown m_globalBestObjectiveBreakdown; // 真实 gbest 对应的误差分解。
 	QVector<QVector<double> > m_lastEvaluatedLogLogData; // 最近一次真实求解得到的 result log-log 曲线。
 	QVector<QVector<double> > m_globalBestLogLogData;    // 当前全局最优对应的 result log-log 曲线。
 	mutable AutoFitObjectiveBreakdown m_lastObjectiveBreakdown; // 最近一次损失评价的误差分解。
 	QVector<QVector<double> > 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<double> m_convergenceHistory; // 每代全局最优误差历史，用于收敛判断。
 
@@ -393,7 +406,7 @@ private:
 	// ===== 精英保护 =====  
 	QVector<double> 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_<run_id>.csv 完整路径。
+	QString m_traceFilePath;             // 本次自动拟合 trace CSV 的完整路径。
-	QString m_traceMetaFilePath;         // pso_baseline_trace_<run_id>.meta.json 完整路径。
+	QString m_traceMetaFilePath;         // 与 trace 匹配的 meta JSON 完整路径。
 	QFile m_traceFile;                   // trace CSV 文件句柄。
 	bool m_surrogateScreeningEnabled;    // 用户配置中的 PSO acceleration 开关。
 	unsigned int m_psoRandomSeed;        // PSO 随机种子，也用于可复现 random audit。

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();

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 lowerFactor = 0.1;
-		const double upperFactor = afterFit ? 2.0 : 10.0;
+		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
+		if((minValue < physicalMin && !nmAutoFitUiNearlyEqual(minValue, physicalMin))
-		   || initialValue < physicalMin || initialValue > physicalMax
+		   || (minValue > physicalMax && !nmAutoFitUiNearlyEqual(minValue, physicalMax))
-		   || maxValue < physicalMin || maxValue > 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(1, skinValue);
-			updateRangeForParameter(2, wellboreStorageValue, false);
+			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);
 		}
 	}