|
|
#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 <QApplication>
|
|
|
#include <QDebug>
|
|
|
#include <QTime>
|
|
|
#include <QDir>
|
|
|
#include <QTextStream>
|
|
|
#include <QFileInfo>
|
|
|
#include <QDateTime>
|
|
|
#include <QtCore/qmath.h>
|
|
|
#include <cmath>
|
|
|
#include <algorithm>
|
|
|
#include <limits>
|
|
|
|
|
|
#ifdef Q_OS_WIN
|
|
|
#include <windows.h>
|
|
|
#include <float.h>
|
|
|
#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<double>& 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<double>& 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<int>& 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<bool>& 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<double> parameters;
|
|
|
QVector<double> coordinates;
|
|
|
AutoFitObjectiveBreakdownLM breakdown;
|
|
|
QVector<QVector<double> > 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<double>& 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<double>& left,
|
|
|
const QVector<double>& 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<QVector<double> > matrix,
|
|
|
QVector<double> rightHandSide,
|
|
|
QVector<double>* 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<QVector<double> >& 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<QVector<double> >* jacobian,
|
|
|
const QVector<double>& oldResidual,
|
|
|
const QVector<double>& newResidual,
|
|
|
const QVector<double>& 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<double>* gradient,
|
|
|
double oldValue,
|
|
|
double newValue,
|
|
|
const QVector<double>& 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_<pid>_<timestamp> 目录。退出时只删除本类创建的目录，启动时顺便清理
|
|
|
// 旧进程遗留的 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<QVector<double> >& 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<double> 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(QDir::tempPath()).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";
|
|
|
|
|
|
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<double> initialFullParams = buildTraceParameterVector(m_initialValues);
|
|
|
QVector<double> targetTime = m_targetLogLogData.size() > 0
|
|
|
? m_targetLogLogData[0] : QVector<double>();
|
|
|
QVector<double> targetPressure = m_targetLogLogData.size() > 1
|
|
|
? m_targetLogLogData[1] : QVector<double>();
|
|
|
QVector<double> targetDerivative = m_targetLogLogData.size() > 2
|
|
|
? m_targetLogLogData[2] : QVector<double>();
|
|
|
|
|
|
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<double>& parameters,
|
|
|
double solverObjective,
|
|
|
bool solverSuccess,
|
|
|
int elapsedMs,
|
|
|
const QString& decision,
|
|
|
const AutoFitObjectiveBreakdownLM* objectiveBreakdown)
|
|
|
{
|
|
|
if(!m_traceFile.isOpen()) {
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
QVector<double> 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);
|
|
|
} else {
|
|
|
for(int i = 0; i < 13; ++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<double> nmCalculationAutoFitLM::buildTraceParameterVector(const QVector<double>& selectedParameters) const
|
|
|
{
|
|
|
// 将 LM 内部使用的“启用参数向量”还原成完整 10 维参数向量。
|
|
|
// 未启用的参数从当前 DataManager 读取，启用的参数用 selectedParameters 覆盖。
|
|
|
// trace CSV 和 meta 使用该完整向量记录一次候选评价。
|
|
|
QVector<double> 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<nmDataVerticalFracturedWell*>(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<nmDataHorizontalFracturedWell*>(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<double> 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<nmDataWellBase*> 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<nmDataVerticalFracturedWell*>(pTargetWell);
|
|
|
if(fracturedWell) {
|
|
|
initialValue = fracturedWell->getDfc().getValue().toDouble();
|
|
|
}
|
|
|
} else if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
|
|
|
nmDataHorizontalFracturedWell* fracturedWell =
|
|
|
dynamic_cast<nmDataHorizontalFracturedWell*>(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<nmDataVerticalFracturedWell*>(pTargetWell);
|
|
|
if(fracturedWell) {
|
|
|
initialValue = fracturedWell->getFractureHalfLength().getValue().toDouble();
|
|
|
}
|
|
|
} else if(pTargetWell && pTargetWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
|
|
|
nmDataHorizontalFracturedWell* fracturedWell =
|
|
|
dynamic_cast<nmDataHorizontalFracturedWell*>(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<double>& parameters,
|
|
|
double* fitness,
|
|
|
AutoFitObjectiveBreakdownLM* breakdown,
|
|
|
QVector<QVector<double> >* 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<QVector<double> > jacobian;
|
|
|
QVector<double> verticalGradient(dimensions, 0.0);
|
|
|
QVector<double> horizontalGradient(dimensions, 0.0);
|
|
|
QVector<double> shapeGradient(dimensions, 0.0);
|
|
|
QVector<bool> jacobianColumnValid(dimensions, false);
|
|
|
|
|
|
// 参数向量的顺序始终与 m_enabledParamIndices 一致，不能按完整参数索引
|
|
|
// 直接访问；下面两个转换函数集中维护这层映射关系。
|
|
|
auto coordinatesFromParameters = [&](const QVector<double>& parameters)
|
|
|
-> QVector<double> {
|
|
|
QVector<double> 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<double>& coordinates)
|
|
|
-> QVector<double> {
|
|
|
QVector<double> 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<QVector<double> >(
|
|
|
residualCount, QVector<double>(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<double> 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<double>* 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<double> 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<int> selectedColumns;
|
|
|
QVector<bool> 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<QVector<double> > normalMatrix(
|
|
|
selectedCount, QVector<double>(selectedCount, 0.0));
|
|
|
QVector<double> 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<double> selectedStep;
|
|
|
bool solved = solveTrustRegionLinearSystem(
|
|
|
normalMatrix, rightHandSide, &selectedStep);
|
|
|
QVector<double> 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<double> 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<double> 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<double>& 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<QVector<double>> 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<QVector<double>> 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<double>& parameters)
|
|
|
{
|
|
|
// 将粒子的“启用参数向量”写回项目数据。
|
|
|
// parameters 的维度必须等于用户勾选的参数数量，顺序由 m_enabledParamIndices 决定。
|
|
|
// 这里不直接跑求解器，只负责把 DataManager 调整到该粒子对应的模型状态。
|
|
|
if(parameters.size() != getEnabledParameterCount()) {
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
updateReservoirParameters(parameters);
|
|
|
updateWellParameters(parameters);
|
|
|
}
|
|
|
|
|
|
void nmCalculationAutoFitLM::updateReservoirParameters(const QVector<double>& 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<double>& parameters)
|
|
|
{
|
|
|
// 更新目标井上的拟合参数。目前井级可拟合参数主要是：
|
|
|
// - skin：写入第一个 perforation；
|
|
|
// - wellboreC：写入井筒储集系数；
|
|
|
// - Dfc/裂缝半长：只写入垂直压裂井或多段压裂水平井。
|
|
|
// 如果目标井不存在或没有射孔数据，这里只记录 debug，不抛异常。
|
|
|
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
|
|
|
|
|
|
// 只更新当前目标井的井参数，避免多井项目中误改其他井。
|
|
|
//QVector<nmDataWellBase*> 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<nmDataVerticalFracturedWell*>(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<nmDataHorizontalFracturedWell*>(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<nmDataVerticalFracturedWell*>(pWell);
|
|
|
if(fracturedWell) {
|
|
|
fracturedWell->getFractureHalfLength().setValue(value);
|
|
|
}
|
|
|
} else if(pWell->getWellType() == NM_WELL_MODEL::Horizontal_Fractured_Well) {
|
|
|
nmDataHorizontalFracturedWell* fracturedWell =
|
|
|
dynamic_cast<nmDataHorizontalFracturedWell*>(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<nmDataVerticalWell*>(pWell);
|
|
|
|
|
|
if(pVerticalWell) {
|
|
|
QVector<nmDataVerticalWell> wells;
|
|
|
wells.append(*pVerticalWell);
|
|
|
dataManager->updateVerticalWells(wells);
|
|
|
}
|
|
|
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
case NM_WELL_MODEL::Vertical_Fractured_Well: {
|
|
|
nmDataVerticalFracturedWell* pVFracturedWell = dynamic_cast<nmDataVerticalFracturedWell*>(pWell);
|
|
|
|
|
|
if(pVFracturedWell) {
|
|
|
QVector<nmDataVerticalFracturedWell> wells;
|
|
|
wells.append(*pVFracturedWell);
|
|
|
dataManager->updateVerticalFracturedWells(wells);
|
|
|
}
|
|
|
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
case NM_WELL_MODEL::Horizontal_Fractured_Well: {
|
|
|
nmDataHorizontalFracturedWell* pHFracturedWell = dynamic_cast<nmDataHorizontalFracturedWell*>(pWell);
|
|
|
|
|
|
if(pHFracturedWell) {
|
|
|
QVector<nmDataHorizontalFracturedWell> wells;
|
|
|
wells.append(*pHFracturedWell);
|
|
|
dataManager->updateHorizontalFracturedWells(wells);
|
|
|
}
|
|
|
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
default:
|
|
|
break;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
// ==================== 求解器相关方法 ====================
|
|
|
|
|
|
QVector<QVector<double>> 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<double>& 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<QVector<double>>& 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<QVector<double>>& 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<QVector<double> >& target,
|
|
|
const QVector<QVector<double> >& result) const
|
|
|
{
|
|
|
// 主目标在目标与模拟曲线的公共时间范围内比较压力和导数残差；上下、左右
|
|
|
// 和形状只负责诊断误差来源和选择参数，避免同一残差在 total 中被重复计算。
|
|
|
// 整个计算过程均位于 log(time)-log(value) 坐标。
|
|
|
const double invalidLoss = 1.0e10;
|
|
|
const double valueFloor = 1.0e-12;
|
|
|
const int numPoints = 80;
|
|
|
m_lastObjectiveBreakdown = AutoFitObjectiveBreakdownLM();
|
|
|
|
|
|
if(!validateLogLogData(target) || !validateLogLogData(result)) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
try {
|
|
|
// 无法比较的采样行先跳过；有限但非正的导数无法进入双对数空间，
|
|
|
// 当前数据又没有逐点有效掩码，因此遇到这种导数时判本次评价无效。
|
|
|
auto prepareCurve = [valueFloor](const QVector<QVector<double> >& data,
|
|
|
int firstIndex,
|
|
|
QVector<QPointF>* pressure,
|
|
|
QVector<QPointF>* 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<QPointF>* curve) {
|
|
|
std::stable_sort(
|
|
|
curve->begin(), curve->end(),
|
|
|
[](const QPointF& left, const QPointF& right) {
|
|
|
return left.x() < right.x();
|
|
|
});
|
|
|
|
|
|
QVector<QPointF> 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<QPointF> targetPressure;
|
|
|
QVector<QPointF> targetDerivative;
|
|
|
QVector<QPointF> resultPressure;
|
|
|
QVector<QPointF> resultDerivative;
|
|
|
// 模拟结果已跳过 DLL 首点，因此误差计算同步忽略目标曲线首点。
|
|
|
if(!prepareCurve(target, 1, &targetPressure, &targetDerivative) ||
|
|
|
!prepareCurve(result, 0, &resultPressure, &resultDerivative)) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
// 在双对数坐标中插值。二分定位用于后面的多次水平配准试算。
|
|
|
auto interpolateLogValue = [valueFloor](
|
|
|
const QVector<QPointF>& 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);
|
|
|
};
|
|
|
|
|
|
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;
|
|
|
}
|
|
|
|
|
|
// 与 PSO 保持一致：只在目标与模拟曲线的时间交集内比较，不再设置
|
|
|
// 覆盖率门槛，也不对交集之外的首尾数据做外推。
|
|
|
const double overlapMinX = qMax(targetMinX, resultMinX);
|
|
|
const double overlapMaxX = qMin(targetMaxX, resultMaxX);
|
|
|
if(overlapMinX >= overlapMaxX) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
QVector<double> commonX(numPoints);
|
|
|
QVector<double> commonLogX(numPoints);
|
|
|
QVector<double> targetLogPressure(numPoints);
|
|
|
QVector<double> targetLogDerivative(numPoints);
|
|
|
const double comparisonLogMinX = qLn(overlapMinX);
|
|
|
const double comparisonLogMaxX = qLn(overlapMaxX);
|
|
|
|
|
|
// 在公共时间范围内生成固定维度的 log-time 网格，保持 LM 残差向量为 160 维。
|
|
|
for(int i = 0; i < numPoints; ++i) {
|
|
|
double logX = comparisonLogMinX +
|
|
|
static_cast<double>(i) *
|
|
|
(comparisonLogMaxX - comparisonLogMinX) /
|
|
|
(numPoints - 1);
|
|
|
commonLogX[i] = logX;
|
|
|
// 首尾直接使用原始端点，避免 exp(log(t)) 的舍入误差越过严格插值边界。
|
|
|
if(i == 0) {
|
|
|
commonX[i] = overlapMinX;
|
|
|
} else if(i == numPoints - 1) {
|
|
|
commonX[i] = overlapMaxX;
|
|
|
} else {
|
|
|
commonX[i] = qExp(logX);
|
|
|
}
|
|
|
|
|
|
if(!interpolateLogValue(
|
|
|
targetPressure, commonX[i],
|
|
|
&targetLogPressure[i]) ||
|
|
|
!interpolateLogValue(
|
|
|
targetDerivative, commonX[i],
|
|
|
&targetLogDerivative[i])) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
QVector<double> pressureResidual(
|
|
|
numPoints, std::numeric_limits<double>::quiet_NaN());
|
|
|
QVector<double> derivativeResidual(
|
|
|
numPoints, std::numeric_limits<double>::quiet_NaN());
|
|
|
|
|
|
// 残差定义为“模拟减目标”：正值表示模拟曲线偏高，负值表示偏低。
|
|
|
for(int i = 0; i < numPoints; ++i) {
|
|
|
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];
|
|
|
}
|
|
|
|
|
|
AutoFitObjectiveBreakdownLM breakdown;
|
|
|
|
|
|
// 在指定中心附近计算普通均方根误差。
|
|
|
auto rmseAround = [](
|
|
|
const QVector<double>& 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<int>(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<double>::quiet_NaN();
|
|
|
};
|
|
|
|
|
|
auto rmse = [&rmseAround](
|
|
|
const QVector<double>& values,
|
|
|
int begin,
|
|
|
int end) -> double {
|
|
|
return rmseAround(values, begin, end, 0.0);
|
|
|
};
|
|
|
|
|
|
// 普通算术平均中心保留上下偏差的符号。
|
|
|
auto meanCenterRange = [](
|
|
|
const QVector<double>& values,
|
|
|
int begin,
|
|
|
int end) -> double {
|
|
|
int validBegin = qMax(0, begin);
|
|
|
int validEnd =
|
|
|
qMin(end, static_cast<int>(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<double>::quiet_NaN();
|
|
|
}
|
|
|
return center / count;
|
|
|
};
|
|
|
|
|
|
// 压力和导数合并后只求一个公共中心，表示两条曲线共同的上下位移。
|
|
|
// 分别去中心会把压力与导数之间真实的相对形状差异一并消除。
|
|
|
auto commonMeanCenterRange = [&meanCenterRange](
|
|
|
const QVector<double>& pressureValues,
|
|
|
const QVector<double>& derivativeValues,
|
|
|
int begin,
|
|
|
int end) -> double {
|
|
|
QVector<double> combined;
|
|
|
int validBegin = qMax(0, begin);
|
|
|
int validEnd = qMin(
|
|
|
end,
|
|
|
qMin(static_cast<int>(pressureValues.size()),
|
|
|
static_cast<int>(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<double>& pressureValues,
|
|
|
const QVector<double>& 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<double>::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 =
|
|
|
(comparisonLogMaxX - comparisonLogMinX) /
|
|
|
(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<double>* shiftedPressure,
|
|
|
QVector<double>* shiftedDerivative) -> bool {
|
|
|
if(!shiftedPressure || !shiftedDerivative) {
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
shiftedPressure->fill(
|
|
|
std::numeric_limits<double>::quiet_NaN(),
|
|
|
numPoints);
|
|
|
shiftedDerivative->fill(
|
|
|
std::numeric_limits<double>::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<double>::quiet_NaN();
|
|
|
double bestCenteredLoss =
|
|
|
std::numeric_limits<double>::infinity();
|
|
|
double bestPhysicalShift = 0.0;
|
|
|
int bestShiftStep = 0;
|
|
|
double bestPressureLoss =
|
|
|
std::numeric_limits<double>::infinity();
|
|
|
double bestPressureShift = 0.0;
|
|
|
double bestDerivativeLoss =
|
|
|
std::numeric_limits<double>::infinity();
|
|
|
double bestDerivativeShift = 0.0;
|
|
|
QVector<double> profileLosses(
|
|
|
2 * halfShiftStepCount + 1,
|
|
|
std::numeric_limits<double>::quiet_NaN());
|
|
|
QVector<double> profileCommonBiases(
|
|
|
2 * halfShiftStepCount + 1,
|
|
|
std::numeric_limits<double>::quiet_NaN());
|
|
|
QVector<double> shiftedPressureResidual;
|
|
|
QVector<double> shiftedDerivativeResidual;
|
|
|
|
|
|
// 位移和公共上下偏移联合求解，避免“先扣上下还是先扣左右”的顺序依赖。
|
|
|
for(int shiftStep = -halfShiftStepCount;
|
|
|
shiftStep <= halfShiftStepCount;
|
|
|
++shiftStep) {
|
|
|
double physicalShift =
|
|
|
shiftStep * physicalShiftStep;
|
|
|
if(!buildShiftResidual(
|
|
|
physicalShift,
|
|
|
registrationBegin,
|
|
|
registrationEnd,
|
|
|
&shiftedPressureResidual,
|
|
|
&shiftedDerivativeResidual)) {
|
|
|
continue;
|
|
|
}
|
|
|
|
|
|
double commonBias = commonMeanCenterRange(
|
|
|
shiftedPressureResidual,
|
|
|
shiftedDerivativeResidual,
|
|
|
registrationBegin,
|
|
|
registrationEnd);
|
|
|
double centeredLoss = jointRmseAround(
|
|
|
shiftedPressureResidual,
|
|
|
shiftedDerivativeResidual,
|
|
|
registrationBegin,
|
|
|
registrationEnd,
|
|
|
commonBias);
|
|
|
double pressureBias = meanCenterRange(
|
|
|
shiftedPressureResidual,
|
|
|
registrationBegin,
|
|
|
registrationEnd);
|
|
|
double derivativeBias = meanCenterRange(
|
|
|
shiftedDerivativeResidual,
|
|
|
registrationBegin,
|
|
|
registrationEnd);
|
|
|
double pressureLoss = rmseAround(
|
|
|
shiftedPressureResidual,
|
|
|
registrationBegin,
|
|
|
registrationEnd,
|
|
|
pressureBias);
|
|
|
double derivativeLoss = rmseAround(
|
|
|
shiftedDerivativeResidual,
|
|
|
registrationBegin,
|
|
|
registrationEnd,
|
|
|
derivativeBias);
|
|
|
|
|
|
if(!isFiniteNumber(centeredLoss) ||
|
|
|
!isFiniteNumber(pressureLoss) ||
|
|
|
!isFiniteNumber(derivativeLoss)) {
|
|
|
continue;
|
|
|
}
|
|
|
int profileIndex = shiftStep + halfShiftStepCount;
|
|
|
profileLosses[profileIndex] = centeredLoss;
|
|
|
profileCommonBiases[profileIndex] = commonBias;
|
|
|
if(shiftStep == 0) {
|
|
|
zeroShiftCenteredLoss = centeredLoss;
|
|
|
}
|
|
|
if(isBetterProfileValue(
|
|
|
centeredLoss,
|
|
|
physicalShift,
|
|
|
bestCenteredLoss,
|
|
|
bestPhysicalShift)) {
|
|
|
bestCenteredLoss = centeredLoss;
|
|
|
bestPhysicalShift = physicalShift;
|
|
|
bestShiftStep = shiftStep;
|
|
|
}
|
|
|
if(isBetterProfileValue(
|
|
|
pressureLoss,
|
|
|
physicalShift,
|
|
|
bestPressureLoss,
|
|
|
bestPressureShift)) {
|
|
|
bestPressureLoss = pressureLoss;
|
|
|
bestPressureShift = physicalShift;
|
|
|
}
|
|
|
if(isBetterProfileValue(
|
|
|
derivativeLoss,
|
|
|
physicalShift,
|
|
|
bestDerivativeLoss,
|
|
|
bestDerivativeShift)) {
|
|
|
bestDerivativeLoss = derivativeLoss;
|
|
|
bestDerivativeShift = physicalShift;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
if(!isFiniteNumber(zeroShiftCenteredLoss) ||
|
|
|
!isFiniteNumber(bestCenteredLoss) ||
|
|
|
!isFiniteNumber(bestPressureLoss) ||
|
|
|
!isFiniteNumber(bestDerivativeLoss)) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
// horizontalGain 是“允许水平位移”相对“固定零位移”减少的均方能量。
|
|
|
// 只有改善足够明显且最优点不是边界，才把位移解释为可靠左右偏差。
|
|
|
double horizontalGain = nestedRmsContribution(
|
|
|
zeroShiftCenteredLoss, bestCenteredLoss);
|
|
|
double horizontalSignalThreshold =
|
|
|
qMax(1.0e-5, zeroShiftCenteredLoss * 0.02);
|
|
|
int bestProfileIndex = bestShiftStep + halfShiftStepCount;
|
|
|
double nearbyProfileLoss =
|
|
|
std::numeric_limits<double>::infinity();
|
|
|
int leftProfileIndex =
|
|
|
bestProfileIndex - shiftSubdivisions;
|
|
|
int rightProfileIndex =
|
|
|
bestProfileIndex + shiftSubdivisions;
|
|
|
if(leftProfileIndex >= 0 &&
|
|
|
leftProfileIndex < profileLosses.size() &&
|
|
|
isFiniteNumber(profileLosses[leftProfileIndex])) {
|
|
|
nearbyProfileLoss = qMin(
|
|
|
nearbyProfileLoss,
|
|
|
profileLosses[leftProfileIndex]);
|
|
|
}
|
|
|
if(rightProfileIndex >= 0 &&
|
|
|
rightProfileIndex < profileLosses.size() &&
|
|
|
isFiniteNumber(profileLosses[rightProfileIndex])) {
|
|
|
nearbyProfileLoss = qMin(
|
|
|
nearbyProfileLoss,
|
|
|
profileLosses[rightProfileIndex]);
|
|
|
}
|
|
|
double profileContrast = isFiniteNumber(nearbyProfileLoss)
|
|
|
? nestedRmsContribution(
|
|
|
nearbyProfileLoss,
|
|
|
bestCenteredLoss)
|
|
|
: 0.0;
|
|
|
bool flatRegistrationProfile =
|
|
|
profileContrast <= horizontalSignalThreshold;
|
|
|
|
|
|
// 平台曲线的 profile 也可能很平，但公共 bias 在各个位移下保持稳定，
|
|
|
// 此时仍能可靠判断上下。只有近优位移会明显改变 bias 才说明上下/左右不可辨识。
|
|
|
double minimumNearOptimalBias =
|
|
|
std::numeric_limits<double>::infinity();
|
|
|
double maximumNearOptimalBias =
|
|
|
-std::numeric_limits<double>::infinity();
|
|
|
for(int i = 0; i < profileLosses.size(); ++i) {
|
|
|
if(isFiniteNumber(profileLosses[i]) &&
|
|
|
isFiniteNumber(profileCommonBiases[i]) &&
|
|
|
profileLosses[i] <=
|
|
|
bestCenteredLoss + horizontalSignalThreshold) {
|
|
|
minimumNearOptimalBias = qMin(
|
|
|
minimumNearOptimalBias,
|
|
|
profileCommonBiases[i]);
|
|
|
maximumNearOptimalBias = qMax(
|
|
|
maximumNearOptimalBias,
|
|
|
profileCommonBiases[i]);
|
|
|
}
|
|
|
}
|
|
|
double nearOptimalBiasSpread =
|
|
|
isFiniteNumber(minimumNearOptimalBias) &&
|
|
|
isFiniteNumber(maximumNearOptimalBias)
|
|
|
? maximumNearOptimalBias - minimumNearOptimalBias
|
|
|
: std::numeric_limits<double>::infinity();
|
|
|
bool commonBiasStable = nearOptimalBiasSpread <= 1.0e-2;
|
|
|
bool horizontalAtBoundary =
|
|
|
halfShiftStepCount > 0 &&
|
|
|
qAbs(bestShiftStep) == halfShiftStepCount;
|
|
|
// 压力和导数通道分别求出的最佳位移若方向相反或相差过大，说明一个
|
|
|
// 单一水平平移无法解释两条曲线，此时标记配准歧义并禁用左右引导。
|
|
|
bool pressureShiftDetected =
|
|
|
qAbs(bestPressureShift) >=
|
|
|
0.5 * physicalShiftStep;
|
|
|
bool derivativeShiftDetected =
|
|
|
qAbs(bestDerivativeShift) >=
|
|
|
0.5 * physicalShiftStep;
|
|
|
bool channelShiftConflict =
|
|
|
pressureShiftDetected &&
|
|
|
derivativeShiftDetected &&
|
|
|
(bestPressureShift * bestDerivativeShift < 0.0 ||
|
|
|
qAbs(bestPressureShift - bestDerivativeShift) >
|
|
|
2.0 * logGridStep);
|
|
|
|
|
|
breakdown.horizontalLoss = horizontalGain;
|
|
|
breakdown.horizontalReliable =
|
|
|
maximumShiftIntervals > 0 &&
|
|
|
!horizontalAtBoundary &&
|
|
|
!channelShiftConflict &&
|
|
|
!flatRegistrationProfile &&
|
|
|
horizontalGain > horizontalSignalThreshold &&
|
|
|
qAbs(bestPhysicalShift) >=
|
|
|
0.5 * physicalShiftStep;
|
|
|
breakdown.registrationAmbiguous =
|
|
|
channelShiftConflict ||
|
|
|
(qAbs(bestPhysicalShift) >=
|
|
|
0.5 * physicalShiftStep &&
|
|
|
!breakdown.horizontalReliable) ||
|
|
|
(flatRegistrationProfile && !commonBiasStable);
|
|
|
breakdown.horizontalPhysicalShift =
|
|
|
breakdown.horizontalReliable
|
|
|
? bestPhysicalShift
|
|
|
: 0.0;
|
|
|
|
|
|
// 可信水平位移确定后，在该位移实际覆盖的最大区间重新计算上下和形状。
|
|
|
int diagnosticBegin = 0;
|
|
|
int diagnosticEnd = numPoints;
|
|
|
while(diagnosticBegin < diagnosticEnd &&
|
|
|
commonLogX[diagnosticBegin] +
|
|
|
breakdown.horizontalPhysicalShift <
|
|
|
resultLogMinX - 1.0e-12) {
|
|
|
++diagnosticBegin;
|
|
|
}
|
|
|
while(diagnosticEnd > diagnosticBegin &&
|
|
|
commonLogX[diagnosticEnd - 1] +
|
|
|
breakdown.horizontalPhysicalShift >
|
|
|
resultLogMaxX + 1.0e-12) {
|
|
|
--diagnosticEnd;
|
|
|
}
|
|
|
if(diagnosticEnd - diagnosticBegin <
|
|
|
minimumRegistrationPoints ||
|
|
|
!buildShiftResidual(
|
|
|
breakdown.horizontalPhysicalShift,
|
|
|
diagnosticBegin,
|
|
|
diagnosticEnd,
|
|
|
&shiftedPressureResidual,
|
|
|
&shiftedDerivativeResidual)) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
double commonBias = commonMeanCenterRange(
|
|
|
shiftedPressureResidual,
|
|
|
shiftedDerivativeResidual,
|
|
|
diagnosticBegin,
|
|
|
diagnosticEnd);
|
|
|
double rawAlignedLoss = jointRmseAround(
|
|
|
shiftedPressureResidual,
|
|
|
shiftedDerivativeResidual,
|
|
|
diagnosticBegin,
|
|
|
diagnosticEnd,
|
|
|
0.0);
|
|
|
double centeredAlignedLoss = jointRmseAround(
|
|
|
shiftedPressureResidual,
|
|
|
shiftedDerivativeResidual,
|
|
|
diagnosticBegin,
|
|
|
diagnosticEnd,
|
|
|
commonBias);
|
|
|
if(!isFiniteNumber(commonBias) ||
|
|
|
!isFiniteNumber(rawAlignedLoss) ||
|
|
|
!isFiniteNumber(centeredAlignedLoss)) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
// 原始对齐误差减去公共中心后的能量差定义为上下误差贡献。只有它相对
|
|
|
// 当前对齐误差足够明显，且配准无歧义时，公共 bias 才可用于有符号选参。
|
|
|
breakdown.verticalCommonBias = commonBias;
|
|
|
breakdown.verticalLoss = nestedRmsContribution(
|
|
|
rawAlignedLoss, centeredAlignedLoss);
|
|
|
breakdown.verticalReliable =
|
|
|
!breakdown.registrationAmbiguous &&
|
|
|
breakdown.verticalLoss >
|
|
|
qMax(1.0e-5, rawAlignedLoss * 0.02);
|
|
|
|
|
|
QVector<double> shapePressure(
|
|
|
numPoints, std::numeric_limits<double>::quiet_NaN());
|
|
|
QVector<double> shapeDerivative(
|
|
|
numPoints, std::numeric_limits<double>::quiet_NaN());
|
|
|
for(int i = diagnosticBegin; i < diagnosticEnd; ++i) {
|
|
|
if(isFiniteNumber(shiftedPressureResidual[i])) {
|
|
|
shapePressure[i] =
|
|
|
shiftedPressureResidual[i] - commonBias;
|
|
|
}
|
|
|
if(isFiniteNumber(shiftedDerivativeResidual[i])) {
|
|
|
shapeDerivative[i] =
|
|
|
shiftedDerivativeResidual[i] - commonBias;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
// shapeLoss 是去除可信左右位移和公共均值中心后的剩余误差。
|
|
|
double shapePressureLoss = rmse(
|
|
|
shapePressure, diagnosticBegin, diagnosticEnd);
|
|
|
double shapeDerivativeLoss = rmse(
|
|
|
shapeDerivative, diagnosticBegin, diagnosticEnd);
|
|
|
if(!isFiniteNumber(shapePressureLoss) ||
|
|
|
!isFiniteNumber(shapeDerivativeLoss)) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
breakdown.shapeLoss = qSqrt(
|
|
|
0.5 *
|
|
|
(shapePressureLoss * shapePressureLoss +
|
|
|
shapeDerivativeLoss * shapeDerivativeLoss));
|
|
|
|
|
|
// 现阶段不识别或单独调度晚期流动段；保留字段只为了维持现有 trace 列。
|
|
|
breakdown.lateDerivativeSlopeBias = 0.0;
|
|
|
breakdown.lateDerivativeTrendLoss = 0.0;
|
|
|
breakdown.lateDerivativeTrendReliable = false;
|
|
|
|
|
|
if(!isFiniteNumber(breakdown.pressureLoss) ||
|
|
|
!isFiniteNumber(breakdown.derivativeLoss) ||
|
|
|
!isFiniteNumber(breakdown.verticalCommonBias) ||
|
|
|
!isFiniteNumber(breakdown.verticalLoss) ||
|
|
|
!isFiniteNumber(breakdown.horizontalPhysicalShift) ||
|
|
|
!isFiniteNumber(breakdown.horizontalLoss) ||
|
|
|
!isFiniteNumber(breakdown.shapeLoss) ||
|
|
|
breakdown.residualVector.size() != 2 * numPoints) {
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
|
|
|
// LM 总目标等于固定残差向量的二范数；压力和导数各占一半能量。
|
|
|
// 上下、左右和形状分量不参与候选排序与接受。
|
|
|
breakdown.total = qSqrt(
|
|
|
0.5 * breakdown.pressureLoss * breakdown.pressureLoss +
|
|
|
0.5 * breakdown.derivativeLoss * breakdown.derivativeLoss);
|
|
|
breakdown.valid =
|
|
|
isFiniteNumber(breakdown.total) &&
|
|
|
breakdown.total >= 0.0;
|
|
|
m_lastObjectiveBreakdown = breakdown;
|
|
|
|
|
|
DEBUG_OUT(
|
|
|
QString("LogLog objective: pressure=%1, derivative=%2, vertical=%3, horizontal=%4, shape=%5, ambiguous=%6, shift=%7, total=%8")
|
|
|
.arg(breakdown.pressureLoss, 0, 'e', 4)
|
|
|
.arg(breakdown.derivativeLoss, 0, 'e', 4)
|
|
|
.arg(breakdown.verticalLoss, 0, 'e', 4)
|
|
|
.arg(breakdown.horizontalLoss, 0, 'e', 4)
|
|
|
.arg(breakdown.shapeLoss, 0, 'e', 4)
|
|
|
.arg(breakdown.registrationAmbiguous)
|
|
|
.arg(breakdown.horizontalPhysicalShift, 0, 'e', 4)
|
|
|
.arg(breakdown.total, 0, 'e', 4));
|
|
|
|
|
|
return breakdown.valid
|
|
|
? qMin(1.0e9, breakdown.total)
|
|
|
: invalidLoss;
|
|
|
} catch(const std::exception& e) {
|
|
|
DEBUG_OUT(
|
|
|
QString("Exception in LogLog error calculation: %1")
|
|
|
.arg(e.what()));
|
|
|
return invalidLoss;
|
|
|
} catch(...) {
|
|
|
DEBUG_OUT("Unknown exception in LogLog error calculation");
|
|
|
return invalidLoss;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
QVector<QVector<double>> nmCalculationAutoFitLM::runSolverDll()
|
|
|
{
|
|
|
// DLL 求解器路径。
|
|
|
// 这个函数负责把当前 DataManager 中的项目状态交给底层数值求解器，
|
|
|
// 数值求解仍包含全部计算井，以保留井间干扰；后处理只提取目标井曲线。
|
|
|
//
|
|
|
// 如果这里失败，通常需要优先检查：HX_NWTM.dll、license、网格/井数据是否完整、
|
|
|
// 目标井是否存在，以及 DataManager 中刚写入的参数是否导致求解器异常。
|
|
|
DEBUG_OUT("SOLVER DLL START");
|
|
|
|
|
|
// 创建任务时绑定当前分析，并在当前线程冻结本次求解所需的全部输入。
|
|
|
nmDataAnalyzeManager* dataManager =
|
|
|
nmDataAnalyzeManager::getCurrentInstance();
|
|
|
if(!dataManager || m_targetWellName.isEmpty()) {
|
|
|
DEBUG_OUT("Data manager or target well is unavailable");
|
|
|
return QVector<QVector<double> >();
|
|
|
}
|
|
|
|
|
|
if(m_evaluationInProgress > 0) {
|
|
|
DEBUG_OUT("DLL Solver already running, skipping");
|
|
|
return QVector<QVector<double>>();
|
|
|
}
|
|
|
|
|
|
++m_evaluationInProgress;
|
|
|
m_lastEvaluatedLogLogData.clear();
|
|
|
QVector<QVector<double>> result;
|
|
|
nmCalculationDllPebiSolverTask* dllTask = nullptr;
|
|
|
|
|
|
try {
|
|
|
DEBUG_OUT("Creating DLL solver task");
|
|
|
dllTask = new nmCalculationDllPebiSolverTask(
|
|
|
m_tempDirectory,
|
|
|
dataManager,
|
|
|
m_targetWellName);
|
|
|
|
|
|
if(m_shouldStop) {
|
|
|
DEBUG_OUT("Should stop - cleaning up and returning empty result");
|
|
|
delete dllTask;
|
|
|
--m_evaluationInProgress;
|
|
|
return result;
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT("Starting DLL solver execution...");
|
|
|
|
|
|
// 异步执行 DLL 任务。循环等待期间持续 processEvents，保证界面不会完全卡死。
|
|
|
dllTask->start();
|
|
|
|
|
|
// 等待完成。最大等待 1 小时，适配大模型慢算；用户停止时会 terminate。
|
|
|
const int maxWait = 3600000; // 1h超时
|
|
|
const int checkInterval = 50;
|
|
|
QTime waitTimer;
|
|
|
waitTimer.start();
|
|
|
|
|
|
while(waitTimer.elapsed() < maxWait) {
|
|
|
// wait(timeout) 会在线程一完成时立即返回，避免原来固定 msleep(500)
|
|
|
// 带来的每次 0～500ms 额外等待；50ms 间隔仍可及时处理停止请求和界面事件。
|
|
|
if(dllTask->wait(checkInterval)) {
|
|
|
DEBUG_OUT("DLL solver task completed");
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
QApplication::processEvents(QEventLoop::ExcludeUserInputEvents, checkInterval);
|
|
|
|
|
|
if(m_shouldStop) {
|
|
|
DEBUG_OUT("DLL solver task terminated by user");
|
|
|
dllTask->terminate();
|
|
|
break;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
// 超时处理
|
|
|
if(dllTask->isRunning()) {
|
|
|
DEBUG_OUT("DLL solver task timeout, terminating...");
|
|
|
dllTask->terminate();
|
|
|
dllTask->wait(2000);
|
|
|
|
|
|
delete dllTask;
|
|
|
dllTask = nullptr;
|
|
|
--m_evaluationInProgress;
|
|
|
m_consecutiveFailures++;
|
|
|
return result;
|
|
|
}
|
|
|
|
|
|
// 线程结束后检查真实执行结果，防止失败时复用上一粒子的旧曲线。
|
|
|
dllTask->wait();
|
|
|
if(!dllTask->wasSuccessful()) {
|
|
|
DEBUG_OUT("DLL solver task reported failure");
|
|
|
delete dllTask;
|
|
|
dllTask = nullptr;
|
|
|
--m_evaluationInProgress;
|
|
|
m_consecutiveFailures++;
|
|
|
return result;
|
|
|
}
|
|
|
|
|
|
// 任务结束后复制其局部结果，删除任务前不再持有任务内部引用。
|
|
|
QVector<QVector<double>> pressureResult = dllTask->getAutoFitResultPressure();
|
|
|
QVector<QVector<double>> logLogResult = dllTask->getAutoFitResultLogLog();
|
|
|
|
|
|
DEBUG_OUT(QString("DLL result verification - Pressure arrays: %1, LogLog arrays: %2")
|
|
|
.arg(pressureResult.size()).arg(logLogResult.size()));
|
|
|
|
|
|
if(pressureResult.size() >= 2) {
|
|
|
DEBUG_OUT(QString("Pressure result - Time points: %1, Pressure points: %2")
|
|
|
.arg(pressureResult[0].size()).arg(pressureResult[1].size()));
|
|
|
|
|
|
if(pressureResult[0].size() > 0) {
|
|
|
DEBUG_OUT(QString("Sample pressure data - Time[0]: %1, Time[last]: %2, P[0]: %3, P[last]: %4")
|
|
|
.arg(pressureResult[0][0])
|
|
|
.arg(pressureResult[0][pressureResult[0].size() - 1])
|
|
|
.arg(pressureResult[1][0])
|
|
|
.arg(pressureResult[1][pressureResult[1].size() - 1]));
|
|
|
}
|
|
|
}
|
|
|
|
|
|
// 数据有效性检查
|
|
|
if(pressureResult.size() >= 2
|
|
|
&& pressureResult[0].size() > 0
|
|
|
&& pressureResult[1].size() > 0
|
|
|
&& validateLogLogData(logLogResult)) {
|
|
|
result = pressureResult;
|
|
|
m_lastEvaluatedLogLogData = logLogResult;
|
|
|
DEBUG_OUT(QString("Got DLL solver result: %1 points").arg(result[0].size()));
|
|
|
m_consecutiveFailures = 0;
|
|
|
|
|
|
// 检查结果是否与之前不同。若连续粒子得到完全相同的压力曲线，
|
|
|
// 可能说明参数没有正确写入 DataManager，或求解器缓存/状态没有刷新。
|
|
|
static QVector<double> lastPressureResult;
|
|
|
bool isDifferentFromLast = false;
|
|
|
|
|
|
if(lastPressureResult.isEmpty() || lastPressureResult.size() != pressureResult[1].size()) {
|
|
|
isDifferentFromLast = true;
|
|
|
} else {
|
|
|
for(int i = 0; i < qMin(5, pressureResult[1].size()); ++i) {
|
|
|
if(qAbs(lastPressureResult[i] - pressureResult[1][i]) > 1e-12) {
|
|
|
isDifferentFromLast = true;
|
|
|
break;
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
|
|
|
if(isDifferentFromLast) {
|
|
|
DEBUG_OUT("RESULT VERIFICATION: Got NEW result data from DLL");
|
|
|
lastPressureResult = pressureResult[1];
|
|
|
} else {
|
|
|
DEBUG_OUT("!!! WARNING: Result data appears to be identical to previous run !!!");
|
|
|
}
|
|
|
|
|
|
} else {
|
|
|
DEBUG_OUT("DLL solver result is empty or invalid");
|
|
|
DEBUG_OUT(QString("Pressure result size: %1, Array sizes: %2, %3")
|
|
|
.arg(pressureResult.size())
|
|
|
.arg(pressureResult.size() > 0 ? pressureResult[0].size() : 0)
|
|
|
.arg(pressureResult.size() > 1 ? pressureResult[1].size() : 0));
|
|
|
m_consecutiveFailures++;
|
|
|
}
|
|
|
|
|
|
} catch(const std::bad_alloc& e) {
|
|
|
DEBUG_OUT(QString("Memory allocation failed in DLL solver: %1").arg(e.what()));
|
|
|
m_consecutiveFailures++;
|
|
|
} catch(const std::exception& e) {
|
|
|
DEBUG_OUT(QString("Exception in DLL solver: %1").arg(e.what()));
|
|
|
m_consecutiveFailures++;
|
|
|
} catch(...) {
|
|
|
DEBUG_OUT("Unknown exception in DLL solver");
|
|
|
m_consecutiveFailures++;
|
|
|
}
|
|
|
|
|
|
// 清理DLL任务
|
|
|
if(dllTask) {
|
|
|
if(dllTask->isRunning()) {
|
|
|
dllTask->terminate();
|
|
|
dllTask->wait(3000);
|
|
|
}
|
|
|
|
|
|
DEBUG_OUT("Cleaning up DLL solver task...");
|
|
|
delete dllTask;
|
|
|
dllTask = nullptr;
|
|
|
}
|
|
|
|
|
|
QApplication::processEvents(QEventLoop::AllEvents, 100);
|
|
|
--m_evaluationInProgress;
|
|
|
|
|
|
DEBUG_OUT(QString("SOLVER DLL END - ResultPoints: %1")
|
|
|
.arg(result.isEmpty() ? 0 : result[0].size()));
|
|
|
|
|
|
return result;
|
|
|
}
|
|
|
|
|
|
bool nmCalculationAutoFitLM::runFinalFullSolver()
|
|
|
{
|
|
|
// 不设置目标井名，任务按原完整模式保存全部井和网格结果。
|
|
|
nmDataAnalyzeManager* dataManager =
|
|
|
nmDataAnalyzeManager::getCurrentInstance();
|
|
|
if(!dataManager) {
|
|
|
DEBUG_OUT("Cannot start final full-field solver without a data manager");
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
if(m_evaluationInProgress > 0) {
|
|
|
DEBUG_OUT("Cannot start final full-field solver while another evaluation is running");
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
++m_evaluationInProgress;
|
|
|
nmCalculationDllPebiSolverTask dllTask(m_tempDirectory, dataManager);
|
|
|
dllTask.start();
|
|
|
|
|
|
const int maxWait = 3600000;
|
|
|
const int checkInterval = 50;
|
|
|
QTime waitTimer;
|
|
|
waitTimer.start();
|
|
|
|
|
|
while(waitTimer.elapsed() < maxWait) {
|
|
|
if(dllTask.wait(checkInterval)) {
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
// 最终完整计算可能持续较长时间，此处需处理停止按钮事件。
|
|
|
QApplication::processEvents(QEventLoop::AllEvents, checkInterval);
|
|
|
|
|
|
if(!dllTask.isRunning()) {
|
|
|
break;
|
|
|
}
|
|
|
|
|
|
if(m_shouldStop) {
|
|
|
DEBUG_OUT("Final full-field solver terminated by user");
|
|
|
dllTask.terminate();
|
|
|
dllTask.wait(2000);
|
|
|
--m_evaluationInProgress;
|
|
|
return false;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
if(dllTask.isRunning()) {
|
|
|
DEBUG_OUT("Final full-field solver timeout, terminating task");
|
|
|
dllTask.terminate();
|
|
|
dllTask.wait(2000);
|
|
|
--m_evaluationInProgress;
|
|
|
return false;
|
|
|
}
|
|
|
|
|
|
dllTask.wait();
|
|
|
// 后台只生成局部结果快照；确认求解和输入版本均有效后再一次性写回当前分析。
|
|
|
const bool succeeded =
|
|
|
dllTask.wasSuccessful() && dllTask.commitResult(dataManager);
|
|
|
--m_evaluationInProgress;
|
|
|
return succeeded;
|
|
|
}
|
|
|
|
|
|
QString nmCalculationAutoFitLM::getStopReasonDescription(StopReasonLM reason) const
|
|
|
{
|
|
|
// 将停止枚举转成人类可读文本，用于日志和运行摘要。
|
|
|
switch(reason) {
|
|
|
case LM_TARGET_ACHIEVED:
|
|
|
return tr("Target error achieved");
|
|
|
|
|
|
case LM_TRUE_CONVERGENCE:
|
|
|
return tr("Algorithm converged to stable solution");
|
|
|
|
|
|
case LM_LOCAL_OPTIMUM:
|
|
|
return tr("Local optimum detected");
|
|
|
|
|
|
case LM_MAX_ITERATIONS:
|
|
|
return tr("Maximum iterations reached");
|
|
|
|
|
|
case LM_USER_STOPPED:
|
|
|
return tr("Stopped by user request");
|
|
|
|
|
|
case LM_CONSECUTIVE_FAILURES:
|
|
|
return tr("Too many consecutive failures");
|
|
|
|
|
|
case LM_OPTIMIZATION_FAILED:
|
|
|
return tr("Optimization failed");
|
|
|
|
|
|
default:
|
|
|
return tr("Unknown reason");
|
|
|
}
|
|
|
}
|