diff --git a/Bin/Config/Lang/cn/nmNum_cn.qm b/Bin/Config/Lang/cn/nmNum_cn.qm
index 349e595..ca1b12f 100644
Binary files a/Bin/Config/Lang/cn/nmNum_cn.qm and b/Bin/Config/Lang/cn/nmNum_cn.qm differ
diff --git a/Bin/Config/Lang/cn/nmNum_cn.ts b/Bin/Config/Lang/cn/nmNum_cn.ts
index 42ebb3f..3d81ec7 100644
--- a/Bin/Config/Lang/cn/nmNum_cn.ts
+++ b/Bin/Config/Lang/cn/nmNum_cn.ts
@@ -638,6 +638,389 @@ Reason: %1
代理评分服务器标准错误输出(stderr):%1
+
+ nmCalculationAutoFitLM
+
+ === User Stop Request Received ===
+ === 用户停止请求已接收 ===
+
+
+ Gracefully stopping LM automatic fitting...
+ 正在平稳停止 LM 自动拟合...
+
+
+ Waiting for current solver evaluation to stop...
+ 正在等待当前求解器评估停止...
+
+
+ LM automatic fitting stop request processed
+ LM 自动拟合停止请求已处理
+
+
+ Stop request received but optimization is not running
+ 收到停止请求,但优化未运行
+
+
+ LM fitting trace: %1
+ LM 拟合跟踪:%1
+
+
+ LM fitting trace metadata: %1
+ LM 拟合跟踪元数据:%1
+
+
+ === LM Run Summary ===
+ === LM 运行摘要 ===
+
+
+ Stop reason: %1
+ 停止原因:%1
+
+
+ Result: %1, final error=%2, iterations=%3, evaluations=%4 (successful=%5, failed=%6)
+ 结果:%1,最终误差=%2,迭代次数=%3,评估次数=%4(成功=%5,失败=%6)
+
+
+ SUCCESS
+ 成功
+
+
+ FAILED
+ 失败
+
+
+ Artifacts: trace=%1
+ 产物:跟踪=%1
+
+
+ Artifacts: trace_meta=%1
+ 产物:跟踪元数据=%1
+
+
+ ERROR: Failed to load configuration from data manager
+ 错误:从数据管理器加载配置失败
+
+
+ Algorithm: LM
+ 算法:LM
+
+
+ Enabled parameters count: %1
+ 启用参数数量:%1
+
+
+ ERROR: No parameters enabled for optimization
+ 错误:没有启用优化参数
+
+
+ ERROR: Target LogLog data is empty or insufficient
+ 错误:目标双对数数据为空或不足
+
+
+ ERROR: Target LogLog data arrays have inconsistent sizes
+ 错误:目标双对数数据数组大小不一致
+
+
+ ERROR: Target well name is empty
+ 错误:目标井名称为空
+
+
+ Target data validation passed (%1 data points)
+ 目标数据验证通过(%1 个数据点)
+
+
+ === Evaluating Initial Solution (Elite Protection) ===
+ === 评估初始解 ===
+
+
+ Initial parameters:
+ 初始参数:
+
+
+ Starting initial solution evaluation...
+ 开始初始解评估...
+
+
+ Initial solution evaluation successful
+ 初始解评估成功
+
+
+ Initial Error: %1
+ 初始误差:%1
+
+
+ Initial solution evaluation failed
+ 初始解评估失败
+
+
+ Exception during initial solution evaluation
+ 初始解评估期间出现异常
+
+
+ Critical exception in automatic fitting: %1
+ 自动拟合发生严重异常:%1
+
+
+ CRITICAL ERROR: %1
+ 严重错误:%1
+
+
+ Unknown critical exception in automatic fitting
+ 自动拟合发生未知严重异常
+
+
+ CRITICAL ERROR: Unknown exception in automatic fitting
+ 严重错误:自动拟合发生未知异常
+
+
+ Applying optimized parameters to model...
+ 正在将优化参数应用到模型...
+
+
+ Final full-field calculation skipped after user stop
+ 用户停止后已跳过最终全场计算
+
+
+ Running final full-field calculation with optimized parameters...
+ 正在使用优化参数运行最终全场计算...
+
+
+ Final full-field calculation completed successfully
+ 最终全场计算成功完成
+
+
+ Final full-field calculation stopped by user
+ 最终全场计算已被用户停止
+
+
+ Optimized parameters were found, but the final full-field calculation failed
+ 已找到优化参数,但最终全场计算失败
+
+
+ ERROR: Final full-field calculation failed
+ 错误:最终全场计算失败
+
+
+ === Optimization Results ===
+ === 优化结果 ===
+
+
+ Final error: %1
+ 最终误差:%1
+
+
+ Total iterations: %1
+ 总迭代次数:%1
+
+
+ Total evaluations: %1 (successful: %2)
+ 总评估次数:%1(成功:%2)
+
+
+ Optimized parameters:
+ 优化参数:
+
+
+ Parameters and full-field results applied successfully to data manager
+ 参数和全场计算结果已成功应用到数据管理器
+
+
+ Optimized parameters applied to data manager
+ 优化参数已应用到数据管理器
+
+
+ ERROR: Failed to apply final parameters: %1
+ 错误:应用最终参数失败:%1
+
+
+ ERROR: Unknown error applying final parameters
+ 错误:应用最终参数时出现未知错误
+
+
+ Target achieved. Best error: %1, Iterations: %2
+ 达到目标。最佳误差:%1,迭代次数:%2
+
+
+ === LM AUTOMATIC FITTING SUCCESSFUL ===
+ === LM 自动拟合成功 ===
+
+
+ Automatic fitting converged to a stable solution. Best error: %1, Iterations: %2
+ 自动拟合已收敛到稳定解。最佳误差:%1,迭代次数:%2
+
+
+ === LM AUTOMATIC FITTING CONVERGED ===
+ === LM 自动拟合已收敛 ===
+
+
+ Automatic fitting reached a local optimum. Best error: %1, Iterations: %2
+ 自动拟合达到局部最优。最佳误差:%1,迭代次数:%2
+
+
+ === LM AUTOMATIC FITTING - LOCAL OPTIMUM ===
+ === LM 自动拟合 - 局部最优 ===
+
+
+ Max iterations reached. Best error: %1, Iterations: %2
+ 达到最大迭代次数。最佳误差:%1,迭代次数:%2
+
+
+ === LM AUTOMATIC FITTING - MAX ITERATIONS ===
+ === LM 自动拟合 - 达到最大迭代次数 ===
+
+
+ Best error: %1, Iterations: %2
+ 最佳误差:%1,迭代次数:%2
+
+
+ === LM AUTOMATIC FITTING STOPPED BY USER ===
+ === LM 自动拟合已被用户停止 ===
+
+
+ Automatic fitting failed due to consecutive failures. Best error: %1, Iterations: %2
+ 自动拟合因连续失败而终止。最佳误差:%1,迭代次数:%2
+
+
+ === LM AUTOMATIC FITTING FAILED ===
+ === LM 自动拟合失败 ===
+
+
+ Automatic fitting ended unexpectedly. Best error: %1, Iterations: %2
+ 自动拟合意外结束。最佳误差:%1,迭代次数:%2
+
+
+ === LM AUTOMATIC FITTING - UNKNOWN END ===
+ === LM 自动拟合 - 未知结束 ===
+
+
+ No valid parameters are available for trust-region fitting
+ 没有可用于 LM 拟合的有效参数
+
+
+ The initial solution and parameter-range midpoint are both invalid
+ 初始解和参数范围中点均无效
+
+
+ === Starting LM Main Loop ===
+ === 开始 LM 主循环 ===
+
+
+ LM starting point error: %1; evaluation budget: %2
+ LM 起点误差:%1;最大评估次数:%2
+
+
+ No effective improvement for %1 consecutive steps; rebuilding sensitivity model for confirmation
+ 连续 %1 次无有效改善,正在重建灵敏度模型进行确认
+
+
+ Effective improvement threshold: max(%1, %2% of baseline error); %3 consecutive ineffective steps trigger convergence confirmation
+ 有效改善阈值:取 %1 与基准误差的 %2% 中较大值;连续 %3 次无有效改善后进行收敛确认
+
+
+ Sensitivity probe accepted: error reduced to %1
+ 已接受灵敏度试算点,误差降至 %1
+
+
+ Sensitivity model rebuilt: %1/%2 parameter columns valid
+ 灵敏度模型重建完成:%1/%2 个参数列有效
+
+
+ Sensitivity rebuild produced no effective improvement; local convergence detected
+ 灵敏度重建后仍无有效改善,判定为局部收敛
+
+
+ vertical deviation
+ 上下偏差
+
+
+ horizontal deviation
+ 左右偏差
+
+
+ shape deviation
+ 形状偏差
+
+
+ total error
+ 总误差
+
+
+ Iteration %1: focus=%2, parameters=%3, error=%4, result=%5
+ 迭代 %1:关注=%2,参数个数=%3,误差=%4,结果=%5
+
+
+ accepted
+ 接受
+
+
+ rejected
+ 拒绝
+
+
+ No initial solution for elite protection
+ 没有可用于精英保护的初始解
+
+
+ === Final Result Validation (Elite Protection) ===
+ === 最终结果验证 ===
+
+
+ Comparing results: Initial=%1, Final=%2
+ 比较结果:初始=%1,最终=%2
+
+
+ Improvement: %1 (%2%)
+ 改善:%1(%2%)
+
+
+ Elite protection triggered: final result is invalid or worse than initial
+ 触发精英保护:最终结果无效或差于初始解
+
+
+ Restoring initial solution as final result
+ 正在恢复初始解作为最终结果
+
+
+ Initial solution restored successfully
+ 初始解恢复成功
+
+
+ Final result validated - solution is not worse than initial
+ 最终结果验证通过,拟合结果不差于初始解
+
+
+ Target error achieved
+ 达到目标误差
+
+
+ Algorithm converged to stable solution
+ 算法已收敛到稳定解
+
+
+ Local optimum detected
+ 检测到局部最优
+
+
+ Maximum iterations reached
+ 达到最大迭代次数
+
+
+ Stopped by user request
+ 根据用户请求停止
+
+
+ Too many consecutive failures
+ 连续失败次数过多
+
+
+ Optimization failed
+ 优化失败
+
+
+ Unknown reason
+ 未知原因
+
+
nmCalculationSolver
@@ -3454,6 +3837,12 @@ Supported types: Vertical, Vertical Fractured, and Horizontal Multi-Fractured We
PSO求解完成:
+
+ %1 Optimization completed:
+
+ %1 拟合完成:
+
+
Best Error: %1
@@ -3484,6 +3873,12 @@ Supported types: Vertical, Vertical Fractured, and Horizontal Multi-Fractured We
PSO拟合被用户强行终止:
+
+ %1 Optimization stopped by user:
+
+ %1 拟合已被用户停止:
+
+
Current parameters have been applied to the model.
当前参数已应用到模型。
@@ -3686,6 +4081,10 @@ Supported types: Vertical, Vertical Fractured, and Horizontal Multi-Fractured We
GA auto fitting started
遗传法自动拟合开始
+
+ %1 auto fitting started
+ %1 自动拟合开始
+
%1 fitting parameters set: MaxIterations=%2, TargetAccuracy=%3, TargetWell=%4
%1 拟合参数设置:最大迭代次数=%2,目标精度=%3,目标井=%4
@@ -3779,8 +4178,8 @@ Supported types: Vertical, Vertical Fractured, and Horizontal Multi-Fractured We
=== %1 拟合会话开始于 %2 ===
- Total fitting time: %1 s (%2 min %3 s)
- 拟合总耗时:%1 s (%2 min %3 s)
+ Total fitting time: %1 s (%2 min %3 s), total solver calls: %4
+ 拟合总耗时:%1 s (%2 min %3 s),共调用求解器 %4 次
Fitting Curve
diff --git a/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h b/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h
new file mode 100644
index 0000000..71e22de
--- /dev/null
+++ b/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h
@@ -0,0 +1,201 @@
+#ifndef NMCALCULATIONAUTOFITLM_H
+#define NMCALCULATIONAUTOFITLM_H
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include "nmCalculation_global.h"
+
+class nmDataWellBase;
+
+// 双对数曲线误差分解。total 是 LM 候选接受和排序的唯一依据,
+// 其余诊断量用于有限差分灵敏度分析和信赖域选参。
+struct AutoFitObjectiveBreakdownLM {
+ bool valid;
+ double total;
+ double pressureLoss;
+ double derivativeLoss;
+ QVector residualVector;
+ double verticalCommonBias;
+ double verticalLoss;
+ bool verticalReliable;
+ double horizontalPhysicalShift;
+ double horizontalLoss;
+ bool horizontalReliable;
+ bool registrationAmbiguous;
+ double shapeLoss;
+ double lateDerivativeSlopeBias;
+ double lateDerivativeTrendLoss;
+ bool lateDerivativeTrendReliable;
+ double coverage;
+
+ AutoFitObjectiveBreakdownLM()
+ : valid(false)
+ , total(1.0e10)
+ , pressureLoss(std::numeric_limits::quiet_NaN())
+ , derivativeLoss(std::numeric_limits::quiet_NaN())
+ , verticalCommonBias(std::numeric_limits::quiet_NaN())
+ , verticalLoss(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())
+ , lateDerivativeSlopeBias(std::numeric_limits::quiet_NaN())
+ , lateDerivativeTrendLoss(std::numeric_limits::quiet_NaN())
+ , lateDerivativeTrendReliable(false)
+ , coverage(std::numeric_limits::quiet_NaN())
+ {}
+};
+
+// 有限差分 + LM/信赖域拟合的停止原因。
+enum StopReasonLM {
+ LM_CONTINUE_OPTIMIZATION = 0,
+ LM_TARGET_ACHIEVED = 1,
+ LM_TRUE_CONVERGENCE = 2,
+ LM_LOCAL_OPTIMUM = 3,
+ LM_MAX_ITERATIONS = 4,
+ LM_USER_STOPPED = 5,
+ LM_CONSECUTIVE_FAILURES = 6,
+ LM_OPTIMIZATION_FAILED = 7
+};
+
+class NMCALCULATION_EXPORT nmCalculationAutoFitLM : public QObject
+{
+ Q_OBJECT
+
+public:
+ explicit nmCalculationAutoFitLM(QObject* parent = 0);
+ ~nmCalculationAutoFitLM();
+
+ void setTargetLogLogData(const QVector >& targetData);
+ bool startAutoFitting();
+ void stopFitting();
+ QVector getBestSolution() const;
+ double getBestFitness() const;
+ AutoFitObjectiveBreakdownLM getLastObjectiveBreakdown() const;
+ QString getLastError() const;
+ bool isRunning() const;
+ int getCurrentIteration() const;
+ int getTotalEvaluations() const;
+ void resetOptimizer();
+ void setTargetWellName(const QString& wellName);
+
+signals:
+ void progressUpdated(int iteration, double bestFitness);
+ void fittingFinished(bool success, const QString& message);
+ void bestCurveUpdated(QVector > targetData,
+ QVector > bestData,
+ int iteration,
+ double fitness);
+ void logMessageGenerated(const QString& message);
+
+private:
+ void initializeTemporaryDirectory();
+ void cleanupTemporaryDirectory();
+ bool removeDirectoryRecursively(const QString& path);
+ void cleanupOldTemporaryDirectories();
+
+ bool loadAllConfigFromDataManager();
+ void loadOptimizationConfig();
+ void loadParameterBounds();
+ void extractUserInitialValues();
+
+ // 有限差分 + LM/信赖域核心算法。
+ StopReasonLM runTrustRegionFitting();
+ bool evaluateTrustRegionPoint(const QVector& parameters,
+ double* fitness,
+ AutoFitObjectiveBreakdownLM* breakdown,
+ QVector >* curve,
+ int* elapsedMs);
+ double evaluateFitness(const QVector& parameters);
+
+ void applyParametersToDataManager(const QVector& parameters);
+ void updateReservoirParameters(const QVector& parameters);
+ void updateWellParameters(const QVector& parameters);
+ void updateWellToDataManager(nmDataWellBase* pWell);
+
+ QVector > runSolver();
+ QVector > runSolverDll();
+ bool runFinalFullSolver();
+ void saveOptimizationResult();
+ void validateAndProtectFinalResult();
+ QString getStopReasonDescription(StopReasonLM reason) const;
+ int getEnabledParameterCount() const;
+
+ // LM 运行轨迹。
+ void initializeTraceFile();
+ void closeTraceFile();
+ void writeTraceHeader();
+ void writeTraceMetaFile();
+ void writeTraceRow(int iteration,
+ int parameterIndex,
+ const QString& phase,
+ const QVector& parameters,
+ double solverObjective,
+ bool solverSuccess,
+ int elapsedMs,
+ const QString& decision,
+ const AutoFitObjectiveBreakdownLM* objectiveBreakdown = nullptr);
+ QVector buildTraceParameterVector(const QVector& selectedParameters) const;
+ void emitRunSummary(bool success, StopReasonLM finalReason);
+
+ bool validateParameters(const QVector& parameters) const;
+ bool validateLogLogData(const QVector >& logLogData) const;
+ bool validateInitialValues() const;
+ bool validateSolverResult(const QVector >& result) const;
+ double calculateLogLogCurveError(const QVector >& target,
+ const QVector >& result) const;
+
+private:
+ bool m_isRunning;
+ bool m_shouldStop;
+ int m_currentIteration;
+ QString m_lastError;
+
+ QVector m_initialValues;
+ QVector m_globalBestPosition;
+ double m_globalBestFitness;
+ AutoFitObjectiveBreakdownLM m_globalBestObjectiveBreakdown;
+ QVector > m_lastEvaluatedLogLogData;
+ QVector > m_globalBestLogLogData;
+ mutable AutoFitObjectiveBreakdownLM m_lastObjectiveBreakdown;
+ QVector > m_userInitialLogLogData;
+ AutoFitObjectiveBreakdownLM m_userInitialObjectiveBreakdown;
+
+ // 参数索引:0 k,1 skin,2 wellboreC,3 phi,4 h,5 Ct,
+ // 6 Cf,7 Swi,8 Dfc,9 fractureHalfLength。
+ QVector m_parameterSelected;
+ QVector m_parameterLower;
+ QVector m_parameterUpper;
+ QVector m_enabledParamIndices;
+ QVector > m_targetLogLogData;
+ QString m_targetWellName;
+
+ int m_maxIterations;
+ double m_targetError;
+ int m_totalEvaluations;
+ int m_successfulEvaluations;
+
+ volatile int m_evaluationInProgress;
+ int m_consecutiveFailures;
+ int m_maxConsecutiveFailures;
+
+ QVector m_userInitialSolution;
+ double m_userInitialFitness;
+ bool m_hasValidUserSolution;
+
+ QString m_tempDirectory;
+ QString m_traceRunId;
+ QString m_traceFilePath;
+ QString m_traceMetaFilePath;
+ QFile m_traceFile;
+};
+
+#endif // NMCALCULATIONAUTOFITLM_H
diff --git a/Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h b/Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h
index 9731259..cc794be 100644
--- a/Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h
+++ b/Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h
@@ -9,7 +9,6 @@
#include
#include
#include
-#include
#include "nmCalculation_global.h"
@@ -20,69 +19,6 @@ class nmDataWellBase;
class QTimer;
class QProcess;
-// 双对数曲线误差分解。该结构同时保存用于候选排序的主目标,以及用于判断
-// 曲线上下、左右和形状偏差的诊断量。total 是唯一的接受和排序依据,诊断量
-// 只参与信赖域选参,不能再次叠加到 total,否则会重复计算同一批曲线残差。
-struct AutoFitObjectiveBreakdown {
- // valid 表示本次曲线评价完整有效;无效评价统一保留 total=1e10。
- // pressureLoss 和 derivativeLoss 均在 log(value) 空间按固定网格计算。
- bool valid;
- double total;
- double pressureLoss;
- double derivativeLoss;
- // 固定目标网格上的普通对数残差。非代理搜索使用它建立完整 Jacobian,
- // 向量平方和与 total 的平方一致。
- QVector residualVector;
-
- // 上下偏差使用压力和导数残差共享的算术平均中心。
- // verticalCommonBias 为正表示模拟曲线整体偏高,为负表示整体偏低;
- // verticalReliable=false 时仍保留数值,但不能据此确定参数调整方向。
- double verticalCommonBias;
- double verticalLoss;
- bool verticalReliable;
-
- // 水平偏差在 log(time) 坐标中计算。physicalShift 为正表示模拟曲线相对
- // 目标偏右,即相同曲线特征在模拟结果中出现得更晚。
- double horizontalPhysicalShift;
- double horizontalLoss;
- bool horizontalReliable;
- // true 表示当前曲线无法可靠区分上下和左右误差;此时禁止使用两类有符号
- // 诊断量选参,但去除公共中心后的 shapeLoss 仍可用于局部选参。
- bool registrationAmbiguous;
-
- // 去除公共均值中心和可信左右偏差后剩余的整体形状误差;verticalReliable
- // 只控制能否把公共中心解释为上下参数方向,不改变 shape 的中心化公式。
- double shapeLoss;
-
- // 兼容现有 trace 列。当前非代理搜索不再单独识别或调度晚期分量。
- double lateDerivativeSlopeBias;
- double lateDerivativeTrendLoss;
- bool lateDerivativeTrendReliable;
-
- // 模拟曲线对目标固定网格的有效覆盖率,取覆盖点比例与连续 log-time
- // 跨度比例中的较小值。低于损失函数门槛时本次评价直接无效。
- double coverage;
-
- AutoFitObjectiveBreakdown()
- : valid(false)
- , total(1.0e10)
- , pressureLoss(std::numeric_limits::quiet_NaN())
- , derivativeLoss(std::numeric_limits::quiet_NaN())
- , verticalCommonBias(std::numeric_limits::quiet_NaN())
- , verticalLoss(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())
- , lateDerivativeSlopeBias(std::numeric_limits::quiet_NaN())
- , lateDerivativeTrendLoss(std::numeric_limits::quiet_NaN())
- , lateDerivativeTrendReliable(false)
- , coverage(std::numeric_limits::quiet_NaN())
- {}
-};
-
// PSO粒子结构
// 这里的 position / velocity / bestPosition 只保存“用户勾选参与拟合的参数”,
// 不是完整的 11 个储层/井筒参数。完整参数向量会在写 trace 或调用代理模型时
@@ -97,8 +33,6 @@ struct AutoFitParticle {
QVector velocity; // 速度
QVector bestPosition; // 真实求解器确认的个体最优位置
QVector guideBestPosition; // 仅用于速度更新的引导位置;不会参与真实 gbest/最终结果
- AutoFitObjectiveBreakdown currentObjectiveBreakdown; // 当前真实评价对应的误差分解
- AutoFitObjectiveBreakdown bestObjectiveBreakdown; // pbest 对应的误差分解
double fitness; // 当前适应度
double bestFitness; // 真实求解器确认的个体最优适应度
double guideBestObjective; // guideBestPosition 对应的真实或代理目标值
@@ -162,10 +96,10 @@ public:
void stopFitting();
QVector getBestSolution() const;
double getBestFitness() const;
- AutoFitObjectiveBreakdown getLastObjectiveBreakdown() const;
QString getLastError() const;
bool isRunning() const;
int getCurrentIteration() const;
+ int getTotalEvaluations() const;
void resetOptimizer();
void setPSOTargetWellName(const QString& wellName);
@@ -204,26 +138,20 @@ private:
void loadOptimizationConfig();
void loadParameterBounds();
- // ===== 自动拟合核心算法 =====
+ // ===== PSO核心算法 =====
//
- // 代理开启时保留原 PSO 筛选流程;代理关闭时使用真实求解器驱动的
- // 诊断灵敏度信赖域搜索,不依赖 pbest/gbest 速度公式。
+ // 主流程:
+ // 1. extractUserInitialValues(): 从当前项目数据中取用户已有初始解;
+ // 2. initializeSwarm(): 根据初始解和上下界生成粒子群;
+ // 3. updateParticle(): 对单个粒子跑真实求解器并计算误差;
+ // 4. updateGlobalBest(): 只用真实求解器误差更新全局最优;
+ // 5. updateVelocityAndPosition(): 按 PSO 公式推进下一代粒子。
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);
// ===== 参数应用方法 =====
//
@@ -239,7 +167,6 @@ private:
// ===== 求解器相关 =====
QVector > runSolver();
QVector> runSolverDll();
- bool runFinalFullSolver();
QVector> runSolverExe();
// ===== 数据处理 =====
@@ -294,8 +221,7 @@ private:
double surrogateObjective,
const QString& screeningDecision,
const QVector& pbestPosition,
- double pbestObjective,
- const AutoFitObjectiveBreakdown* objectiveBreakdown = nullptr);
+ double pbestObjective);
void writeIterationTraceRows();
QVector buildTraceParameterVector(const QVector& selectedParameters) const;
void resetRunSummary();
@@ -351,34 +277,30 @@ private:
bool m_isRunning; // 当前是否有一次自动拟合正在运行。
bool m_shouldStop; // 用户停止标志;主循环和求解器等待循环会定期检查它。
bool m_isPaused; // 预留暂停标志;主循环中有暂停等待逻辑。
- int m_currentIteration; // 当前自动拟合迭代序号,从 0 开始。
+ int m_currentIteration; // 当前 PSO 迭代序号,从 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; // 用户初始解对应的误差分解。
// ===== 优化配置 =====
//
// 参数索引约定:
// 0 k 渗透率;1 skin 表皮系数;2 wellboreC 井筒储集;
// 3 phi 孔隙度;4 h 储层厚度;5 Ct 综合压缩系数;
- // 6 Cf 岩石压缩系数;7 Swi 初始含水饱和度;
- // 8 Dfc 裂缝导流能力;9 fractureHalfLength 裂缝半长。
+ // 6 Cf 岩石压缩系数;7 Swi 初始含水饱和度。
// m_enabledParamIndices 保存被用户勾选的参数索引,粒子的 position 维度与它一致。
- QVector m_parameterSelected; // 完整 10 个参数是否被用户勾选参与拟合。
- QVector m_parameterLower; // 完整 10 个参数的搜索下界。
- QVector m_parameterUpper; // 完整 10 个参数的搜索上界。
- QVector m_enabledParamIndices; // 被勾选参数在完整 10 维体系中的索引。
+ QVector m_parameterSelected; // 完整 8 个参数是否被用户勾选参与拟合。
+ QVector m_parameterLower; // 完整 8 个参数的搜索下界。
+ QVector m_parameterUpper; // 完整 8 个参数的搜索上界。
+ QVector m_enabledParamIndices; // 被勾选参数在完整 8 维体系中的索引。
QVector > m_targetLogLogData; // 目标井 history log-log 曲线:time/pressure/derivative。
QString m_targetWellName; // 目标井名称;读写井参数和读取模拟曲线都依赖它。
@@ -391,7 +313,7 @@ private:
double m_socialParam; // 群体学习因子,控制粒子靠近全局 gbest 的程度。
// ===== 统计信息 =====
- int m_totalEvaluations; // 真实求解器评价总次数,包含粒子评价和方向试算。
+ int m_totalEvaluations; // 已调用真实求解器评价的粒子总数。
int m_successfulEvaluations; // 真实求解器成功且误差有效的评价次数。
QVector m_convergenceHistory; // 每代全局最优误差历史,用于收敛判断。
@@ -407,7 +329,7 @@ private:
// ===== 精英保护 =====
QVector m_userInitialSolution; // 用户初始解参数,若最终改进不足会恢复它。
double m_userInitialFitness; // 用户初始解真实误差。
- double m_improvementThreshold; // 仅用于日志区分显著改进和微小改进。
+ double m_improvementThreshold; // 最终结果相对初始解至少需要达到的改进阈值。
bool m_hasValidUserSolution; // 初始解是否成功跑过真实求解器。
int m_consecutiveFailedIterations; // 连续失败迭代次数
@@ -436,8 +358,8 @@ private:
// 这些字段只描述代理筛选和运行复盘,不参与 PSO 数学更新。
bool m_traceEnabled; // 是否写出 trace CSV/meta 文件。
QString m_traceRunId; // 本次运行 ID,作为 trace/candidate/score 文件名的一部分。
- QString m_traceFilePath; // 本次自动拟合 trace CSV 的完整路径。
- QString m_traceMetaFilePath; // 与 trace 匹配的 meta JSON 完整路径。
+ QString m_traceFilePath; // pso_baseline_trace_.csv 完整路径。
+ QString m_traceMetaFilePath; // pso_baseline_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 a8c4f76..8d10dc8 100644
--- a/Include/nmNum/nmSubWxs/nmWxAutomaticFitting.h
+++ b/Include/nmNum/nmSubWxs/nmWxAutomaticFitting.h
@@ -21,6 +21,7 @@
#include "nmDataWellBase.h"
#include "nmDataAutomaticFitting.h"
#include "nmCalculationAutoFitPSO.h"
+#include "nmCalculationAutoFitLM.h"
#include "nmWxAutomaticfittingStart.h"
#include "nmSubWxs_global.h"
@@ -109,6 +110,7 @@ private:
// 自动拟合相关成员
nmCalculationAutoFitPSO* m_autoFitterPSO;
+ nmCalculationAutoFitLM* m_autoFitterLM;
QProgressDialog* m_progressDialog;
QTimer* m_progressTimer;
bool m_autoParameterRanges;
diff --git a/Include/nmNum/nmSubWxs/nmWxAutomaticFittingStart.h b/Include/nmNum/nmSubWxs/nmWxAutomaticFittingStart.h
index 3e9067e..45aa745 100644
--- a/Include/nmNum/nmSubWxs/nmWxAutomaticFittingStart.h
+++ b/Include/nmNum/nmSubWxs/nmWxAutomaticFittingStart.h
@@ -23,9 +23,11 @@
#include
#include "nmCalculationAutoFitPSO.h"
+#include "nmCalculationAutoFitLM.h"
// 前向声明
class nmCalculationAutoFitPSO;
+class nmCalculationAutoFitLM;
class QPainter;
class QColor;
class QPaintEvent;
@@ -74,8 +76,9 @@ public:
explicit nmWxAutomaticfittingStart(QWidget *parent = 0);
~nmWxAutomaticfittingStart();
- // PSO算法接口
+ // 算法接口
void setAutoFitter(nmCalculationAutoFitPSO* autoFitter);
+ void setAutoFitter(nmCalculationAutoFitLM* autoFitter);
// 通用设置接口
@@ -151,6 +154,8 @@ private:
// 算法实例
nmCalculationAutoFitPSO* m_autoFitterPSO;
+ nmCalculationAutoFitLM* m_autoFitterLM;
+ QString m_algorithmName;
// 拟合参数
int m_maxIterations;
diff --git a/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp b/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp
new file mode 100644
index 0000000..a3bac99
--- /dev/null
+++ b/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp
@@ -0,0 +1,4396 @@
+#include "nmCalculationAutoFitLM.h"
+#include "nmCalculationDllPebiSolverTask.h"
+#include "nmDataAnalyzeManager.h"
+#include "nmDataWellBase.h"
+#include "nmDataVerticalWell.h"
+#include "nmDataVerticalFracturedWell.h"
+#include "nmDataHorizontalFracturedWell.h"
+#include "nmDataReservoir.h"
+#include "nmDataAutomaticFitting.h"
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#ifdef Q_OS_WIN
+#include
+#include
+#define DEBUG_OUT(msg) OutputDebugStringA(QString("[AutoFitLM] %1\n").arg(msg).toLocal8Bit().data())
+#endif
+
+static inline bool isFiniteNumber(double value)
+{
+#ifdef Q_OS_WIN
+ return _finite(value) != 0;
+#else
+ return std::isfinite(value);
+#endif
+}
+
+// 从嵌套 RMSE 中提取被消除的独立误差贡献。
+static double nestedRmsContribution(double reducedModelLoss,
+ double fullModelLoss)
+{
+ return qSqrt(qMax(0.0,
+ reducedModelLoss * reducedModelLoss -
+ fullModelLoss * fullModelLoss));
+}
+
+static inline void msleep(int ms)
+{
+#ifdef Q_OS_WIN
+ Sleep(ms);
+#else
+ Q_UNUSED(ms);
+#endif
+}
+
+static QString csvEscape(const QString& text)
+{
+ QString escaped = text;
+ escaped.replace("\"", "\"\"");
+ return QString("\"%1\"").arg(escaped);
+}
+
+static QString traceNumber(double value)
+{
+ return isFiniteNumber(value) ? QString::number(value, 'g', 17) : QString();
+}
+
+static QString traceParamAt(const QVector& params, int index)
+{
+ return (index >= 0 && index < params.size())
+ ? traceNumber(params[index])
+ : QString();
+}
+
+static QString jsonEscape(const QString& text)
+{
+ QString escaped = text;
+ escaped.replace("\\", "\\\\");
+ escaped.replace("\"", "\\\"");
+ escaped.replace("\b", "\\b");
+ escaped.replace("\f", "\\f");
+ escaped.replace("\n", "\\n");
+ escaped.replace("\r", "\\r");
+ escaped.replace("\t", "\\t");
+ return QString("\"%1\"").arg(escaped);
+}
+
+static QString jsonNumber(double value)
+{
+ return isFiniteNumber(value) ? QString::number(value, 'g', 17) : QString("null");
+}
+
+static QString jsonDoubleArray(const QVector& values)
+{
+ QStringList items;
+ for(int i = 0; i < values.size(); ++i) {
+ items << jsonNumber(values[i]);
+ }
+ return QString("[%1]").arg(items.join(","));
+}
+
+static QString jsonIntArray(const QVector& values)
+{
+ QStringList items;
+ for(int i = 0; i < values.size(); ++i) {
+ items << QString::number(values[i]);
+ }
+ return QString("[%1]").arg(items.join(","));
+}
+
+static QString jsonBoolArray(const QVector& values)
+{
+ QStringList items;
+ for(int i = 0; i < values.size(); ++i) {
+ items << (values[i] ? "true" : "false");
+ }
+ return QString("[%1]").arg(items.join(","));
+}
+
+static QString jsonStringArray(const QStringList& values)
+{
+ QStringList items;
+ for(int i = 0; i < values.size(); ++i) {
+ items << jsonEscape(values[i]);
+ }
+ return QString("[%1]").arg(items.join(","));
+}
+
+// 信赖域搜索统一在 [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;
+ AutoFitObjectiveBreakdownLM breakdown;
+ QVector > curve;
+ double fitness;
+ int elapsedMs;
+ bool valid;
+
+ TrustRegionEvaluation()
+ : fitness(1.0e10)
+ , elapsedMs(-1)
+ , valid(false)
+ {}
+};
+
+// LM 只使用固定长度、全部有限的普通残差。
+static bool trustRegionResidualsValid(
+ const AutoFitObjectiveBreakdownLM& breakdown)
+{
+ // 损失函数固定使用 80 个压力点和 80 个导数点。严格校验长度,避免
+ // Jacobian 沿用旧维度后访问另一候选的短残差向量。
+ if(!breakdown.valid || breakdown.residualVector.size() != 160) {
+ 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 AutoFitObjectiveBreakdownLM& breakdown,
+ double diagnosisThreshold)
+{
+ int component = TRUST_REGION_TOTAL_COMPONENT;
+ double largestLoss = diagnosisThreshold;
+
+ if(breakdown.verticalReliable &&
+ isFiniteNumber(breakdown.verticalLoss) &&
+ breakdown.verticalLoss > largestLoss) {
+ component = TRUST_REGION_VERTICAL_COMPONENT;
+ largestLoss = breakdown.verticalLoss;
+ }
+ if(breakdown.horizontalReliable &&
+ !breakdown.registrationAmbiguous &&
+ isFiniteNumber(breakdown.horizontalLoss) &&
+ breakdown.horizontalLoss > largestLoss) {
+ component = TRUST_REGION_HORIZONTAL_COMPONENT;
+ largestLoss = breakdown.horizontalLoss;
+ }
+ if(isFiniteNumber(breakdown.shapeLoss) &&
+ breakdown.shapeLoss > largestLoss) {
+ component = TRUST_REGION_SHAPE_COMPONENT;
+ }
+
+ return component;
+}
+
+// 求解选中参数对应的阻尼正规方程。上下和左右诊断量保留方向;形状没有
+// 天然正负,因此使用 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];
+ }
+}
+
+static QStringList traceParameterNames()
+{
+ QStringList names;
+ names << "k"
+ << "skin"
+ << "wellboreC"
+ << "phi"
+ << "h"
+ << "Ct"
+ << "Cf"
+ << "Swi"
+ << "Dfc"
+ << "fractureHalfLength";
+ return names;
+}
+
+nmCalculationAutoFitLM::nmCalculationAutoFitLM(QObject* parent)
+ : QObject(parent)
+ , m_isRunning(false)
+ , m_shouldStop(false)
+ , m_currentIteration(0)
+ , m_globalBestFitness(1e10)
+ , m_maxIterations(100)
+ , m_targetError(0.001)
+ , m_totalEvaluations(0)
+ , m_successfulEvaluations(0)
+ , m_evaluationInProgress(0)
+ , m_consecutiveFailures(0)
+ , m_userInitialFitness(1e10)
+ , m_maxConsecutiveFailures(3)
+ , m_hasValidUserSolution(false)
+ , m_targetWellName("")
+ , m_traceRunId("")
+ , m_traceFilePath("")
+ , m_traceMetaFilePath("")
+{
+ // LM 对象只初始化信赖域运行状态和求解器临时目录。
+ initializeTemporaryDirectory();
+ DEBUG_OUT("LM automatic fitting calculator initialized");
+}
+
+// 析构函数:停止仍在进行的拟合、关闭 trace 文件并清理临时目录。
+// 自动拟合可能在 UI 线程中被窗口关闭打断,因此析构时要尽量温和地等待当前评价结束;
+// 如果等待超时,再强制清除运行标志,避免对象销毁后还有信号回调访问成员变量。
+nmCalculationAutoFitLM::~nmCalculationAutoFitLM()
+{
+ if(m_isRunning) {
+ m_shouldStop = true;
+ int waitCount = 0;
+ while(m_isRunning && waitCount < 100) {
+ QApplication::processEvents(QEventLoop::ExcludeUserInputEvents, 50);
+ msleep(50);
+ waitCount++;
+ }
+ if(m_isRunning) {
+ DEBUG_OUT("Force stopping LM fitting after timeout");
+ m_isRunning = false;
+ }
+ }
+
+ closeTraceFile();
+ cleanupTemporaryDirectory();
+ disconnect(this, nullptr, nullptr, nullptr);
+}
+
+// ===== 临时目录工具 =====
+//
+// 真实求解器 DLL 和自动拟合过程会产生中间文件,因此每次创建独立的
+// autofit_temp__ 目录。退出时只删除本类创建的目录,启动时顺便清理
+// 旧进程遗留的 autofit_temp_*,避免长期调试后应用目录被临时文件堆满。
+
+void nmCalculationAutoFitLM::initializeTemporaryDirectory()
+{
+ // 先清理历史遗留目录,再为本次对象创建唯一目录。
+ // 目录名包含进程 ID 和毫秒时间戳,通常足够唯一;counter 是极端重名时的兜底。
+ cleanupOldTemporaryDirectories();
+
+ QString timestamp = QDateTime::currentDateTime().toString("yyyyMMdd_hhmmss_zzz");
+ QString processId = QString::number(QCoreApplication::applicationPid());
+
+ m_tempDirectory = QApplication::applicationDirPath() +
+ "/autofit_temp_" + processId + "_" + timestamp;
+
+ int counter = 0;
+ QString originalPath = m_tempDirectory;
+
+ while(QDir(m_tempDirectory).exists() && counter < 100) {
+ m_tempDirectory = originalPath + "_" + QString::number(counter);
+ counter++;
+ }
+
+ if(QDir().mkpath(m_tempDirectory)) {
+ DEBUG_OUT(QString("Initialized temp directory: %1").arg(m_tempDirectory));
+ } else {
+ DEBUG_OUT(QString("Warning: Failed to create temp directory: %1").arg(m_tempDirectory));
+ m_tempDirectory = QApplication::applicationDirPath();
+ }
+}
+
+void nmCalculationAutoFitLM::cleanupTemporaryDirectory()
+{
+ // 析构或用户停止时调用。删除失败一般是文件仍被 DLL/系统占用,
+ // 这里只记录 debug 信息,不让清理失败影响 UI 退出。
+ if(QDir(m_tempDirectory).exists()) {
+ if(removeDirectoryRecursively(m_tempDirectory)) {
+ DEBUG_OUT("Temp directory cleaned up successfully");
+ } else {
+ DEBUG_OUT("Warning: Failed to clean up temp directory completely");
+ }
+ }
+}
+
+bool nmCalculationAutoFitLM::removeDirectoryRecursively(const QString& path)
+{
+ // Qt 旧版本没有统一可用的 removeRecursively 行为时,用本函数递归删除。
+ // 调用方传入的是本类创建的临时目录或旧 autofit_temp_* 目录。
+ QDir dir(path);
+
+ if(!dir.exists()) {
+ return true;
+ }
+
+ // 递归删除子目录和文件。包含 Hidden,避免隐藏中间文件阻塞目录删除。
+ QFileInfoList entries = dir.entryInfoList(QDir::NoDotAndDotDot | QDir::AllEntries | QDir::Hidden);
+ bool allRemoved = true;
+
+ for(int i = 0; i < entries.size(); ++i) {
+ const QFileInfo& entry = entries[i];
+
+ if(entry.isDir()) {
+ if(!removeDirectoryRecursively(entry.absoluteFilePath())) {
+ allRemoved = false;
+ }
+ } else {
+ QFile file(entry.absoluteFilePath());
+
+ // 处理只读文件。某些求解器输出可能带只读属性,删除前先补写权限。
+ if(!file.permissions().testFlag(QFile::WriteUser)) {
+ file.setPermissions(file.permissions() | QFile::WriteUser);
+ }
+
+ if(!file.remove()) {
+ DEBUG_OUT(QString("Failed to remove file: %1").arg(entry.absoluteFilePath()));
+ allRemoved = false;
+ }
+ }
+ }
+
+ // 删除目录本身
+ if(allRemoved) {
+ return dir.rmdir(path);
+ }
+
+ return false;
+}
+
+void nmCalculationAutoFitLM::cleanupOldTemporaryDirectories()
+{
+ // 应用启动或新建自动拟合对象时清理旧目录。
+ // 不删除当前进程 ID 对应目录,防止同进程内多个拟合对象或正在运行的求解器被误删。
+ QString appDir = QApplication::applicationDirPath();
+ QDir dir(appDir);
+
+ // 获取当前进程ID,避免删除当前进程可能使用的目录
+ QString currentProcessId = QString::number(QCoreApplication::applicationPid());
+
+ // 查找所有以 "autofit_temp_" 开头的目录
+ QStringList filters;
+ filters << "autofit_temp_*";
+ QFileInfoList tempDirs = dir.entryInfoList(filters, QDir::Dirs | QDir::NoDotAndDotDot);
+
+ if(tempDirs.isEmpty()) {
+ DEBUG_OUT("No old temporary directories found");
+ return;
+ }
+
+ DEBUG_OUT(QString("Found %1 potential old temporary directories to clean up").arg(tempDirs.size()));
+
+ int cleanedCount = 0;
+ int failedCount = 0;
+ int skippedCount = 0;
+
+ for(int i = 0; i < tempDirs.size(); ++i) {
+ const QFileInfo& dirInfo = tempDirs[i];
+ QString dirPath = dirInfo.absoluteFilePath();
+ QString dirName = dirInfo.fileName();
+
+ // 检查是否是当前进程的目录(虽然理论上不应该存在,但为了安全起见)
+ if(dirName.contains("_" + currentProcessId + "_")) {
+ DEBUG_OUT(QString("Skipping current process directory: %1").arg(dirName));
+ skippedCount++;
+ continue;
+ }
+
+ // 尝试删除目录
+ DEBUG_OUT(QString("Attempting to remove old temp directory: %1").arg(dirName));
+
+ if(removeDirectoryRecursively(dirPath)) {
+ DEBUG_OUT(QString("Successfully cleaned up: %1").arg(dirName));
+ cleanedCount++;
+ } else {
+ DEBUG_OUT(QString("Failed to clean up: %1 (may be in use by another process)").arg(dirName));
+ failedCount++;
+ }
+ }
+
+ // 输出清理统计信息
+ DEBUG_OUT(QString("Old temp directories cleanup summary: %1 removed, %2 failed, %3 skipped")
+ .arg(cleanedCount).arg(failedCount).arg(skippedCount));
+}
+
+
+// ==================== 公共接口方法 ====================
+
+void nmCalculationAutoFitLM::setTargetLogLogData(const QVector >& targetData)
+{
+ // 目标曲线由界面层从目标井 history log-log 传入。
+ // 约定 targetData[0]=time,targetData[1]=pressure,targetData[2]=pressure derivative。
+ m_targetLogLogData = targetData;
+ DEBUG_OUT(QString("Target LogLog data set: %1 arrays").arg(targetData.size()));
+
+ if(targetData.size() >= 3) {
+ DEBUG_OUT(QString("LogLog data points: X=%1, Y1=%2, Y2=%3")
+ .arg(targetData[0].size())
+ .arg(targetData[1].size())
+ .arg(targetData[2].size()));
+ }
+}
+
+
+void nmCalculationAutoFitLM::stopFitting()
+{
+ // 用户点击停止时只设置请求标志,让 LM 主循环和求解器等待逻辑自然退出。
+ if(m_isRunning) {
+ emit logMessageGenerated(tr("=== User Stop Request Received ==="));
+ emit logMessageGenerated(tr("Gracefully stopping LM automatic fitting..."));
+ m_shouldStop = true;
+
+ // 给当前评价一个短暂的自然退出时间。若仍在运行,
+ // runSolverDll() 会在下一个等待周期检查 m_shouldStop 并结束任务。
+ int waitCount = 0;
+
+ while(m_evaluationInProgress > 0 && waitCount < 30) {
+ QApplication::processEvents(QEventLoop::ExcludeUserInputEvents, 50);
+ msleep(50);
+ waitCount++;
+ }
+
+ if(m_evaluationInProgress > 0) {
+ emit logMessageGenerated(tr("Waiting for current solver evaluation to stop..."));
+ }
+
+ emit logMessageGenerated(tr("LM automatic fitting stop request processed"));
+ } else {
+ emit logMessageGenerated(tr("Stop request received but optimization is not running"));
+ }
+}
+
+bool nmCalculationAutoFitLM::isRunning() const
+{
+ // UI 查询当前是否处于自动拟合运行状态。
+ return m_isRunning;
+}
+
+int nmCalculationAutoFitLM::getCurrentIteration() const
+{
+ // UI 进度条和日志展示用的当前迭代序号。
+ return m_currentIteration;
+}
+
+int nmCalculationAutoFitLM::getTotalEvaluations() const
+{
+ // 返回本轮拟合实际调用真实求解器的总次数,供完成日志展示。
+ return m_totalEvaluations;
+}
+
+QVector nmCalculationAutoFitLM::getBestSolution() const
+{
+ // 返回紧凑的“启用参数向量”,顺序与 m_enabledParamIndices 一致。
+ return m_globalBestPosition;
+}
+
+double nmCalculationAutoFitLM::getBestFitness() const
+{
+ // 当前全局最优真实误差。越小越好,1e10 附近通常表示尚无有效解。
+ return m_globalBestFitness;
+}
+
+AutoFitObjectiveBreakdownLM nmCalculationAutoFitLM::getLastObjectiveBreakdown() const
+{
+ // 返回最近一次损失评价的误差分解,供界面或后续优化逻辑读取。
+ return m_lastObjectiveBreakdown;
+}
+
+QString nmCalculationAutoFitLM::getLastError() const
+{
+ // 上一次失败的人类可读错误信息,主要给 UI 层弹窗或日志使用。
+ return m_lastError;
+}
+
+void nmCalculationAutoFitLM::resetOptimizer()
+{
+ // 清空一次运行产生的状态,但不销毁对象。
+ // 配置字段会在 startAutoFitting() 中重新从 DataManager 读取;
+ // trace 文件先关闭,避免新一轮 run 继续写到旧 CSV。
+ closeTraceFile();
+ m_globalBestPosition.clear();
+ m_globalBestFitness = 1e10;
+ m_globalBestObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
+ m_lastEvaluatedLogLogData.clear();
+ m_globalBestLogLogData.clear();
+ m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
+ m_userInitialLogLogData.clear();
+ m_userInitialObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
+ m_currentIteration = 0;
+ m_totalEvaluations = 0;
+ m_successfulEvaluations = 0;
+ m_lastError.clear();
+ m_initialValues.clear();
+ m_userInitialSolution.clear();
+ m_userInitialFitness = 1e10;
+ m_hasValidUserSolution = false;
+ m_traceMetaFilePath.clear();
+
+ DEBUG_OUT("LM optimizer reset");
+}
+
+void nmCalculationAutoFitLM::setTargetWellName(const QString& wellName)
+{
+ // 目标井名是贯穿拟合流程的关键索引:
+ // 读目标曲线、写 skin/wellboreC、求解后取 resultLogLog 都依赖这个名字。
+ m_targetWellName = wellName;
+}
+
+void nmCalculationAutoFitLM::initializeTraceFile()
+{
+ // LM 轨迹独立保存在应用目录,不依赖机器学习目录或外部评分进程。
+ closeTraceFile();
+ m_traceRunId = QDateTime::currentDateTime().toString("yyyyMMdd_hhmmss_zzz");
+ QDir traceDir(QDir(QApplication::applicationDirPath()).absoluteFilePath("autofit_lm_trace"));
+
+ if(!traceDir.exists() && !QDir().mkpath(traceDir.absolutePath())) {
+ DEBUG_OUT(QString("Failed to create LM trace directory: %1").arg(traceDir.absolutePath()));
+ m_traceFilePath.clear();
+ m_traceMetaFilePath.clear();
+ return;
+ }
+
+ m_traceFilePath = traceDir.absoluteFilePath(
+ QString("lm_trust_region_trace_%1.csv").arg(m_traceRunId));
+ m_traceMetaFilePath = traceDir.absoluteFilePath(
+ QString("lm_trust_region_trace_%1.meta.json").arg(m_traceRunId));
+ m_traceFile.setFileName(m_traceFilePath);
+
+ if(!m_traceFile.open(QIODevice::WriteOnly | QIODevice::Text)) {
+ DEBUG_OUT(QString("Failed to open LM trace file: %1").arg(m_traceFilePath));
+ m_traceFilePath.clear();
+ m_traceMetaFilePath.clear();
+ return;
+ }
+
+ writeTraceHeader();
+ writeTraceMetaFile();
+ emit logMessageGenerated(tr("LM fitting trace: %1").arg(m_traceFilePath));
+ if(!m_traceMetaFilePath.isEmpty()) {
+ emit logMessageGenerated(
+ tr("LM fitting trace metadata: %1").arg(m_traceMetaFilePath));
+ }
+}
+
+void nmCalculationAutoFitLM::closeTraceFile()
+{
+ if(m_traceFile.isOpen()) {
+ m_traceFile.flush();
+ m_traceFile.close();
+ }
+}
+
+void nmCalculationAutoFitLM::writeTraceHeader()
+{
+ if(!m_traceFile.isOpen()) {
+ return;
+ }
+
+ QStringList cols;
+ cols << "run_id"
+ << "iteration"
+ << "parameter_index"
+ << "phase"
+ << "k"
+ << "skin"
+ << "wellboreC"
+ << "phi"
+ << "h"
+ << "Ct"
+ << "Cf"
+ << "Swi"
+ << "Dfc"
+ << "fractureHalfLength"
+ << "solver_objective"
+ << "solver_success"
+ << "elapsed_ms"
+ << "decision"
+ << "enabled_param_indices"
+ << "pressure_loss"
+ << "derivative_loss"
+ << "vertical_common_bias"
+ << "vertical_loss"
+ << "vertical_reliable"
+ << "horizontal_physical_shift"
+ << "horizontal_loss"
+ << "horizontal_reliable"
+ << "shape_loss"
+ << "late_trend_loss"
+ << "late_slope_bias"
+ << "late_trend_reliable"
+ << "registration_ambiguous"
+ << "coverage";
+
+ QTextStream out(&m_traceFile);
+ out << cols.join(",") << "\n";
+}
+
+void nmCalculationAutoFitLM::writeTraceMetaFile()
+{
+ if(m_traceMetaFilePath.isEmpty()) {
+ return;
+ }
+
+ QFile metaFile(m_traceMetaFilePath);
+ if(!metaFile.open(QIODevice::WriteOnly | QIODevice::Text)) {
+ DEBUG_OUT(QString("Failed to open LM trace meta file: %1").arg(m_traceMetaFilePath));
+ m_traceMetaFilePath.clear();
+ return;
+ }
+
+ QStringList parameterNames = traceParameterNames();
+ QStringList enabledNames;
+ for(int i = 0; i < m_enabledParamIndices.size(); ++i) {
+ int paramIndex = m_enabledParamIndices[i];
+ enabledNames << ((paramIndex >= 0 && paramIndex < parameterNames.size())
+ ? parameterNames[paramIndex]
+ : QString::number(paramIndex));
+ }
+
+ QVector initialFullParams = buildTraceParameterVector(m_initialValues);
+ QVector targetTime = m_targetLogLogData.size() > 0
+ ? m_targetLogLogData[0] : QVector();
+ QVector targetPressure = m_targetLogLogData.size() > 1
+ ? m_targetLogLogData[1] : QVector();
+ QVector targetDerivative = m_targetLogLogData.size() > 2
+ ? m_targetLogLogData[2] : QVector();
+
+ QTextStream out(&metaFile);
+ out << "{\n";
+ out << " \"schema_version\": 1,\n";
+ out << " \"trace_type\": \"finite_difference_lm_trust_region\",\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";
+ out << " \"target\": {\n";
+ out << " \"well_name\": " << jsonEscape(m_targetWellName) << ",\n";
+ out << " \"time\": " << jsonDoubleArray(targetTime) << ",\n";
+ out << " \"pressure\": " << jsonDoubleArray(targetPressure) << ",\n";
+ out << " \"derivative\": " << jsonDoubleArray(targetDerivative) << "\n";
+ out << " },\n";
+ out << " \"lm\": {\n";
+ out << " \"max_iterations\": " << m_maxIterations << ",\n";
+ out << " \"target_error\": " << jsonNumber(m_targetError) << "\n";
+ out << " },\n";
+ out << " \"parameters\": {\n";
+ out << " \"names\": " << jsonStringArray(parameterNames) << ",\n";
+ out << " \"enabled_indices\": " << jsonIntArray(m_enabledParamIndices) << ",\n";
+ out << " \"enabled_names\": " << jsonStringArray(enabledNames) << ",\n";
+ out << " \"selected_flags\": " << jsonBoolArray(m_parameterSelected) << ",\n";
+ out << " \"lower\": " << jsonDoubleArray(m_parameterLower) << ",\n";
+ out << " \"upper\": " << jsonDoubleArray(m_parameterUpper) << ",\n";
+ out << " \"initial_selected\": " << jsonDoubleArray(m_initialValues) << ",\n";
+ out << " \"initial_full\": " << jsonDoubleArray(initialFullParams) << "\n";
+ out << " }\n";
+ out << "}\n";
+ metaFile.flush();
+ metaFile.close();
+}
+
+void nmCalculationAutoFitLM::writeTraceRow(
+ int iteration,
+ int parameterIndex,
+ const QString& phase,
+ const QVector& parameters,
+ double solverObjective,
+ bool solverSuccess,
+ int elapsedMs,
+ const QString& decision,
+ const AutoFitObjectiveBreakdownLM* objectiveBreakdown)
+{
+ if(!m_traceFile.isOpen()) {
+ return;
+ }
+
+ QVector fullParams = buildTraceParameterVector(parameters);
+ QStringList enabledIndices;
+ for(int i = 0; i < m_enabledParamIndices.size(); ++i) {
+ enabledIndices << QString::number(m_enabledParamIndices[i]);
+ }
+
+ QStringList cols;
+ cols << csvEscape(m_traceRunId)
+ << QString::number(iteration)
+ << QString::number(parameterIndex)
+ << csvEscape(phase);
+ for(int i = 0; i < 10; ++i) {
+ cols << traceParamAt(fullParams, i);
+ }
+ cols << traceNumber(solverObjective)
+ << QString::number(solverSuccess ? 1 : 0)
+ << QString::number(elapsedMs)
+ << csvEscape(decision)
+ << csvEscape(enabledIndices.join(";"));
+
+ if(objectiveBreakdown && objectiveBreakdown->valid) {
+ cols << traceNumber(objectiveBreakdown->pressureLoss)
+ << traceNumber(objectiveBreakdown->derivativeLoss)
+ << traceNumber(objectiveBreakdown->verticalCommonBias)
+ << traceNumber(objectiveBreakdown->verticalLoss)
+ << QString::number(objectiveBreakdown->verticalReliable ? 1 : 0)
+ << traceNumber(objectiveBreakdown->horizontalPhysicalShift)
+ << traceNumber(objectiveBreakdown->horizontalLoss)
+ << QString::number(objectiveBreakdown->horizontalReliable ? 1 : 0)
+ << traceNumber(objectiveBreakdown->shapeLoss)
+ << traceNumber(objectiveBreakdown->lateDerivativeTrendLoss)
+ << traceNumber(objectiveBreakdown->lateDerivativeSlopeBias)
+ << QString::number(objectiveBreakdown->lateDerivativeTrendReliable ? 1 : 0)
+ << QString::number(objectiveBreakdown->registrationAmbiguous ? 1 : 0)
+ << traceNumber(objectiveBreakdown->coverage);
+ } else {
+ for(int i = 0; i < 14; ++i) {
+ cols << QString();
+ }
+ }
+
+ QTextStream out(&m_traceFile);
+ out << cols.join(",") << "\n";
+ m_traceFile.flush();
+}
+
+void nmCalculationAutoFitLM::emitRunSummary(bool success, StopReasonLM finalReason)
+{
+ // 汇总仅描述 LM 迭代和真实求解器评价。
+ emit logMessageGenerated(tr("=== LM Run Summary ==="));
+ emit logMessageGenerated(tr("Stop reason: %1").arg(getStopReasonDescription(finalReason)));
+ emit logMessageGenerated(
+ tr("Result: %1, final error=%2, iterations=%3, evaluations=%4 (successful=%5, failed=%6)")
+ .arg(success ? tr("SUCCESS") : tr("FAILED"))
+ .arg(m_globalBestFitness, 0, 'e', 4)
+ .arg(m_currentIteration + 1)
+ .arg(m_totalEvaluations)
+ .arg(m_successfulEvaluations)
+ .arg(m_totalEvaluations - m_successfulEvaluations));
+
+ if(!m_traceFilePath.isEmpty()) {
+ emit logMessageGenerated(tr("Artifacts: trace=%1").arg(m_traceFilePath));
+ }
+ if(!m_traceMetaFilePath.isEmpty()) {
+ emit logMessageGenerated(tr("Artifacts: trace_meta=%1").arg(m_traceMetaFilePath));
+ }
+}
+
+QVector nmCalculationAutoFitLM::buildTraceParameterVector(const QVector& selectedParameters) const
+{
+ // 将 LM 内部使用的“启用参数向量”还原成完整 10 维参数向量。
+ // 未启用的参数从当前 DataManager 读取,启用的参数用 selectedParameters 覆盖。
+ // trace CSV 和 meta 使用该完整向量记录一次候选评价。
+ QVector fullParams(10, 0.0);
+
+ nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
+
+ if(dataManager) {
+ nmDataReservoir reservoirData = dataManager->getReservoirDataCopy();
+ fullParams[0] = reservoirData.getPermeability().getValue().toDouble();
+ fullParams[3] = reservoirData.getPorosity().getValue().toDouble();
+ fullParams[4] = reservoirData.getThickness().getValue().toDouble();
+ fullParams[5] = reservoirData.getCt().getValue().toDouble();
+ fullParams[6] = reservoirData.getCf().getValue().toDouble();
+ fullParams[7] = reservoirData.getSwi().getValue().toDouble();
+
+ nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName);
+
+ if(pTargetWell) {
+ nmDataPerforation* perforation = pTargetWell->getPerforation(0);
+ if(perforation) {
+ fullParams[1] = perforation->getSkin().getValue().toDouble();
+ }
+ fullParams[2] = pTargetWell->getWellboreStorage().getValue().toDouble();
+
+ // Dfc 只存在于两类压裂井,普通井在完整向量中保持为 0。
+ if(pTargetWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
+ nmDataVerticalFracturedWell* fracturedWell =
+ dynamic_cast(pTargetWell);
+ if(fracturedWell) {
+ fullParams[8] = fracturedWell->getDfc().getValue().toDouble();
+ fullParams[9] = fracturedWell->getFractureHalfLength().getValue().toDouble();
+ }
+ } else if(pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
+ nmDataHorizontalFracturedWell* fracturedWell =
+ dynamic_cast(pTargetWell);
+ if(fracturedWell) {
+ fullParams[8] = fracturedWell->getDfc().getValue().toDouble();
+ fullParams[9] = fracturedWell->getFractureHalfLength().getValue().toDouble();
+ }
+ }
+ }
+ }
+
+ for(int i = 0; i < selectedParameters.size() && i < m_enabledParamIndices.size(); ++i) {
+ int paramIndex = m_enabledParamIndices[i];
+
+ if(paramIndex >= 0 && paramIndex < fullParams.size()) {
+ fullParams[paramIndex] = selectedParameters[i];
+ }
+ }
+
+ return fullParams;
+}
+
+// ==================== 数据加载方法 ====================
+
+bool nmCalculationAutoFitLM::loadAllConfigFromDataManager()
+{
+ // 统一从 DataManager 加载本次运行所需配置。
+ // UI 层只负责把用户选择保存到 nmDataAutomaticFitting,本类从这里开始完全数据驱动。
+ try {
+ loadOptimizationConfig();
+ loadParameterBounds();
+ extractUserInitialValues(); // 直接提取初始值,无需条件判断
+ return true;
+ } catch(...) {
+ m_lastError = "Failed to load configuration from data manager";
+ return false;
+ }
+}
+
+void nmCalculationAutoFitLM::loadOptimizationConfig()
+{
+ // LM 只读取迭代次数和目标误差,其他数值控制保留在现有算法实现中。
+ nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
+ nmDataAutomaticFitting fittingData = dataManager->getAutomaticFittingDataCopy();
+
+ m_maxIterations = fittingData.getIterationCount().getValue().toInt();
+ m_targetError = fittingData.getErrorTolerance().getValue().toDouble();
+ DEBUG_OUT(QString("Loaded LM config: iterations=%1, error=%2")
+ .arg(m_maxIterations).arg(m_targetError));
+}
+
+void nmCalculationAutoFitLM::loadParameterBounds()
+{
+ // 读取用户勾选的拟合参数及上下界。
+ //
+ // 这里构建三个核心数组:
+ // - m_parameterSelected[10]:完整参数体系中每个参数是否参与拟合;
+ // - m_parameterLower/Upper[10]:完整参数体系的搜索上下界;
+ // - m_enabledParamIndices:把粒子内部紧凑向量映射回完整参数索引。
+ nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
+ nmDataAutomaticFitting fittingData = dataManager->getAutomaticFittingDataCopy();
+
+ // 获取参数选择状态
+ m_parameterSelected.resize(10);
+ m_parameterSelected[0] = fittingData.getPermeabilitySelected();
+ m_parameterSelected[1] = fittingData.getSkinSelected();
+ m_parameterSelected[2] = fittingData.getWellboreStorageSelected();
+ m_parameterSelected[3] = fittingData.getPorositySelected();
+ m_parameterSelected[4] = fittingData.getThicknessSelected();
+ m_parameterSelected[5] = fittingData.getCtSelected();
+ m_parameterSelected[6] = fittingData.getCfSelected();
+ m_parameterSelected[7] = fittingData.getSwiSelected();
+ m_parameterSelected[8] = fittingData.getFractureConductivitySelected();
+ m_parameterSelected[9] = fittingData.getFractureHalfLengthSelected();
+
+ // 获取参数边界
+ m_parameterLower.resize(10);
+ m_parameterUpper.resize(10);
+
+ m_parameterLower[0] = fittingData.getPermeabilityMin().getValue().toDouble();
+ m_parameterUpper[0] = fittingData.getPermeabilityMax().getValue().toDouble();
+
+ m_parameterLower[1] = fittingData.getSkinMin().getValue().toDouble();
+ m_parameterUpper[1] = fittingData.getSkinMax().getValue().toDouble();
+
+ m_parameterLower[2] = fittingData.getWellboreStorageMin().getValue().toDouble();
+ m_parameterUpper[2] = fittingData.getWellboreStorageMax().getValue().toDouble();
+
+ m_parameterLower[3] = fittingData.getPorosityMin().getValue().toDouble();
+ m_parameterUpper[3] = fittingData.getPorosityMax().getValue().toDouble();
+
+ m_parameterLower[4] = fittingData.getThicknessMin().getValue().toDouble();
+ m_parameterUpper[4] = fittingData.getThicknessMax().getValue().toDouble();
+
+ m_parameterLower[5] = fittingData.getCtMin().getValue().toDouble();
+ m_parameterUpper[5] = fittingData.getCtMax().getValue().toDouble();
+
+ m_parameterLower[6] = fittingData.getCfMin().getValue().toDouble();
+ m_parameterUpper[6] = fittingData.getCfMax().getValue().toDouble();
+
+ m_parameterLower[7] = fittingData.getSwiMin().getValue().toDouble();
+ m_parameterUpper[7] = fittingData.getSwiMax().getValue().toDouble();
+
+ m_parameterLower[8] = fittingData.getFractureConductivityMin().getValue().toDouble();
+ m_parameterUpper[8] = fittingData.getFractureConductivityMax().getValue().toDouble();
+
+ m_parameterLower[9] = fittingData.getFractureHalfLengthMin().getValue().toDouble();
+ m_parameterUpper[9] = fittingData.getFractureHalfLengthMax().getValue().toDouble();
+
+ // 更新启用参数索引
+ m_enabledParamIndices.clear();
+
+ for(int i = 0; i < m_parameterSelected.size(); ++i) {
+ if(m_parameterSelected[i]) {
+ m_enabledParamIndices.append(i);
+ }
+ }
+
+ DEBUG_OUT(QString("Loaded parameter bounds: %1 enabled parameters")
+ .arg(m_enabledParamIndices.size()));
+}
+
+// ==================== 自动拟合核心方法 ====================
+bool nmCalculationAutoFitLM::startAutoFitting()
+{
+ // 总入口只负责准备数据、调用有限差分 + LM/信赖域,并写回最终结果。
+ StopReasonLM finalReason = LM_CONTINUE_OPTIMIZATION;
+
+ if(m_isRunning) {
+ m_lastError = "Auto fitting is already running";
+ return false;
+ }
+
+ try {
+ if(!loadAllConfigFromDataManager()) {
+ emit logMessageGenerated(tr("ERROR: Failed to load configuration from data manager"));
+ return false;
+ }
+
+ emit logMessageGenerated(tr("Algorithm: LM"));
+ const int enabledParams = getEnabledParameterCount();
+ emit logMessageGenerated(tr("Enabled parameters count: %1").arg(enabledParams));
+
+ if(enabledParams == 0) {
+ m_lastError = "No parameters enabled for optimization";
+ emit logMessageGenerated(tr("ERROR: No parameters enabled for optimization"));
+ return false;
+ }
+
+ if(m_targetLogLogData.size() < 3) {
+ m_lastError = "Target LogLog data is empty or insufficient";
+ emit logMessageGenerated(tr("ERROR: Target LogLog data is empty or insufficient"));
+ return false;
+ }
+
+ if(m_targetLogLogData[0].size() != m_targetLogLogData[1].size() ||
+ m_targetLogLogData[0].size() != m_targetLogLogData[2].size()) {
+ m_lastError = "Target LogLog data arrays have inconsistent sizes";
+ emit logMessageGenerated(tr("ERROR: Target LogLog data arrays have inconsistent sizes"));
+ return false;
+ }
+
+ if(m_targetWellName.isEmpty()) {
+ m_lastError = "Target well name is empty";
+ emit logMessageGenerated(tr("ERROR: Target well name is empty"));
+ return false;
+ }
+
+ emit logMessageGenerated(
+ tr("Target data validation passed (%1 data points)")
+ .arg(m_targetLogLogData[0].size()));
+
+ // resetOptimizer() 会清空运行状态,因此先保存从 DataManager 提取的初始值。
+ QVector savedInitialValues = m_initialValues;
+ resetOptimizer();
+ m_isRunning = true;
+ m_shouldStop = false;
+ m_currentIteration = 0;
+ m_consecutiveFailures = 0;
+ m_initialValues = savedInitialValues;
+ initializeTraceFile();
+
+ // 先用真实求解器评价用户当前模型,供最终精英保护使用。
+ if(!savedInitialValues.isEmpty()) {
+ m_userInitialSolution = savedInitialValues;
+ emit logMessageGenerated(tr("=== Evaluating Initial Solution (Elite Protection) ==="));
+
+ QString paramStr = tr("Initial parameters: ");
+ for(int i = 0; i < m_userInitialSolution.size(); ++i) {
+ paramStr += QString("[%1]=%2 ")
+ .arg(i).arg(m_userInitialSolution[i], 0, 'f', 6);
+ }
+ emit logMessageGenerated(paramStr);
+
+ try {
+ emit logMessageGenerated(
+ tr("Starting initial solution evaluation..."));
+ QTime initialEvalTimer;
+ initialEvalTimer.start();
+ m_totalEvaluations++;
+ m_userInitialFitness = evaluateFitness(m_userInitialSolution);
+ const int initialEvalElapsedMs = initialEvalTimer.elapsed();
+
+ if(m_userInitialFitness < 1e9) {
+ m_successfulEvaluations++;
+ m_hasValidUserSolution = true;
+ m_globalBestFitness = m_userInitialFitness;
+ 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));
+ emit bestCurveUpdated(m_targetLogLogData,
+ m_globalBestLogLogData,
+ 0,
+ m_globalBestFitness);
+ } else {
+ m_hasValidUserSolution = false;
+ emit logMessageGenerated(tr("Initial solution evaluation failed"));
+ }
+
+ writeTraceRow(-1,
+ -1,
+ "initial_solution",
+ m_userInitialSolution,
+ m_userInitialFitness,
+ m_userInitialFitness < 1e9,
+ initialEvalElapsedMs,
+ m_hasValidUserSolution ? "valid" : "invalid",
+ m_hasValidUserSolution
+ ? &m_userInitialObjectiveBreakdown : nullptr);
+ } catch(...) {
+ m_hasValidUserSolution = false;
+ emit logMessageGenerated(tr("Exception during initial solution evaluation"));
+ }
+
+ m_initialValues = savedInitialValues;
+ }
+
+ finalReason = runTrustRegionFitting();
+ validateAndProtectFinalResult();
+
+ if(!m_globalBestPosition.isEmpty() && m_globalBestObjectiveBreakdown.valid) {
+ // 精英保护之后记录最终行,保证轨迹与实际写回参数一致。
+ writeTraceRow(m_currentIteration,
+ -1,
+ "trust_region_final",
+ m_globalBestPosition,
+ m_globalBestFitness,
+ m_globalBestFitness < 1.0e9,
+ -1,
+ "final_result",
+ &m_globalBestObjectiveBreakdown);
+ }
+
+ if(m_globalBestFitness < m_targetError) {
+ finalReason = LM_TARGET_ACHIEVED;
+ }
+ } catch(const std::exception& e) {
+ m_lastError = QString(tr("Critical exception in automatic fitting: %1")).arg(e.what());
+ emit logMessageGenerated(tr("CRITICAL ERROR: %1").arg(e.what()));
+ closeTraceFile();
+ cleanupTemporaryDirectory();
+ m_isRunning = false;
+ emit fittingFinished(false, m_lastError);
+ return false;
+ } catch(...) {
+ m_lastError = tr("Unknown critical exception in automatic fitting");
+ emit logMessageGenerated(tr("CRITICAL ERROR: Unknown exception in automatic fitting"));
+ closeTraceFile();
+ cleanupTemporaryDirectory();
+ m_isRunning = false;
+ emit fittingFinished(false, m_lastError);
+ return false;
+ }
+
+ bool finalFullSolverSucceeded = true;
+ bool finalFullSolverExecuted = false;
+
+ if(!m_globalBestPosition.isEmpty()) {
+ try {
+ emit logMessageGenerated(tr("Applying optimized parameters to model..."));
+ applyParametersToDataManager(m_globalBestPosition);
+
+ // 裂缝参数会改变网格输入;标记失效后,最终求解任务会基于新快照重建网格。
+ const bool fractureGridParameterSelected =
+ (m_parameterSelected.size() > 8 && m_parameterSelected[8]) ||
+ (m_parameterSelected.size() > 9 && m_parameterSelected[9]);
+ if(fractureGridParameterSelected) {
+ nmDataAnalyzeManager* dataManager =
+ nmDataAnalyzeManager::getCurrentInstance();
+ if(!dataManager) {
+ throw std::runtime_error("Data manager is unavailable");
+ }
+ dataManager->invalidatePebiGrid();
+ }
+
+ if(m_shouldStop) {
+ emit logMessageGenerated(
+ tr("Final full-field calculation skipped after user stop"));
+ } else {
+ emit logMessageGenerated(
+ tr("Running final full-field calculation with optimized parameters..."));
+ finalFullSolverExecuted = true;
+ finalFullSolverSucceeded = runFinalFullSolver();
+
+ if(finalFullSolverSucceeded) {
+ emit logMessageGenerated(
+ tr("Final full-field calculation completed successfully"));
+ } else if(m_shouldStop) {
+ finalFullSolverExecuted = false;
+ finalFullSolverSucceeded = true;
+ emit logMessageGenerated(
+ tr("Final full-field calculation stopped by user"));
+ } else {
+ m_lastError =
+ tr("Optimized parameters were found, but the final full-field calculation failed");
+ emit logMessageGenerated(
+ tr("ERROR: Final full-field calculation failed"));
+ }
+ }
+
+ saveOptimizationResult();
+ emit logMessageGenerated(tr("=== Optimization Results ==="));
+ emit logMessageGenerated(
+ tr("Final error: %1").arg(m_globalBestFitness, 0, 'e', 4));
+ emit logMessageGenerated(
+ tr("Total iterations: %1").arg(m_currentIteration + 1));
+ emit logMessageGenerated(
+ tr("Total evaluations: %1 (successful: %2)")
+ .arg(m_totalEvaluations).arg(m_successfulEvaluations));
+
+ QString finalParams = tr("Optimized parameters: ");
+ for(int i = 0; i < m_globalBestPosition.size(); ++i) {
+ finalParams += QString("[%1]=%2 ")
+ .arg(i).arg(m_globalBestPosition[i], 0, 'f', 6);
+ }
+ emit logMessageGenerated(finalParams);
+
+ if(finalFullSolverExecuted && finalFullSolverSucceeded) {
+ emit logMessageGenerated(
+ tr("Parameters and full-field results applied successfully to data manager"));
+ } else if(!finalFullSolverExecuted) {
+ emit logMessageGenerated(
+ tr("Optimized parameters applied to data manager"));
+ }
+ } catch(const std::exception& e) {
+ finalFullSolverSucceeded = false;
+ m_lastError = QString("Failed to apply final parameters: %1").arg(e.what());
+ emit logMessageGenerated(
+ tr("ERROR: Failed to apply final parameters: %1").arg(e.what()));
+ } catch(...) {
+ finalFullSolverSucceeded = false;
+ m_lastError = "Failed to apply final parameters due to unknown error";
+ emit logMessageGenerated(
+ tr("ERROR: Unknown error applying final parameters"));
+ }
+ }
+
+ m_isRunning = false;
+ bool success = false;
+ QString message;
+
+ if(finalReason == LM_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("=== LM AUTOMATIC FITTING SUCCESSFUL ==="));
+ } else if(finalReason == LM_TRUE_CONVERGENCE) {
+ success = true;
+ message = 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(tr("=== LM AUTOMATIC FITTING CONVERGED ==="));
+ } else if(finalReason == LM_LOCAL_OPTIMUM) {
+ success = true;
+ message = 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(tr("=== LM AUTOMATIC FITTING - LOCAL OPTIMUM ==="));
+ } else if(finalReason == LM_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("=== LM AUTOMATIC FITTING - MAX ITERATIONS ==="));
+ } else if(finalReason == LM_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("=== LM AUTOMATIC FITTING STOPPED BY USER ==="));
+ } else if(finalReason == LM_CONSECUTIVE_FAILURES) {
+ message = 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(tr("=== LM AUTOMATIC FITTING FAILED ==="));
+ } else {
+ message = QString(
+ tr("Automatic fitting ended unexpectedly. Best error: %1, Iterations: %2"))
+ .arg(m_globalBestFitness, 0, 'e', 4)
+ .arg(m_currentIteration + 1);
+ emit logMessageGenerated(tr("=== LM AUTOMATIC FITTING - UNKNOWN END ==="));
+ }
+
+ if(!finalFullSolverSucceeded) {
+ success = false;
+ message = m_lastError;
+ }
+
+ emitRunSummary(success, finalReason);
+ emit progressUpdated(m_maxIterations, m_globalBestFitness);
+ QApplication::processEvents();
+ msleep(200);
+ QApplication::processEvents();
+ closeTraceFile();
+ emit fittingFinished(success, message);
+ cleanupTemporaryDirectory();
+ return success;
+}
+
+void nmCalculationAutoFitLM::extractUserInitialValues()
+{
+ // 从当前项目模型读取用户已有初始参数。
+ // 只提取用户勾选的参数,并按 m_enabledParamIndices 的顺序写入 m_initialValues。
+ // 这些值用于初始解真实评价、LM 起点和最终精英保护。
+ nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
+ nmDataReservoir reservoirData = dataManager->getReservoirDataCopy();
+ //QVector wells = dataManager->getWellDataList();
+
+ nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName);
+
+ m_initialValues.clear();
+
+ // 按照启用参数的顺序提取初始值。井参数来自目标井,储层参数来自 reservoirData。
+ for(int i = 0; i < m_enabledParamIndices.size(); ++i) {
+ int paramIndex = m_enabledParamIndices[i];
+ double initialValue = 0.0;
+
+ switch(paramIndex) {
+ case 0: // 渗透率
+ initialValue = reservoirData.getPermeability().getValue().toDouble();
+ break;
+
+ case 1: // 表皮系数
+ if(pTargetWell) {
+ initialValue = pTargetWell->getPerforation(0)->getSkin().getValue().toDouble();
+ }
+
+ break;
+
+ case 2: // 井筒储集系数
+ if(pTargetWell) {
+ initialValue = pTargetWell->getWellboreStorage().getValue().toDouble();
+ }
+
+ break;
+
+ case 3: // 孔隙度
+ initialValue = reservoirData.getPorosity().getValue().toDouble();
+ break;
+
+ case 4: // 储层厚度
+ initialValue = reservoirData.getThickness().getValue().toDouble();
+ break;
+
+ case 5: // 综合压缩系数
+ initialValue = reservoirData.getCt().getValue().toDouble();
+ break;
+
+ case 6: // 岩石压缩系数
+ initialValue = reservoirData.getCf().getValue().toDouble();
+ break;
+
+ case 7: // 初始含水饱和度
+ initialValue = reservoirData.getSwi().getValue().toDouble();
+ break;
+
+ case 8: // 裂缝导流能力
+ if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
+ nmDataVerticalFracturedWell* fracturedWell =
+ dynamic_cast(pTargetWell);
+ if(fracturedWell) {
+ initialValue = fracturedWell->getDfc().getValue().toDouble();
+ }
+ } else if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
+ nmDataHorizontalFracturedWell* fracturedWell =
+ dynamic_cast(pTargetWell);
+ if(fracturedWell) {
+ initialValue = fracturedWell->getDfc().getValue().toDouble();
+ }
+ }
+ break;
+
+ case 9: // 裂缝半长
+ if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
+ nmDataVerticalFracturedWell* fracturedWell =
+ dynamic_cast(pTargetWell);
+ if(fracturedWell) {
+ initialValue = fracturedWell->getFractureHalfLength().getValue().toDouble();
+ }
+ } else if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
+ nmDataHorizontalFracturedWell* fracturedWell =
+ dynamic_cast(pTargetWell);
+ if(fracturedWell) {
+ initialValue = fracturedWell->getFractureHalfLength().getValue().toDouble();
+ }
+ }
+ break;
+ }
+
+ m_initialValues.append(initialValue);
+ }
+
+ DEBUG_OUT(QString("Extracted %1 user initial values").arg(m_initialValues.size()));
+
+ for(int i = 0; i < m_initialValues.size(); ++i) {
+ DEBUG_OUT(QString(" Initial[%1] = %2").arg(i).arg(m_initialValues[i], 0, 'e', 3));
+ }
+
+ // 验证初始值
+ if(!validateInitialValues()) {
+ DEBUG_OUT("Warning: Some initial values are outside parameter bounds");
+ }
+}
+
+bool nmCalculationAutoFitLM::evaluateTrustRegionPoint(
+ const QVector& parameters,
+ double* fitness,
+ AutoFitObjectiveBreakdownLM* 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;
+}
+
+StopReasonLM nmCalculationAutoFitLM::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 LM_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;
+ // 误差下降至少达到绝对 1e-5 且相对当前有效基准 0.2% 才算有效改善。
+ // 更小的下降仍保留为最佳解,但不能反复清除停滞状态、延长拟合时间。
+ const double effectiveRelativeImprovement = 2.0e-3;
+ const double effectiveAbsoluteImprovement = 1.0e-5;
+ const int maximumIneffectiveSteps = 3;
+
+ // damping 是 LM 阻尼;拒绝或预测失准时增大,真实下降与预测一致时减小。
+ // 两组累计量控制 Jacobian 重建,避免长期使用已偏离当前工作点的局部模型。
+ double trustRadius = 0.12;
+ double damping = 1.0e-2;
+ int consecutiveRejectedSteps = 0;
+ int consecutiveSolverFailures = 0;
+ int acceptedSinceRebuild = 0;
+ int consecutiveIneffectiveSteps = 0;
+ double movementSinceRebuild = 0.0;
+ bool rebuildRequested = true;
+ bool modelRebuiltAtMinimumRadius = false;
+ bool stagnationConfirmationRequested = false;
+ StopReasonLM stopReason = LM_MAX_ITERATIONS;
+
+ // jacobian 的行对应固定 160 维残差,列对应用户勾选的参数。
+ // 三个 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_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 {
+ 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,
+ current.valid ? "midpoint_valid" : "midpoint_invalid",
+ current.valid ? ¤t.breakdown : nullptr);
+ if(!current.valid) {
+ m_lastError = tr("The initial solution and parameter-range midpoint are both invalid");
+ return m_shouldStop
+ ? LM_USER_STOPPED
+ : LM_OPTIMIZATION_FAILED;
+ }
+ publishAcceptedPoint(current);
+ }
+
+ restoreEvaluationState(current);
+ emit logMessageGenerated(tr("=== Starting LM Main Loop ==="));
+ emit logMessageGenerated(
+ tr("LM starting point error: %1; evaluation budget: %2")
+ .arg(current.fitness, 0, 'e', 4)
+ .arg(maximumEvaluations));
+
+ // 有效改善始终相对“上一次有效改善后的误差”累计判断,避免一连串微小
+ // 下降每次都清零计数;累计达到门槛后才开始新的有效改善基准。
+ double effectiveImprovementBaseline = current.fitness;
+ auto registerEffectiveImprovement = [&](double fitness) -> bool {
+ const double requiredImprovement = qMax(
+ effectiveAbsoluteImprovement,
+ qAbs(effectiveImprovementBaseline) *
+ effectiveRelativeImprovement);
+ const double improvement = effectiveImprovementBaseline - fitness;
+ if(improvement < requiredImprovement) {
+ return false;
+ }
+
+ effectiveImprovementBaseline = fitness;
+ consecutiveIneffectiveSteps = 0;
+ stagnationConfirmationRequested = false;
+ return true;
+ };
+
+ // 连续三次没有有效改善时只请求一次灵敏度重建。重建完成后由主循环
+ // 直接检查累计改善,仍达不到门槛就判定局部收敛,不再继续微小试探。
+ auto recordIneffectiveStep = [&]() -> bool {
+ ++consecutiveIneffectiveSteps;
+ if(consecutiveIneffectiveSteps < maximumIneffectiveSteps) {
+ return false;
+ }
+
+ if(stagnationConfirmationRequested) {
+ return true;
+ }
+
+ consecutiveIneffectiveSteps = 0;
+ stagnationConfirmationRequested = true;
+ rebuildRequested = true;
+ emit logMessageGenerated(
+ tr("No effective improvement for %1 consecutive steps; "
+ "rebuilding sensitivity model for confirmation")
+ .arg(maximumIneffectiveSteps));
+ return false;
+ };
+
+ emit logMessageGenerated(
+ tr("Effective improvement threshold: max(%1, %2% of baseline error); "
+ "%3 consecutive ineffective steps trigger convergence confirmation")
+ .arg(effectiveAbsoluteImprovement, 0, 'e', 2)
+ .arg(effectiveRelativeImprovement * 100.0, 0, 'f', 2)
+ .arg(maximumIneffectiveSteps));
+
+ if(current.fitness < m_targetError) {
+ return LM_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,
+ decision,
+ 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);
+ writeTraceRow(m_currentIteration,
+ bestProbeColumn,
+ "trust_region_sensitivity_accept",
+ current.parameters,
+ current.fitness,
+ true,
+ 0,
+ "accepted_cached_probe",
+ ¤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) {
+ const bool confirmingStagnation =
+ stagnationConfirmationRequested;
+ if(!rebuildSensitivity()) {
+ stopReason = m_shouldStop
+ ? LM_USER_STOPPED
+ : LM_LOCAL_OPTIMUM;
+ break;
+ }
+ if(current.fitness < m_targetError) {
+ stopReason = LM_TARGET_ACHIEVED;
+ break;
+ }
+ if(m_totalEvaluations >= maximumEvaluations) {
+ stopReason = LM_MAX_ITERATIONS;
+ break;
+ }
+ const bool rebuildEffective =
+ registerEffectiveImprovement(current.fitness);
+ if(confirmingStagnation && !rebuildEffective) {
+ emit logMessageGenerated(
+ tr("Sensitivity rebuild produced no effective improvement; "
+ "local convergence detected"));
+ stopReason = LM_LOCAL_OPTIMUM;
+ 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 = LM_LOCAL_OPTIMUM;
+ break;
+ }
+ trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5);
+ damping = qMin(1.0e8, damping * 4.0);
+ rebuildRequested = true;
+ if(recordIneffectiveStep()) {
+ stopReason = LM_LOCAL_OPTIMUM;
+ break;
+ }
+ 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 = LM_LOCAL_OPTIMUM;
+ break;
+ }
+ rebuildRequested = true;
+ }
+ if(recordIneffectiveStep()) {
+ stopReason = LM_LOCAL_OPTIMUM;
+ break;
+ }
+ 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,
+ "solver_invalid",
+ nullptr);
+ restoreEvaluationState(current);
+ if(consecutiveRejectedSteps >= 2) {
+ rebuildRequested = true;
+ }
+ if(consecutiveSolverFailures >= m_maxConsecutiveFailures) {
+ stopReason = LM_CONSECUTIVE_FAILURES;
+ break;
+ }
+ if(recordIneffectiveStep()) {
+ stopReason = LM_LOCAL_OPTIMUM;
+ break;
+ }
+ continue;
+ }
+
+ consecutiveSolverFailures = 0;
+ // 有效候选即使最终被拒绝,也提供了一条真实割线,可用于修正下一轮
+ // 局部模型;是否成为新工作点仍只由下面的 total 严格比较决定。
+ const AutoFitObjectiveBreakdownLM 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;
+
+ 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;
+ }
+ }
+
+ // 候选只要更优就继续作为 current 保存;是否足以解除停滞,则统一
+ // 相对上一次有效改善基准判断。拒绝和微小改善都会累计无效次数。
+ const bool effectiveImprovement =
+ registerEffectiveImprovement(current.fitness);
+ if(!effectiveImprovement && recordIneffectiveStep()) {
+ stopReason = LM_LOCAL_OPTIMUM;
+ }
+
+ writeTraceRow(m_currentIteration,
+ -1,
+ "trust_region_candidate",
+ candidate.parameters,
+ candidate.fitness,
+ true,
+ candidate.elapsedMs,
+ accepted
+ ? QString("accepted_%1").arg(componentName)
+ : QString("rejected_%1").arg(componentName),
+ &candidate.breakdown);
+
+ QString componentDisplayName = componentName;
+ if(componentName == "vertical") {
+ componentDisplayName = tr("vertical deviation");
+ } else if(componentName == "horizontal") {
+ componentDisplayName = tr("horizontal deviation");
+ } else if(componentName == "shape") {
+ componentDisplayName = tr("shape deviation");
+ } else if(componentName == "total") {
+ componentDisplayName = tr("total error");
+ }
+
+ emit logMessageGenerated(
+ tr("Iteration %1: focus=%2, parameters=%3, error=%4, result=%5")
+ .arg(iteration + 1)
+ .arg(componentDisplayName)
+ .arg(selectedColumns.size())
+ .arg(candidate.fitness, 0, 'e', 4)
+ .arg(accepted ? tr("accepted") : tr("rejected")));
+ emit progressUpdated(iteration + 1, m_globalBestFitness);
+
+ if(stopReason == LM_LOCAL_OPTIMUM) {
+ break;
+ }
+
+ if(current.fitness < m_targetError) {
+ stopReason = LM_TARGET_ACHIEVED;
+ break;
+ }
+ if(trustRadius <= minimumTrustRadius * 1.01 &&
+ consecutiveRejectedSteps >= 2) {
+ if(modelRebuiltAtMinimumRadius) {
+ stopReason = LM_LOCAL_OPTIMUM;
+ break;
+ }
+ rebuildRequested = true;
+ }
+ }
+
+ if(completedIterations > 0) {
+ m_currentIteration = completedIterations - 1;
+ }
+ restoreEvaluationState(current);
+
+ if(m_shouldStop) {
+ return LM_USER_STOPPED;
+ }
+ if(current.fitness < m_targetError) {
+ return LM_TARGET_ACHIEVED;
+ }
+ if(stopReason == LM_CONSECUTIVE_FAILURES ||
+ stopReason == LM_LOCAL_OPTIMUM ||
+ stopReason == LM_OPTIMIZATION_FAILED) {
+ return stopReason;
+ }
+ return LM_MAX_ITERATIONS;
+}
+
+double nmCalculationAutoFitLM::evaluateFitness(const QVector& parameters)
+{
+ // LM 候选评价函数,也是自动拟合最核心的闭环:
+ // 1. 校验候选参数是否在用户设置的上下界和基本物理范围内;
+ // 2. 将参数写入 DataManager 的储层/目标井对象;
+ // 3. 调用真实数值求解器,生成模拟结果;
+ // 4. 从本次求解任务读取目标井 result log-log 曲线;
+ // 5. 与目标 history log-log 曲线计算误差,误差越小代表拟合越好。
+ //
+ // 返回 1e10 表示候选评价失败或结果不可用。
+ const QString funcName = QString("evaluateError[%1]").arg(m_currentIteration);
+ static int callCount = 0;
+ callCount++;
+ m_lastEvaluatedLogLogData.clear();
+ m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
+
+ try {
+ DEBUG_OUT(QString("%1: Call #%2 - Starting evaluation with %3 parameters")
+ .arg(funcName).arg(callCount).arg(parameters.size()));
+
+ // 打印参数值
+ QString paramStr = "Parameters: ";
+
+ for(int i = 0; i < parameters.size(); ++i) {
+ paramStr += QString("[%1]=%2 ").arg(i).arg(parameters[i], 0, 'f', 6);
+ }
+
+ DEBUG_OUT(QString("%1: %2").arg(funcName).arg(paramStr));
+
+ // 1. 参数有效性检查。这里先拦截明显非法的候选,
+ // 避免把非有限数、越界值或极端危险值传给求解器。
+ if(!validateParameters(parameters)) {
+ DEBUG_OUT(QString("%1: Call #%2 - VALIDATION FAILED").arg(funcName).arg(callCount));
+
+ // 详细检查每个参数
+ for(int i = 0; i < parameters.size(); ++i) {
+ if(!isFiniteNumber(parameters[i])) {
+ DEBUG_OUT(QString(" -> Param[%1] is NOT finite: %2").arg(i).arg(parameters[i]));
+ }
+
+ if(i < m_enabledParamIndices.size()) {
+ int paramIndex = m_enabledParamIndices[i];
+
+ if(paramIndex >= 0 && paramIndex < m_parameterLower.size()) {
+ double lower = m_parameterLower[paramIndex];
+ double upper = m_parameterUpper[paramIndex];
+
+ if(parameters[i] < lower) {
+ DEBUG_OUT(QString(" -> Param[%1]=%2 < lower bound %3")
+ .arg(i).arg(parameters[i]).arg(lower));
+ }
+
+ if(parameters[i] > upper) {
+ DEBUG_OUT(QString(" -> Param[%1]=%2 > upper bound %3")
+ .arg(i).arg(parameters[i]).arg(upper));
+ }
+
+ // 检查危险值。这些条件不是严格物理模型定义,
+ // 而是工程保护:避免求解器在明显异常输入下崩溃或返回无意义曲线。
+ switch(paramIndex) {
+ case 0: // 渗透率
+ if(parameters[i] <= 1e-6) {
+ DEBUG_OUT(QString(" -> REJECTED: Permeability too small: %1").arg(parameters[i]));
+ }
+
+ break;
+
+ case 2: // 井筒储集系数
+ if(parameters[i] <= 1e-8) {
+ DEBUG_OUT(QString(" -> REJECTED: Wellbore storage too small: %1").arg(parameters[i]));
+ }
+
+ break;
+
+ case 3: // 孔隙度
+ if(parameters[i] <= 1e-4 || parameters[i] >= 0.95) {
+ DEBUG_OUT(QString(" -> REJECTED: Unrealistic porosity: %1").arg(parameters[i]));
+ }
+
+ break;
+ }
+ }
+ }
+ }
+
+ return 1e10;
+ }
+
+ DEBUG_OUT(QString("%1: Call #%2 - Parameters validated OK").arg(funcName).arg(callCount));
+
+ // 2. 数据管理器检查。后续参数写回和求解器组装都依赖当前 DataManager。
+ nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
+
+ if(!dataManager) {
+ DEBUG_OUT(QString("%1: Call #%2 - DataManager is NULL").arg(funcName).arg(callCount));
+ return 1e10;
+ }
+
+ DEBUG_OUT(QString("%1: Call #%2 - DataManager OK").arg(funcName).arg(callCount));
+
+ // 3. 应用参数。parameters 的顺序与 m_enabledParamIndices 对齐,
+ // applyParametersToDataManager() 会把它们拆分写入储层参数和目标井参数。
+ try {
+ DEBUG_OUT(QString("%1: Call #%2 - Applying parameters...").arg(funcName).arg(callCount));
+ applyParametersToDataManager(parameters);
+ DEBUG_OUT(QString("%1: Call #%2 - Parameters applied successfully").arg(funcName).arg(callCount));
+ } catch(const std::exception& e) {
+ DEBUG_OUT(QString("%1: Call #%2 - FAILED to apply parameters: %3")
+ .arg(funcName).arg(callCount).arg(e.what()));
+ return 1e10;
+ }
+
+ // Dfc 和裂缝半长属于网格输入。标记网格失效,使下一次任务基于当前参数快照重建。
+ const bool fractureGridParameterSelected =
+ (m_parameterSelected.size() > 8 && m_parameterSelected[8]) ||
+ (m_parameterSelected.size() > 9 && m_parameterSelected[9]);
+ if(fractureGridParameterSelected) {
+ dataManager->invalidatePebiGrid();
+ }
+
+ // 4. 运行求解器。真实求解器偶发失败时允许重试,避免一次 DLL 调用异常
+ // 直接让整个粒子评价失败。
+ QVector> solverResult;
+ const int maxRetries = 2;
+ bool solverSuccess = false;
+
+ for(int retry = 0; retry <= maxRetries; ++retry) {
+ if(m_shouldStop) {
+ DEBUG_OUT(QString("%1: Call #%2 - User stop requested").arg(funcName).arg(callCount));
+ return 1e10;
+ }
+
+ try {
+ DEBUG_OUT(QString("%1: Call #%2 - Solver attempt %3/%4")
+ .arg(funcName).arg(callCount).arg(retry + 1).arg(maxRetries + 1));
+
+ if(retry > 0) {
+ DEBUG_OUT(QString("%1: Call #%2 - Retry delay...").arg(funcName).arg(callCount));
+ msleep(1000);
+ }
+
+ solverResult = runSolver();
+
+ // 详细检查求解器结果
+ if(solverResult.isEmpty()) {
+ DEBUG_OUT(QString("%1: Call #%2 - Solver returned EMPTY result").arg(funcName).arg(callCount));
+ } else {
+ DEBUG_OUT(QString("%1: Call #%2 - Solver returned %3 arrays")
+ .arg(funcName).arg(callCount).arg(solverResult.size()));
+
+ for(int i = 0; i < solverResult.size(); ++i) {
+ DEBUG_OUT(QString(" -> Array[%1] size: %2").arg(i).arg(solverResult[i].size()));
+ }
+
+ if(validateSolverResult(solverResult)) {
+ DEBUG_OUT(QString("%1: Call #%2 - Solver result VALIDATED on attempt %3")
+ .arg(funcName).arg(callCount).arg(retry + 1));
+ solverSuccess = true;
+ break;
+ } else {
+ DEBUG_OUT(QString("%1: Call #%2 - Solver result VALIDATION FAILED on attempt %3")
+ .arg(funcName).arg(callCount).arg(retry + 1));
+ }
+ }
+
+ } catch(const std::exception& e) {
+ DEBUG_OUT(QString("%1: Call #%2 - Solver EXCEPTION on attempt %3: %4")
+ .arg(funcName).arg(callCount).arg(retry + 1).arg(e.what()));
+ }
+ }
+
+ if(!solverSuccess) {
+ DEBUG_OUT(QString("%1: Call #%2 - ALL SOLVER ATTEMPTS FAILED").arg(funcName).arg(callCount));
+ return 1e10;
+ }
+
+ // 5. 获取 LogLog 数据。runSolverDll() 直接从求解任务复制目标井曲线,
+ // 不再依赖 DataManager 中可能被其它井或上一粒子改写的共享结果。
+ QVector> resultLogLogData = m_lastEvaluatedLogLogData;
+
+ try {
+ if(!validateLogLogData(resultLogLogData)) {
+ DEBUG_OUT(QString("%1: Call #%2 - LogLog data VALIDATION FAILED")
+ .arg(funcName).arg(callCount));
+
+ // 详细输出LogLog数据问题
+ if(resultLogLogData.size() < 3) {
+ DEBUG_OUT(QString(" -> LogLog arrays count: %1 (need 3)")
+ .arg(resultLogLogData.size()));
+ } else {
+ DEBUG_OUT(QString(" -> LogLog array sizes: X=%1, Y1=%2, Y2=%3")
+ .arg(resultLogLogData[0].size())
+ .arg(resultLogLogData[1].size())
+ .arg(resultLogLogData[2].size()));
+
+ // 检查数据有效性
+ for(int i = 0; i < qMin(5, resultLogLogData[0].size()); ++i) {
+ if(!isFiniteNumber(resultLogLogData[0][i]) ||
+ !isFiniteNumber(resultLogLogData[1][i]) ||
+ !isFiniteNumber(resultLogLogData[2][i])) {
+ DEBUG_OUT(QString(" -> Invalid data at index %1: X=%2, Y1=%3, Y2=%4")
+ .arg(i)
+ .arg(resultLogLogData[0][i])
+ .arg(resultLogLogData[1][i])
+ .arg(resultLogLogData[2][i]));
+ }
+ }
+ }
+
+ return 1e10;
+ }
+
+ DEBUG_OUT(QString("%1: Call #%2 - LogLog data validated, size: %3")
+ .arg(funcName).arg(callCount).arg(resultLogLogData[0].size()));
+
+ } catch(const std::exception& e) {
+ DEBUG_OUT(QString("%1: Call #%2 - Error getting LogLog data: %3")
+ .arg(funcName).arg(callCount).arg(e.what()));
+ return 1e10;
+ }
+
+ // 6. 计算误差。这里比较的是目标井 history log-log 与当前模拟 result log-log。
+ double error;
+
+ try {
+ DEBUG_OUT(QString("%1: Call #%2 - Calculating error...")
+ .arg(funcName).arg(callCount));
+
+ error = calculateLogLogCurveError(m_targetLogLogData, resultLogLogData);
+
+ if(!isFiniteNumber(error) || error < 0) {
+ DEBUG_OUT(QString("%1: Call #%2 - INVALID error value: %3")
+ .arg(funcName).arg(callCount).arg(error));
+ return 1e10;
+ }
+
+ DEBUG_OUT(QString("%1: Call #%2 - SUCCESS! Error = %3")
+ .arg(funcName).arg(callCount).arg(error, 0, 'e', 6));
+ // 保存最后一次有效曲线,供 LM 候选评价和精英保护复用。
+ m_lastEvaluatedLogLogData = resultLogLogData;
+
+ } catch(const std::exception& e) {
+ DEBUG_OUT(QString("%1: Call #%2 - Error calculation FAILED: %3")
+ .arg(funcName).arg(callCount).arg(e.what()));
+ return 1e10;
+ }
+
+ return error;
+
+ } catch(const std::exception& e) {
+ DEBUG_OUT(QString("%1: Call #%2 - TOP-LEVEL EXCEPTION: %3")
+ .arg(funcName).arg(callCount).arg(e.what()));
+ return 1e10;
+ } catch(...) {
+ DEBUG_OUT(QString("%1: Call #%2 - UNKNOWN TOP-LEVEL EXCEPTION")
+ .arg(funcName).arg(callCount));
+ return 1e10;
+ }
+}
+
+// ==================== 参数应用方法 ====================
+
+void nmCalculationAutoFitLM::applyParametersToDataManager(const QVector& parameters)
+{
+ // 将粒子的“启用参数向量”写回项目数据。
+ // parameters 的维度必须等于用户勾选的参数数量,顺序由 m_enabledParamIndices 决定。
+ // 这里不直接跑求解器,只负责把 DataManager 调整到该粒子对应的模型状态。
+ if(parameters.size() != getEnabledParameterCount()) {
+ return;
+ }
+
+ updateReservoirParameters(parameters);
+ updateWellParameters(parameters);
+}
+
+void nmCalculationAutoFitLM::updateReservoirParameters(const QVector& parameters)
+{
+ // 更新储层级参数。井级参数 skin/wellboreC 不在这里改,由 updateWellParameters() 负责。
+ // 这里先取 DataManager 中 reservoir 的副本,修改后再整体写回 DataManager。
+ nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
+ nmDataReservoir reservoirData = dataManager->getReservoirDataCopy();
+
+ // paramIndex 是粒子 position 中的索引;i 是完整 10 个参数体系中的索引。
+ // 只有 m_parameterSelected[i] 为 true 时,才从 parameters 中消费一个值。
+ int paramIndex = 0;
+
+ for(int i = 0; i < m_parameterSelected.size(); ++i) {
+ if(m_parameterSelected[i] && paramIndex < parameters.size()) {
+ double value = parameters[paramIndex];
+
+ switch(i) {
+ case 0: // 渗透率
+ reservoirData.getPermeability().setValue(value);
+ break;
+
+ case 3: // 孔隙度
+ reservoirData.getPorosity().setValue(value);
+ break;
+
+ case 4: // 储层厚度
+ reservoirData.getThickness().setValue(value);
+ break;
+
+ case 5: // 综合压缩系数
+ reservoirData.getCt().setValue(value);
+ break;
+
+ case 6: // 岩石压缩系数
+ reservoirData.getCf().setValue(value);
+ break;
+
+ case 7: // 初始含水饱和度
+ reservoirData.getSwi().setValue(value);
+ break;
+ }
+
+ paramIndex++;
+ }
+ }
+
+ // 更新数据管理器
+ dataManager->updateReservoirData(reservoirData);
+}
+
+void nmCalculationAutoFitLM::updateWellParameters(const QVector& parameters)
+{
+ // 更新目标井上的拟合参数。目前井级可拟合参数主要是:
+ // - skin:写入第一个 perforation;
+ // - wellboreC:写入井筒储集系数;
+ // - Dfc/裂缝半长:只写入垂直压裂井或多段压裂水平井。
+ // 如果目标井不存在或没有射孔数据,这里只记录 debug,不抛异常。
+ nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
+
+ // 只更新当前目标井的井参数,避免多井项目中误改其他井。
+ //QVector wells = dataManager->getWellDataList();
+ nmDataWellBase* pWell = dataManager->findWellByName(m_targetWellName);
+
+ if(!pWell) return;
+
+ //nmDataWellBase* pWell = wells[0]; // 使用第一口井
+
+ int paramIndex = 0;
+
+ for(int i = 0; i < m_parameterSelected.size(); ++i) {
+ if(m_parameterSelected[i] && paramIndex < parameters.size()) {
+ double value = parameters[paramIndex];
+
+ switch(i) {
+ case 1: { // 表皮系数
+ nmDataPerforation* perf = pWell->getPerforation(0);
+
+ if(perf) {
+ nmDataAttribute skinAttr = perf->getSkin();
+ skinAttr.setValue(value);
+ perf->setSkin(skinAttr);
+ }
+ }
+ break;
+
+ case 2: { // 井筒储集系数
+ nmDataAttribute wellboreAttr = pWell->getWellboreStorage();
+ wellboreAttr.setValue(value);
+ pWell->setWellboreStorage(wellboreAttr);
+ }
+ break;
+
+ case 8: { // 裂缝导流能力
+ if(pWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
+ nmDataVerticalFracturedWell* fracturedWell =
+ dynamic_cast(pWell);
+ if(fracturedWell) {
+ nmDataAttribute dfc = fracturedWell->getDfc();
+ dfc.setValue(value);
+ fracturedWell->setDfc(dfc);
+ }
+ } else if(pWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
+ nmDataHorizontalFracturedWell* fracturedWell =
+ dynamic_cast(pWell);
+ if(fracturedWell) {
+ nmDataAttribute dfc = fracturedWell->getDfc();
+ dfc.setValue(value);
+ fracturedWell->setDfc(dfc);
+ }
+ }
+ }
+ break;
+
+ case 9: { // 裂缝半长
+ // 直接修改井对象中的属性,复用已有信号重算裂缝端点。
+ if(pWell->getWellType() == NM_WELL_MODEL::Vertical_Fractured_Well) {
+ nmDataVerticalFracturedWell* fracturedWell =
+ dynamic_cast(pWell);
+ if(fracturedWell) {
+ fracturedWell->getFractureHalfLength().setValue(value);
+ }
+ } else if(pWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
+ nmDataHorizontalFracturedWell* fracturedWell =
+ dynamic_cast(pWell);
+ if(fracturedWell) {
+ fracturedWell->getFractureHalfLength().setValue(value);
+ }
+ }
+ }
+ break;
+ }
+
+ paramIndex++;
+ }
+ }
+
+ // 根据井类型更新到数据管理器
+ updateWellToDataManager(pWell);
+}
+
+void nmCalculationAutoFitLM::updateWellToDataManager(nmDataWellBase* pWell)
+{
+ if(!pWell) return;
+
+ nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
+ NM_WELL_MODEL wellType = pWell->getWellType();
+
+ // DataManager 内部按井型维护不同容器。修改基类指针后,需要根据实际井型
+ // 调用对应 update 接口,才能让后续求解器组装读到最新 skin / wellboreC。
+ switch(wellType) {
+ case NM_WELL_MODEL::Vertical_Well: {
+ nmDataVerticalWell* pVerticalWell = dynamic_cast(pWell);
+
+ if(pVerticalWell) {
+ QVector wells;
+ wells.append(*pVerticalWell);
+ dataManager->updateVerticalWells(wells);
+ }
+
+ break;
+ }
+
+ case NM_WELL_MODEL::Vertical_Fractured_Well: {
+ nmDataVerticalFracturedWell* pVFracturedWell = dynamic_cast(pWell);
+
+ if(pVFracturedWell) {
+ QVector wells;
+ wells.append(*pVFracturedWell);
+ dataManager->updateVerticalFracturedWells(wells);
+ }
+
+ break;
+ }
+
+ case NM_WELL_MODEL::Horizontal_Fractured_Well: {
+ nmDataHorizontalFracturedWell* pHFracturedWell = dynamic_cast(pWell);
+
+ if(pHFracturedWell) {
+ QVector wells;
+ wells.append(*pHFracturedWell);
+ dataManager->updateHorizontalFracturedWells(wells);
+ }
+
+ break;
+ }
+
+ default:
+ break;
+ }
+}
+
+// ==================== 求解器相关方法 ====================
+
+QVector> nmCalculationAutoFitLM::runSolver()
+{
+ // 真实求解器统一走 DLL 方式,返回值由 evaluateFitness() 继续校验。
+ return runSolverDll();
+}
+
+// ==================== 数据处理方法 ====================
+// ==================== 算法辅助方法 ====================
+
+void nmCalculationAutoFitLM::saveOptimizationResult()
+{
+ // 当前函数只做日志记录。真正把最优参数写回项目数据的是
+ // startAutoFitting() 结束阶段的 applyParametersToDataManager(m_globalBestPosition)。
+ DEBUG_OUT(QString("Optimization result: error=%1, evaluations=%2/%3")
+ .arg(m_globalBestFitness, 0, 'e', 4)
+ .arg(m_successfulEvaluations)
+ .arg(m_totalEvaluations));
+}
+
+void nmCalculationAutoFitLM::validateAndProtectFinalResult()
+{
+ // 最终精英保护只阻止无效结果或真正变差的结果。任何真实误差下降都应保留,
+ // 不能再用固定百分比门槛把已经找到的更优解恢复成初始值。
+ if(!m_hasValidUserSolution) {
+ emit logMessageGenerated(tr("No initial solution for elite protection"));
+ return;
+ }
+
+ emit logMessageGenerated(tr("=== Final Result Validation (Elite Protection) ==="));
+
+ // 使用已有的评估结果
+ double finalFitness = m_globalBestFitness;
+ double initialFitness = m_userInitialFitness;
+
+ emit logMessageGenerated(tr("Comparing results: Initial=%1, Final=%2")
+ .arg(initialFitness, 0, 'e', 4).arg(finalFitness, 0, 'e', 4));
+
+ 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(!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 - solution is not worse than initial"));
+ }
+}
+
+// ==================== 工具方法 ====================
+
+int nmCalculationAutoFitLM::getEnabledParameterCount() const
+{
+ // 返回粒子维度,即用户勾选参与拟合的参数数量。
+ int count = 0;
+
+ for(int i = 0; i < m_parameterSelected.size(); ++i) {
+ if(m_parameterSelected[i]) count++;
+ }
+
+ return count;
+}
+
+// ==================== 验证和处理方法 ====================
+
+bool nmCalculationAutoFitLM::validateParameters(const QVector& parameters) const
+{
+ // 参数物理范围已由拟合窗口统一校验;候选评价只检查维度、有限数和
+ // 用户设置的上下界,避免另一套硬编码阈值与实际搜索范围冲突。
+ if(parameters.size() != getEnabledParameterCount()) {
+ return false;
+ }
+
+ for(int i = 0; i < parameters.size(); ++i) {
+ if(!isFiniteNumber(parameters[i])) {
+ return false;
+ }
+
+ // 检查参数范围
+ if(i < m_enabledParamIndices.size()) {
+ int paramIndex = m_enabledParamIndices[i];
+
+ if(paramIndex >= 0 && paramIndex < m_parameterLower.size()) {
+ if(parameters[i] < m_parameterLower[paramIndex] ||
+ parameters[i] > m_parameterUpper[paramIndex]) {
+ return false;
+ }
+ }
+ }
+ }
+
+ return true;
+}
+
+bool nmCalculationAutoFitLM::validateLogLogData(const QVector>& logLogData) const
+{
+ // 校验双对数曲线结构。约定:
+ // logLogData[0]=time,logLogData[1]=pressure,logLogData[2]=pressure derivative。
+ // 三列必须长度一致,且至少有足够点数用于插值和误差计算。
+ if(logLogData.size() < 3) {
+ DEBUG_OUT("LogLog data has less than 3 arrays");
+ return false;
+ }
+
+ // 检查数组大小一致性
+ int size = logLogData[0].size();
+
+ if(size == 0) {
+ DEBUG_OUT("Empty LogLog data");
+ return false;
+ }
+
+ if(logLogData[1].size() != size || logLogData[2].size() != size) {
+ DEBUG_OUT(QString("LogLog data size mismatch: X=%1, Y1=%2, Y2=%3")
+ .arg(logLogData[0].size())
+ .arg(logLogData[1].size())
+ .arg(logLogData[2].size()));
+ return false;
+ }
+
+ // 检查最小数据点数
+ if(size < 5) {
+ DEBUG_OUT(QString("Too few LogLog data points: %1").arg(size));
+ return false;
+ }
+
+ // 数据有效性检查
+ for(int i = 0; i < size; ++i) {
+ if(!isFiniteNumber(logLogData[0][i]) ||
+ !isFiniteNumber(logLogData[1][i]) ||
+ !isFiniteNumber(logLogData[2][i])) {
+ DEBUG_OUT(QString("Invalid LogLog data at index %1").arg(i));
+ return false;
+ }
+ }
+
+ return true;
+}
+
+bool nmCalculationAutoFitLM::validateInitialValues() const
+{
+ // 检查当前模型读取出的初始参数是否和用户勾选维度一致,并且在上下界内。
+ // 如果初始值越界,算法仍可继续,但日志会提示,因为精英保护可能不可用或效果变差。
+ if(m_initialValues.size() != m_enabledParamIndices.size()) {
+ DEBUG_OUT("Initial values count mismatch with enabled parameters");
+ return false;
+ }
+
+ bool allValid = true;
+
+ for(int i = 0; i < m_initialValues.size(); ++i) {
+ int paramIndex = m_enabledParamIndices[i];
+ double value = m_initialValues[i];
+
+ if(!isFiniteNumber(value)) {
+ DEBUG_OUT(QString("Initial value[%1] is not finite: %2").arg(i).arg(value));
+ allValid = false;
+ continue;
+ }
+
+ if(paramIndex < m_parameterLower.size() && paramIndex < m_parameterUpper.size()) {
+ double minVal = m_parameterLower[paramIndex];
+ double maxVal = m_parameterUpper[paramIndex];
+
+ if(value < minVal || value > maxVal) {
+ DEBUG_OUT(QString("Initial value[%1] = %2 is outside bounds [%3, %4]")
+ .arg(i).arg(value).arg(minVal).arg(maxVal));
+ allValid = false;
+ }
+ }
+ }
+
+ return allValid;
+}
+
+bool nmCalculationAutoFitLM::validateSolverResult(const QVector>& result) const
+{
+ // 校验求解器压力结果。这里检查的是 pressure result,至少需要 time 和 pressure 两列。
+ // result log-log 的结构会在 validateLogLogData() 中另行检查。
+ if(result.size() < 2) {
+ DEBUG_OUT("Solver result has less than 2 arrays");
+ return false;
+ }
+
+ if(result[0].size() != result[1].size()) {
+ DEBUG_OUT(QString("Size mismatch: X=%1, Y=%2").arg(result[0].size()).arg(result[1].size()));
+ return false;
+ }
+
+ if(result[0].size() == 0) {
+ DEBUG_OUT("Empty solver result");
+ return false;
+ }
+
+ // 检查最小数据点数
+ if(result[0].size() < 10) {
+ DEBUG_OUT(QString("Too few data points: %1").arg(result[0].size()));
+ return false;
+ }
+
+ // 数据有效性检查
+ for(int i = 0; i < result[0].size(); ++i) {
+ if(!isFiniteNumber(result[0][i]) || !isFiniteNumber(result[1][i])) {
+ DEBUG_OUT(QString("Invalid data at index %1: X=%2, Y=%3")
+ .arg(i).arg(result[0][i]).arg(result[1][i]));
+ return false;
+ }
+ }
+
+ // 检查X值单调性
+ bool isMonotonic = true;
+
+ for(int i = 1; i < result[0].size(); ++i) {
+ if(result[0][i] <= result[0][i - 1]) {
+ isMonotonic = false;
+ break;
+ }
+ }
+
+ if(!isMonotonic) {
+ DEBUG_OUT("X values are not monotonically increasing");
+ }
+
+ return true;
+}
+
+double nmCalculationAutoFitLM::calculateLogLogCurveError(
+ const QVector >& target,
+ const QVector >& result) const
+{
+ // 主目标只比较固定网格上的压力和导数残差;上下、左右和形状只负责诊断
+ // 误差来源和选择参数,避免同一残差在 total 中被重复计算。整个计算过程均
+ // 位于 log(time)-log(value) 坐标,因此得到的是相对尺度偏差而非原始压力量纲。
+ const double invalidLoss = 1.0e10;
+ const double valueFloor = 1.0e-12;
+ const double minimumCoverage = 0.95;
+ const int numPoints = 80;
+ m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
+
+ if(!validateLogLogData(target) || !validateLogLogData(result)) {
+ return invalidLoss;
+ }
+
+ try {
+ // 无法比较的采样行先跳过;有限但非正的导数无法进入双对数空间,
+ // 当前数据又没有逐点有效掩码,因此遇到这种导数时判本次评价无效。
+ auto prepareCurve = [valueFloor](const QVector >& data,
+ int firstIndex,
+ QVector* pressure,
+ QVector* derivative) -> bool {
+ if(!pressure || !derivative || data.size() < 3 ||
+ data[0].size() != data[1].size() ||
+ data[0].size() != data[2].size() ||
+ firstIndex < 0 || firstIndex >= data[0].size()) {
+ return false;
+ }
+
+ for(int i = firstIndex; i < data[0].size(); ++i) {
+ 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(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();
+ });
+
+ QVector unique;
+ unique.reserve(curve->size());
+ for(int i = 0; i < curve->size(); ++i) {
+ 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;
+ };
+
+ sortAndUnique(pressure);
+ sortAndUnique(derivative);
+ return pressure->size() >= 3 && derivative->size() >= 3;
+ };
+
+ QVector targetPressure;
+ QVector targetDerivative;
+ QVector resultPressure;
+ QVector resultDerivative;
+ // 模拟结果已跳过 DLL 首点,因此误差计算同步忽略目标曲线首点。
+ if(!prepareCurve(target, 1, &targetPressure, &targetDerivative) ||
+ !prepareCurve(result, 0, &resultPressure, &resultDerivative)) {
+ return invalidLoss;
+ }
+
+ // 在双对数坐标中插值。二分定位用于后面的多次水平配准试算。
+ 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 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;
+ }
+ }
+
+ 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));
+
+ if(qAbs(denominator) <= 1.0e-12) {
+ *value = leftLogY;
+ } else {
+ double ratio = (qLn(x) - leftLogX) / denominator;
+ *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);
+ };
+
+ const double targetMinX = targetPressure.first().x();
+ 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;
+ }
+
+ QVector commonX(numPoints);
+ QVector commonLogX(numPoints);
+ QVector targetLogPressure(numPoints);
+ QVector targetLogDerivative(numPoints);
+ 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);
+ commonLogX[i] = logX;
+ // 首尾直接使用原始端点,避免 exp(log(t)) 的舍入误差越过严格插值边界。
+ if(i == 0) {
+ commonX[i] = targetMinX;
+ } else if(i == numPoints - 1) {
+ commonX[i] = targetMaxX;
+ } else {
+ commonX[i] = qExp(logX);
+ }
+
+ if(!interpolateLogValue(
+ targetPressure, commonX[i],
+ &targetLogPressure[i]) ||
+ !interpolateLogValue(
+ targetDerivative, commonX[i],
+ &targetLogDerivative[i])) {
+ return invalidLoss;
+ }
+ }
+
+ 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;
+ }
+
+ double resultLogPressure = 0.0;
+ double resultLogDerivative = 0.0;
+ if(!interpolateLogValue(
+ resultPressure, commonX[i],
+ &resultLogPressure) ||
+ !interpolateLogValue(
+ resultDerivative, commonX[i],
+ &resultLogDerivative)) {
+ // 已位于结果时间范围内却无法插值说明数据存在内部断点,不能补线。
+ return invalidLoss;
+ }
+
+ pressureResidual[i] =
+ resultLogPressure - targetLogPressure[i];
+ derivativeResidual[i] =
+ resultLogDerivative - targetLogDerivative[i];
+ ++supportedCount;
+ if(firstSupported < 0) {
+ firstSupported = i;
+ }
+ lastSupported = i;
+ }
+
+ // 覆盖率同时约束“有效点数量”和“连续时间跨度”。取两者较小值可避免
+ // 点数很多但只集中在局部时段的候选被误认为覆盖充分。
+ AutoFitObjectiveBreakdownLM breakdown;
+ breakdown.coverage = supportedCount > 0
+ ? static_cast(supportedCount) /
+ numPoints
+ : 0.0;
+ if(firstSupported >= 0 && lastSupported >= firstSupported) {
+ double targetSpan =
+ qMax(1.0e-12, targetLogMaxX - targetLogMinX);
+ double coveredSpan =
+ commonLogX[lastSupported] - commonLogX[firstSupported];
+ breakdown.coverage = qMin(
+ breakdown.coverage,
+ qMax(0.0, coveredSpan / targetSpan));
+ }
+
+ // 覆盖率只判断候选是否有效,不再加入固定惩罚,避免所有误差被整体抬高。
+ if(breakdown.coverage < minimumCoverage) {
+ breakdown.total = invalidLoss;
+ m_lastObjectiveBreakdown = breakdown;
+ return invalidLoss;
+ }
+
+ // 通过门槛后最多只缺少首尾少量目标点。按模拟曲线端点趋势补齐后,
+ // 每个候选仍在固定 80 点上计算均方根误差,不能靠少算难拟合端点获益。
+ for(int i = 0; i < numPoints; ++i) {
+ 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];
+ }
+
+ // 在指定中心附近计算普通均方根误差。
+ auto rmseAround = [](
+ 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 = validBegin; i < validEnd; ++i) {
+ if(!isFiniteNumber(values[i])) {
+ continue;
+ }
+
+ double difference = values[i] - center;
+ sum += difference * difference;
+ ++count;
+ }
+
+ return count > 0
+ ? qSqrt(sum / count)
+ : std::numeric_limits::quiet_NaN();
+ };
+
+ auto rmse = [&rmseAround](
+ const QVector& values,
+ int begin,
+ int end) -> double {
+ return rmseAround(values, begin, end, 0.0);
+ };
+
+ // 普通算术平均中心保留上下偏差的符号。
+ auto meanCenterRange = [](
+ 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 = validBegin; i < validEnd; ++i) {
+ if(isFiniteNumber(values[i])) {
+ center += values[i];
+ ++count;
+ }
+ }
+ if(count == 0) {
+ return std::numeric_limits::quiet_NaN();
+ }
+ return center / count;
+ };
+
+ // 压力和导数合并后只求一个公共中心,表示两条曲线共同的上下位移。
+ // 分别去中心会把压力与导数之间真实的相对形状差异一并消除。
+ auto commonMeanCenterRange = [&meanCenterRange](
+ 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 meanCenterRange(combined, 0, combined.size());
+ };
+
+ // 两个通道按能量等权合并,返回值与单通道 RMSE 保持同一量纲。
+ auto jointRmseAround = [&rmseAround](
+ const QVector& pressureValues,
+ const QVector& derivativeValues,
+ int begin,
+ int end,
+ double center) -> double {
+ double pressureLoss = rmseAround(
+ pressureValues, begin, end, center);
+ double derivativeLoss = rmseAround(
+ 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 =
+ rmse(pressureResidual, 0, numPoints);
+ breakdown.derivativeLoss =
+ rmse(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) {
+ breakdown.residualVector.append(
+ residualScale *
+ pressureResidual[i]);
+ }
+ for(int i = 0; i < numPoints; ++i) {
+ breakdown.residualVector.append(
+ residualScale *
+ 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;
+ }
+
+ // 对数时间还原后再次限制到原始端点,避免 exp(log(t)) 的
+ // 舍入误差越过严格插值边界。
+ double shiftedX = qBound(
+ resultMinX,
+ qExp(qBound(resultLogMinX,
+ shiftedLogX,
+ resultLogMaxX)),
+ resultMaxX);
+ 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