You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
nmWTAI-Platform/Src/nmNum/nmCalculation/nmCalculationAutoFitPSO.cpp

6005 lines
237 KiB
C++

This file contains ambiguous Unicode characters!

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

#include "nmCalculationAutoFitPSO.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 "nmCalculationPebiGrid.h"
#include <QApplication>
#include <QDebug>
#include <QTime>
#include <QDir>
#include <QTextStream>
#include <QFileInfo>
#include <QProcess>
#include <QProcessEnvironment>
#include <QPair>
#include <QtCore/qmath.h>
#include <cmath>
#include <algorithm>
#include <limits>
#include <QTimer>
#ifdef Q_OS_WIN
#include <windows.h>
#include <float.h>
#define DEBUG_OUT(msg) OutputDebugStringA(QString("[AutoFit] %1\n").arg(msg).toLocal8Bit().data())
#endif
// ===== PSO基础常量 =====
// MIN_FITNESS_IMPROVEMENT判断适应度是否有有效改善的最小阈值。
// VELOCITY_LIMIT_FACTOR单次速度最大值占参数搜索区间的比例防止粒子一步跨太远。
// CONVERGENCE_CHECK_INTERVAL收敛检查的基础间隔部分历史逻辑和日志会使用它。
const double nmCalculationAutoFitPSO::MIN_FITNESS_IMPROVEMENT = 1e-8;
const double nmCalculationAutoFitPSO::VELOCITY_LIMIT_FACTOR = 0.2;
const int nmCalculationAutoFitPSO::CONVERGENCE_CHECK_INTERVAL = 5;
// ===== 代理模型筛选参数 =====
//
// PSO acceleration 的原则是“少跑一部分真实求解器,但不能完全信任代理模型”:
// - keep_fraction代理评分排名靠前的粒子进入真实求解器
// - audit_fraction从被筛掉的粒子中随机抽查一部分避免代理模型系统性漏掉好解
// - min_solver_fraction每代至少有一定比例粒子跑真实求解器保证 pbest/gbest 仍由真实误差支撑;
// - warmup_iterations前几代全量真实评价让每个粒子先拿到可靠的个体最优
// - full_solver_interval周期性全量审计用于纠偏代理模型排序误差。
static const double kSurrogateKeepFraction = 0.50; // 常规工况:代理排序前 50% 粒子进入真实求解器。
static const double kSurrogateHighRiskKeepFraction = 0.70; // 高风险流量制度:保留更多 top 粒子。
static const double kSurrogateAuditFraction = 0.05; // 从被筛掉候选中随机审计 5%。
static const double kSurrogateMinSolverFraction = 0.55; // 常规工况:每代至少 55% 粒子跑真实求解器。
static const double kSurrogateHighRiskMinSolverFraction = 0.80; // 高风险流量制度:每代至少 80% 粒子跑真实求解器。
static const int kSurrogateWarmupIterations = 2; // 前 2 代全量真实评价,给每个粒子建立 pbest。
static const int kSurrogateFullSolverInterval = 10; // 每 10 代做一次全量真实求解审计。
static const int kSurrogateScoringMaxAttempts = 2; // Python 评分失败时最多重试 2 次。
static const bool kUseFixedPsoSeed = true; // 固定随机种子便于复现 trace 和交接演示。
static const unsigned int kFixedPsoSeed = 1792008679u; // 当前固定种子值。
// ===== 代理模型训练域约束 =====
//
// 这些范围必须与 ML/nmWTAI-ML/configs/data_gen_family_random_v2_q.yaml 中的
// 数据生成范围保持一致。超出范围时C++ 侧会强制退回真实求解器,避免代理模型
// 在训练域外给出不可靠排序。
static const double kSurrogateKMin = 1.0e-3; // 代理模型训练域:渗透率下界。
static const double kSurrogateKMax = 10.0; // 代理模型训练域:渗透率上界。
static const double kSurrogateSkinMin = -10.0; // 代理模型训练域skin 下界。
static const double kSurrogateSkinMax = 10.0; // 代理模型训练域skin 上界。
static const double kSurrogateWellboreCMin = 1.0e-4; // 代理模型训练域:井筒储集下界。
static const double kSurrogateWellboreCMax = 2.0; // 代理模型训练域:井筒储集上界。
static const double kSurrogatePhiMin = 1.0e-2; // 代理模型训练域:孔隙度下界。
static const double kSurrogatePhiMax = 0.50; // 代理模型训练域:孔隙度上界。
static const double kSurrogateHMin = 2.0; // 代理模型训练域:储层厚度下界。
static const double kSurrogateHMax = 50.0; // 代理模型训练域:储层厚度上界。
static const double kSurrogateCfFixed = 4.315e-4; // 当前代理模型按近似固定 Cf 训练。
static const double kSurrogateCfRelTol = 0.25; // Cf 允许相对偏差,超过则不使用代理筛选。
static const int kSurrogateMinProdSections = 3; // 代理训练流量制度:生产段数量下界。
static const int kSurrogateMaxProdSections = 5; // 代理训练流量制度:生产段数量上界。
static const double kSurrogateProdTimeMin = 24.0; // 代理训练流量制度:生产总时长下界。
static const double kSurrogateProdTimeMax = 220.0; // 代理训练流量制度:生产总时长上界。
static const double kSurrogateShutinTimeMin = 12.0; // 代理训练流量制度:关井段时长下界。
static const double kSurrogateShutinTimeMax = 140.0; // 代理训练流量制度:关井段时长上界。
static const double kSurrogateQMin = 50.0; // 代理训练流量制度:产量下界。
static const double kSurrogateQMax = 500.0; // 代理训练流量制度:产量上界。
static const double kSurrogateQZeroThreshold = 1.0e-6; // 判断最后一段是否关井的零产量阈值。
static const double kSurrogateMaxRateStepRatio = 2.2; // 相邻生产段产量最大跳变比。
static const double kSurrogateStrongDeclineEndStartRatio = 0.70; // 末段/首段产量小于该值视为强递减。
static const double kSurrogateRateTrendRelTol = 1.0e-9; // 判断产量趋势时的相对容差。
// 判断 double 是否为有限数。PSO、误差函数、trace 写出都依赖它过滤 NaN/Inf。
static inline bool isFiniteNumber(double value)
{
#ifdef Q_OS_WIN
return _finite(value) != 0;
#else
return std::isfinite(value);
#endif
}
static inline bool isInClosedRange(double value, double lower, double upper)
{
return isFiniteNumber(value) && value >= lower && value <= upper;
}
// 判断 value 是否接近 reference。用于 Cf 这类“近似固定”的训练域检查。
static inline bool isNearRelative(double value, double reference, double relTol)
{
if(!isFiniteNumber(value) || !isFiniteNumber(reference)) {
return false;
}
return qAbs(value - reference) <= qMax(1e-12, qAbs(reference) * relTol);
}
// Windows 下的毫秒级睡眠封装。用于重试间隔、UI 事件循环间隔和模拟模式节奏控制。
static inline void msleep(int ms)
{
#ifdef Q_OS_WIN
Sleep(ms);
#endif
}
static QString csvEscape(const QString& text)
{
// CSV 字段转义。trace 中的 run_id、phase、screeningDecision 等字符串都走这里,
// 防止逗号或引号破坏表格结构。
QString escaped = text;
escaped.replace("\"", "\"\"");
return QString("\"%1\"").arg(escaped);
}
static QString traceNumber(double value)
{
// trace CSV 中的数值格式化。不可用数值写空字段,便于后续 pandas/Excel 读取。
if(!isFiniteNumber(value)) {
return QString();
}
return QString::number(value, 'g', 17);
}
static QString traceParamAt(const QVector<double>& params, int index)
{
// 从完整参数向量中安全取值并格式化。越界时写空字段,避免 trace 写出崩溃。
if(index < 0 || index >= params.size()) {
return QString();
}
return traceNumber(params[index]);
}
static QString jsonEscape(const QString& text)
{
// 手写 JSON 字符串转义。trace meta 不依赖第三方 JSON writer
// 所以所有字符串字段写入前必须统一经过这里。
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)
{
// JSON 数值格式化。NaN/Inf 在 JSON 中没有合法表示,因此写 null。
if(!isFiniteNumber(value)) {
return "null";
}
return QString::number(value, 'g', 17);
}
static QString jsonDoubleArray(const QVector<double>& values)
{
// 将 double 数组写成 JSON 数组,用于 trace meta 的目标曲线、参数上下界等。
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)
{
// 将 int 数组写成 JSON 数组,用于 enabled_param_indices 等字段。
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)
{
// 将 bool 数组写成 JSON 数组,用于参数是否启用等字段。
QStringList items;
for(int i = 0; i < values.size(); ++i) {
items << (values[i] ? "true" : "false");
}
return QString("[%1]").arg(items.join(","));
}
static double deterministicAuditRandom01(unsigned int seed, int generation, int particleIndex)
{
// 代理筛选的随机审计需要“随机但可复现”:
// 同一个 PSO 种子、代号、粒子号会得到同一个 0-1 值,方便复盘 trace。
unsigned int x = seed;
x ^= 0x9e3779b9u + static_cast<unsigned int>(generation + 1) + (x << 6) + (x >> 2);
x ^= 0x85ebca6bu + static_cast<unsigned int>(particleIndex + 1) + (x << 6) + (x >> 2);
x ^= x >> 16;
x *= 0x7feb352du;
x ^= x >> 15;
x *= 0x846ca68bu;
x ^= x >> 16;
return static_cast<double>(x) / 4294967295.0;
}
static QString jsonStringArray(const QStringList& values)
{
// 将字符串列表写成 JSON 数组,用于参数名、启用参数名等字段。
QStringList items;
for(int i = 0; i < values.size(); ++i) {
items << jsonEscape(values[i]);
}
return QString("[%1]").arg(items.join(","));
}
static QString jsonPointArray(const QVector<QPointF>& points)
{
// 将 QPointF 数组写成 [{"x":...,"y":...}],用于 trace meta 保存流量制度。
QStringList items;
for(int i = 0; i < points.size(); ++i) {
items << QString("{\"x\":%1,\"y\":%2}")
.arg(jsonNumber(points[i].x()))
.arg(jsonNumber(points[i].y()));
}
return QString("[%1]").arg(items.join(","));
}
static bool isValidMlRootDirectory(const QString& mlRootPath)
{
// 判断一个目录是不是可用的 ML/nmWTAI-ML 根目录。
// 目前以评分脚本是否存在作为最小可用条件。
if(mlRootPath.isEmpty()) {
return false;
}
QFileInfo scriptInfo(QDir(mlRootPath).absoluteFilePath("scripts/score_pso_candidates.py"));
return scriptInfo.exists() && scriptInfo.isFile();
}
static QString findExecutableInPath(const QString& executableName)
{
// 在系统 PATH 中查找可执行文件。找不到 .venv 里的 python.exe 时会用它兜底。
if(executableName.isEmpty()) {
return QString();
}
QString pathValue = QProcessEnvironment::systemEnvironment().value("PATH");
if(pathValue.isEmpty()) {
pathValue = QString::fromLocal8Bit(qgetenv("PATH"));
}
if(pathValue.isEmpty()) {
return QString();
}
QStringList executableCandidates;
executableCandidates << executableName;
#ifdef Q_OS_WIN
if(!executableName.contains('.')) {
executableCandidates << (executableName + ".exe")
<< (executableName + ".bat")
<< (executableName + ".cmd");
}
#endif
#ifdef Q_OS_WIN
QChar pathSeparator(';');
#else
QChar pathSeparator(':');
#endif
QStringList pathEntries = pathValue.split(pathSeparator, QString::SkipEmptyParts);
for(int i = 0; i < pathEntries.size(); ++i) {
QDir dir(pathEntries[i].trimmed());
if(!dir.exists()) {
continue;
}
for(int j = 0; j < executableCandidates.size(); ++j) {
QFileInfo executableInfo(dir.absoluteFilePath(executableCandidates[j]));
if(executableInfo.exists() && executableInfo.isFile()) {
return QDir::cleanPath(executableInfo.absoluteFilePath());
}
}
}
return QString();
}
static QStringList traceParameterNames()
{
// trace 和 trace meta 使用的完整参数名顺序。
// 这个顺序必须与 buildTraceParameterVector() 和 m_parameterSelected 的 0-10 索引一致。
QStringList names;
names << "k"
<< "skin"
<< "wellboreC"
<< "phi"
<< "initialPressure"
<< "h"
<< "Ct"
<< "Cf"
<< "Soi"
<< "Swi"
<< "Sgi";
return names;
}
// 构造函数:初始化 PSO 参数、收敛判断阈值、trace/代理统计和模拟模式状态。
// 这里设置的是默认值,真正运行前会由 loadOptimizationConfig() 和
// loadParameterBounds() 从 DataManager 读取界面配置覆盖一部分字段。
nmCalculationAutoFitPSO::nmCalculationAutoFitPSO(QObject* parent)
: QObject(parent)
, m_isRunning(false)
, m_shouldStop(false)
, m_isPaused(false)
, m_currentIteration(0)
, m_globalBestFitness(1e10)
, m_previousBestFitness(1e10)
, m_swarmSize(20)
, m_maxIterations(100)
, m_targetError(0.001)
, m_inertiaWeight(0.4)
, m_cognitiveParam(2.0)
, m_socialParam(1.2)
, m_totalEvaluations(0)
, m_successfulEvaluations(0)
, m_evaluationInProgress(0)
, m_consecutiveFailures(0)
, m_improvementThreshold(0.05)
, m_userInitialFitness(1e10)
, m_consecutiveFailedIterations(0)
, m_maxConsecutiveFailures(3)
, m_hasValidUserSolution(false)
, m_diversityThreshold(0.05)
, m_convergenceVarianceThreshold(1e-8)
, m_velocityConvergenceThreshold(0.01)
, m_trueConvergenceWindow(15)
, m_localOptimumWindow(8)
, m_nearTargetFactor(2.0)
, m_farTargetFactor(10.0)
, m_targetWellName("")
, m_traceEnabled(true)
, m_traceRunId("")
, m_traceFilePath("")
, m_traceMetaFilePath("")
, m_surrogateScreeningEnabled(false)
, m_psoRandomSeed(0)
, m_surrogateRunContextSummary("")
, m_surrogateWarmupIterationCount(0)
, m_surrogatePeriodicAuditIterationCount(0)
, m_surrogateContextBlockedIterationCount(0)
, m_surrogateActiveIterationCount(0)
, m_surrogateSelectedParticleCount(0)
, m_surrogateScreenedParticleCount(0)
, m_surrogateTopKParticleCount(0)
, m_surrogateAuditParticleCount(0)
, m_surrogateFallbackParticleCount(0)
, m_surrogateDomainParticleCount(0)
, m_surrogateMinFloorParticleCount(0)
, m_surrogateScoringProcess(nullptr)
, m_surrogateScoringServerKey("")
, m_simulationMode(false)
, m_simulationTimer(nullptr)
, m_simulationIteration(0)
, m_simulationMaxIterations(100)
, m_simulationTargetError(0.001)
, m_simulationStartError(0.5)
, m_simulationCurrentError(0.5)
, m_simulationLogInterval(2)
, m_simulationCurrentParticle(0)
, m_simulationSuccessCount(0)
, m_simulationFailCount(0)
, m_simulationIterationStarted(false)
, m_simulationPreviousIterError(0.5)
{
DEBUG_OUT(QString("PSO Constructor: this=0x%1").arg((quintptr)this, 0, 16));
initializeTemporaryDirectory();
DEBUG_OUT("AutoFit calculator initialized (data-driven mode)");
DEBUG_OUT("PSO Constructor completed");
}
// 析构函数:停止仍在进行的拟合、关闭代理评分服务、关闭 trace 文件并清理临时目录。
// 自动拟合可能在 UI 线程中被窗口关闭打断,因此析构时要尽量温和地等待当前评价结束;
// 如果等待超时,再强制清除运行标志,避免对象销毁后还有信号回调或后台进程访问成员变量。
nmCalculationAutoFitPSO::~nmCalculationAutoFitPSO()
{
DEBUG_OUT(QString("PSO Destructor: this=0x%1").arg((quintptr)this, 0, 16));
if(m_simulationTimer) {
m_simulationTimer->stop();
delete m_simulationTimer;
m_simulationTimer = nullptr;
}
// Stop algorithm.
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 PSO - timeout reached");
m_isRunning = false;
}
}
// Clean temp directory.
stopSurrogateScoringServer();
closeTraceFile();
cleanupTemporaryDirectory();
// Disconnect all signals to avoid callbacks after destruction.
disconnect(this, nullptr, nullptr, nullptr);
DEBUG_OUT("PSO Destructor completed");
}
// ===== 临时目录工具 =====
//
// 真实求解器 DLL 和自动拟合过程会产生中间文件,因此每次创建独立的
// autofit_temp_<pid>_<timestamp> 目录。退出时只删除本类创建的目录,启动时顺便清理
// 旧进程遗留的 autofit_temp_*,避免长期调试后应用目录被临时文件堆满。
void nmCalculationAutoFitPSO::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 nmCalculationAutoFitPSO::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 nmCalculationAutoFitPSO::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 nmCalculationAutoFitPSO::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 nmCalculationAutoFitPSO::setTargetLogLogData(const QVector<QVector<double> >& targetData)
{
// 目标曲线由界面层从目标井 history log-log 传入。
// 约定 targetData[0]=timetargetData[1]=pressuretargetData[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 nmCalculationAutoFitPSO::stopFitting()
{
// 用户点击停止时走这里。停止策略是“请求式停止”:
// 先置 m_shouldStop让主循环/求解器等待逻辑自然退出;短时间内还在评价时再重置计数。
// 这样可以减少 DLL 任务被硬中断导致的数据状态残留。
if(m_simulationMode && m_simulationTimer) {
m_simulationTimer->stop();
}
if(m_isRunning) {
DEBUG_OUT("Stop request received, setting stop flag...");
// 添加停止日志
emit logMessageGenerated(tr("=== User Stop Request Received ==="));
emit logMessageGenerated(tr("Gracefully stopping PSO optimization..."));
m_shouldStop = true;
// 等待当前评估完成,缩短超时时间
int waitCount = 0;
while(m_evaluationInProgress > 0 && waitCount < 30) { // 减少等待时间
QApplication::processEvents(QEventLoop::ExcludeUserInputEvents, 50);
msleep(50);
waitCount++;
}
// 超时时强制重置
if(m_evaluationInProgress > 0) {
DEBUG_OUT("Force resetting evaluation counter");
emit logMessageGenerated(tr("Force stopping current evaluation..."));
m_evaluationInProgress = 0;
}
// 确保运行标志被清除
m_isRunning = false;
emit logMessageGenerated(tr("PSO optimization stop request processed"));
DEBUG_OUT("Stop request processed");
} else {
DEBUG_OUT("Stop request received but PSO is not running");
emit logMessageGenerated(tr("Stop request received but optimization is not running"));
}
closeTraceFile();
cleanupTemporaryDirectory();
}
bool nmCalculationAutoFitPSO::isRunning() const
{
// UI 查询当前是否处于自动拟合运行状态。
return m_isRunning;
}
int nmCalculationAutoFitPSO::getCurrentIteration() const
{
// UI 进度条和日志展示用的当前迭代序号。
return m_currentIteration;
}
QVector<double> nmCalculationAutoFitPSO::getBestSolution() const
{
// 返回紧凑的“启用参数向量”,顺序与 m_enabledParamIndices 一致。
return m_globalBestPosition;
}
double nmCalculationAutoFitPSO::getBestFitness() const
{
// 当前全局最优真实误差。越小越好1e10 附近通常表示尚无有效解。
return m_globalBestFitness;
}
QString nmCalculationAutoFitPSO::getLastError() const
{
// 上一次失败的人类可读错误信息,主要给 UI 层弹窗或日志使用。
return m_lastError;
}
void nmCalculationAutoFitPSO::resetOptimizer()
{
// 清空一次运行产生的状态,但不销毁对象。
// 配置字段会在 startAutoFitting() 中重新从 DataManager 读取;
// trace 文件先关闭,避免新一轮 run 继续写到旧 CSV。
closeTraceFile();
m_swarm.clear();
m_globalBestPosition.clear();
m_globalBestFitness = 1e10;
m_previousBestFitness = 1e10;
m_lastEvaluatedLogLogData.clear();
m_globalBestLogLogData.clear();
m_userInitialLogLogData.clear();
m_currentIteration = 0;
m_totalEvaluations = 0;
m_successfulEvaluations = 0;
m_convergenceHistory.clear();
m_lastError.clear();
m_initialValues.clear();
m_userInitialSolution.clear();
m_userInitialFitness = 1e10;
m_hasValidUserSolution = false;
m_diversityHistory.clear();
m_velocityHistory.clear();
m_particleStagnationHistory.clear();
m_traceMetaFilePath.clear();
resetRunSummary();
DEBUG_OUT("Optimizer reset");
}
void nmCalculationAutoFitPSO::setPSOTargetWellName(const QString& wellName)
{
// 目标井名是贯穿拟合流程的关键索引:
// 读目标曲线、写 skin/wellboreC、求解后取 resultLogLog 都依赖这个名字。
m_targetWellName = wellName;
}
void nmCalculationAutoFitPSO::initializeTraceFile()
{
// 创建本次 PSO 的可复盘文件:
// - pso_baseline_trace_<run_id>.csv逐代逐粒子的参数、真实误差、代理误差和筛选决策
// - pso_baseline_trace_<run_id>.meta.json目标曲线、流量制度、参数上下界、PSO/代理配置。
//
// 代理模型评分脚本也会读取 meta.json因此 trace meta 不是单纯日志而是C++ 与 Python 代理模型之间的运行上下文契约。
if(!m_traceEnabled) {
return;
}
closeTraceFile();
captureSurrogateRunContextSummary();
m_traceRunId = QDateTime::currentDateTime().toString("yyyyMMdd_hhmmss_zzz");
QDir traceDir(QDir(getMlRootPath()).absoluteFilePath("data/temp"));
if(!traceDir.exists() && !QDir().mkpath(traceDir.absolutePath())) {
DEBUG_OUT(QString("Failed to create PSO trace directory: %1").arg(traceDir.absolutePath()));
m_traceFilePath.clear();
return;
}
m_traceFilePath = traceDir.absoluteFilePath(QString("pso_baseline_trace_%1.csv").arg(m_traceRunId));
m_traceMetaFilePath = traceDir.absoluteFilePath(QString("pso_baseline_trace_%1.meta.json").arg(m_traceRunId));
m_traceFile.setFileName(m_traceFilePath);
if(!m_traceFile.open(QIODevice::WriteOnly | QIODevice::Text)) {
DEBUG_OUT(QString("Failed to open PSO trace file: %1").arg(m_traceFilePath));
m_traceFilePath.clear();
m_traceMetaFilePath.clear();
return;
}
writeTraceHeader();
writeTraceMetaFile();
DEBUG_OUT(QString("PSO baseline trace initialized: %1").arg(m_traceFilePath));
emit logMessageGenerated(tr("PSO baseline trace: %1").arg(m_traceFilePath));
if(!m_traceMetaFilePath.isEmpty()) {
emit logMessageGenerated(tr("PSO baseline trace meta: %1").arg(m_traceMetaFilePath));
}
emit logMessageGenerated(tr("PSO surrogate screening: %1, model=%2, keep=%3, audit=%4, warmup=%5, min_solver=%6")
.arg(isSurrogateScreeningEnabled() ? tr("enabled") : tr("disabled"))
.arg(getSurrogateTag())
.arg(getSurrogateKeepFraction(), 0, 'f', 2)
.arg(getSurrogateAuditFraction(), 0, 'f', 2)
.arg(getSurrogateWarmupIterations())
.arg(getSurrogateMinSolverFraction(), 0, 'f', 2));
if(!isSurrogateScreeningEnabled()) {
emit logMessageGenerated(tr("Surrogate screening is turned off in the automatic fitting configuration"));
} else {
if(m_surrogateRunContextSummary == "passed") {
emit logMessageGenerated(tr("Surrogate run context check: passed"));
} else {
emit logMessageGenerated(tr("Surrogate run context check: %1").arg(m_surrogateRunContextSummary));
}
}
}
void nmCalculationAutoFitPSO::resetRunSummary()
{
// 重置本次运行的代理筛选统计。这里只记录统计口径,不参与 PSO 更新。
m_surrogateRunContextSummary.clear();
m_surrogateWarmupIterationCount = 0;
m_surrogatePeriodicAuditIterationCount = 0;
m_surrogateContextBlockedIterationCount = 0;
m_surrogateActiveIterationCount = 0;
m_surrogateSelectedParticleCount = 0;
m_surrogateScreenedParticleCount = 0;
m_surrogateTopKParticleCount = 0;
m_surrogateAuditParticleCount = 0;
m_surrogateFallbackParticleCount = 0;
m_surrogateDomainParticleCount = 0;
m_surrogateMinFloorParticleCount = 0;
}
void nmCalculationAutoFitPSO::captureSurrogateRunContextSummary()
{
// 在 trace 初始化阶段记录一次“本项目是否支持代理模型”的结论。
// 后续每代仍会再次检查,防止运行中 DataManager 状态发生变化。
if(!isSurrogateScreeningEnabled()) {
m_surrogateRunContextSummary = "disabled_by_config";
return;
}
QString contextReason;
if(isSurrogateRunContextSupported(&contextReason)) {
m_surrogateRunContextSummary = "passed";
} else {
m_surrogateRunContextSummary = QString("failed - %1").arg(contextReason);
}
}
void nmCalculationAutoFitPSO::emitRunSummary(bool success, StopReasonPSO finalReason)
{
// 运行结束时输出一组交接/排障最有用的摘要:
// 真实求解次数、成功失败数、代理模型启用情况、各类筛选决策计数,以及 trace 文件路径。
emit logMessageGenerated(tr("=== PSO 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 ? "SUCCESS" : "FAILED")
.arg(m_globalBestFitness, 0, 'e', 4)
.arg(m_currentIteration + 1)
.arg(m_totalEvaluations)
.arg(m_successfulEvaluations)
.arg(m_totalEvaluations - m_successfulEvaluations));
emit logMessageGenerated(tr("Surrogate config: %1, context=%2")
.arg(isSurrogateScreeningEnabled() ? "enabled" : "disabled")
.arg(m_surrogateRunContextSummary.isEmpty() ? QString("unknown") : m_surrogateRunContextSummary));
emit logMessageGenerated(tr("Surrogate iterations: active=%1, warmup=%2, periodic_audit=%3, context_blocked=%4")
.arg(m_surrogateActiveIterationCount)
.arg(m_surrogateWarmupIterationCount)
.arg(m_surrogatePeriodicAuditIterationCount)
.arg(m_surrogateContextBlockedIterationCount));
emit logMessageGenerated(tr("Surrogate particles: selected=%1, screened_out=%2, topk=%3, random_audit=%4, fallback=%5, domain_gate=%6, min_solver_floor=%7")
.arg(m_surrogateSelectedParticleCount)
.arg(m_surrogateScreenedParticleCount)
.arg(m_surrogateTopKParticleCount)
.arg(m_surrogateAuditParticleCount)
.arg(m_surrogateFallbackParticleCount)
.arg(m_surrogateDomainParticleCount)
.arg(m_surrogateMinFloorParticleCount));
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));
}
}
void nmCalculationAutoFitPSO::closeTraceFile()
{
// trace 与代理评分 server 生命周期绑定在一次 PSO run 内。
// 关闭 trace 时同时停止 server避免 Python 进程继续占用模型或 stdout 管道。
stopSurrogateScoringServer();
if(m_traceFile.isOpen()) {
m_traceFile.flush();
m_traceFile.close();
}
}
void nmCalculationAutoFitPSO::writeTraceHeader()
{
// trace CSV 字段说明:
// - 当前粒子参数只记录代理模型关心的 k/skin/wellboreC/phi/h/Cf
// - solver_objective 是真实求解器误差;
// - surrogate_objective 是 Python 代理评分;
// - screening_decision 说明该粒子为什么跑/不跑真实求解器;
// - pbest/gbest 字段用于离线复盘 PSO 更新是否只依赖真实误差。
if(!m_traceFile.isOpen()) {
return;
}
QStringList cols;
cols << "run_id"
<< "generation"
<< "particle_id"
<< "phase"
<< "k"
<< "skin"
<< "wellboreC"
<< "phi"
<< "h"
<< "Cf"
<< "solver_objective"
<< "solver_success"
<< "elapsed_ms"
<< "surrogate_objective"
<< "screening_decision"
<< "pbest_objective"
<< "pbest_k"
<< "pbest_skin"
<< "pbest_wellboreC"
<< "pbest_phi"
<< "pbest_h"
<< "pbest_Cf"
<< "gbest_objective"
<< "gbest_k"
<< "gbest_skin"
<< "gbest_wellboreC"
<< "gbest_phi"
<< "gbest_h"
<< "gbest_Cf"
<< "enabled_param_indices";
QTextStream out(&m_traceFile);
out << cols.join(",") << "\n";
}
void nmCalculationAutoFitPSO::writeTraceMetaFile()
{
// 写 trace 的配套 JSON。它把“如何解释 CSV”的上下文固化下来
// 目标井、目标 log-log 曲线、流量制度、参数名和上下界、初始参数、PSO 配置、
// 代理模型 tag/stage/筛选比例等。Python 评分脚本会读取这个文件来构造模型输入。
if(m_traceMetaFilePath.isEmpty()) {
return;
}
QFile metaFile(m_traceMetaFilePath);
if(!metaFile.open(QIODevice::WriteOnly | QIODevice::Text)) {
DEBUG_OUT(QString("Failed to open PSO trace meta file: %1").arg(m_traceMetaFilePath));
m_traceMetaFilePath.clear();
return;
}
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
nmDataWellBase* pTargetWell = dataManager ? dataManager->findWellByName(m_targetWellName) : nullptr;
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>();
QVector<QPointF> flowPoints = pTargetWell ? pTargetWell->getFlowPoints() : QVector<QPointF>();
QTextStream out(&metaFile);
out << "{\n";
out << " \"schema_version\": 1,\n";
out << " \"trace_type\": \"pso_baseline_replay_meta\",\n";
out << " \"run_id\": " << jsonEscape(m_traceRunId) << ",\n";
out << " \"created_at\": " << jsonEscape(QDateTime::currentDateTime().toString(Qt::ISODate)) << ",\n";
out << " \"trace_csv\": " << jsonEscape(QFileInfo(m_traceFilePath).fileName()) << ",\n";
out << " \"target\": {\n";
out << " \"well_name\": " << jsonEscape(m_targetWellName) << ",\n";
out << " \"well_type\": " << (pTargetWell ? QString::number(static_cast<int>(pTargetWell->getWellType())) : "null") << ",\n";
out << " \"section_policy\": \"current_target_well_indexF\",\n";
out << " \"section_index\": " << (pTargetWell ? QString::number(pTargetWell->getIndexF()) : "null") << ",\n";
out << " \"target_loglog\": {\n";
out << " \"time\": " << jsonDoubleArray(targetTime) << ",\n";
out << " \"pressure\": " << jsonDoubleArray(targetPressure) << ",\n";
out << " \"derivative\": " << jsonDoubleArray(targetDerivative) << "\n";
out << " },\n";
out << " \"flow_points\": " << jsonPointArray(flowPoints) << "\n";
out << " },\n";
out << " \"pso\": {\n";
out << " \"swarm_size\": " << m_swarmSize << ",\n";
out << " \"max_iterations\": " << m_maxIterations << ",\n";
out << " \"target_error\": " << jsonNumber(m_targetError) << ",\n";
out << " \"inertia_weight\": " << jsonNumber(m_inertiaWeight) << ",\n";
out << " \"cognitive_param\": " << jsonNumber(m_cognitiveParam) << ",\n";
out << " \"social_param\": " << jsonNumber(m_socialParam) << ",\n";
out << " \"random_seed\": " << m_psoRandomSeed << "\n";
out << " },\n";
out << " \"surrogate_screening\": {\n";
out << " \"enabled\": " << (isSurrogateScreeningEnabled() ? "true" : "false") << ",\n";
out << " \"model_tag\": " << jsonEscape(getSurrogateTag()) << ",\n";
out << " \"keep_fraction\": " << jsonNumber(getSurrogateKeepFraction()) << ",\n";
out << " \"audit_fraction\": " << jsonNumber(getSurrogateAuditFraction()) << ",\n";
out << " \"min_solver_fraction\": " << jsonNumber(getSurrogateMinSolverFraction()) << ",\n";
out << " \"warmup_iterations\": " << getSurrogateWarmupIterations() << ",\n";
out << " \"full_solver_interval\": " << getSurrogateFullSolverInterval() << ",\n";
out << " \"audit_random_source\": \"deterministic_hash\"\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 << " \"replay_notes\": [\n";
out << " \"CSV rows contain true solver objectives for evaluated particles.\",\n";
out << " \"This meta file contains the target curve and run context needed for surrogate replay.\",\n";
out << " \"pbest and gbest must be updated only from true solver objectives.\"\n";
out << " ]\n";
out << "}\n";
metaFile.flush();
metaFile.close();
}
QVector<double> nmCalculationAutoFitPSO::buildTraceParameterVector(const QVector<double>& selectedParameters) const
{
// 将粒子内部使用的“启用参数向量”还原成完整 11 维参数向量。
// 未启用的参数从当前 DataManager 读取,启用的参数用 selectedParameters 覆盖。
// trace CSV、候选 CSV、代理训练域检查都需要这个完整向量。
QVector<double> fullParams(11, 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.getInitialPressure().getValue().toDouble();
fullParams[5] = reservoirData.getThickness().getValue().toDouble();
fullParams[6] = reservoirData.getCt().getValue().toDouble();
fullParams[7] = reservoirData.getCf().getValue().toDouble();
fullParams[8] = reservoirData.getSoi().getValue().toDouble();
fullParams[9] = reservoirData.getSwi().getValue().toDouble();
fullParams[10] = reservoirData.getSgi().getValue().toDouble();
nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName);
if(pTargetWell) {
fullParams[1] = pTargetWell->getPerforation(0)->getSkin().getValue().toDouble();
fullParams[2] = pTargetWell->getWellboreStorage().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;
}
void nmCalculationAutoFitPSO::writeTraceRow(int generation,
int particleIndex,
const QString& phase,
const QVector<double>& parameters,
double solverObjective,
bool solverSuccess,
int elapsedMs,
double surrogateObjective,
const QString& screeningDecision,
const QVector<double>& pbestPosition,
double pbestObjective)
{
// 写一行 trace。generation=-1/particleIndex=-1 表示用户初始解;
// 普通粒子行的 phase 为 particle_solver 或 particle_not_evaluated。
// 被代理筛掉的粒子 solver_objective 写空值,但仍保留 surrogate_objective 和 decision。
if(!m_traceFile.isOpen()) {
return;
}
QVector<double> currentParams = buildTraceParameterVector(parameters);
QVector<double> pbestParams = pbestPosition.isEmpty() ? QVector<double>() : buildTraceParameterVector(pbestPosition);
QVector<double> gbestParams = m_globalBestPosition.isEmpty() ? QVector<double>() : buildTraceParameterVector(m_globalBestPosition);
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(generation)
<< QString::number(particleIndex)
<< csvEscape(phase)
<< traceParamAt(currentParams, 0)
<< traceParamAt(currentParams, 1)
<< traceParamAt(currentParams, 2)
<< traceParamAt(currentParams, 3)
<< traceParamAt(currentParams, 5)
<< traceParamAt(currentParams, 7)
<< traceNumber(solverObjective)
<< QString::number(solverSuccess ? 1 : 0)
<< QString::number(elapsedMs)
<< traceNumber(surrogateObjective)
<< csvEscape(screeningDecision)
<< traceNumber(pbestObjective)
<< traceParamAt(pbestParams, 0)
<< traceParamAt(pbestParams, 1)
<< traceParamAt(pbestParams, 2)
<< traceParamAt(pbestParams, 3)
<< traceParamAt(pbestParams, 5)
<< traceParamAt(pbestParams, 7)
<< traceNumber(m_globalBestFitness)
<< traceParamAt(gbestParams, 0)
<< traceParamAt(gbestParams, 1)
<< traceParamAt(gbestParams, 2)
<< traceParamAt(gbestParams, 3)
<< traceParamAt(gbestParams, 5)
<< traceParamAt(gbestParams, 7)
<< csvEscape(enabledIndices.join(";"));
QTextStream out(&m_traceFile);
out << cols.join(",") << "\n";
m_traceFile.flush();
}
void nmCalculationAutoFitPSO::writeIterationTraceRows()
{
// 每代真实评价和全局最优更新后写出所有粒子。
// 这样 trace 中能同时看到“这一代谁被筛掉”和“当前 pbest/gbest 是什么”。
if(!m_traceFile.isOpen()) {
return;
}
for(int i = 0; i < m_swarm.size(); ++i) {
const AutoFitParticle& particle = m_swarm[i];
double solverObjective = particle.evaluatedThisIteration
? particle.fitness
: std::numeric_limits<double>::quiet_NaN();
writeTraceRow(m_currentIteration,
i,
particle.evaluatedThisIteration ? "particle_solver" : "particle_not_evaluated",
particle.position,
solverObjective,
particle.lastEvaluationSuccess,
particle.lastEvaluationElapsedMs,
particle.surrogateObjective,
particle.screeningDecision,
particle.bestPosition,
particle.bestFitness);
}
}
bool nmCalculationAutoFitPSO::isSurrogateScreeningEnabled() const
{
// 对应界面上的 PSO acceleration / surrogate screening 开关。
return m_surrogateScreeningEnabled;
}
QString nmCalculationAutoFitPSO::getMlRootPath() const
{
// 查找 ML/nmWTAI-ML 根目录。程序可能从 build/bin/release 等不同目录启动,
// 因此按多个相对路径尝试。isValidMlRootDirectory() 用评分脚本是否存在来确认。
QString appDir = QApplication::applicationDirPath();
QString cwd = QDir::currentPath();
QStringList candidateRoots;
candidateRoots << QDir::cleanPath(appDir + "/../../ML/nmWTAI-ML")
<< QDir::cleanPath(appDir + "/../ML/nmWTAI-ML")
<< QDir::cleanPath(cwd + "/ML/nmWTAI-ML")
<< QDir::cleanPath(cwd + "/../ML/nmWTAI-ML")
<< QDir::cleanPath(cwd + "/../../ML/nmWTAI-ML");
for(int i = 0; i < candidateRoots.size(); ++i) {
if(isValidMlRootDirectory(candidateRoots[i])) {
return candidateRoots[i];
}
}
return candidateRoots.isEmpty()
? QString()
: candidateRoots.first();
}
QString nmCalculationAutoFitPSO::getPythonExecutablePath() const
{
// 优先使用 ML 项目内的虚拟环境 python其次使用 PATH 中的 python。
// 这样开发机上只要 ML/nmWTAI-ML/.venv 已配置,就不依赖全局 Python 环境。
QString mlRoot = getMlRootPath();
QStringList candidateExecutables;
if(!mlRoot.isEmpty()) {
candidateExecutables << QDir(mlRoot).absoluteFilePath(".venv/Scripts/python.exe");
candidateExecutables << QDir(mlRoot).absoluteFilePath(".venv/bin/python");
}
for(int i = 0; i < candidateExecutables.size(); ++i) {
QFileInfo pythonInfo(candidateExecutables[i]);
if(pythonInfo.exists() && pythonInfo.isFile()) {
return QDir::cleanPath(pythonInfo.absoluteFilePath());
}
}
QString resolvedPython = findExecutableInPath("python");
if(!resolvedPython.isEmpty()) {
return QDir::cleanPath(resolvedPython);
}
return "python";
}
QString nmCalculationAutoFitPSO::getSurrogateTag() const
{
// 当前 C++ 绑定的代理模型实验标识,对应 ML 侧训练输出目录/checkpoint 命名。
// 如果以后换模型,需要同步更新这里和 ML 脚本可识别的 tag。
return "family_random_v2_q_50k_noslope";
}
QString nmCalculationAutoFitPSO::getSurrogateStage() const
{
// 当前代理模型的数据生成阶段。v2_q 表示训练时包含流量制度 q 的变化。
return "family_random_v2_q";
}
double nmCalculationAutoFitPSO::getSurrogateKeepFraction() const
{
// 强递减生产制度下代理排序风险更高,因此保留更多 top 粒子进入真实求解器。
if(getSurrogateStage().contains("v2_q") && isStrongDecliningProductionSchedule()) {
return kSurrogateHighRiskKeepFraction;
}
return kSurrogateKeepFraction;
}
double nmCalculationAutoFitPSO::getSurrogateAuditFraction() const
{
// 随机审计比例,审计的是代理模型筛掉的候选。
return kSurrogateAuditFraction;
}
double nmCalculationAutoFitPSO::getSurrogateMinSolverFraction() const
{
// 每代真实求解器最低覆盖比例。强递减生产制度下提高到高风险配置。
if(getSurrogateStage().contains("v2_q") && isStrongDecliningProductionSchedule()) {
return kSurrogateHighRiskMinSolverFraction;
}
return kSurrogateMinSolverFraction;
}
int nmCalculationAutoFitPSO::getSurrogateWarmupIterations() const
{
// 前若干代不筛选,确保粒子先建立真实 pbest。
return kSurrogateWarmupIterations;
}
int nmCalculationAutoFitPSO::getSurrogateFullSolverInterval() const
{
// 周期性全量真实评价间隔,用于审计代理模型排序。
return kSurrogateFullSolverInterval;
}
bool nmCalculationAutoFitPSO::writeSurrogateCandidateCsv(const QString& candidatePath) const
{
// 写给 Python 代理模型的候选粒子 CSV。
//
// 文件位置通常是 ML/nmWTAI-ML/data/temp/pso_screen_candidates_<run_id>_gen<N>.csv。
// 字段顺序必须与 ML 脚本 score_pso_candidates*.py 保持一致:
// particle_id,k,skin,wellboreC,phi,h,Cf
//
// 这里写的是“完整参数体系”中的关键字段,不是粒子内部紧凑向量。
// 例如用户没有勾选 Cf 时Cf 会从当前 DataManager 取值写入 CSV。
QFile file(candidatePath);
if(!file.open(QIODevice::WriteOnly | QIODevice::Text)) {
DEBUG_OUT(QString("Failed to open surrogate candidate file: %1").arg(candidatePath));
return false;
}
QTextStream out(&file);
// 代理模型只需要这些输入字段;其它参数当前不在训练输入集中。
out << "particle_id,k,skin,wellboreC,phi,h,Cf\n";
for(int i = 0; i < m_swarm.size(); ++i) {
QVector<double> params = buildTraceParameterVector(m_swarm[i].position);
QStringList cols;
cols << QString::number(i)
<< traceParamAt(params, 0)
<< traceParamAt(params, 1)
<< traceParamAt(params, 2)
<< traceParamAt(params, 3)
<< traceParamAt(params, 5)
<< traceParamAt(params, 7);
out << cols.join(",") << "\n";
}
file.flush();
file.close();
return true;
}
bool nmCalculationAutoFitPSO::runSurrogateScoringProcess(const QString& candidatePath, const QString& scorePath, QString* failureReason)
{
// 代理评分统一入口。优先走常驻 Python server失败后降级为一次性脚本。
// 常驻 server 可以避免每一代重复加载模型,速度更快;一次性脚本更稳,适合作为兜底。
QString serverFailure;
if(requestSurrogateScoresFromServer(candidatePath, scorePath, &serverFailure)) {
return true;
}
emit logMessageGenerated(tr("Surrogate scoring server unavailable; falling back to one-shot Python scoring"));
if(!serverFailure.isEmpty()) {
emit logMessageGenerated(tr("Surrogate scoring server detail: %1").arg(serverFailure));
}
stopSurrogateScoringServer();
return runSurrogateScoringScriptOnce(candidatePath, scorePath, failureReason);
}
bool nmCalculationAutoFitPSO::ensureSurrogateScoringServer(QString* failureReason)
{
// 确保常驻代理评分 server 已启动且与当前 run 配置匹配。
//
// serverKey 由 python 路径、脚本路径、trace meta、tag、stage 组成。
// 只要其中任何一个变化,就停止旧进程并启动新进程,避免拿旧模型或旧目标曲线评分。
if(m_traceMetaFilePath.isEmpty() || !QFileInfo(m_traceMetaFilePath).exists()) {
QString reason = "trace meta file is missing";
DEBUG_OUT(QString("Surrogate scoring server skipped: %1").arg(reason));
if(failureReason) {
*failureReason = reason;
}
return false;
}
QString mlRoot = getMlRootPath();
QString scriptPath = QDir(mlRoot).absoluteFilePath("scripts/score_pso_candidates_server.py");
if(!QFileInfo(scriptPath).exists()) {
QString reason = QString("surrogate scoring server script not found: %1").arg(scriptPath);
DEBUG_OUT(reason);
if(failureReason) {
*failureReason = reason;
}
return false;
}
QString pythonPath = getPythonExecutablePath();
QString serverKey = QString("%1|%2|%3|%4|%5")
.arg(pythonPath)
.arg(scriptPath)
.arg(m_traceMetaFilePath)
.arg(getSurrogateTag())
.arg(getSurrogateStage());
if(m_surrogateScoringProcess &&
m_surrogateScoringProcess->state() == QProcess::Running &&
m_surrogateScoringServerKey == serverKey) {
return true;
}
stopSurrogateScoringServer();
QStringList args;
// -u 让 Python stdout/stderr 不缓冲C++ 可以及时读到 ready/ok 行。
args << "-u"
<< scriptPath
<< "--trace-meta" << m_traceMetaFilePath
<< "--tag" << getSurrogateTag()
<< "--stage" << getSurrogateStage();
m_surrogateScoringProcess = new QProcess(this);
m_surrogateScoringProcess->setWorkingDirectory(mlRoot);
emit logMessageGenerated(tr("Starting surrogate scoring server: python=%1, mlRoot=%2")
.arg(pythonPath)
.arg(mlRoot));
m_surrogateScoringProcess->start(pythonPath, args);
if(!m_surrogateScoringProcess->waitForStarted(15000)) {
QString reason = QString("failed to start scoring server: python=%1, script=%2, error=%3")
.arg(pythonPath)
.arg(scriptPath)
.arg(m_surrogateScoringProcess->errorString());
DEBUG_OUT(reason);
stopSurrogateScoringServer();
if(failureReason) {
*failureReason = reason;
}
return false;
}
QTime timer;
timer.start();
QByteArray readyLine;
while(timer.elapsed() < 120000) {
// server 启动协议Python 首先输出一行 JSON包含 ready=true 或 ok=false。
// C++ 不做完整 JSON 解析,只检查关键字段,降低 Qt 版本依赖。
if(m_surrogateScoringProcess->state() != QProcess::Running) {
QString stdOut = QString::fromLocal8Bit(m_surrogateScoringProcess->readAllStandardOutput()).trimmed();
QString stdErr = QString::fromLocal8Bit(m_surrogateScoringProcess->readAllStandardError()).trimmed();
QString reason = QString("scoring server exited before ready: stdout=%1, stderr=%2")
.arg(stdOut.isEmpty() ? QString("<empty>") : stdOut)
.arg(stdErr.isEmpty() ? QString("<empty>") : stdErr);
DEBUG_OUT(reason);
stopSurrogateScoringServer();
if(failureReason) {
*failureReason = reason;
}
return false;
}
if(!m_surrogateScoringProcess->canReadLine()) {
m_surrogateScoringProcess->waitForReadyRead(500);
QApplication::processEvents();
continue;
}
readyLine = m_surrogateScoringProcess->readLine().trimmed();
if(readyLine.contains("\"ready\"") && readyLine.contains("true")) {
m_surrogateScoringServerKey = serverKey;
QString stdErr = QString::fromLocal8Bit(m_surrogateScoringProcess->readAllStandardError()).trimmed();
if(!stdErr.isEmpty()) {
emit logMessageGenerated(tr("Surrogate scoring server stderr: %1").arg(stdErr));
}
return true;
}
if(readyLine.contains("\"ok\"") && readyLine.contains("false")) {
QString reason = QString("scoring server failed during initialization: %1").arg(QString::fromLocal8Bit(readyLine));
DEBUG_OUT(reason);
stopSurrogateScoringServer();
if(failureReason) {
*failureReason = reason;
}
return false;
}
}
QString stdErr = QString::fromLocal8Bit(m_surrogateScoringProcess->readAllStandardError()).trimmed();
QString reason = QString("scoring server ready timeout: stderr=%1")
.arg(stdErr.isEmpty() ? QString("<empty>") : stdErr);
DEBUG_OUT(reason);
stopSurrogateScoringServer();
if(failureReason) {
*failureReason = reason;
}
return false;
}
bool nmCalculationAutoFitPSO::requestSurrogateScoresFromServer(const QString& candidatePath, const QString& scorePath, QString* failureReason)
{
// 向常驻 Python server 发送一次评分请求。
// 请求通过 stdin 传 JSON 行server 读取候选 CSV 后把评分写入 scorePath
// 再通过 stdout 返回 ok=true/false。
if(!ensureSurrogateScoringServer(failureReason)) {
return false;
}
QFile::remove(scorePath);
QString command = QString("{\"cmd\":\"score\",\"candidates\":%1,\"output\":%2}\n")
.arg(jsonEscape(candidatePath))
.arg(jsonEscape(scorePath));
QByteArray bytes = command.toUtf8();
qint64 written = m_surrogateScoringProcess->write(bytes);
if(written != bytes.size() || !m_surrogateScoringProcess->waitForBytesWritten(10000)) {
QString reason = QString("failed to send scoring request to server: written=%1 expected=%2 error=%3")
.arg(written)
.arg(bytes.size())
.arg(m_surrogateScoringProcess->errorString());
DEBUG_OUT(reason);
if(failureReason) {
*failureReason = reason;
}
return false;
}
QTime timer;
timer.start();
while(timer.elapsed() < 120000) {
// 每次请求最多等待 120 秒。超时或 server 退出都会让本代代理筛选降级为全量真实求解器。
if(m_surrogateScoringProcess->state() != QProcess::Running) {
QString stdOut = QString::fromLocal8Bit(m_surrogateScoringProcess->readAllStandardOutput()).trimmed();
QString stdErr = QString::fromLocal8Bit(m_surrogateScoringProcess->readAllStandardError()).trimmed();
QString reason = QString("scoring server exited during request: stdout=%1, stderr=%2")
.arg(stdOut.isEmpty() ? QString("<empty>") : stdOut)
.arg(stdErr.isEmpty() ? QString("<empty>") : stdErr);
DEBUG_OUT(reason);
if(failureReason) {
*failureReason = reason;
}
return false;
}
if(!m_surrogateScoringProcess->canReadLine()) {
m_surrogateScoringProcess->waitForReadyRead(500);
QApplication::processEvents();
continue;
}
QByteArray lineBytes = m_surrogateScoringProcess->readLine().trimmed();
QString line = QString::fromLocal8Bit(lineBytes);
if(line.contains("\"ok\"") && line.contains("true")) {
QFileInfo scoreInfo(scorePath);
if(scoreInfo.exists() && scoreInfo.size() > 0) {
QString stdErr = QString::fromLocal8Bit(m_surrogateScoringProcess->readAllStandardError()).trimmed();
if(!stdErr.isEmpty()) {
emit logMessageGenerated(tr("Surrogate scoring server stderr: %1").arg(stdErr));
}
return true;
}
QString reason = QString("scoring server reported success but score file is missing: %1").arg(scorePath);
DEBUG_OUT(reason);
if(failureReason) {
*failureReason = reason;
}
return false;
}
if(line.contains("\"ok\"") && line.contains("false")) {
QString reason = QString("scoring server request failed: %1").arg(line);
DEBUG_OUT(reason);
if(failureReason) {
*failureReason = reason;
}
return false;
}
}
QString reason = "scoring server request timed out";
DEBUG_OUT(reason);
if(failureReason) {
*failureReason = reason;
}
return false;
}
void nmCalculationAutoFitPSO::stopSurrogateScoringServer()
{
// 停止常驻 Python 评分进程。
// 先发送 quit 指令,短时间内不退出再 kill避免应用关闭时遗留 python.exe。
if(!m_surrogateScoringProcess) {
return;
}
if(m_surrogateScoringProcess->state() == QProcess::Running) {
m_surrogateScoringProcess->write("{\"cmd\":\"quit\"}\n");
m_surrogateScoringProcess->waitForBytesWritten(1000);
m_surrogateScoringProcess->waitForFinished(3000);
if(m_surrogateScoringProcess->state() == QProcess::Running) {
m_surrogateScoringProcess->kill();
m_surrogateScoringProcess->waitForFinished(3000);
}
}
delete m_surrogateScoringProcess;
m_surrogateScoringProcess = nullptr;
m_surrogateScoringServerKey.clear();
}
bool nmCalculationAutoFitPSO::runSurrogateScoringScriptOnce(const QString& candidatePath, const QString& scorePath, QString* failureReason)
{
// 一次性代理评分脚本兜底路径。
//
// 输入:
// - candidateswriteSurrogateCandidateCsv() 写出的粒子参数;
// - trace-meta目标曲线、流量制度、参数上下界等上下文
// - tag/stage选择 ML 侧模型。
//
// 输出:
// - scorePathPython 写出的评分 CSV随后由 readSurrogateScores() 读取。
if(m_traceMetaFilePath.isEmpty() || !QFileInfo(m_traceMetaFilePath).exists()) {
QString reason = "trace meta file is missing";
DEBUG_OUT(QString("Surrogate scoring skipped: %1").arg(reason));
if(failureReason) {
*failureReason = reason;
}
return false;
}
QString mlRoot = getMlRootPath();
QString scriptPath = QDir(mlRoot).absoluteFilePath("scripts/score_pso_candidates.py");
if(!QFileInfo(scriptPath).exists()) {
QString reason = QString("surrogate scoring script not found: %1").arg(scriptPath);
DEBUG_OUT(reason);
if(failureReason) {
*failureReason = reason;
}
return false;
}
QString pythonPath = getPythonExecutablePath();
QStringList args;
args << scriptPath
<< "--candidates" << candidatePath
<< "--trace-meta" << m_traceMetaFilePath
<< "--output" << scorePath
<< "--tag" << getSurrogateTag()
<< "--stage" << getSurrogateStage();
QFile::remove(scorePath);
QString lastFailure;
for(int attempt = 1; attempt <= kSurrogateScoringMaxAttempts; ++attempt) {
// 评分脚本可能因首次加载模型、文件占用或环境抖动失败,因此允许有限重试。
emit logMessageGenerated(tr("Surrogate scoring attempt %1/%2: python=%3, mlRoot=%4")
.arg(attempt)
.arg(kSurrogateScoringMaxAttempts)
.arg(pythonPath)
.arg(mlRoot));
QProcess process;
process.setWorkingDirectory(mlRoot);
process.start(pythonPath, args);
if(!process.waitForStarted(10000)) {
lastFailure = QString("failed to start scoring process (attempt %1/%2): python=%3, script=%4, error=%5")
.arg(attempt)
.arg(kSurrogateScoringMaxAttempts)
.arg(pythonPath)
.arg(scriptPath)
.arg(process.errorString());
DEBUG_OUT(lastFailure);
if(attempt < kSurrogateScoringMaxAttempts) {
emit logMessageGenerated(tr("Surrogate scoring start failed, retrying once: %1").arg(process.errorString()));
msleep(500);
continue;
}
if(failureReason) {
*failureReason = lastFailure;
}
return false;
}
if(!process.waitForFinished(120000)) {
process.kill();
process.waitForFinished(3000);
QString stdOut = QString::fromLocal8Bit(process.readAllStandardOutput()).trimmed();
QString stdErr = QString::fromLocal8Bit(process.readAllStandardError()).trimmed();
lastFailure = QString("scoring process timed out (attempt %1/%2): python=%3, script=%4, stdout=%5, stderr=%6")
.arg(attempt)
.arg(kSurrogateScoringMaxAttempts)
.arg(pythonPath)
.arg(scriptPath)
.arg(stdOut.isEmpty() ? QString("<empty>") : stdOut)
.arg(stdErr.isEmpty() ? QString("<empty>") : stdErr);
DEBUG_OUT(lastFailure);
if(attempt < kSurrogateScoringMaxAttempts) {
emit logMessageGenerated(tr("Surrogate scoring timed out, retrying once"));
msleep(500);
continue;
}
if(failureReason) {
*failureReason = lastFailure;
}
return false;
}
QString stdOut = QString::fromLocal8Bit(process.readAllStandardOutput()).trimmed();
QString stdErr = QString::fromLocal8Bit(process.readAllStandardError()).trimmed();
if(process.exitStatus() != QProcess::NormalExit || process.exitCode() != 0) {
lastFailure = QString("scoring process failed (attempt %1/%2): exitStatus=%3, exitCode=%4, python=%5, script=%6, stdout=%7, stderr=%8")
.arg(attempt)
.arg(kSurrogateScoringMaxAttempts)
.arg(process.exitStatus() == QProcess::NormalExit ? "NormalExit" : "CrashExit")
.arg(process.exitCode())
.arg(pythonPath)
.arg(scriptPath)
.arg(stdOut.isEmpty() ? QString("<empty>") : stdOut)
.arg(stdErr.isEmpty() ? QString("<empty>") : stdErr);
DEBUG_OUT(lastFailure);
if(attempt < kSurrogateScoringMaxAttempts) {
emit logMessageGenerated(tr("Surrogate scoring process failed, retrying once: exitCode=%1").arg(process.exitCode()));
msleep(500);
continue;
}
if(failureReason) {
*failureReason = lastFailure;
}
return false;
}
QFileInfo scoreInfo(scorePath);
if(!scoreInfo.exists() || scoreInfo.size() <= 0) {
lastFailure = QString("scoring process finished without a usable score file (attempt %1/%2): output=%3, stdout=%4, stderr=%5")
.arg(attempt)
.arg(kSurrogateScoringMaxAttempts)
.arg(scorePath)
.arg(stdOut.isEmpty() ? QString("<empty>") : stdOut)
.arg(stdErr.isEmpty() ? QString("<empty>") : stdErr);
DEBUG_OUT(lastFailure);
if(attempt < kSurrogateScoringMaxAttempts) {
emit logMessageGenerated(tr("Surrogate scoring produced no score file, retrying once"));
msleep(500);
continue;
}
if(failureReason) {
*failureReason = lastFailure;
}
return false;
}
if(!stdErr.isEmpty()) {
emit logMessageGenerated(tr("Surrogate scoring stderr: %1").arg(stdErr));
}
if(!stdOut.isEmpty()) {
emit logMessageGenerated(tr("Surrogate scoring stdout: %1").arg(stdOut));
}
return true;
}
if(failureReason) {
*failureReason = lastFailure;
}
return false;
}
QVector<double> nmCalculationAutoFitPSO::readSurrogateScores(const QString& scorePath) const
{
// 读取 Python 输出的评分 CSV。
//
// 当前脚本约定至少包含:
// particle_id,surrogate_objective,...,success
// 本函数只读取 particle_id、score 和 success 标志。失败或缺失的粒子保持 NaN
// 后续 buildSurrogateEvaluationMask() 会通过 fallback_gate 强制跑真实求解器。
QVector<double> scores(m_swarm.size(), std::numeric_limits<double>::quiet_NaN());
QFile file(scorePath);
if(!file.open(QIODevice::ReadOnly | QIODevice::Text)) {
DEBUG_OUT(QString("Failed to open surrogate score file: %1").arg(scorePath));
return scores;
}
QTextStream in(&file);
bool firstLine = true;
while(!in.atEnd()) {
QString line = in.readLine().trimmed();
if(line.isEmpty()) {
continue;
}
if(firstLine) {
firstLine = false;
continue;
}
QStringList parts = line.split(",");
if(parts.size() < 5) {
continue;
}
bool idOk = false;
bool scoreOk = false;
int particleId = parts[0].toInt(&idOk);
double score = parts[1].toDouble(&scoreOk);
bool success = parts[4].trimmed() == "1";
if(idOk && scoreOk && success && particleId >= 0 && particleId < scores.size()) {
scores[particleId] = score;
}
}
return scores;
}
bool nmCalculationAutoFitPSO::isSurrogateRunContextSupported(QString* reason) const
{
// 判断“当前项目工况”是否适合使用代理模型。
// 这是运行级别的 gate只要项目整体不满足训练域条件本次 PSO 全程退回真实求解器。
// 这样做比盲目调用代理模型更稳,因为代理模型只在特定 PEBI/单井/垂直井/无复杂几何场景训练过。
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
if(!dataManager) {
if(reason) *reason = "data manager is unavailable";
return false;
}
// 网格类型 gate当前代理模型只按 PEBI 网格样本训练。
if(dataManager->getGridType() != NM_Grid_PEBI) {
if(reason) *reason = "surrogate model is trained only for PEBI grid cases";
return false;
}
nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName);
if(!pTargetWell) {
if(reason) *reason = QString("target well '%1' is missing").arg(m_targetWellName);
return false;
}
// 井型 gate当前只支持垂直井。水平井/压裂井等几何特征不在训练集中。
if(pTargetWell->getWellType() != NM_WELL_MODEL::Vertical_Well) {
if(reason) *reason = "surrogate model is currently limited to vertical-well cases";
return false;
}
int validCalculationWells = 0;
QVector<QPair<NM_WELL_MODEL, QString> > calculationWells = dataManager->getCalculationWells();
for(int i = 0; i < calculationWells.size(); ++i) {
if(calculationWells[i].first != NM_WELL_MODEL::Unknow_Well) {
validCalculationWells++;
}
}
// 单井 gate训练数据假设只有一口有效计算井多井干扰会改变压力响应形态。
if(validCalculationWells != 1) {
if(reason) *reason = QString("surrogate model expects one active calculation well, got %1").arg(validCalculationWells);
return false;
}
// 几何 gate断层、裂缝、区域变化会引入训练集中没有的边界/非均质影响。
if(!dataManager->getFaultDataList().isEmpty() ||
!dataManager->getFractureDataList().isEmpty() ||
!dataManager->getRegionDataList().isEmpty()) {
if(reason) *reason = "surrogate model is not enabled for changed fault/fracture/region geometry";
return false;
}
// 参数 gate代理输入目前只覆盖 k/skin/wellboreC/phi/h。
// 用户勾选其它参数时,代理无法可靠反映这些参数变化,直接禁用代理筛选。
QVector<int> allowedParamIndices;
allowedParamIndices << 0 << 1 << 2 << 3 << 5;
for(int i = 0; i < m_enabledParamIndices.size(); ++i) {
if(!allowedParamIndices.contains(m_enabledParamIndices[i])) {
if(reason) {
QStringList names = traceParameterNames();
int idx = m_enabledParamIndices[i];
QString name = (idx >= 0 && idx < names.size()) ? names[idx] : QString::number(idx);
*reason = QString("parameter '%1' is outside surrogate input set").arg(name);
}
return false;
}
}
nmDataReservoir reservoirData = dataManager->getReservoirDataCopy();
double cf = reservoirData.getCf().getValue().toDouble();
// Cf gate当前模型按近似固定 Cf 训练,允许小范围相对偏差。
if(!isNearRelative(cf, kSurrogateCfFixed, kSurrogateCfRelTol)) {
if(reason) {
*reason = QString("Cf=%1 is outside surrogate fixed-Cf domain").arg(cf, 0, 'g', 10);
}
return false;
}
QVector<QPointF> flowPoints = pTargetWell->getFlowPoints();
// flowPoints 有时第一段是 (0,0) 占位点,不代表真实生产段,先移除。
if(!flowPoints.isEmpty() && qAbs(flowPoints[0].x()) <= 1e-12 && qAbs(flowPoints[0].y()) <= 1e-12) {
flowPoints.remove(0);
}
if(flowPoints.size() < 2) {
if(reason) *reason = "flow schedule is too short for surrogate screening";
return false;
}
// 约定最后一段是关井段,因此生产段数量 = 总段数 - 1。
int productionSections = flowPoints.size() - 1;
if(productionSections < kSurrogateMinProdSections || productionSections > kSurrogateMaxProdSections) {
if(reason) {
*reason = QString("production section count %1 is outside surrogate range [%2,%3]")
.arg(productionSections)
.arg(kSurrogateMinProdSections)
.arg(kSurrogateMaxProdSections);
}
return false;
}
double productionTime = 0.0;
for(int i = 0; i < productionSections; ++i) {
double dt = flowPoints[i].x();
double q = flowPoints[i].y();
if(!isFiniteNumber(dt) || dt <= 0.0 || !isInClosedRange(q, kSurrogateQMin, kSurrogateQMax)) {
if(reason) {
*reason = QString("production schedule section %1 is outside surrogate dt/q domain").arg(i + 1);
}
return false;
}
if(i > 0) {
// 相邻生产段产量跳变不能太大,否则代理模型可能没见过这种制度。
double previousQ = flowPoints[i - 1].y();
double stepRatio = qMax(q, previousQ) / qMax(1e-12, qMin(q, previousQ));
if(stepRatio > kSurrogateMaxRateStepRatio) {
if(reason) {
*reason = QString("rate step ratio %1 exceeds surrogate domain").arg(stepRatio, 0, 'g', 10);
}
return false;
}
}
productionTime += dt;
}
// 旧阶段模型对强递减制度不稳v2_q 阶段允许但会提高真实求解比例。
if(!getSurrogateStage().contains("v2_q")) {
double endStartRatio = std::numeric_limits<double>::quiet_NaN();
if(isStrongDecliningProductionSchedule(&endStartRatio)) {
if(reason) {
*reason = QString("strongly decreasing production schedule ratio %1 is outside surrogate domain")
.arg(endStartRatio, 0, 'g', 10);
}
return false;
}
}
if(!isInClosedRange(productionTime, kSurrogateProdTimeMin, kSurrogateProdTimeMax)) {
if(reason) {
*reason = QString("production total time %1 is outside surrogate range").arg(productionTime, 0, 'g', 10);
}
return false;
}
double shutinTime = flowPoints.last().x();
double shutinRate = flowPoints.last().y();
if(!isInClosedRange(shutinTime, kSurrogateShutinTimeMin, kSurrogateShutinTimeMax) ||
qAbs(shutinRate) > kSurrogateQZeroThreshold) {
if(reason) *reason = "last flow section is not a trained shut-in section";
return false;
}
// 当前目标段必须是最后的关井段,对应试井解释的双对数分析段。
if(pTargetWell->getIndexF() != flowPoints.size()) {
if(reason) {
*reason = QString("target section index %1 is not the last section %2")
.arg(pTargetWell->getIndexF())
.arg(flowPoints.size());
}
return false;
}
return true;
}
bool nmCalculationAutoFitPSO::isSurrogateCandidateInDomain(const QVector<double>& selectedParameters, QString* reason) const
{
// 判断“单个粒子候选参数”是否在代理模型训练域内。
// 与 isSurrogateRunContextSupported() 不同,这里检查的是粒子参数值本身。
// 超出范围的粒子会被 domain_gate 强制跑真实求解器。
QVector<double> params = buildTraceParameterVector(selectedParameters);
double k = params.size() > 0 ? params[0] : std::numeric_limits<double>::quiet_NaN();
double skin = params.size() > 1 ? params[1] : std::numeric_limits<double>::quiet_NaN();
double wellboreC = params.size() > 2 ? params[2] : std::numeric_limits<double>::quiet_NaN();
double phi = params.size() > 3 ? params[3] : std::numeric_limits<double>::quiet_NaN();
double h = params.size() > 5 ? params[5] : std::numeric_limits<double>::quiet_NaN();
double cf = params.size() > 7 ? params[7] : std::numeric_limits<double>::quiet_NaN();
if(!isInClosedRange(k, kSurrogateKMin, kSurrogateKMax)) {
if(reason) *reason = QString("k=%1").arg(k, 0, 'g', 10);
return false;
}
if(!isInClosedRange(skin, kSurrogateSkinMin, kSurrogateSkinMax)) {
if(reason) *reason = QString("skin=%1").arg(skin, 0, 'g', 10);
return false;
}
if(!isInClosedRange(wellboreC, kSurrogateWellboreCMin, kSurrogateWellboreCMax)) {
if(reason) *reason = QString("wellboreC=%1").arg(wellboreC, 0, 'g', 10);
return false;
}
if(!isInClosedRange(phi, kSurrogatePhiMin, kSurrogatePhiMax)) {
if(reason) *reason = QString("phi=%1").arg(phi, 0, 'g', 10);
return false;
}
if(!isInClosedRange(h, kSurrogateHMin, kSurrogateHMax)) {
if(reason) *reason = QString("h=%1").arg(h, 0, 'g', 10);
return false;
}
if(!isNearRelative(cf, kSurrogateCfFixed, kSurrogateCfRelTol)) {
if(reason) *reason = QString("Cf=%1").arg(cf, 0, 'g', 10);
return false;
}
return true;
}
bool nmCalculationAutoFitPSO::forceSolverByFallbackGate(const QVector<double>& selectedParameters) const
{
// 经验 fallback gate。
// 有些点虽然还在训练域边界内,但组合起来属于代理模型容易排序失真的区域,
// 例如极低渗透率配大井筒储集、极端 skin 等。遇到这些候选时直接跑真实求解器。
QVector<double> params = buildTraceParameterVector(selectedParameters);
double k = params.size() > 0 ? params[0] : std::numeric_limits<double>::quiet_NaN();
double skin = params.size() > 1 ? params[1] : std::numeric_limits<double>::quiet_NaN();
double wellboreC = params.size() > 2 ? params[2] : std::numeric_limits<double>::quiet_NaN();
double phi = params.size() > 3 ? params[3] : std::numeric_limits<double>::quiet_NaN();
double h = params.size() > 5 ? params[5] : std::numeric_limits<double>::quiet_NaN();
if(!isFiniteNumber(k) || !isFiniteNumber(skin) || !isFiniteNumber(wellboreC) || !isFiniteNumber(phi) || !isFiniteNumber(h)) {
return true;
}
if(k <= 1e-3 && wellboreC > 0.5) {
return true;
}
if(phi < 0.03 && h < 4.0) {
return true;
}
if(k < 0.1 && h >= 50.0) {
return true;
}
if(qAbs(skin) > 6.0) {
return true;
}
return false;
}
bool nmCalculationAutoFitPSO::isStrongDecliningProductionSchedule(double* endStartRatio) const
{
// 判断生产段产量是否整体非增且首末下降明显。
// 这类流量制度代理模型风险更高,所以 keep/min_solver 会自动提高。
if(endStartRatio) {
*endStartRatio = std::numeric_limits<double>::quiet_NaN();
}
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
if(!dataManager) {
return false;
}
nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName);
if(!pTargetWell) {
return false;
}
QVector<QPointF> flowPoints = pTargetWell->getFlowPoints();
if(!flowPoints.isEmpty() && qAbs(flowPoints[0].x()) <= 1e-12 && qAbs(flowPoints[0].y()) <= 1e-12) {
flowPoints.remove(0);
}
if(flowPoints.size() < 2) {
return false;
}
int productionSections = flowPoints.size() - 1;
bool productionRatesNonIncreasing = true;
bool productionRatesHaveMeaningfulDrop = false;
for(int i = 0; i < productionSections; ++i) {
double q = flowPoints[i].y();
if(!isFiniteNumber(q)) {
return false;
}
if(i > 0) {
double previousQ = flowPoints[i - 1].y();
if(q > previousQ * (1.0 + kSurrogateRateTrendRelTol)) {
productionRatesNonIncreasing = false;
}
if(q < previousQ * (1.0 - kSurrogateRateTrendRelTol)) {
productionRatesHaveMeaningfulDrop = true;
}
}
}
if(!productionRatesNonIncreasing || !productionRatesHaveMeaningfulDrop) {
return false;
}
double firstRate = flowPoints.first().y();
double lastProductionRate = flowPoints[productionSections - 1].y();
double ratio = lastProductionRate / qMax(1e-12, firstRate);
if(endStartRatio) {
*endStartRatio = ratio;
}
return ratio < kSurrogateStrongDeclineEndStartRatio;
}
QVector<bool> nmCalculationAutoFitPSO::buildSurrogateEvaluationMask()
{
// 返回值 evaluateMask 与 m_swarm 一一对应:
// true 表示该粒子本代需要继续跑真实数值求解器;
// false 表示该粒子本代被代理模型筛掉,只记录 screened_out不更新 pbest/gbest。
//
// 这个函数是 PSO acceleration 的唯一入口。为了保证结果可信,它包含多层保护:
// 1. 工况不支持代理模型时,全量真实求解器;
// 2. 前几代 warmup全量真实求解器
// 3. 周期性审计,全量真实求解器;
// 4. 代理评分失败时,全量真实求解器;
// 5. top-k + random audit + fallback + min solver floor 共同决定最终 mask。
QVector<bool> evaluateMask(m_swarm.size(), true);
for(int i = 0; i < m_swarm.size(); ++i) {
m_swarm[i].surrogateObjective = std::numeric_limits<double>::quiet_NaN();
m_swarm[i].selectedForSolver = true;
m_swarm[i].selectedByAudit = false;
m_swarm[i].screeningDecision = "full_solver";
}
if(!isSurrogateScreeningEnabled()) {
// 用户关闭 PSO acceleration保持原始 PSO 行为,所有粒子都跑真实求解器。
return evaluateMask;
}
QString contextReason;
if(!isSurrogateRunContextSupported(&contextReason)) {
// 当前项目工况不在代理模型训练域内。这里不报错,只降级为全量真实求解器,
// 这样用户仍能完成自动拟合,只是没有代理加速收益。
m_surrogateContextBlockedIterationCount++;
for(int i = 0; i < m_swarm.size(); ++i) {
m_swarm[i].screeningDecision = "context_gate_full_solver";
}
emit logMessageGenerated(tr("Surrogate screening disabled for this run context: %1").arg(contextReason));
return evaluateMask;
}
// 前几代全量真实评价,让每个粒子都先获得可信的个体最优 pbest。
// 否则代理模型一开始就筛掉粒子,可能导致某些粒子从未有真实误差基准。
if(m_currentIteration < getSurrogateWarmupIterations()) {
m_surrogateWarmupIterationCount++;
for(int i = 0; i < m_swarm.size(); ++i) {
m_swarm[i].screeningDecision = "warmup_full_solver";
}
return evaluateMask;
}
if(getSurrogateFullSolverInterval() > 0 &&
m_currentIteration > 0 &&
((m_currentIteration + 1) % getSurrogateFullSolverInterval()) == 0) {
// 周期性全量审计。它牺牲少量速度,换取对代理排序偏差的长期纠偏能力。
m_surrogatePeriodicAuditIterationCount++;
for(int i = 0; i < m_swarm.size(); ++i) {
m_swarm[i].screeningDecision = "periodic_full_solver";
}
emit logMessageGenerated(tr("Surrogate audit iteration %1: running full solver for all particles")
.arg(m_currentIteration + 1));
return evaluateMask;
}
QString mlRoot = getMlRootPath();
QDir tempDir(QDir(mlRoot).absoluteFilePath("data/temp"));
if(!tempDir.exists()) {
QDir().mkpath(tempDir.absolutePath());
}
QString candidatePath = tempDir.absoluteFilePath(QString("pso_screen_candidates_%1_gen%2.csv")
.arg(m_traceRunId)
.arg(m_currentIteration));// C++ 写给 Python 的候选粒子参数
QString scorePath = tempDir.absoluteFilePath(QString("pso_screen_scores_%1_gen%2.csv")
.arg(m_traceRunId)
.arg(m_currentIteration));// Python 写回给 C++ 的代理评分结果
QString scoringFailureReason;
if(!writeSurrogateCandidateCsv(candidatePath) || !runSurrogateScoringProcess(candidatePath, scorePath, &scoringFailureReason)) {
// Python 环境、checkpoint、processed 数据或脚本异常时都走这里。
// 降级策略是全量真实求解器,保证代理模块故障不会中断 PSO 拟合。
emit logMessageGenerated(tr("Surrogate screening unavailable; falling back to full solver for iteration %1")
.arg(m_currentIteration + 1));
if(!scoringFailureReason.isEmpty()) {
emit logMessageGenerated(tr("Surrogate scoring failure detail: %1").arg(scoringFailureReason));
}
return evaluateMask;
}
QVector<double> scores = readSurrogateScores(scorePath);
QVector<int> finiteScoreIndices;
// 只接受成功写出的有限代理目标函数。失败候选后面会通过 fallback 强制真实评价。
for(int i = 0; i < scores.size(); ++i) {
m_swarm[i].surrogateObjective = scores[i];
if(isFiniteNumber(scores[i])) {
finiteScoreIndices.append(i);
}
}
if(finiteScoreIndices.isEmpty()) {
emit logMessageGenerated(tr("Surrogate scoring returned no valid scores; falling back to full solver"));
return evaluateMask;
}
std::sort(finiteScoreIndices.begin(), finiteScoreIndices.end(), [&](int a, int b) {
return scores[a] < scores[b];
});
// 代理目标函数越小越好。先将所有粒子设为 screened_out再逐步放回需要真实评价的粒子。
evaluateMask.fill(false);
for(int i = 0; i < m_swarm.size(); ++i) {
m_swarm[i].selectedForSolver = false;
m_swarm[i].screeningDecision = "screened_out";
}
int keepN = qMax(1, qCeil(m_swarm.size() * getSurrogateKeepFraction()));
keepN = qMin(keepN, finiteScoreIndices.size());
// 代理评分 top-k这些粒子最可能给出好解优先跑真实求解器。
for(int rank = 0; rank < keepN; ++rank) {
int idx = finiteScoreIndices[rank];
evaluateMask[idx] = true;
m_swarm[idx].selectedForSolver = true;
m_swarm[idx].screeningDecision = "surrogate_topk";
}
QVector<QPair<double, int> > auditCandidates;
// 随机审计:从代理模型认为“不好”的候选中抽查一部分。
// 使用 deterministicAuditRandom01 保证同一随机种子下 trace 可复现。
for(int i = 0; i < finiteScoreIndices.size(); ++i) {
int idx = finiteScoreIndices[i];
if(!evaluateMask[idx]) {
auditCandidates.append(qMakePair(deterministicAuditRandom01(m_psoRandomSeed, m_currentIteration, idx), idx));
}
}
std::sort(auditCandidates.begin(), auditCandidates.end(), [](const QPair<double, int>& a, const QPair<double, int>& b) {
return a.first < b.first;
});
int auditN = qMin(auditCandidates.size(), qCeil(m_swarm.size() * getSurrogateAuditFraction()));
for(int i = 0; i < auditN; ++i) {
int idx = auditCandidates[i].second;
evaluateMask[idx] = true;
m_swarm[idx].selectedForSolver = true;
m_swarm[idx].selectedByAudit = true;
m_swarm[idx].screeningDecision = "random_audit";
}
int fallbackCount = 0;
int domainFallbackCount = 0;
// fallback gate只要候选点超出代理训练域、代理评分无效、粒子还没有可靠 pbest
// 或落入经验高风险区域,就强制真实求解。这样做是为了保护 PSO 收敛质量。
for(int i = 0; i < m_swarm.size(); ++i) {
QString domainReason;
bool domainFallback = !isSurrogateCandidateInDomain(m_swarm[i].position, &domainReason);
bool fallback = !isFiniteNumber(scores[i]) ||
m_swarm[i].bestFitness >= 1e9 ||
forceSolverByFallbackGate(m_swarm[i].position) ||
domainFallback;
if(fallback) {
if(!evaluateMask[i]) {
fallbackCount++;
}
if(domainFallback) {
domainFallbackCount++;
}
evaluateMask[i] = true;
m_swarm[i].selectedForSolver = true;
m_swarm[i].screeningDecision = domainFallback ? "domain_gate" : "fallback_gate";
}
}
int selectedCount = 0;
for(int i = 0; i < evaluateMask.size(); ++i) {
if(evaluateMask[i]) {
selectedCount++;
}
}
int minSolverN = qMax(1, qCeil(m_swarm.size() * getSurrogateMinSolverFraction()));
if(selectedCount < minSolverN) {
// 最小真实求解比例:防止代理模型过度筛选,导致一代中真实误差信息太少。
QVector<QPair<double, int> > extraCandidates;
for(int i = 0; i < m_swarm.size(); ++i) {
if(!evaluateMask[i] && isFiniteNumber(scores[i])) {
extraCandidates.append(qMakePair(scores[i], i));
}
}
std::sort(extraCandidates.begin(), extraCandidates.end(), [](const QPair<double, int>& a, const QPair<double, int>& b) {
return a.first < b.first;
});
for(int i = 0; i < extraCandidates.size() && selectedCount < minSolverN; ++i) {
int idx = extraCandidates[i].second;
evaluateMask[idx] = true;
m_swarm[idx].selectedForSolver = true;
m_swarm[idx].screeningDecision = "min_solver_floor";
selectedCount++;
}
}
emit logMessageGenerated(tr("Surrogate screening iteration %1: selected %2/%3 particles (keep=%4, audit=%5, min_solver=%6, fallback=%7, domain=%8)")
.arg(m_currentIteration + 1)
.arg(selectedCount)
.arg(m_swarm.size())
.arg(getSurrogateKeepFraction(), 0, 'f', 2)
.arg(getSurrogateAuditFraction(), 0, 'f', 2)
.arg(getSurrogateMinSolverFraction(), 0, 'f', 2)
.arg(fallbackCount)
.arg(domainFallbackCount));
m_surrogateActiveIterationCount++;
for(int i = 0; i < m_swarm.size(); ++i) {
const QString& decision = m_swarm[i].screeningDecision;
if(m_swarm[i].selectedForSolver) {
m_surrogateSelectedParticleCount++;
}
if(decision == "screened_out") {
m_surrogateScreenedParticleCount++;
} else if(decision == "surrogate_topk") {
m_surrogateTopKParticleCount++;
} else if(decision == "random_audit") {
m_surrogateAuditParticleCount++;
} else if(decision == "fallback_gate") {
m_surrogateFallbackParticleCount++;
} else if(decision == "domain_gate") {
m_surrogateDomainParticleCount++;
} else if(decision == "min_solver_floor") {
m_surrogateMinFloorParticleCount++;
}
}
return evaluateMask;
}
//QString nmCalculationAutoFitPSO::getTargetWellName() const
//{
// return m_targetWellName;
//}
// ==================== 数据加载方法 ====================
bool nmCalculationAutoFitPSO::loadAllConfigFromDataManager()
{
// 统一从 DataManager 加载本次运行所需配置。
// UI 层只负责把用户选择保存到 nmDataAutomaticFitting本类从这里开始完全数据驱动。
try {
loadOptimizationConfig();
loadParameterBounds();
extractUserInitialValues(); // 直接提取初始值,无需条件判断
return true;
} catch(...) {
m_lastError = "Failed to load configuration from data manager";
return false;
}
}
void nmCalculationAutoFitPSO::loadOptimizationConfig()
{
// 读取“优化控制”配置:迭代次数、目标误差、是否启用代理筛选。
// PSO 的 swarm/inertia/cognitive/social 目前是类内固定默认值,
// 不是界面输入项;后续若要开放高级设置,可以从这里扩展。
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
nmDataAutomaticFitting fittingData = dataManager->getAutomaticFittingDataCopy();
// 加载基本配置
m_maxIterations = fittingData.getIterationCount().getValue().toInt();
m_targetError = fittingData.getErrorTolerance().getValue().toDouble();
m_surrogateScreeningEnabled = fittingData.getSurrogateScreeningEnabled();
// 设置 PSO 默认参数:
// swarmSize=20 控制每代粒子数inertia/cognitive/social 控制搜索惯性、
// 个体经验和群体经验的权重。adaptiveParameterUpdate() 会在运行中微调这些值。
m_swarmSize = 20;
m_inertiaWeight = 0.4;
m_cognitiveParam = 2.0;
m_socialParam = 1.2;
DEBUG_OUT(QString("Loaded optimization config: iterations=%1, error=%2, surrogate=%3")
.arg(m_maxIterations).arg(m_targetError).arg(m_surrogateScreeningEnabled ? "enabled" : "disabled"));
}
void nmCalculationAutoFitPSO::loadParameterBounds()
{
// 读取用户勾选的拟合参数及上下界。
//
// 这里构建三个核心数组:
// - m_parameterSelected[11]:完整参数体系中每个参数是否参与拟合;
// - m_parameterLower/Upper[11]:完整参数体系的搜索上下界;
// - m_enabledParamIndices把粒子内部紧凑向量映射回完整参数索引。
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
nmDataAutomaticFitting fittingData = dataManager->getAutomaticFittingDataCopy();
// 获取参数选择状态
m_parameterSelected.resize(11);
m_parameterSelected[0] = fittingData.getPermeabilitySelected();
m_parameterSelected[1] = fittingData.getSkinSelected();
m_parameterSelected[2] = fittingData.getWellboreStorageSelected();
m_parameterSelected[3] = fittingData.getPorositySelected();
m_parameterSelected[4] = fittingData.getInitialPressureSelected();
m_parameterSelected[5] = fittingData.getThicknessSelected();
m_parameterSelected[6] = fittingData.getCtSelected();
m_parameterSelected[7] = fittingData.getCfSelected();
m_parameterSelected[8] = fittingData.getSoiSelected();
m_parameterSelected[9] = fittingData.getSwiSelected();
m_parameterSelected[10] = fittingData.getSgiSelected();
// 获取参数边界
m_parameterLower.resize(11);
m_parameterUpper.resize(11);
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.getInitialPressureMin().getValue().toDouble();
m_parameterUpper[4] = fittingData.getInitialPressureMax().getValue().toDouble();
m_parameterLower[5] = fittingData.getThicknessMin().getValue().toDouble();
m_parameterUpper[5] = fittingData.getThicknessMax().getValue().toDouble();
m_parameterLower[6] = fittingData.getCtMin().getValue().toDouble();
m_parameterUpper[6] = fittingData.getCtMax().getValue().toDouble();
m_parameterLower[7] = fittingData.getCfMin().getValue().toDouble();
m_parameterUpper[7] = fittingData.getCfMax().getValue().toDouble();
m_parameterLower[8] = fittingData.getSoiMin().getValue().toDouble();
m_parameterUpper[8] = fittingData.getSoiMax().getValue().toDouble();
m_parameterLower[9] = fittingData.getSwiMin().getValue().toDouble();
m_parameterUpper[9] = fittingData.getSwiMax().getValue().toDouble();
m_parameterLower[10] = fittingData.getSgiMin().getValue().toDouble();
m_parameterUpper[10] = fittingData.getSgiMax().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()));
}
// ==================== PSO算法核心方法 ====================
bool nmCalculationAutoFitPSO::startAutoFitting()
{
// 自动拟合的总入口。可以把这个函数当成 PSO 的“运行剧本”:
// 读取配置 -> 校验输入 -> 评价用户初始解 -> 初始化粒子群 ->
// 按代循环评价粒子 -> 更新全局最优 -> 判断停止 -> 保存结果。
if(m_isRunning) {
m_lastError = "Auto fitting is already running";
return false;
}
// 发送初始化日志
//emit logMessageGenerated(tr("=== PSO Automatic Fitting Started ==="));
emit logMessageGenerated(tr("Algorithm: Particle Swarm Optimization"));
try {
// 从 DataManager 读取界面保存的自动拟合配置。
// 本类不直接依赖 UI 控件,便于后续从脚本或其他入口复用。
if(!loadAllConfigFromDataManager()) {
emit logMessageGenerated(tr("ERROR: Failed to load configuration from data manager"));
return false;
}
if(m_simulationMode) {
// 调试/演示用快速路径,不调用真实求解器。正式工况通常不走这里。
DEBUG_OUT("=== SIMULATION MODE ACTIVATED ===");
runSimulatedFitting();
return true;
}
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.isEmpty() || 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;
}
// 目标曲线必须至少包含三列time、pressure、pressure derivative
// 且三列长度一致。后续误差函数会按 time 做插值对齐。
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;
}
emit logMessageGenerated(tr("Target data validation passed (%1 data points)").arg(m_targetLogLogData[0].size()));
// 使用保存的初始值进行精英保护。resetOptimizer() 会清空部分运行状态,
// 所以先把用户当前模型参数缓存下来,后面再恢复用于初始解评价和粒子初始化。
QVector<double> savedInitialValues = m_initialValues;
// 重置状态
resetOptimizer();
m_isRunning = true;
m_shouldStop = false;
m_isPaused = false;
m_currentIteration = 0;
m_consecutiveFailures = 0;
if(kUseFixedPsoSeed) {
// 固定随机种子便于复现 trace 和代理筛选效果。若做生产随机搜索,
// 可以关闭 kUseFixedPsoSeed改用当前时间生成种子。
m_psoRandomSeed = kFixedPsoSeed;
} else {
QTime seedTime = QTime::currentTime();
m_psoRandomSeed = static_cast<unsigned int>(QDateTime::currentDateTime().toTime_t());
m_psoRandomSeed ^= static_cast<unsigned int>((seedTime.hour() << 22) ^ (seedTime.minute() << 16) ^ (seedTime.second() << 10) ^ seedTime.msec());
}
if(m_psoRandomSeed == 0u) {
m_psoRandomSeed = 1u;
}
qsrand(m_psoRandomSeed);
m_consecutiveFailedIterations = 0; // 重置连续失败迭代计数
// 精英保护:先评价用户当前模型。
// 如果后续 PSO 没找到更好解,最终结果不会比用户已有模型更差。
m_initialValues = savedInitialValues;
initializeTraceFile();
if(!savedInitialValues.isEmpty()) {
m_userInitialSolution = savedInitialValues;
emit logMessageGenerated(tr("=== Evaluating Initial Solution (Elite Protection) ==="));
// 输出初始参数值,方便在 trace / 日志中确认初始解对应的参数顺序。
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..."));
DEBUG_OUT(QString("Before evaluateError: m_userInitialSolution size = %1").arg(m_userInitialSolution.size()));
if(m_userInitialSolution.isEmpty()) {
emit logMessageGenerated(tr("ERROR: m_userInitialSolution is empty!"));
return false;
}
QTime initialEvalTimer;
initialEvalTimer.start();
m_userInitialFitness = evaluateFitness(m_userInitialSolution);
int initialEvalElapsedMs = initialEvalTimer.elapsed();
if(m_userInitialFitness < 1e9) {
m_hasValidUserSolution = true;
m_globalBestFitness = m_userInitialFitness;
m_globalBestPosition = m_userInitialSolution;
m_userInitialLogLogData = m_lastEvaluatedLogLogData;
m_globalBestLogLogData = m_userInitialLogLogData;
emit logMessageGenerated(tr("Initial solution evaluation successful"));
emit logMessageGenerated(tr("Initial Error: %1").arg(m_userInitialFitness, 0, 'e', 4));
emit logMessageGenerated(tr("Elite protection activated - initial solution will be preserved if no significant improvement found"));
emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, 0, m_globalBestFitness);
} else {
m_hasValidUserSolution = false;
emit logMessageGenerated(tr("Initial solution evaluation failed - starting with random initialization"));
}
writeTraceRow(-1,
-1,
"initial_solution",
m_userInitialSolution,
m_userInitialFitness,
m_userInitialFitness < 1e9,
initialEvalElapsedMs,
std::numeric_limits<double>::quiet_NaN(),
"initial_solution",
m_hasValidUserSolution ? m_userInitialSolution : QVector<double>(),
m_userInitialFitness);
} catch(...) {
m_hasValidUserSolution = false;
emit logMessageGenerated(tr("Exception during initial solution evaluation"));
}
// 恢复初始值供粒子群初始化使用
m_initialValues = savedInitialValues;
}
// 初始化粒子群。粒子维度等于用户勾选的参数数量,而不是固定 11 维。
if(kUseFixedPsoSeed) {
emit logMessageGenerated(tr("PSO random seed: %1 ").arg(m_psoRandomSeed));
} else {
emit logMessageGenerated(tr("PSO random seed: %1 ").arg(m_psoRandomSeed));
}
initializeSwarm();
emit logMessageGenerated(tr("Swarm initialized: %1 particles, %2 dimensions").arg(m_swarmSize).arg(getEnabledParameterCount()));
// PSO主循环每一代先决定哪些粒子需要真实评价再更新全局最优和粒子位置。
emit logMessageGenerated(tr("=== Starting PSO Main Loop ==="));
for(m_currentIteration = 0; m_currentIteration < m_maxIterations && !m_shouldStop; ++m_currentIteration) {
// 每次迭代都输出标题,或者只在重要迭代输出详细信息
bool shouldOutputDetail = (m_currentIteration % qMax(1, m_maxIterations / 10) == 0) ||
(m_currentIteration < 5) ||
(m_currentIteration >= m_maxIterations - 2); // 总是输出前5次和最后2次的详细信息
// 每次迭代都输出标题
emit logMessageGenerated(tr("--- Iteration %1/%2 ---").arg(m_currentIteration + 1).arg(m_maxIterations));
// 只在特定迭代输出详细统计信息
if(shouldOutputDetail) {
emit logMessageGenerated(tr("Current best error: %1").arg(m_globalBestFitness, 0, 'e', 4));
emit logMessageGenerated(tr("Total evaluations: %1 (successful: %2, failures: %3)")
.arg(m_totalEvaluations).arg(m_successfulEvaluations).arg(m_totalEvaluations - m_successfulEvaluations));
}
// 检查暂停状态
while(m_isPaused && !m_shouldStop) {
QApplication::processEvents();
msleep(100);
}
if(m_shouldStop) {
emit logMessageGenerated(tr("Optimization stopped by user request"));
break;
}
int currentIterationSuccessful = 0;
int currentIterationTotal = 0;
int currentIterationFailed = 0;
int currentIterationSkipped = 0;
QVector<bool> evaluateMask = buildSurrogateEvaluationMask();
for(int i = 0; i < m_swarmSize && !m_shouldStop; ++i) {
try {
if(i < evaluateMask.size() && !evaluateMask[i]) {
// 被代理模型筛掉的粒子本代不跑真实求解器。
// 这里不会用 surrogateObjective 更新 pbest避免代理误差污染真实最优。
AutoFitParticle& particle = m_swarm[i];
particle.evaluatedThisIteration = false;
particle.lastEvaluationSuccess = false;
particle.lastEvaluationElapsedMs = -1;
particle.selectedForSolver = false;
if(particle.screeningDecision.isEmpty()) {
particle.screeningDecision = "screened_out";
}
currentIterationSkipped++;
continue;
}
double previousFitness = m_swarm[i].fitness;
// 对单个粒子执行真实评价:写参数 -> 跑求解器 -> 算曲线误差 -> 更新 pbest。
updateParticle(i);
currentIterationTotal++;
if(m_swarm[i].fitness < 1e9) {
currentIterationSuccessful++;
if(shouldOutputDetail && m_swarm[i].fitness < previousFitness) {
// 计算改进幅度
double improvement = (previousFitness - m_swarm[i].fitness) / qMax(1e-10, previousFitness);
if(improvement > 0.01) { // 1%以上改进才输出
emit logMessageGenerated(tr(" Particle %1 improved: %2 -> %3 (%4% improvement)")
.arg(i + 1).arg(previousFitness, 0, 'e', 3).arg(m_swarm[i].fitness, 0, 'e', 3)
.arg(improvement * 100, 0, 'f', 2));
}
}
} else {
currentIterationFailed++;
if(shouldOutputDetail) {
emit logMessageGenerated(tr(" Particle %1: evaluation failed").arg(i + 1));
}
}
// 每评估2个粒子处理一次事件
if(i % 2 == 0) {
QApplication::processEvents();
}
} catch(const std::exception& e) {
if(shouldOutputDetail) {
emit logMessageGenerated(tr(" Particle %1: Exception: %2").arg(i + 1).arg(e.what()));
}
m_swarm[i].fitness = 1e10;
currentIterationTotal++;
currentIterationFailed++;
m_totalEvaluations++;
} catch(...) {
if(shouldOutputDetail) {
emit logMessageGenerated(tr(" Particle %1: Unknown exception").arg(i + 1));
}
m_swarm[i].fitness = 1e10;
currentIterationTotal++;
currentIterationFailed++;
m_totalEvaluations++;
}
}
if(m_shouldStop) break;
double currentSuccessRate = currentIterationTotal > 0 ?
(double)currentIterationSuccessful / currentIterationTotal : 0.0;
// 更新连续失败迭代计数
if(currentIterationSuccessful == 0) {
m_consecutiveFailedIterations++;
emit logMessageGenerated(tr("WARNING: No successful evaluations in iteration %1 (consecutive failures: %2)")
.arg(m_currentIteration + 1).arg(m_consecutiveFailedIterations));
} else {
m_consecutiveFailedIterations = 0; // 重置连续失败计数
}
// 输出当前迭代统计
if(shouldOutputDetail) {
emit logMessageGenerated(tr("Current iteration stats: %1 successful, %2 failed, %3 skipped out of %4 particles (solver success rate: %5%)")
.arg(currentIterationSuccessful).arg(currentIterationFailed)
.arg(currentIterationSkipped)
.arg(m_swarmSize).arg(currentSuccessRate * 100, 0, 'f', 1));
}
// 检查是否需要停止优化
if(m_consecutiveFailedIterations >= m_maxConsecutiveFailures) {
m_lastError = QString("Too many consecutive failed iterations (%1)").arg(m_consecutiveFailedIterations);
emit logMessageGenerated(tr("ERROR: Too many consecutive failed iterations (%1/%2) - stopping optimization")
.arg(m_consecutiveFailedIterations).arg(m_maxConsecutiveFailures));
break;
}
// 警告低成功率但不立即停止
if(currentSuccessRate < 0.5 && m_currentIteration > 3) {
emit logMessageGenerated(tr("WARNING: Low success rate (%1%) in iteration %2, but continuing optimization")
.arg(currentSuccessRate * 100, 0, 'f', 1).arg(m_currentIteration + 1));
}
// 更新全局最优。updateGlobalBest() 只读取粒子的真实 bestFitness
// 不使用代理模型的 surrogateObjective。
//double previousGlobalBest = m_globalBestFitness;
updateGlobalBest();
// 记录本代所有粒子的真实/代理误差和筛选决策,用于复盘和排障。
writeIterationTraceRows();
// 自适应参数更新
adaptiveParameterUpdate(m_currentIteration);
emit progressUpdated(m_currentIteration + 1, m_globalBestFitness);
// 记录收敛历史
m_convergenceHistory.append(m_globalBestFitness);
// 每5次迭代输出收敛状态
if((m_currentIteration + 1) % 5 == 0) {
double recentImprovement = 0.0;
if(m_convergenceHistory.size() >= 5) {
double oldBest = m_convergenceHistory[m_convergenceHistory.size() - 5];
recentImprovement = (oldBest - m_globalBestFitness) / qMax(1e-10, oldBest);
}
emit logMessageGenerated(tr(" Convergence check: recent improvement = %1% over last 5 iterations")
.arg(recentImprovement * 100, 0, 'f', 4));
}
// 更新收敛指标
updateConvergenceMetrics();
// 智能收敛判断:目标达成、稳定收敛、局部最优、连续失败等都会在这里决策。
StopReasonPSO stopReason = analyzeOptimizationStatus();
if(stopReason == PSO_TARGET_ACHIEVED) {
emit logMessageGenerated(tr("=== TARGET ACHIEVED ==="));
emit logMessageGenerated(tr("Target error achieved! Current error: %1 < Target: %2")
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_targetError, 0, 'e', 4));
emit logMessageGenerated(tr("Optimization completed successfully after %1 iterations")
.arg(m_currentIteration + 1));
break;
} else if(stopReason == PSO_TRUE_CONVERGENCE) {
emit logMessageGenerated(tr("=== TRUE CONVERGENCE DETECTED ==="));
emit logMessageGenerated(tr("Algorithm has converged to a stable solution"));
emit logMessageGenerated(tr("Final error: %1 after %2 iterations")
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1));
emit logMessageGenerated(tr("Solution quality: %1 (1.0 = target achieved)")
.arg(m_targetError / qMax(1e-10, m_globalBestFitness), 0, 'f', 3));
break;
} else if(stopReason == PSO_LOCAL_OPTIMUM) {
emit logMessageGenerated(tr("=== LOCAL OPTIMUM DETECTED ==="));
emit logMessageGenerated(tr("Algorithm appears to be trapped in local optimum"));
emit logMessageGenerated(tr("Current error: %1 after %2 iterations")
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1));
emit logMessageGenerated(tr("Suggestion: Try restarting with different parameters or larger search space"));
break;
} else if(stopReason == PSO_CONSECUTIVE_FAILURES) {
emit logMessageGenerated(tr("=== CONSECUTIVE FAILURES ==="));
emit logMessageGenerated(tr("Too many consecutive failed iterations (%1/%2)")
.arg(m_consecutiveFailedIterations).arg(m_maxConsecutiveFailures));
break;
} else if(stopReason == PSO_CONTINUE_OPTIMIZATION) {
// 继续优化每10次迭代输出一次状态
if(m_currentIteration > 0 && m_currentIteration % 10 == 0) {
double diversity = calculateSwarmDiversity();
double avgVel = calculateAverageVelocity();
double stagnation = calculateParticleStagnationRate();
emit logMessageGenerated(tr(" Optimization status: diversity=%1, avgVel=%2, stagnation=%3%")
.arg(diversity, 0, 'f', 4).arg(avgVel, 0, 'f', 4).arg(stagnation * 100, 0, 'f', 1));
}
}
// 更新粒子位置和速度
updateVelocityAndPosition();
// 输出迭代结束标记
if(shouldOutputDetail) {
emit logMessageGenerated(tr(" Iteration %1 completed - Current best: %2")
.arg(m_currentIteration + 1).arg(m_globalBestFitness, 0, 'e', 4));
}
// 强制处理事件,保持界面响应
QApplication::processEvents();
// 在迭代间稍作停顿,减少系统负载
if(m_currentIteration % 3 == 2) {
msleep(200);
}
}
// 最终结果验证和保护
validateAndProtectFinalResult();
} catch(const std::exception& e) {
m_lastError = QString(tr("Critical exception in PSO main loop: %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 = QString(tr("Unknown critical exception in PSO main loop"));
emit logMessageGenerated(tr("CRITICAL ERROR: Unknown exception in PSO main loop"));
closeTraceFile();
cleanupTemporaryDirectory();
m_isRunning = false;
emit fittingFinished(false, m_lastError);
return false;
}
m_isRunning = false;
// 应用最终参数
if(!m_globalBestPosition.isEmpty()) {
try {
emit logMessageGenerated(tr("Applying optimized parameters to model..."));
applyParametersToDataManager(m_globalBestPosition);
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);
emit logMessageGenerated(tr("Parameters applied successfully to data manager"));
} catch(const std::exception& e) {
emit logMessageGenerated(tr("ERROR: Failed to apply final parameters: %1").arg(e.what()));
m_lastError = QString("Failed to apply final parameters: %1").arg(e.what());
} catch(...) {
emit logMessageGenerated(tr("ERROR: Unknown error applying final parameters"));
m_lastError = "Failed to apply final parameters due to unknown error";
}
}
// 判断系统确定最终结果
bool success;
QString message;
StopReasonPSO finalReason = analyzeOptimizationStatus();
if(finalReason == PSO_TARGET_ACHIEVED) {
success = true;
message = QString(tr("Target achieved. Best error: %1, Iterations: %2"))
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1);
emit logMessageGenerated(tr("=== PSO OPTIMIZATION SUCCESSFUL ==="));
} else if(finalReason == PSO_TRUE_CONVERGENCE) {
success = true;
message = QString(tr("PSO optimization converged to stable solution. Best error: %1, Iterations: %2"))
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1);
emit logMessageGenerated(tr("=== PSO OPTIMIZATION CONVERGED ==="));
} else if(finalReason == PSO_LOCAL_OPTIMUM) {
success = true;
message = QString(tr("PSO optimization trapped in local optimum. Best error: %1, Iterations: %2"))
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1);
emit logMessageGenerated(tr("=== PSO OPTIMIZATION - LOCAL OPTIMUM ==="));
} else if(finalReason == PSO_MAX_ITERATIONS) {
success = true;
message = QString(tr("Max iterations reached. Best error: %1, Iterations: %2"))
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1);
emit logMessageGenerated(tr("=== PSO OPTIMIZATION - MAX ITERATIONS ==="));
} else if(finalReason == PSO_USER_STOPPED) {
success = true;
message = QString(tr("Best error: %1, Iterations: %2"))
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1);
emit logMessageGenerated(tr("=== PSO OPTIMIZATION STOPPED BY USER ==="));
} else if(finalReason == PSO_CONSECUTIVE_FAILURES) {
success = false;
message = QString(tr("PSO optimization failed due to consecutive failures. Best error: %1, Iterations: %2"))
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1);
emit logMessageGenerated(tr("=== PSO OPTIMIZATION FAILED ==="));
} else {
success = false;
message = QString(tr("PSO optimization ended unexpectedly. Best error: %1, Iterations: %2"))
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_currentIteration + 1);
emit logMessageGenerated(tr("=== PSO OPTIMIZATION - UNKNOWN END ==="));
}
emitRunSummary(success, finalReason);
// 先发送最终进度更新确保进度条达到100%
emit progressUpdated(m_maxIterations, m_globalBestFitness);
// 确保界面更新完成
QApplication::processEvents();
// 短暂延迟,让进度条动画完成
msleep(200);
// 再次确保界面更新
QApplication::processEvents();
// 现在才发送完成信号,这会触发消息框
closeTraceFile();
emit fittingFinished(success, message);
cleanupTemporaryDirectory();
return success;
}
void nmCalculationAutoFitPSO::extractUserInitialValues()
{
// 从当前项目模型读取用户已有初始参数。
// 只提取用户勾选的参数,并按 m_enabledParamIndices 的顺序写入 m_initialValues。
// 这些值用于:
// 1. 初始解真实评价;
// 2. 粒子群引导初始化;
// 3. 最终精英保护。
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.getInitialPressure().getValue().toDouble();
break;
case 5: // 储层厚度
initialValue = reservoirData.getThickness().getValue().toDouble();
break;
case 6: // 综合压缩系数
initialValue = reservoirData.getCt().getValue().toDouble();
break;
case 7: // 岩石压缩系数
initialValue = reservoirData.getCf().getValue().toDouble();
break;
case 8: // 初始含油饱和度
initialValue = reservoirData.getSoi().getValue().toDouble();
break;
case 9: // 初始含水饱和度
initialValue = reservoirData.getSwi().getValue().toDouble();
break;
case 10: // 初始含气饱和度
initialValue = reservoirData.getSgi().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");
}
}
void nmCalculationAutoFitPSO::initializeSwarm()
{
// 初始化粒子群。
// 若用户已有初始模型有效,则一部分粒子围绕初始值做局部搜索;
// 另一部分粒子在全局上下界内随机分布,保留跳出局部区域的能力。
int dimensions = getEnabledParameterCount();
if(dimensions == 0) return;
m_swarm.resize(m_swarmSize);
m_globalBestPosition.resize(dimensions);
bool hasValidInitials = !m_initialValues.isEmpty() && m_initialValues.size() >= dimensions;
// 引导搜索比例。guided 粒子围绕用户初始值,非 guided 粒子做全局随机搜索。
//int guidedCount = hasValidInitials ? qMax(2, (m_swarmSize * 3) / 4) : qMax(1, m_swarmSize / 3);
int guidedCount = hasValidInitials ? qMax(2, m_swarmSize / 2) : qMax(1, m_swarmSize / 3);
DEBUG_OUT(QString("Enhanced swarm initialization: %1 particles, %2 guided, %3 random")
.arg(m_swarmSize).arg(guidedCount).arg(m_swarmSize - guidedCount));
for(int i = 0; i < m_swarmSize; ++i) {
AutoFitParticle& particle = m_swarm[i];
particle.position.resize(dimensions);
particle.velocity.resize(dimensions);
particle.bestPosition.resize(dimensions);
particle.bestFitness = 1e10;
particle.fitness = 1e10;
particle.evaluatedThisIteration = false;
particle.lastEvaluationSuccess = false;
particle.lastEvaluationElapsedMs = -1;
particle.surrogateObjective = 1e10;
particle.selectedForSolver = true;
particle.selectedByAudit = false;
particle.screeningDecision = "full_solver";
for(int j = 0; j < dimensions; ++j) {
int paramIndex = m_enabledParamIndices[j];
double range = m_parameterUpper[paramIndex] - m_parameterLower[paramIndex];
if(i == 0 && hasValidInitials) {
// 第一个粒子:完全使用用户初始值,确保初始模型本身参与粒子群。
particle.position[j] = m_initialValues[j];
} else if(i < guidedCount && hasValidInitials) {
// 引导粒子采用分层搜索半径:
// 近邻精细搜索 + 中等扰动 + 较大扰动,覆盖初始模型周围不同尺度。
double searchRadius;
if(i <= guidedCount / 3) {
searchRadius = range * 0.03; // 3%范围内精细搜索
} else if(i <= guidedCount * 2 / 3) {
searchRadius = range * 0.08; // 8%范围内搜索
} else {
searchRadius = range * 0.15; // 15%范围内搜索
}
double offset = (random01() - 0.5) * searchRadius;
particle.position[j] = m_initialValues[j] + offset;
} else {
// 少数粒子进行全局搜索,防止所有粒子都被初始模型附近的局部最优限制。
particle.position[j] = m_parameterLower[paramIndex] + random01() * range;
}
}
// 边界检查
clampToLimits(particle.position);
// 速度初始化。引导粒子速度更小,避免刚开始就离开初始模型邻域;
// 全局粒子速度稍大,便于探索更宽的搜索空间。
for(int j = 0; j < dimensions; ++j) {
int paramIndex = m_enabledParamIndices[j];
double range = m_parameterUpper[paramIndex] - m_parameterLower[paramIndex];
double velocityFactor;
if(i < guidedCount) {
velocityFactor = VELOCITY_LIMIT_FACTOR * 0.15; // 更小的速度
} else {
velocityFactor = VELOCITY_LIMIT_FACTOR * 0.8; // 适中速度
}
particle.velocity[j] = (random01() - 0.5) * range * velocityFactor;
}
particle.bestPosition = particle.position;
}
}
void nmCalculationAutoFitPSO::updateParticle(int particleIndex)
{
if(particleIndex < 0 || particleIndex >= m_swarm.size()) return;
AutoFitParticle& particle = m_swarm[particleIndex];
// 单粒子真实评价入口。
// 这里调用 evaluateFitness(),因此会真实写 DataManager、调用求解器、计算误差。
// 被代理模型筛掉的粒子不会进入这个函数。
particle.evaluatedThisIteration = false;
particle.lastEvaluationSuccess = false;
particle.lastEvaluationElapsedMs = -1;
particle.selectedForSolver = true;
if(particle.screeningDecision.isEmpty()) {
particle.screeningDecision = "full_solver";
}
QTime evalTimer;
evalTimer.start();
particle.fitness = evaluateFitness(particle.position);
particle.currentLogLogData = m_lastEvaluatedLogLogData;
particle.lastEvaluationElapsedMs = evalTimer.elapsed();
particle.evaluatedThisIteration = true;
particle.lastEvaluationSuccess = (particle.fitness < 1e9);
m_totalEvaluations++;
if(particle.fitness < 1e9) {
m_successfulEvaluations++;
}
// 更新个体最优 pbest。这里使用的是真实求解器误差 particle.fitness
// 不是代理模型给出的 surrogateObjective。
if(particle.fitness < particle.bestFitness) {
particle.bestFitness = particle.fitness;
particle.bestPosition = particle.position;
particle.bestLogLogData = particle.currentLogLogData;
DEBUG_OUT(QString("Particle %1 improved: error = %2")
.arg(particleIndex).arg(particle.fitness, 0, 'e', 4));
}
}
void nmCalculationAutoFitPSO::updateGlobalBest()
{
// 保存上一轮的全局最优,用于后续自适应参数调整等
m_previousBestFitness = m_globalBestFitness;
bool globalBestUpdated = false;
// 遍历所有粒子,寻找比当前 global best 更好的个体最优。
// 注意particle.bestFitness 只有在真实求解器评价成功后才会更新。
for(int i = 0; i < m_swarm.size(); ++i) {
const AutoFitParticle& particle = m_swarm[i];
// 只要个体最优比当前全局最优小,就认为是更好的解
if(particle.bestFitness < m_globalBestFitness) {
double improvement = m_globalBestFitness - particle.bestFitness;
double relativeImprovement =
improvement / qMax(1e-10, qAbs(m_globalBestFitness));
// 区分显著改进和微小改进,但无论如何都会更新全局最优
if(relativeImprovement > m_improvementThreshold) {
DEBUG_OUT(QString(tr("Global best updated with %1% improvement: %2"))
.arg(relativeImprovement * 100.0, 0, 'f', 3)
.arg(particle.bestFitness, 0, 'e', 4));
} else {
DEBUG_OUT(QString(tr("Global best updated (minor improvement %1% < %2%) to %3"))
.arg(relativeImprovement * 100.0, 0, 'f', 3)
.arg(m_improvementThreshold * 100.0, 0, 'f', 2)
.arg(particle.bestFitness, 0, 'e', 4));
}
// 无论相对改进是否超过阈值,都要更新全局最优
m_globalBestFitness = particle.bestFitness;
m_globalBestPosition = particle.bestPosition;
m_globalBestLogLogData = particle.bestLogLogData;
globalBestUpdated = true;
emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, m_currentIteration + 1, m_globalBestFitness);
}
}
// 只有在本轮没有找到任何更好的解时才考虑恢复用户初始解
if(!globalBestUpdated && m_hasValidUserSolution && m_userInitialFitness < m_globalBestFitness) {
double improvement = m_globalBestFitness - m_userInitialFitness;
double relativeImprovement =
improvement / qMax(1e-10, qAbs(m_globalBestFitness));
DEBUG_OUT("Elite protection: Restoring user initial solution as global best");
DEBUG_OUT(QString("Elite solution is better than current best by %1%")
.arg(relativeImprovement * 100.0, 0, 'f', 3));
m_globalBestFitness = m_userInitialFitness;
m_globalBestPosition = m_userInitialSolution;
m_globalBestLogLogData = m_userInitialLogLogData;
emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, m_currentIteration + 1, m_globalBestFitness);
}
}
//void nmCalculationAutoFitPSO::updateGlobalBest()
//{
// m_previousBestFitness = m_globalBestFitness;
// bool globalBestUpdated = false;
//
// for(int i = 0; i < m_swarm.size(); ++i) {
// const AutoFitParticle& particle = m_swarm[i];
//
// if(particle.bestFitness < m_globalBestFitness) {
// // 只有显著改进时才更新
// double improvement = (m_globalBestFitness - particle.bestFitness);
// double relativeImprovement = improvement / qMax(1e-10, qAbs(m_globalBestFitness));
//
// if(relativeImprovement > m_improvementThreshold) {
// m_globalBestFitness = particle.bestFitness;
// m_globalBestPosition = particle.bestPosition;
// globalBestUpdated = true;
//
// DEBUG_OUT(QString("Global best updated with %1% improvement: %2")
// .arg(relativeImprovement * 100, 0, 'f', 3)
// .arg(m_globalBestFitness, 0, 'e', 4));
// } else {
// DEBUG_OUT(QString("Rejecting minor improvement: %1% (threshold: %2%)")
// .arg(relativeImprovement * 100, 0, 'f', 4)
// .arg(m_improvementThreshold * 100, 0, 'f', 2));
// }
// }
// }
//
// // 在没有找到更好解时才检查
// if(!globalBestUpdated && m_hasValidUserSolution && m_userInitialFitness < m_globalBestFitness) {
// DEBUG_OUT("Elite protection: No significant improvement found, checking initial solution");
//
// // 检查初始解是否仍然是最优的
// double improvement = m_globalBestFitness - m_userInitialFitness;
// double relativeImprovement = improvement / qMax(1e-10, qAbs(m_globalBestFitness));
//
// if(relativeImprovement > m_improvementThreshold * 0.5) { // 使用更宽松的阈值
// DEBUG_OUT("Elite protection: Restoring user initial solution");
// m_globalBestFitness = m_userInitialFitness;
// m_globalBestPosition = m_userInitialSolution;
// }
// }
//}
void nmCalculationAutoFitPSO::updateVelocityAndPosition()
{
for(int i = 0; i < m_swarm.size(); ++i) {
AutoFitParticle& particle = m_swarm[i];
for(int j = 0; j < particle.position.size(); ++j) {
// PSO速度更新公式
// inertia 保留上一轮搜索方向;
// cognitive 项让粒子靠近自己的历史最优;
// social 项让粒子靠近群体全局最优。
double r1 = random01();
double r2 = random01();
particle.velocity[j] = m_inertiaWeight * particle.velocity[j] +
m_cognitiveParam * r1 * (particle.bestPosition[j] - particle.position[j]) +
m_socialParam * r2 * (m_globalBestPosition[j] - particle.position[j]);
// 速度限制。不同参数量纲差异很大,所以按该参数搜索区间的一定比例限速。
int paramIndex = m_enabledParamIndices[j];
double range = m_parameterUpper[paramIndex] - m_parameterLower[paramIndex];
double maxVel = range * VELOCITY_LIMIT_FACTOR;
particle.velocity[j] = qMax(-maxVel, qMin(maxVel, particle.velocity[j]));
// 更新位置
particle.position[j] += particle.velocity[j];
}
// 边界约束。PSO 位置更新后如果越界,直接夹回用户设置的上下界。
clampToLimits(particle.position);
}
}
//double nmCalculationAutoFitPSO::evaluateFitness(const QVector<double>& parameters)
//{
// const QString funcName = QString("evaluateError[%1]").arg(m_currentIteration);
// static int callCount = 0;
// callCount++;
//
// 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 - Invalid parameters").arg(funcName).arg(callCount));
// return 1e10;
// }
//
// // 2. 检查数据管理器状态
// nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
//
// if(!dataManager) {
// DEBUG_OUT(QString("%1: Call #%2 - DataManager is null").arg(funcName).arg(callCount));
// return 1e10;
// }
//
// // 3. 应用参数到数据管理器
// try {
// DEBUG_OUT(QString("%1: Call #%2 - Applying parameters to DataManager").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;
// } catch(...) {
// DEBUG_OUT(QString("%1: Call #%2 - Unknown error applying parameters").arg(funcName).arg(callCount));
// return 1e10;
// }
//
// // 4. 运行求解器
// QVector<QVector<double> > solverResult;
// const int maxRetries = 2;
// bool solverSuccess = false;
//
// for(int retry = 0; retry <= maxRetries; ++retry) {
// if(m_shouldStop) 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 before solver attempt")
// .arg(funcName).arg(callCount));
// msleep(1000); // 增加延迟时间
// }
//
// solverResult = runSolver();
//
// if(!solverResult.isEmpty() && validateSolverResult(solverResult)) {
// DEBUG_OUT(QString("%1: Call #%2 - Solver successful on attempt %3, result size: %4")
// .arg(funcName).arg(callCount).arg(retry + 1).arg(solverResult[0].size()));
// solverSuccess = true;
// break;
// } else {
// DEBUG_OUT(QString("%1: Call #%2 - Solver failed on attempt %3 - empty or invalid result")
// .arg(funcName).arg(callCount).arg(retry + 1));
//
// if(retry < maxRetries) {
// DEBUG_OUT(QString("%1: Call #%2 - Will retry solver").arg(funcName).arg(callCount));
// }
// }
//
// } 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()));
// } catch(...) {
// DEBUG_OUT(QString("%1: Call #%2 - Unknown solver exception on attempt %3")
// .arg(funcName).arg(callCount).arg(retry + 1));
// }
// }
//
// if(!solverSuccess) {
// DEBUG_OUT(QString("%1: Call #%2 - All solver attempts failed").arg(funcName).arg(callCount));
// return 1e10;
// }
//
// // 5. 获取双对数结果数据
// QVector<QVector<double> > resultLogLogData;
//
// try {
// //wells = dataManager->getWellDataList(); // 重新获取,确保是最新的
// nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName);
//
// if(pTargetWell) {
// resultLogLogData = pTargetWell->getResultLogLog();
//
// if(!validateLogLogData(resultLogLogData)) {
// DEBUG_OUT(QString("%1: Call #%2 - Invalid result LogLog data").arg(funcName).arg(callCount));
// return 1e10;
// }
//
// DEBUG_OUT(QString("%1: Call #%2 - LogLog result data obtained, size: %3")
// .arg(funcName).arg(callCount).arg(resultLogLogData[0].size()));
// } else {
// DEBUG_OUT(QString("%1: Call #%2 - No wells found after solver").arg(funcName).arg(callCount));
// return 1e10;
// }
// } catch(const std::exception& e) {
// DEBUG_OUT(QString("%1: Call #%2 - Error getting LogLog result: %3").arg(funcName).arg(callCount).arg(e.what()));
// return 1e10;
// } catch(...) {
// DEBUG_OUT(QString("%1: Call #%2 - Unknown error getting LogLog result").arg(funcName).arg(callCount));
// return 1e10;
// }
//
// // 6. 双对数曲线对齐和误差计算
// double error;
//
// try {
// 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 - Evaluation successful, error = %3")
// .arg(funcName).arg(callCount).arg(error, 0, 'e', 6));
//
// } catch(const std::exception& e) {
// DEBUG_OUT(QString("%1: Call #%2 - LogLog error calculation failed: %3").arg(funcName).arg(callCount).arg(e.what()));
// return 1e10;
// } catch(...) {
// DEBUG_OUT(QString("%1: Call #%2 - Unknown error in LogLog error calculation").arg(funcName).arg(callCount));
// 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;
// }
//}
double nmCalculationAutoFitPSO::evaluateFitness(const QVector<double>& parameters)
{
// 粒子适应度评价函数,也是自动拟合最核心的闭环:
// 1. 校验粒子参数是否在用户设置的上下界和基本物理范围内;
// 2. 将参数写入 DataManager 的储层/目标井对象;
// 3. 调用真实数值求解器,生成模拟结果;
// 4. 从目标井读取模拟后的 result log-log 曲线;
// 5. 与目标 history log-log 曲线计算误差,误差越小代表拟合越好。
//
// 返回 1e10 表示该粒子评价失败或结果不可用。PSO 会把它当成很差的解。
const QString funcName = QString("evaluateError[%1]").arg(m_currentIteration);
static int callCount = 0;
callCount++;
m_lastEvaluatedLogLogData.clear();
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. 参数有效性检查。这里先在 PSO 层拦截明显非法的候选,
// 避免把非有限数、越界值或极端危险值传给求解器。
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;
}
// 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 数据。runSolver() 会更新 DataManager 中目标井的计算结果,
// 这里再从目标井读取 resultLogLogData 作为模拟曲线。
QVector<QVector<double>> resultLogLogData;
try {
nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName);
if(!pTargetWell) {
DEBUG_OUT(QString("%1: Call #%2 - Target well '%3' NOT FOUND")
.arg(funcName).arg(callCount).arg(m_targetWellName));
return 1e10;
}
DEBUG_OUT(QString("%1: Call #%2 - Target well found: %3")
.arg(funcName).arg(callCount).arg(m_targetWellName));
resultLogLogData = pTargetWell->getResultLogLog();
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));
// 保存最后一次有效曲线,供 updateParticle() 写入粒子 currentLogLogData。
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 nmCalculationAutoFitPSO::applyParametersToDataManager(const QVector<double>& parameters)
{
// 将粒子的“启用参数向量”写回项目数据。
// parameters 的维度必须等于用户勾选的参数数量,顺序由 m_enabledParamIndices 决定。
// 这里不直接跑求解器,只负责把 DataManager 调整到该粒子对应的模型状态。
if(parameters.size() != getEnabledParameterCount()) {
return;
}
updateReservoirParameters(parameters);
updateWellParameters(parameters);
}
void nmCalculationAutoFitPSO::updateReservoirParameters(const QVector<double>& parameters)
{
// 更新储层级参数。井级参数 skin/wellboreC 不在这里改,由 updateWellParameters() 负责。
// 这里先取 DataManager 中 reservoir 的副本,修改后再整体写回 DataManager。
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
nmDataReservoir reservoirData = dataManager->getReservoirDataCopy();
// paramIndex 是粒子 position 中的索引i 是完整 11 个参数体系中的索引。
// 只有 m_parameterSelected[i] 为 true 时,才从 parameters 中消费一个值。
int paramIndex = 0;
for(int i = 0; i < m_parameterSelected.size() && i < 11; ++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.getInitialPressure().setValue(value);
break;
case 5: // 储层厚度
reservoirData.getThickness().setValue(value);
break;
case 6: // 综合压缩系数
reservoirData.getCt().setValue(value);
break;
case 7: // 岩石压缩系数
reservoirData.getCf().setValue(value);
break;
case 8: // 初始含油饱和度
reservoirData.getSoi().setValue(value);
break;
case 9: // 初始含水饱和度
reservoirData.getSwi().setValue(value);
break;
case 10: // 初始含气饱和度
reservoirData.getSgi().setValue(value);
break;
}
paramIndex++;
}
}
// 更新数据管理器
dataManager->updateReservoirData(reservoirData);
}
void nmCalculationAutoFitPSO::updateWellParameters(const QVector<double>& parameters)
{
// 更新目标井上的拟合参数。目前井级可拟合参数主要是:
// - skin写入第一个 perforation
// - wellboreC写入井筒储集系数。
// 如果目标井不存在或没有射孔数据,这里只记录 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 < 9; ++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;
}
paramIndex++;
}
}
// 根据井类型更新到数据管理器
updateWellToDataManager(pWell);
}
void nmCalculationAutoFitPSO::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>> nmCalculationAutoFitPSO::runSolver()
{
// 真实求解器入口。目前默认走 DLL 方式,保留 runSolverExe() 是历史/备用路径。
// evaluateFitness() 不关心具体求解器形态,只要求返回可验证的双对数结果。
//return runSolverExe();
return runSolverDll();
}
// ==================== 数据处理方法 ====================
QVector<QPointF> nmCalculationAutoFitPSO::interpolateData(
const QVector<QPointF>& source, const QVector<double>& targetX) const
{
// 将一条曲线插值到指定 targetX 网格。
// 双对数误差计算前会分别把目标曲线和模拟曲线插到同一个 commonX 上,
// 这样每个时间点才能逐点比较。
QVector<QPointF> result;
if(source.isEmpty() || targetX.isEmpty()) {
DEBUG_OUT("Warning: Empty data in interpolation");
return result;
}
// 数据清理和验证合并:跳过 NaN/Inf 和极端大值,尽量保留可用点继续计算。
// 这里不立即把整个粒子判失败,是为了提高求解器偶发异常输出下的鲁棒性。
QVector<QPointF> validSource;
const double MAX_REASONABLE_VALUE = 1e12;
for(int i = 0; i < source.size(); ++i) {
const QPointF& point = source[i];
// 检查数值有效性
if(!isFiniteNumber(point.x()) || !isFiniteNumber(point.y())) {
DEBUG_OUT(QString("Skipping invalid data point at index %1: X=%2, Y=%3")
.arg(i).arg(point.x()).arg(point.y()));
continue;
}
// 检查极大值 - 跳过而不是失败
if(qAbs(point.y()) > MAX_REASONABLE_VALUE) {
DEBUG_OUT(QString("Skipping extremely large Y value at index %1: %2")
.arg(i).arg(point.y()));
continue;
}
// 只保留有效的数据点
validSource.append(point);
}
// 检查清理后的数据是否足够
if(validSource.size() < 3) {
DEBUG_OUT(QString("Insufficient valid data points after cleaning: %1")
.arg(validSource.size()));
return result; // 返回空结果,但不算失败
}
DEBUG_OUT(QString("Data cleaning: %1 -> %2 valid points")
.arg(source.size()).arg(validSource.size()));
// 对清理后的数据按 X 排序。后续线性插值要求横坐标单调。
QVector<QPointF> sortedSource = validSource;
for(int i = 0; i < sortedSource.size() - 1; ++i) {
for(int j = 0; j < sortedSource.size() - 1 - i; ++j) {
if(sortedSource[j].x() > sortedSource[j + 1].x()) {
QPointF temp = sortedSource[j];
sortedSource[j] = sortedSource[j + 1];
sortedSource[j + 1] = temp;
}
}
}
double sourceMinX = sortedSource.first().x();
double sourceMaxX = sortedSource.last().x();
double targetMinX = targetX[0];
double targetMaxX = targetX[0];
for(int i = 1; i < targetX.size(); ++i) {
if(targetX[i] < targetMinX) targetMinX = targetX[i];
if(targetX[i] > targetMaxX) targetMaxX = targetX[i];
}
DEBUG_OUT(QString("Source X range: [%1, %2], Target X range: [%3, %4]")
.arg(sourceMinX).arg(sourceMaxX).arg(targetMinX).arg(targetMaxX));
// 插值算法:
// - 在源数据范围内用线性内插;
// - 目标 X 落在源范围外时做有限斜率外推;
// - 找不到区间时使用最近点兜底。
for(int i = 0; i < targetX.size(); ++i) {
double x = targetX[i];
if(!isFiniteNumber(x)) continue;
double y = 0.0;
// 插值逻辑
if(x <= sourceMinX) {
if(sortedSource.size() >= 2) {
double dx = sortedSource[1].x() - sortedSource[0].x();
if(qAbs(dx) > 1e-10) {
double slope = (sortedSource[1].y() - sortedSource[0].y()) / dx;
slope = qMax(-1e6, qMin(1e6, slope));
y = sortedSource[0].y() + slope * (x - sortedSource[0].x());
} else {
y = sortedSource[0].y();
}
} else {
y = sortedSource[0].y();
}
} else if(x >= sourceMaxX) {
if(sortedSource.size() >= 2) {
int lastIdx = sortedSource.size() - 1;
double dx = sortedSource[lastIdx].x() - sortedSource[lastIdx - 1].x();
if(qAbs(dx) > 1e-10) {
double slope = (sortedSource[lastIdx].y() - sortedSource[lastIdx - 1].y()) / dx;
slope = qMax(-1e6, qMin(1e6, slope));
y = sortedSource[lastIdx].y() + slope * (x - sortedSource[lastIdx].x());
} else {
y = sortedSource[lastIdx].y();
}
} else {
y = sortedSource.last().y();
}
} else {
// 内插
bool found = false;
for(int j = 0; j < sortedSource.size() - 1; ++j) {
if(x >= sortedSource[j].x() && x <= sortedSource[j + 1].x()) {
double dx = sortedSource[j + 1].x() - sortedSource[j].x();
if(qAbs(dx) > 1e-10) {
double ratio = (x - sortedSource[j].x()) / dx;
y = sortedSource[j].y() + ratio * (sortedSource[j + 1].y() - sortedSource[j].y());
} else {
y = sortedSource[j].y();
}
found = true;
break;
}
}
if(!found) {
// 使用最近点
double minDist = 1e10;
for(int k = 0; k < sortedSource.size(); ++k) {
double dist = qAbs(sortedSource[k].x() - x);
if(dist < minDist) {
minDist = dist;
y = sortedSource[k].y();
}
}
}
}
// 最终数值检查
if(!isFiniteNumber(y)) {
y = sortedSource.size() > 0 ? sortedSource[0].y() : 1.0;
}
y = qMax(-1e12, qMin(1e12, y));
result.append(QPointF(x, y));
}
if(result.isEmpty()) {
DEBUG_OUT("LogLog interpolation failed");
return result;
}
DEBUG_OUT(QString("Interpolation completed: %1 -> %2 points")
.arg(validSource.size()).arg(result.size()));
return result;
}
// ==================== 算法辅助方法 ====================
void nmCalculationAutoFitPSO::adaptiveParameterUpdate(int iteration)
{
// 自适应调整 PSO 系数。
// 随迭代推进降低惯性权重,让搜索从全局探索逐步转向局部收敛;
// 当近期改进很小时,提高 cognitive、降低 social让粒子更多参考自身历史最优。
double progress = static_cast<double>(iteration) / m_maxIterations;
// 更快降低惯性权重,增强局部搜索
m_inertiaWeight = 0.4 - 0.25 * progress;
if(iteration > 0) {
double improvement = m_previousBestFitness - m_globalBestFitness;
double relativeImprovement = improvement / qMax(1e-10, qAbs(m_previousBestFitness));
if(relativeImprovement < 1e-5) {
// 改进很小时,增强局部搜索
m_cognitiveParam = qMin(2.5, m_cognitiveParam * 1.05);
m_socialParam = qMax(0.8, m_socialParam * 0.95);
} else {
// 有明显改进时,保持平衡
m_cognitiveParam = qMax(1.5, qMin(2.2, m_cognitiveParam));
m_socialParam = qMax(1.0, qMin(1.5, m_socialParam));
}
}
}
void nmCalculationAutoFitPSO::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 nmCalculationAutoFitPSO::validateAndProtectFinalResult()
{
// 最终精英保护。
// PSO 是随机启发式算法,某些工况下可能没有找到比用户初始模型更好的结果。
// 这里用初始解误差与最终全局最优误差做比较,如果改进不足,就恢复初始解,
// 避免“自动拟合”把已有模型调坏。
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));
// 计算改进程度
double improvement = initialFitness - finalFitness;
double relativeImprovement = improvement / qMax(1e-10, qAbs(initialFitness));
emit logMessageGenerated(tr("Improvement: %1 (%2%)")
.arg(improvement, 0, 'e', 4).arg(relativeImprovement * 100, 0, 'f', 2));
if(relativeImprovement < m_improvementThreshold) {
emit logMessageGenerated(tr("Elite protection triggered: insufficient improvement"));
emit logMessageGenerated(tr("Threshold: %1%, Actual: %2%")
.arg(m_improvementThreshold * 100, 0, 'f', 2)
.arg(relativeImprovement * 100, 0, 'f', 4));
emit logMessageGenerated(tr("Restoring initial solution as final result"));
m_globalBestFitness = initialFitness;
m_globalBestPosition = m_userInitialSolution;
m_globalBestLogLogData = m_userInitialLogLogData;
emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, m_currentIteration + 1, m_globalBestFitness);
emit logMessageGenerated(tr("Initial solution restored successfully"));
} else {
emit logMessageGenerated(tr("Final result validated - significant improvement achieved"));
}
}
// ==================== 工具方法 ====================
double nmCalculationAutoFitPSO::random01() const
{
// PSO 使用 Qt 的 qrand/qsrand。随机种子在 startAutoFitting() 中设置,
// 当前默认固定种子,便于复现实验和 trace。
return static_cast<double>(qrand()) / RAND_MAX;
}
void nmCalculationAutoFitPSO::clampToLimits(QVector<double>& parameters) const
{
// 将粒子位置夹回用户给定上下界。parameters 是紧凑的启用参数向量,
// 需要通过 m_enabledParamIndices 找到它在完整参数体系中的上下界。
for(int i = 0; i < parameters.size() && i < m_enabledParamIndices.size(); ++i) {
int paramIndex = m_enabledParamIndices[i];
if(paramIndex >= 0 && paramIndex < m_parameterLower.size()) {
parameters[i] = qMax(m_parameterLower[paramIndex],
qMin(m_parameterUpper[paramIndex], parameters[i]));
}
}
}
int nmCalculationAutoFitPSO::getEnabledParameterCount() const
{
// 返回粒子维度,即用户勾选参与拟合的参数数量。
int count = 0;
for(int i = 0; i < m_parameterSelected.size(); ++i) {
if(m_parameterSelected[i]) count++;
}
return count;
}
// ==================== 验证和处理方法 ====================
bool nmCalculationAutoFitPSO::validateParameters(const QVector<double>& parameters) const
{
// 参数合法性检查分两层:
// 1. 与用户界面设置一致:维度、有限数、上下界;
// 2. 求解器保护:拦截会导致数值崩溃或明显无物理意义的极端值。
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;
}
}
}
}
// 直接拦截会导致求解器数值崩溃的参数值
for(int i = 0; i < parameters.size() && i < m_enabledParamIndices.size(); ++i) {
int paramIndex = m_enabledParamIndices[i];
double value = parameters[i];
switch(paramIndex) {
case 0: // 渗透率:必须大于零
if(value <= 1e-8) {
DEBUG_OUT(QString("Rejecting near-zero permeability: %1").arg(value));
return false;
}
break;
case 2: // 井筒储集系数:必须大于零
if(value <= 1e-10) {
DEBUG_OUT(QString("Rejecting near-zero wellbore storage: %1").arg(value));
return false;
}
break;
case 3: // 孔隙度:必须在合理范围
if(value <= 1e-6 || value >= 0.99) {
DEBUG_OUT(QString("Rejecting unrealistic porosity: %1").arg(value));
return false;
}
break;
case 6: // 综合压缩系数:必须大于零
if(value <= 1e-8) {
DEBUG_OUT(QString("Rejecting near-zero total compressibility: %1").arg(value));
return false;
}
break;
}
}
return true;
}
bool nmCalculationAutoFitPSO::validateLogLogData(const QVector<QVector<double>>& logLogData) const
{
// 校验双对数曲线结构。约定:
// logLogData[0]=timelogLogData[1]=pressurelogLogData[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 nmCalculationAutoFitPSO::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 nmCalculationAutoFitPSO::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 nmCalculationAutoFitPSO::calculateCurveError(
const QVector<QPointF>& curve1, const QVector<QPointF>& curve2) const
{
// 两条已经对齐到同一 X 网格的曲线误差。
// 使用“相对误差为主、绝对误差为辅”的点误差,并对极端值做保护。
if(curve1.size() != curve2.size() || curve1.isEmpty()) {
return 1e10;
}
// 预检查:确保没有无穷大值
for(int i = 0; i < curve1.size(); ++i) {
if(!isFiniteNumber(curve1[i].y()) || !isFiniteNumber(curve2[i].y())) {
DEBUG_OUT(QString("Infinite value detected at index %1: Y1=%2, Y2=%3")
.arg(i).arg(curve1[i].y()).arg(curve2[i].y()));
return 1e10;
}
if(qAbs(curve1[i].y()) > 1e12 || qAbs(curve2[i].y()) > 1e12) {
DEBUG_OUT(QString("Extremely large value detected at index %1")
.arg(i));
return 1e10;
}
}
double totalError = 0.0;
double totalWeight = 0.0;
int validPoints = 0;
for(int i = 0; i < curve1.size(); ++i) {
double y1 = curve1[i].y();
double y2 = curve2[i].y();
// 跳过异常值
if(!isFiniteNumber(y1) || !isFiniteNumber(y2)) {
continue;
}
// 自适应权重Y 值越大权重越小,避免大幅值段完全主导误差。
double weightFactor = qMin(100.0, qAbs(y1) * 0.01);
double weight = 1.0 / (1.0 + weightFactor);
// 相对误差和绝对误差的组合 - 添加数值保护
double yMax = qMax(qAbs(y1), qAbs(y2));
yMax = qMax(1e-12, yMax); // 防止除零
double relativeError = qAbs(y1 - y2) / yMax;
double absoluteError = qAbs(y1 - y2);
// 限制误差值避免爆炸
relativeError = qMin(1e6, relativeError);
absoluteError = qMin(1e6, absoluteError);
// 误差组合:相对误差为主,绝对误差为辅
double pointError = 0.7 * relativeError + 0.3 * absoluteError;
if(isFiniteNumber(pointError) && pointError < 1e10) {
totalError += weight * pointError * pointError;
totalWeight += weight;
validPoints++;
}
}
if(totalWeight > 0 && validPoints > 0) {
double result = sqrt(totalError / totalWeight);
// 最终检查
if(!isFiniteNumber(result)) {
DEBUG_OUT("Final error calculation produced infinite result");
return 1e10;
}
return qMin(1e9, result); // 限制最大误差值
} else {
DEBUG_OUT(QString("No valid points for error calculation: validPoints=%1")
.arg(validPoints));
return 1e10;
}
}
double nmCalculationAutoFitPSO::calculateLogLogCurveError(
const QVector<QVector<double> >& target,
const QVector<QVector<double> >& result) const
{
// 双对数曲线误差计算。
//
// target 通常来自目标井历史曲线result 来自当前粒子参数下的模拟曲线。
// 两条曲线的时间点往往不完全一致,所以这里先取两者时间范围的重叠区间,
// 再在公共时间网格上插值对齐,最后分别计算压力曲线和压力导数曲线误差。
//
// 返回值越小表示拟合越好;返回 1e10 表示曲线无效或无法比较。
// 验证数据
if(!validateLogLogData(target) || !validateLogLogData(result)) {
return 1e10;
}
try {
// 数据对齐:找到目标曲线与模拟曲线 time 轴的重叠区域。
// 不在重叠区域内的点不参与误差,避免外推导致误差失真。
double targetMinX = target[0][0];
double targetMaxX = target[0][0];
for(int i = 1; i < target[0].size(); ++i) {
if(target[0][i] < targetMinX) targetMinX = target[0][i];
if(target[0][i] > targetMaxX) targetMaxX = target[0][i];
}
double resultMinX = result[0][0];
double resultMaxX = result[0][0];
for(int i = 1; i < result[0].size(); ++i) {
if(result[0][i] < resultMinX) resultMinX = result[0][i];
if(result[0][i] > resultMaxX) resultMaxX = result[0][i];
}
double overlapMinX = qMax(targetMinX, resultMinX);
double overlapMaxX = qMin(targetMaxX, resultMaxX);
if(overlapMinX >= overlapMaxX) {
DEBUG_OUT("No overlap between target and result LogLog curves");
return 1e10;
}
// 生成公共 X 网格进行插值。使用对数均匀网格,是为了给早期时间段
// 更多分辨率;试井双对数曲线的早期形态通常对参数识别很敏感。
QVector<double> commonX;
int numPoints = 50;
if(overlapMinX > 0 && overlapMaxX > 0) {
// 对数空间均匀分布
double logMin = qLn(overlapMinX);
double logMax = qLn(overlapMaxX);
for(int i = 0; i < numPoints; ++i) {
double logX = logMin + i * (logMax - logMin) / (numPoints - 1);
double x = qExp(logX);
// 数值保护
if(!isFiniteNumber(x) || x <= 0) {
continue;
}
commonX.append(x);
}
DEBUG_OUT("Using log-uniform grid for better early-time coverage");
}
if(commonX.isEmpty()) {
DEBUG_OUT("Failed to generate common X grid");
return 1e10;
}
// 插值目标曲线。target[1] 是压力target[2] 是压力导数。
QVector<QPointF> targetCurve1, targetCurve2;
for(int i = 0; i < target[0].size(); ++i) {
// 检查数据有效性
if(isFiniteNumber(target[0][i]) && isFiniteNumber(target[1][i]) &&
isFiniteNumber(target[2][i])) {
targetCurve1.append(QPointF(target[0][i], target[1][i]));
targetCurve2.append(QPointF(target[0][i], target[2][i]));
}
}
if(targetCurve1.isEmpty() || targetCurve2.isEmpty()) {
DEBUG_OUT("Target curves are empty after filtering");
return 1e10;
}
QVector<QPointF> alignedTarget1 = interpolateData(targetCurve1, commonX);
QVector<QPointF> alignedTarget2 = interpolateData(targetCurve2, commonX);
// 插值结果曲线。result 与 target 使用同一 commonX保证逐点可比。
QVector<QPointF> resultCurve1, resultCurve2;
for(int i = 0; i < result[0].size(); ++i) {
// 检查数据有效性
if(isFiniteNumber(result[0][i]) && isFiniteNumber(result[1][i]) &&
isFiniteNumber(result[2][i])) {
resultCurve1.append(QPointF(result[0][i], result[1][i]));
resultCurve2.append(QPointF(result[0][i], result[2][i]));
}
}
if(resultCurve1.isEmpty() || resultCurve2.isEmpty()) {
DEBUG_OUT("Result curves are empty after filtering");
return 1e10;
}
QVector<QPointF> alignedResult1 = interpolateData(resultCurve1, commonX);
QVector<QPointF> alignedResult2 = interpolateData(resultCurve2, commonX);
// 检查插值结果
if(alignedTarget1.isEmpty() || alignedTarget2.isEmpty() ||
alignedResult1.isEmpty() || alignedResult2.isEmpty()) {
DEBUG_OUT("LogLog interpolation failed");
return 1e10;
}
if(alignedTarget1.size() != alignedResult1.size() ||
alignedTarget2.size() != alignedResult2.size()) {
DEBUG_OUT("LogLog interpolation size mismatch");
return 1e10;
}
// 计算两条曲线的误差。当前压力和导数各占 50%。
// 如果后续要让导数形态更重要,可以从这里调整权重。
double error1 = calculateCurveError(alignedTarget1, alignedResult1);
double error2 = calculateCurveError(alignedTarget2, alignedResult2);
// 检查个别误差是否有效
if(!isFiniteNumber(error1) || error1 > 1e9) {
DEBUG_OUT(QString("Curve1 error is invalid: %1").arg(error1));
error1 = 1e10;
}
if(!isFiniteNumber(error2) || error2 > 1e9) {
DEBUG_OUT(QString("Curve2 error is invalid: %1").arg(error2));
error2 = 1e10;
}
// 组合误差 - 添加保护
double combinedError;
if(error1 > 1e9 && error2 > 1e9) {
combinedError = 1e10;
} else if(error1 > 1e9) {
combinedError = error2;
} else if(error2 > 1e9) {
combinedError = error1;
} else {
combinedError = 0.5 * error1 + 0.5 * error2;
}
DEBUG_OUT(QString("LogLog errors: Curve1=%1, Curve2=%2, Combined=%3")
.arg(error1, 0, 'e', 4).arg(error2, 0, 'e', 4).arg(combinedError, 0, 'e', 4));
return qMin(1e9, combinedError);
} catch(const std::exception& e) {
DEBUG_OUT(QString("Exception in LogLog error calculation: %1").arg(e.what()));
return 1e10;
} catch(...) {
DEBUG_OUT("Unknown exception in LogLog error calculation");
return 1e10;
}
}
double nmCalculationAutoFitPSO::calculateWeightedPointError(
double target, double result, double timeWeight) const
{
if(!isFiniteNumber(target) || !isFiniteNumber(result)) {
return 1e10;
}
// 对数空间误差:关注数量级差异,适合双对数曲线。
double logTarget = qLn(qMax(1e-10, qAbs(target)));
double logResult = qLn(qMax(1e-10, qAbs(result)));
double logError = qAbs(logTarget - logResult);
// 线性空间相对误差:补充约束实际幅值差异。
double yMax = qMax(qAbs(target), qAbs(result));
double relativeError = qAbs(target - result) / qMax(1e-10, yMax);
// 组合误差:对数误差占 70%,线性误差占 30%。
// timeWeight 当前没有参与最终公式,保留参数是为了后续按时间段加权扩展。
return 0.7 * logError + 0.3 * relativeError;
}
QVector<QVector<double>> nmCalculationAutoFitPSO::runSolverDll()
{
// DLL 求解器路径。
// 这个函数负责把当前 DataManager 中的项目状态交给底层数值求解器,
// 并让目标井生成 result log-log 数据。evaluateFitness() 后续会从目标井读取该结果。
//
// 如果这里失败通常需要优先检查HX_NWTM.dll、license、网格/井数据是否完整、
// 目标井是否存在,以及 DataManager 中刚写入的参数是否导致求解器异常。
DEBUG_OUT("SOLVER DLL START");
if(m_evaluationInProgress > 0) {
DEBUG_OUT("DLL Solver already running, skipping");
return QVector<QVector<double>>();
}
++m_evaluationInProgress;
QVector<QVector<double>> result;
nmCalculationDllPebiSolverTask* dllTask = nullptr;
try {
DEBUG_OUT("Creating DLL solver task");
dllTask = new nmCalculationDllPebiSolverTask(m_tempDirectory);
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。
int waitTime = 0;
const int maxWait = 3600000; // 1h超时
const int checkInterval = 500;
while(waitTime < maxWait) {
QApplication::processEvents(QEventLoop::ExcludeUserInputEvents, 100);
if(!dllTask->isRunning()) {
DEBUG_OUT("DLL solver task completed");
break;
}
if(m_shouldStop) {
DEBUG_OUT("DLL solver task terminated by user");
dllTask->terminate();
break;
}
msleep(checkInterval);
waitTime += checkInterval;
}
// 超时处理
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;
}
// 验证结果数据是否已更新。DLL 任务会把结果写回 DataManager 中的目标井对象。
nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
//QVector<nmDataWellBase*> wells = dataManager->getWellDataList();
nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName);
if(!pTargetWell) {
DEBUG_OUT("No wells found in data manager after DLL execution");
delete dllTask;
--m_evaluationInProgress;
return result;
}
// 验证结果数据。evaluateFitness() 最终用的是 logLogResult
// 但这里返回 pressureResult 给 validateSolverResult() 做基本求解成功判断。
QVector<QVector<double>> pressureResult = pTargetWell->getResultPressure();
QVector<QVector<double>> logLogResult = pTargetWell->getResultLogLog();
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) {
result = pressureResult;
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;
}
//QVector<QVector<double>> nmCalculationAutoFitPSO::runSolverExe()
//{
// DEBUG_OUT("SOLVER EXE START");
//
// if(m_evaluationInProgress > 0) {
// DEBUG_OUT("EXE Solver already running, skipping");
// return QVector<QVector<double>>();
// }
//
// ++m_evaluationInProgress;
// QVector<QVector<double>> result;
// nmCalculationExeSolverTask* exeTask = nullptr;
//
// try {
// DEBUG_OUT("Creating EXE solver task");
// exeTask = new nmCalculationExeSolverTask(QString());
//
// if(m_shouldStop) {
// DEBUG_OUT("Should stop - cleaning up and returning empty result");
// delete exeTask;
// --m_evaluationInProgress;
// return result;
// }
//
// DEBUG_OUT("Starting EXE solver execution...");
//
// // 启动
// bool executeSuccess = exeTask->execute();
//
// // 检查执行状态
// if(!executeSuccess) {
// DEBUG_OUT(QString("EXE solver execution failed: %1").arg(exeTask->getLastError()));
// DEBUG_OUT(QString("EXE solver exit code: %1").arg(exeTask->getExitCode()));
//
// // 清理并返回空结果
// delete exeTask;
// exeTask = nullptr;
// --m_evaluationInProgress;
// m_consecutiveFailures++;
// return result;
// }
//
// DEBUG_OUT("EXE solver execution completed successfully");
//
// // 验证结果数据是否已更新
// nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance();
// QVector<nmDataWellBase*> wells = dataManager->getWellDataList();
//
// if(wells.isEmpty()) {
// DEBUG_OUT("No wells found in data manager after EXE execution");
// delete exeTask;
// --m_evaluationInProgress;
// return result;
// }
//
// // 验证结果数据的时间戳或唯一性
// QVector<QVector<double>> pressureResult = wells[0]->getResultPressure();
// QVector<QVector<double>> logLogResult = wells[0]->getResultLogLog();
//
// DEBUG_OUT(QString("Raw 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(logLogResult.size() >= 3) {
// DEBUG_OUT(QString("LogLog result - X points: %1, Y1 points: %2, Y2 points: %3")
// .arg(logLogResult[0].size()).arg(logLogResult[1].size()).arg(logLogResult[2].size()));
//
// if(logLogResult[0].size() > 0) {
// DEBUG_OUT(QString("Sample LogLog data - X[0]: %1, X[last]: %2, Y1[0]: %3, Y2[0]: %4")
// .arg(logLogResult[0][0])
// .arg(logLogResult[0][logLogResult[0].size()-1])
// .arg(logLogResult[1][0])
// .arg(logLogResult[2][0]));
// }
// }
//
// // 数据有效性检查
// if(pressureResult.size() >= 2 && pressureResult[0].size() > 0 && pressureResult[1].size() > 0) {
// result = pressureResult;
// DEBUG_OUT(QString("Got EXE solver result: %1 points").arg(result[0].size()));
// m_consecutiveFailures = 0;
//
// // 检查结果是否与之前不同
// 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 EXE");
// lastPressureResult = pressureResult[1];
// } else {
// DEBUG_OUT("!!! WARNING: Result data appears to be identical to previous run !!!");
// }
//
// } else {
// DEBUG_OUT("EXE 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::exception& e) {
// DEBUG_OUT(QString("Exception in EXE solver: %1").arg(e.what()));
// m_consecutiveFailures++;
// } catch(...) {
// DEBUG_OUT("Unknown exception in EXE solver");
// m_consecutiveFailures++;
// }
//
// // 清理EXE任务
// if(exeTask) {
// DEBUG_OUT("Cleaning up EXE solver task...");
// delete exeTask;
// exeTask = nullptr;
// }
//
// QApplication::processEvents(QEventLoop::AllEvents, 100);
// --m_evaluationInProgress;
//
// DEBUG_OUT(QString("SOLVER EXE END - ResultPoints: %1")
// .arg(result.isEmpty() ? 0 : result[0].size()));
//
// return result;
//}
StopReasonPSO nmCalculationAutoFitPSO::analyzeOptimizationStatus()
{
// 按优先级判断本轮是否停止。
// 目标达成和最大迭代是硬条件;真收敛/局部最优需要足够历史数据支撑。
// 返回值只描述原因,真正的日志和 UI 收尾在 startAutoFitting() 中处理。
// 1. 检查用户停止
if(m_shouldStop) {
return PSO_USER_STOPPED;
}
// 2. 检查连续失败
if(m_consecutiveFailedIterations >= m_maxConsecutiveFailures) {
return PSO_CONSECUTIVE_FAILURES;
}
// 3. 检查是否达到目标精度
if(m_globalBestFitness < m_targetError) {
return PSO_TARGET_ACHIEVED;
}
// 4. 检查是否达到最大迭代数
if(m_currentIteration >= m_maxIterations - 1) {
return PSO_MAX_ITERATIONS;
}
// 5. 需要足够的历史数据才能判断收敛
if(m_convergenceHistory.size() < m_localOptimumWindow) {
return PSO_CONTINUE_OPTIMIZATION;
}
// 6. 检查真正的收敛
if(m_convergenceHistory.size() >= m_trueConvergenceWindow && checkTrueConvergence()) {
return PSO_TRUE_CONVERGENCE;
}
// 7. 检查局部最优陷阱
if(checkLocalOptimumTrap()) {
return PSO_LOCAL_OPTIMUM;
}
return PSO_CONTINUE_OPTIMIZATION;
}
bool nmCalculationAutoFitPSO::checkTrueConvergence() const
{
// 真收敛判断:不是只看“最近没有改进”,而是同时看解质量、误差稳定性、
// 粒子群多样性、平均速度、长期改进和粒子停滞率。
// 这样可以减少把局部卡住误判成正常收敛的概率。
if(m_convergenceHistory.size() < m_trueConvergenceWindow) {
return false;
}
// 1. 检查解质量 - 如果已经接近目标,小改进可能是真收敛
bool nearTarget = (m_globalBestFitness < m_targetError * m_nearTargetFactor);
// 2. 检查适应度稳定性 - 长期小幅波动
double recentVariance = calculateFitnessVariance(10);
double recentMean = 0.0;
int windowSize = qMin(10, m_convergenceHistory.size());
// 计算最近窗口的均值
for(int i = m_convergenceHistory.size() - windowSize; i < m_convergenceHistory.size(); ++i) {
recentMean += m_convergenceHistory[i];
}
recentMean /= windowSize;
double relativeVariance = recentVariance / qMax(1e-10, recentMean * recentMean);
bool stableError = (relativeVariance < m_convergenceVarianceThreshold);
// 3. 检查粒子群多样性 - 应该收敛到同一区域
double currentDiversity = calculateSwarmDiversity();
bool lowDiversity = (currentDiversity < m_diversityThreshold);
// 4. 检查速度收敛 - 粒子应该几乎停止移动
double avgVelocity = calculateAverageVelocity();
bool velocityConverged = (avgVelocity < m_velocityConvergenceThreshold);
// 5. 检查长期改进趋势
double longTermImprovement = calculateLongTermImprovement(m_trueConvergenceWindow);
bool minimalLongTermImprovement = (longTermImprovement < 1e-4); // 0.01%
// 6. 检查粒子停滞情况
double stagnationRate = calculateParticleStagnationRate();
bool mostParticlesConverged = (stagnationRate > 0.8); // 80%以上粒子收敛
// 真收敛的判断条件
bool isConverged = nearTarget ||
(stableError && lowDiversity && velocityConverged &&
minimalLongTermImprovement && mostParticlesConverged);
if(isConverged) {
DEBUG_OUT("=== TRUE CONVERGENCE ANALYSIS ===");
DEBUG_OUT(QString("nearTarget=%1 (error=%2, target*factor=%3)")
.arg(nearTarget).arg(m_globalBestFitness, 0, 'e', 4)
.arg(m_targetError * m_nearTargetFactor, 0, 'e', 4));
DEBUG_OUT(QString("stableError=%1 (relativeVariance=%2)")
.arg(stableError).arg(relativeVariance, 0, 'e', 6));
DEBUG_OUT(QString("lowDiversity=%1 (diversity=%2, threshold=%3)")
.arg(lowDiversity).arg(currentDiversity, 0, 'f', 6).arg(m_diversityThreshold));
DEBUG_OUT(QString("velocityConverged=%1 (avgVel=%2)")
.arg(velocityConverged).arg(avgVelocity, 0, 'f', 6));
DEBUG_OUT(QString("longTermImprovement=%1%, stagnationRate=%2%")
.arg(longTermImprovement * 100, 0, 'f', 4).arg(stagnationRate * 100, 0, 'f', 1));
}
return isConverged;
}
bool nmCalculationAutoFitPSO::checkLocalOptimumTrap() const
{
// 局部最优判断:当误差离目标还远,但短期几乎不改进,
// 且粒子群过度聚集/异常分散或大部分粒子停滞时,认为可能卡住。
if(m_convergenceHistory.size() < m_localOptimumWindow) {
return false;
}
// 1. 检查解质量 - 如果距离目标还很远,停滞就可能是局部最优
bool farFromTarget = (m_globalBestFitness > m_targetError * m_farTargetFactor);
// 2. 检查短期改进 - 近期改进非常小
double shortTermImprovement = calculateLongTermImprovement(m_localOptimumWindow);
bool poorShortTermImprovement = (shortTermImprovement < 1e-5); // 0.001%
// 3. 检查粒子多样性 - 可能过早聚集或无效分散
double currentDiversity = calculateSwarmDiversity();
bool problematicDiversity = (currentDiversity < m_diversityThreshold * 0.1) ||
(currentDiversity > m_diversityThreshold * 5.0);
// 4. 检查粒子停滞率
double stagnationRate = calculateParticleStagnationRate();
bool mostParticlesStagnant = (stagnationRate > 0.9); // 90%以上停滞但解质量差
// 5. 检查适应度方差 - 可能卡在平坦区域
double fitnessVariance = calculateFitnessVariance(m_localOptimumWindow);
bool flatFitnessLandscape = (fitnessVariance < m_convergenceVarianceThreshold * 0.1);
// 局部最优的判断条件
bool isLocalOptimum = farFromTarget && poorShortTermImprovement &&
(problematicDiversity || (mostParticlesStagnant && flatFitnessLandscape));
if(isLocalOptimum) {
DEBUG_OUT("=== LOCAL OPTIMUM ANALYSIS ===");
DEBUG_OUT(QString("farFromTarget=%1 (error=%2, target*factor=%3)")
.arg(farFromTarget).arg(m_globalBestFitness, 0, 'e', 4)
.arg(m_targetError * m_farTargetFactor, 0, 'e', 4));
DEBUG_OUT(QString("poorShortTermImprovement=%1 (improvement=%2%)")
.arg(poorShortTermImprovement).arg(shortTermImprovement * 100, 0, 'f', 4));
DEBUG_OUT(QString("problematicDiversity=%1 (diversity=%2)")
.arg(problematicDiversity).arg(currentDiversity, 0, 'f', 6));
DEBUG_OUT(QString("stagnationRate=%1%, errorVariance=%2")
.arg(stagnationRate * 100, 0, 'f', 1).arg(fitnessVariance, 0, 'e', 6));
}
return isLocalOptimum;
}
double nmCalculationAutoFitPSO::calculateSwarmDiversity() const
{
// 粒子群多样性:逐维计算粒子位置标准差,再按该参数搜索范围归一化。
// 返回值越小,说明粒子越集中;越大,说明搜索仍然分散。
if(m_swarm.size() < 2) return 0.0;
int dimensions = getEnabledParameterCount();
if(dimensions == 0) return 0.0;
double totalDiversity = 0.0;
for(int dim = 0; dim < dimensions; ++dim) {
// 计算该维度上所有粒子的均值
double mean = 0.0;
for(int i = 0; i < m_swarm.size(); ++i) {
mean += m_swarm[i].position[dim];
}
mean /= m_swarm.size();
// 计算该维度上的方差
double variance = 0.0;
for(int i = 0; i < m_swarm.size(); ++i) {
double diff = m_swarm[i].position[dim] - mean;
variance += diff * diff;
}
variance /= m_swarm.size();
// 归一化到参数范围
int paramIndex = m_enabledParamIndices[dim];
double range = m_parameterUpper[paramIndex] - m_parameterLower[paramIndex];
double normalizedStd = sqrt(variance) / qMax(1e-10, range);
totalDiversity += normalizedStd;
}
return totalDiversity / dimensions;
}
double nmCalculationAutoFitPSO::calculateAverageVelocity() const
{
// 平均速度:逐维速度绝对值按参数范围归一化。
// 收敛阶段速度应逐渐变小。
if(m_swarm.isEmpty()) return 0.0;
int dimensions = getEnabledParameterCount();
if(dimensions == 0) return 0.0;
double totalVelocity = 0.0;
for(int i = 0; i < m_swarm.size(); ++i) {
for(int dim = 0; dim < dimensions; ++dim) {
int paramIndex = m_enabledParamIndices[dim];
double range = m_parameterUpper[paramIndex] - m_parameterLower[paramIndex];
double normalizedVel = qAbs(m_swarm[i].velocity[dim]) / qMax(1e-10, range);
totalVelocity += normalizedVel;
}
}
return totalVelocity / (static_cast<double>(m_swarm.size()) * dimensions);
}
double nmCalculationAutoFitPSO::calculateParticleStagnationRate() const
{
// 粒子停滞率:当前 fitness 与 pbest 很接近时认为该粒子停滞。
// 该指标用于辅助判断真收敛或局部最优。
if(m_swarm.isEmpty()) return 0.0;
int stagnantParticles = 0;
for(int i = 0; i < m_swarm.size(); ++i) {
// 如果当前适应度与个体最优非常接近,认为该粒子停滞
double improvement = (m_swarm[i].bestFitness - m_swarm[i].fitness) /
qMax(1e-10, qAbs(m_swarm[i].bestFitness));
if(qAbs(improvement) < 1e-6) { // 改进小于0.0001%
stagnantParticles++;
}
}
return static_cast<double>(stagnantParticles) / m_swarm.size();
}
double nmCalculationAutoFitPSO::calculateFitnessVariance(int windowSize) const
{
// 最近 windowSize 代全局最优误差的方差。
// 方差越小,说明全局最优曲线越稳定。
if(m_convergenceHistory.size() < windowSize) {
return 1e10; // 数据不足,返回大值
}
// 计算最近windowSize次迭代的方差
double mean = 0.0;
int startIdx = m_convergenceHistory.size() - windowSize;
for(int i = startIdx; i < m_convergenceHistory.size(); ++i) {
mean += m_convergenceHistory[i];
}
mean /= windowSize;
double variance = 0.0;
for(int i = startIdx; i < m_convergenceHistory.size(); ++i) {
double diff = m_convergenceHistory[i] - mean;
variance += diff * diff;
}
variance /= windowSize;
return variance;
}
double nmCalculationAutoFitPSO::calculateLongTermImprovement(int windowSize) const
{
// 最近 windowSize 代相对改进幅度。
// 数据不足时返回 1.0,表示“还不能认为没有改进”。
if(m_convergenceHistory.size() < windowSize) {
return 1.0; // 数据不足,假设有改进
}
double oldFitness = m_convergenceHistory[m_convergenceHistory.size() - windowSize];
double improvement = (oldFitness - m_globalBestFitness) / qMax(1e-10, qAbs(oldFitness));
return improvement;
}
QString nmCalculationAutoFitPSO::getStopReasonDescription(StopReasonPSO reason) const
{
// 将停止枚举转成人类可读文本,用于日志和运行摘要。
switch(reason) {
case PSO_TARGET_ACHIEVED:
return tr("Target error achieved");
case PSO_TRUE_CONVERGENCE:
return tr("Algorithm converged to stable solution");
case PSO_LOCAL_OPTIMUM:
return tr("Local optimum detected");
case PSO_MAX_ITERATIONS:
return tr("Maximum iterations reached");
case PSO_USER_STOPPED:
return tr("Stopped by user request");
case PSO_CONSECUTIVE_FAILURES:
return tr("Too many consecutive failures");
case PSO_OPTIMIZATION_FAILED:
return tr("Optimization failed");
default:
return tr("Unknown reason");
}
}
void nmCalculationAutoFitPSO::updateConvergenceMetrics()
{
// 每代结束后缓存收敛辅助指标,只保留最近 50 个点,避免长时间运行无限增长。
// 更新多样性历史
m_diversityHistory.append(calculateSwarmDiversity());
// 更新速度历史
m_velocityHistory.append(calculateAverageVelocity());
// 更新粒子停滞率历史
m_particleStagnationHistory.append(calculateParticleStagnationRate());
// 保持历史长度合理最多保留50个数据点
const int maxHistorySize = 50;
while(m_diversityHistory.size() > maxHistorySize) {
m_diversityHistory.remove(0);
}
while(m_velocityHistory.size() > maxHistorySize) {
m_velocityHistory.remove(0);
}
while(m_particleStagnationHistory.size() > maxHistorySize) {
m_particleStagnationHistory.remove(0);
}
}
void nmCalculationAutoFitPSO::setSimulationMode(bool enabled)
{
// 仿真模式开关。它用于 UI 演示/调试进度和日志,不调用真实 DLL 求解器。
// 正式自动拟合流程通常保持 false。
m_simulationMode = enabled;
DEBUG_OUT(QString("Simulation mode %1").arg(enabled ? "enabled" : "disabled"));
}
void nmCalculationAutoFitPSO::setSimulationTargetParams(const QVector<double>& targetParams, double targetError)
{
// 设置仿真模式的“目标参数”和“目标误差”。
// 仿真过程会把参数和误差平滑地从初始值过渡到这里。
m_simulationTargetParams = targetParams;
m_simulationTargetError = targetError;
DEBUG_OUT(QString("Simulation target params set: %1 parameters, target error: %2")
.arg(targetParams.size()).arg(targetError, 0, 'e', 4));
}
QVector<QVector<double> > nmCalculationAutoFitPSO::buildSimulatedLogLogData(double progress) const
{
// 为仿真模式生成一条看起来逐渐贴近目标曲线的 log-log 曲线。
// progress=0 表示偏离目标较多progress=1 基本等于目标曲线。
QVector<QVector<double> > simulated;
if(m_targetLogLogData.size() < 3) {
return simulated;
}
int count = qMin(m_targetLogLogData[0].size(), qMin(m_targetLogLogData[1].size(), m_targetLogLogData[2].size()));
simulated.resize(3);
simulated[0].reserve(count);
simulated[1].reserve(count);
simulated[2].reserve(count);
double clampedProgress = qMax(0.0, qMin(1.0, progress));
double pressureOffset = 0.35 * (1.0 - clampedProgress);
double derivativeOffset = 0.28 * (1.0 - clampedProgress);
for(int i = 0; i < count; ++i) {
double x = m_targetLogLogData[0][i];
double y1 = m_targetLogLogData[1][i];
double y2 = m_targetLogLogData[2][i];
if(!isFiniteNumber(x) || !isFiniteNumber(y1) || !isFiniteNumber(y2) || x <= 0.0 || y1 <= 0.0 || y2 <= 0.0) {
continue;
}
double localWave = 1.0 + 0.04 * (1.0 - clampedProgress) * qSin(i * 0.17);
simulated[0].append(x);
simulated[1].append(qMax(1e-12, y1 * (1.0 + pressureOffset) * localWave));
simulated[2].append(qMax(1e-12, y2 * (1.0 - derivativeOffset) / localWave));
}
return simulated;
}
void nmCalculationAutoFitPSO::runSimulatedFitting()
{
// 仿真模式入口:模拟完整 PSO 日志、进度、bestCurveUpdated 信号和最终写回。
// 它不评价粒子、不调用代理模型、不调用 DLL只用于界面联调或演示自动拟合流程。
DEBUG_OUT("=== SIMULATED PSO AUTO FITTING START ===");
m_isRunning = true;
m_shouldStop = false;
m_simulationIteration = 0;
m_simulationMaxIterations = m_maxIterations;
m_simulationCurrentParticle = 0;
m_simulationSuccessCount = 0;
m_simulationFailCount = 0;
m_simulationIterationStarted = false;
// 设置起始误差
m_simulationStartError = 0.7324;
m_simulationCurrentError = m_simulationStartError;
m_simulationPreviousIterError = m_simulationStartError;
int enabledParams = getEnabledParameterCount();
// 从界面提取的初始值
m_simulationStartParams = m_initialValues;
m_globalBestPosition = m_simulationStartParams;
m_globalBestFitness = m_simulationStartError;
m_globalBestLogLogData = buildSimulatedLogLogData(0.0);
emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, 0, m_globalBestFitness);
// 发送初始化日志
emit logMessageGenerated(tr("=== PSO Automatic Fitting Started ==="));
emit logMessageGenerated(tr("Algorithm: Particle Swarm Optimization"));
emit logMessageGenerated(tr("Enabled parameters count: %1").arg(enabledParams));
emit logMessageGenerated(tr("Target data validation passed (%1 data points)")
.arg(m_targetLogLogData.isEmpty() ? 156 : m_targetLogLogData[0].size()));
emit logMessageGenerated(tr("=== Evaluating Initial Solution (Elite Protection) ==="));
QString paramStr = tr("Initial parameters: ");
for(int i = 0; i < m_simulationStartParams.size(); ++i) {
paramStr += QString("[%1]=%2 ").arg(i).arg(m_simulationStartParams[i], 0, 'f', 6);
}
emit logMessageGenerated(paramStr);
emit logMessageGenerated(tr("Starting initial solution evaluation..."));
QApplication::processEvents();
msleep(500);
emit logMessageGenerated(tr("evaluateError returned: %1").arg(m_simulationStartError, 0, 'e', 10));
emit logMessageGenerated(tr("Taking SUCCESS branch (Error < 1e9)"));
emit logMessageGenerated(tr("Initial solution evaluation successful"));
emit logMessageGenerated(tr("Initial Error: %1").arg(m_simulationStartError, 0, 'e', 4));
emit logMessageGenerated(tr("Elite protection activated - initial solution will be preserved if no significant improvement found"));
emit logMessageGenerated(tr("Swarm initialized: %1 particles, %2 dimensions").arg(m_swarmSize).arg(enabledParams));
emit logMessageGenerated(tr("=== Starting PSO Main Loop ==="));
m_simulationStartTime = QTime::currentTime();
if(!m_simulationTimer) {
m_simulationTimer = new QTimer(this);
connect(m_simulationTimer, SIGNAL(timeout()), this, SLOT(onSimulationTimerTick()));
}
int intervalMs = m_simulationLogInterval * 1000;
m_simulationTimer->start(intervalMs);
DEBUG_OUT(QString("Simulation timer started with interval: %1 ms").arg(intervalMs));
}
void nmCalculationAutoFitPSO::onSimulationTimerTick()
{
// 仿真模式用 QTimer 分步输出日志,避免一次性刷完,效果更接近真实长任务。
// 检查是否需要停止
if(m_shouldStop) {
m_simulationTimer->stop();
emit logMessageGenerated(tr("=== User Stop Request Received ==="));
emit logMessageGenerated(tr("Gracefully stopping PSO optimization..."));
emit logMessageGenerated(tr("PSO optimization stop request processed"));
m_globalBestPosition = calculateSimulatedParams(m_simulationIteration);
m_globalBestFitness = m_simulationCurrentError;
m_isRunning = false;
QString message = QString(tr("Stopped by user. Best error: %1, Iterations: %2"))
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_simulationIteration);
emit logMessageGenerated(tr("=== PSO OPTIMIZATION STOPPED BY USER ==="));
emit progressUpdated(m_simulationMaxIterations, m_globalBestFitness);
QApplication::processEvents();
msleep(200);
emit fittingFinished(true, message);
return;
}
// 检查是否完成所有迭代
if(m_simulationIteration >= m_simulationMaxIterations) {
m_simulationTimer->stop();
finishSimulation();
return;
}
// 检查是否需要开始新的迭代
if(!m_simulationIterationStarted) {
startNewSimulationIteration();
return;
}
// 处理当前粒子
if(m_simulationCurrentParticle < m_swarmSize) {
emitParticleLog(m_simulationCurrentParticle);
m_simulationCurrentParticle++;
// 每处理2个粒子处理一次事件
if(m_simulationCurrentParticle % 2 == 0) {
QApplication::processEvents();
}
} else {
// 当前迭代的所有粒子都处理完了,输出迭代统计
finishCurrentIteration();
}
}
void nmCalculationAutoFitPSO::emitSimulationLog(int iteration)
{
// 旧版仿真日志函数:一次性输出某一代的模拟粒子情况。
// 当前 runSimulatedFitting() 主要通过 startNewSimulationIteration()/emitParticleLog()
// 分步输出,但保留本函数便于调试或恢复旧演示节奏。
m_simulationCurrentError = calculateSimulatedError(iteration);
if(m_simulationCurrentError < m_globalBestFitness) {
QVector<double> currentParams = calculateSimulatedParams(iteration);
m_globalBestFitness = m_simulationCurrentError;
m_globalBestPosition = currentParams;
}
bool shouldOutputDetail = (iteration % qMax(1, m_simulationMaxIterations / 10) == 0) ||
(iteration <= 5) ||
(iteration >= m_simulationMaxIterations - 2);
emit logMessageGenerated(tr("--- Iteration %1/%2 ---").arg(iteration).arg(m_simulationMaxIterations));
if(shouldOutputDetail) {
emit logMessageGenerated(tr("Current best error: %1").arg(m_globalBestFitness, 0, 'e', 4));
emit logMessageGenerated(tr("Total evaluations: %1 (successful: %2, failures: %3)")
.arg(iteration * m_swarmSize)
.arg(iteration * m_swarmSize)
.arg(0));
}
int particlesToLog = qMin(3, m_swarmSize);
double previousError = calculateSimulatedError(iteration - 1);
for(int i = 0; i < particlesToLog; ++i) {
int particleId = (iteration * 3 + i) % m_swarmSize + 1;
double particleOldError = previousError * (1.0 + 0.05 * (qrand() % 100) / 100.0);
double particleNewError = m_simulationCurrentError * (1.0 + 0.03 * (qrand() % 100) / 100.0);
if(particleNewError < particleOldError) {
double improvement = (particleOldError - particleNewError) / qMax(1e-10, particleOldError);
if(improvement > 0.01) {
emit logMessageGenerated(tr(" Particle %1 improved: %2 -> %3 (%4% improvement)")
.arg(particleId)
.arg(particleOldError, 0, 'e', 3)
.arg(particleNewError, 0, 'e', 3)
.arg(improvement * 100, 0, 'f', 2));
} else {
emit logMessageGenerated(tr(" Particle %1: minor improvement %2% (< 1% threshold)")
.arg(particleId).arg(improvement * 100, 0, 'f', 3));
}
} else {
emit logMessageGenerated(tr(" Particle %1: no improvement (Error: %2)")
.arg(particleId).arg(particleNewError, 0, 'e', 3));
}
}
int successCount = m_swarmSize - (qrand() % 2);
int failCount = m_swarmSize - successCount;
double successRate = (double)successCount / m_swarmSize * 100;
if(shouldOutputDetail) {
emit logMessageGenerated(tr("Current iteration stats: %1 successful, %2 failed out of %3 particles (success rate: %4%)")
.arg(successCount).arg(failCount).arg(m_swarmSize).arg(successRate, 0, 'f', 1));
}
if(iteration % 5 == 0 && iteration > 0) {
double recentImprovement = (previousError - m_simulationCurrentError) / qMax(1e-10, previousError);
emit logMessageGenerated(tr(" Convergence check: recent improvement = %1% over last 5 iterations")
.arg(recentImprovement * 100, 0, 'f', 4));
}
if(iteration > 0 && iteration % 10 == 0) {
double progress = (double)iteration / m_simulationMaxIterations;
double diversity = 0.15 * (1.0 - progress * 0.8) + 0.01;
double avgVel = 0.08 * (1.0 - progress * 0.9) + 0.005;
double stagnation = 20.0 + progress * 60.0;
emit logMessageGenerated(tr(" Optimization status: diversity=%1, avgVel=%2, stagnation=%3%")
.arg(diversity, 0, 'f', 4)
.arg(avgVel, 0, 'f', 4)
.arg(stagnation, 0, 'f', 1));
}
if(shouldOutputDetail) {
emit logMessageGenerated(tr(" Iteration %1 completed - Current best: %2")
.arg(iteration).arg(m_globalBestFitness, 0, 'e', 4));
}
}
double nmCalculationAutoFitPSO::calculateSimulatedError(int iteration)
{
// 生成随迭代递减的模拟误差,并叠加少量随机扰动。
// 使用非线性 decay让前期下降较快、后期逐渐变慢更像真实优化曲线。
if(iteration <= 0) return m_simulationStartError;
if(iteration >= m_simulationMaxIterations) return m_simulationTargetError;
double progress = (double)iteration / m_simulationMaxIterations;
double decayFactor = 1.0 - pow(progress, 0.6);
double error = m_simulationTargetError +
(m_simulationStartError - m_simulationTargetError) * decayFactor;
double noise = (qrand() % 1000 - 500) / 500000.0;
error = error * (1.0 + noise);
error = qMax(m_simulationTargetError, error);
return error;
}
QVector<double> nmCalculationAutoFitPSO::calculateSimulatedParams(int iteration)
{
// 生成随迭代从初始参数平滑靠近目标参数的模拟参数。
// 返回向量仍然是“启用参数向量”,可直接传给 applyParametersToDataManager()。
QVector<double> currentParams;
if(iteration <= 0) return m_simulationStartParams;
if(iteration >= m_simulationMaxIterations) return m_simulationTargetParams;
double progress = (double)iteration / m_simulationMaxIterations;
double blendFactor = 1.0 - pow(1.0 - progress, 0.6);
for(int i = 0; i < m_simulationStartParams.size() && i < m_simulationTargetParams.size(); ++i) {
double startVal = m_simulationStartParams[i];
double targetVal = m_simulationTargetParams[i];
double currentVal = startVal + (targetVal - startVal) * blendFactor;
double noise = (qrand() % 1000 - 500) / 50000.0;
currentVal = currentVal * (1.0 + noise);
currentParams.append(currentVal);
}
return currentParams;
}
void nmCalculationAutoFitPSO::startNewSimulationIteration()
{
// 开始仿真模式的一代:更新误差/最优曲线,输出该代标题和累计评价数。
m_simulationIteration++;
m_simulationCurrentParticle = 0;
m_simulationSuccessCount = 0;
m_simulationFailCount = 0;
m_simulationIterationStarted = true;
// 保存上一次迭代的误差
m_simulationPreviousIterError = m_simulationCurrentError;
// 计算当前迭代的目标误差
m_simulationCurrentError = calculateSimulatedError(m_simulationIteration);
if(m_simulationCurrentError < m_globalBestFitness) {
QVector<double> currentParams = calculateSimulatedParams(m_simulationIteration);
m_globalBestFitness = m_simulationCurrentError;
m_globalBestPosition = currentParams;
double progress = (double)m_simulationIteration / qMax(1, m_simulationMaxIterations);
m_globalBestLogLogData = buildSimulatedLogLogData(progress);
emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, m_simulationIteration, m_globalBestFitness);
}
bool shouldOutputDetail = (m_simulationIteration % qMax(1, m_simulationMaxIterations / 10) == 0) ||
(m_simulationIteration <= 5) ||
(m_simulationIteration >= m_simulationMaxIterations - 2);
// 输出迭代标题
emit logMessageGenerated(tr("--- Iteration %1/%2 ---").arg(m_simulationIteration).arg(m_simulationMaxIterations));
if(shouldOutputDetail) {
emit logMessageGenerated(tr("Current best error: %1").arg(m_globalBestFitness, 0, 'e', 4));
emit logMessageGenerated(tr("Total evaluations: %1 (successful: %2, failures: %3)")
.arg((m_simulationIteration - 1) * m_swarmSize)
.arg((m_simulationIteration - 1) * m_swarmSize)
.arg(0));
}
}
void nmCalculationAutoFitPSO::emitParticleLog(int particleIndex)
{
// 输出单个仿真粒子的日志。这里随机生成失败、显著改进、微小改进、无改进四类情况,
// 目的是让 UI 日志表现接近真实 PSO而不是用于任何真实优化计算。
int particleId = particleIndex + 1;
// 为每个粒子生成模拟的误差值
double particleNewError = m_simulationCurrentError * (0.97 + 0.06 * (qrand() % 100) / 100.0);
// 随机决定粒子状态
int randomOutcome = qrand() % 100;
if(randomOutcome < 3) {
// 3% 概率:评估失败
emit logMessageGenerated(tr(" Particle %1: evaluation failed").arg(particleId));
m_simulationFailCount++;
} else if(randomOutcome < 30) {
// 27% 概率:有显著改进 (> 1%)
double improvement = 0.01 + 0.12 * (qrand() % 100) / 100.0;
double oldErr = particleNewError / (1.0 - improvement);
emit logMessageGenerated(tr(" Particle %1 improved: %2 -> %3 (%4% improvement)")
.arg(particleId)
.arg(oldErr, 0, 'e', 3)
.arg(particleNewError, 0, 'e', 3)
.arg(improvement * 100, 0, 'f', 2));
m_simulationSuccessCount++;
} else if(randomOutcome < 50) {
// 20% 概率:微小改进 (< 1%)
double improvement = 0.001 + 0.008 * (qrand() % 100) / 100.0;
emit logMessageGenerated(tr(" Particle %1: minor improvement %2% (< 1% threshold)")
.arg(particleId)
.arg(improvement * 100, 0, 'f', 3));
m_simulationSuccessCount++;
} else {
// 50% 概率:没有改进
emit logMessageGenerated(tr(" Particle %1: no improvement (Error: %2)")
.arg(particleId)
.arg(particleNewError, 0, 'e', 3));
m_simulationSuccessCount++;
}
}
void nmCalculationAutoFitPSO::finishCurrentIteration()
{
// 完成仿真模式的一代:输出统计、收敛检查和进度更新,然后等待下一次 timer tick。
m_simulationIterationStarted = false;
bool shouldOutputDetail = (m_simulationIteration % qMax(1, m_simulationMaxIterations / 10) == 0) ||
(m_simulationIteration <= 5) ||
(m_simulationIteration >= m_simulationMaxIterations - 2);
// 每5次迭代输出"发现更好解"
if(m_simulationIteration % 5 == 0 && m_simulationIteration > 0) {
emit logMessageGenerated(tr("Better solution found: %1").arg(m_globalBestFitness, 0, 'e', 4));
}
// 输出当前迭代统计
double successRate = (double)m_simulationSuccessCount / m_swarmSize * 100;
if(shouldOutputDetail) {
emit logMessageGenerated(tr("Current iteration stats: %1 successful, %2 failed out of %3 particles (success rate: %4%)")
.arg(m_simulationSuccessCount).arg(m_simulationFailCount)
.arg(m_swarmSize).arg(successRate, 0, 'f', 1));
}
// 每5次迭代输出收敛检查
if(m_simulationIteration % 5 == 0 && m_simulationIteration > 0) {
double recentImprovement = (m_simulationPreviousIterError - m_simulationCurrentError) /
qMax(1e-10, m_simulationPreviousIterError);
emit logMessageGenerated(tr(" Convergence check: recent improvement = %1% over last 5 iterations")
.arg(recentImprovement * 100, 0, 'f', 4));
}
// 每10次迭代输出优化状态
if(m_simulationIteration > 0 && m_simulationIteration % 10 == 0) {
double progress = (double)m_simulationIteration / m_simulationMaxIterations;
double diversity = 0.15 * (1.0 - progress * 0.8) + 0.01;
double avgVel = 0.08 * (1.0 - progress * 0.9) + 0.005;
double stagnation = 20.0 + progress * 60.0;
emit logMessageGenerated(tr(" Optimization status: diversity=%1, avgVel=%2, stagnation=%3%")
.arg(diversity, 0, 'f', 4)
.arg(avgVel, 0, 'f', 4)
.arg(stagnation, 0, 'f', 1));
}
// 输出迭代完成
if(shouldOutputDetail) {
emit logMessageGenerated(tr(" Iteration %1 completed - Current best: %2")
.arg(m_simulationIteration).arg(m_globalBestFitness, 0, 'e', 4));
}
// 更新进度
emit progressUpdated(m_simulationIteration, m_simulationCurrentError);
QApplication::processEvents();
}
void nmCalculationAutoFitPSO::finishSimulation()
{
// 仿真模式正常结束:直接采用预设目标参数/目标误差作为最终结果,
// 发送与真实拟合相同风格的日志、进度和 fittingFinished 信号。
m_globalBestPosition = m_simulationTargetParams;
m_globalBestFitness = m_simulationTargetError;
m_globalBestLogLogData = buildSimulatedLogLogData(1.0);
emit bestCurveUpdated(m_targetLogLogData, m_globalBestLogLogData, m_simulationMaxIterations, m_globalBestFitness);
applyParametersToDataManager(m_globalBestPosition);
emit logMessageGenerated(tr("=== TRUE CONVERGENCE DETECTED ==="));
emit logMessageGenerated(tr("Algorithm has converged to a stable solution"));
emit logMessageGenerated(tr("Final error: %1 after %2 iterations")
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_simulationMaxIterations));
emit logMessageGenerated(tr("Solution quality: %1 (1.0 = target achieved)")
.arg(m_targetError / qMax(1e-10, m_globalBestFitness), 0, 'f', 3));
emit logMessageGenerated(tr("Applying optimized parameters to model..."));
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_simulationMaxIterations));
emit logMessageGenerated(tr("Total evaluations: %1 (successful: %2)")
.arg(m_simulationMaxIterations * m_swarmSize)
.arg(m_simulationMaxIterations * m_swarmSize));
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);
emit logMessageGenerated(tr("Parameters applied successfully to data manager"));
emit logMessageGenerated(tr("=== PSO OPTIMIZATION CONVERGED ==="));
m_isRunning = false;
emit progressUpdated(m_simulationMaxIterations, m_globalBestFitness);
QApplication::processEvents();
msleep(200);
QString message = QString(tr("PSO optimization converged to stable solution. Best error: %1, Iterations: %2"))
.arg(m_globalBestFitness, 0, 'e', 4).arg(m_simulationMaxIterations);
emit fittingFinished(true, message);
}