#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 #include #include #include #include #include #include #include #include #include #include #include #include #include #ifdef Q_OS_WIN #include #include #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.0; // 相邻生产段产量最大跳变比。 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& 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& 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& 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& 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(generation + 1) + (x << 6) + (x >> 2); x ^= 0x85ebca6bu + static_cast(particleIndex + 1) + (x << 6) + (x >> 2); x ^= x >> 16; x *= 0x7feb352du; x ^= x >> 15; x *= 0x846ca68bu; x ^= x >> 16; return static_cast(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& 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__ 目录。退出时只删除本类创建的目录,启动时顺便清理 // 旧进程遗留的 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 >& targetData) { // 目标曲线由界面层从目标井 history log-log 传入。 // 约定 targetData[0]=time,targetData[1]=pressure,targetData[2]=pressure derivative。 m_targetLogLogData = targetData; DEBUG_OUT(QString("Target LogLog data set: %1 arrays").arg(targetData.size())); if(targetData.size() >= 3) { DEBUG_OUT(QString("LogLog data points: X=%1, Y1=%2, Y2=%3") .arg(targetData[0].size()) .arg(targetData[1].size()) .arg(targetData[2].size())); } } void 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 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_.csv:逐代逐粒子的参数、真实误差、代理误差和筛选决策; // - pso_baseline_trace_.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 initialFullParams = buildTraceParameterVector(m_initialValues); QVector targetTime = m_targetLogLogData.size() > 0 ? m_targetLogLogData[0] : QVector(); QVector targetPressure = m_targetLogLogData.size() > 1 ? m_targetLogLogData[1] : QVector(); QVector targetDerivative = m_targetLogLogData.size() > 2 ? m_targetLogLogData[2] : QVector(); QVector flowPoints = pTargetWell ? pTargetWell->getFlowPoints() : QVector(); 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(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 nmCalculationAutoFitPSO::buildTraceParameterVector(const QVector& selectedParameters) const { // 将粒子内部使用的“启用参数向量”还原成完整 11 维参数向量。 // 未启用的参数从当前 DataManager 读取,启用的参数用 selectedParameters 覆盖。 // trace CSV、候选 CSV、代理训练域检查都需要这个完整向量。 QVector 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& parameters, double solverObjective, bool solverSuccess, int elapsedMs, double surrogateObjective, const QString& screeningDecision, const QVector& 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 currentParams = buildTraceParameterVector(parameters); QVector pbestParams = pbestPosition.isEmpty() ? QVector() : buildTraceParameterVector(pbestPosition); QVector gbestParams = m_globalBestPosition.isEmpty() ? QVector() : 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::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__gen.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 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("") : stdOut) .arg(stdErr.isEmpty() ? QString("") : 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("") : 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("") : stdOut) .arg(stdErr.isEmpty() ? QString("") : 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) { // 一次性代理评分脚本兜底路径。 // // 输入: // - candidates:writeSurrogateCandidateCsv() 写出的粒子参数; // - trace-meta:目标曲线、流量制度、参数上下界等上下文; // - tag/stage:选择 ML 侧模型。 // // 输出: // - scorePath:Python 写出的评分 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("") : stdOut) .arg(stdErr.isEmpty() ? QString("") : 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("") : stdOut) .arg(stdErr.isEmpty() ? QString("") : 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("") : stdOut) .arg(stdErr.isEmpty() ? QString("") : 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 nmCalculationAutoFitPSO::readSurrogateScores(const QString& scorePath) const { // 读取 Python 输出的评分 CSV。 // // 当前脚本约定至少包含: // particle_id,surrogate_objective,...,success // 本函数只读取 particle_id、score 和 success 标志。失败或缺失的粒子保持 NaN, // 后续 buildSurrogateEvaluationMask() 会通过 fallback_gate 强制跑真实求解器。 QVector scores(m_swarm.size(), std::numeric_limits::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 > 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 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 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::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& selectedParameters, QString* reason) const { // 判断“单个粒子候选参数”是否在代理模型训练域内。 // 与 isSurrogateRunContextSupported() 不同,这里检查的是粒子参数值本身。 // 超出范围的粒子会被 domain_gate 强制跑真实求解器。 QVector params = buildTraceParameterVector(selectedParameters); double k = params.size() > 0 ? params[0] : std::numeric_limits::quiet_NaN(); double skin = params.size() > 1 ? params[1] : std::numeric_limits::quiet_NaN(); double wellboreC = params.size() > 2 ? params[2] : std::numeric_limits::quiet_NaN(); double phi = params.size() > 3 ? params[3] : std::numeric_limits::quiet_NaN(); double h = params.size() > 5 ? params[5] : std::numeric_limits::quiet_NaN(); double cf = params.size() > 7 ? params[7] : std::numeric_limits::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& selectedParameters) const { // 经验 fallback gate。 // 有些点虽然还在训练域边界内,但组合起来属于代理模型容易排序失真的区域, // 例如极低渗透率配大井筒储集、极端 skin 等。遇到这些候选时直接跑真实求解器。 QVector params = buildTraceParameterVector(selectedParameters); double k = params.size() > 0 ? params[0] : std::numeric_limits::quiet_NaN(); double skin = params.size() > 1 ? params[1] : std::numeric_limits::quiet_NaN(); double wellboreC = params.size() > 2 ? params[2] : std::numeric_limits::quiet_NaN(); double phi = params.size() > 3 ? params[3] : std::numeric_limits::quiet_NaN(); double h = params.size() > 5 ? params[5] : std::numeric_limits::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::quiet_NaN(); } nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance(); if(!dataManager) { return false; } nmDataWellBase* pTargetWell = dataManager->findWellByName(m_targetWellName); if(!pTargetWell) { return false; } QVector 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 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 evaluateMask(m_swarm.size(), true); for(int i = 0; i < m_swarm.size(); ++i) { m_swarm[i].surrogateObjective = std::numeric_limits::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 scores = readSurrogateScores(scorePath); QVector 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 > 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& a, const QPair& 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 > 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& a, const QPair& 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 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(QDateTime::currentDateTime().toTime_t()); m_psoRandomSeed ^= static_cast((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::quiet_NaN(), "initial_solution", m_hasValidUserSolution ? m_userInitialSolution : QVector(), 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 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 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& 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 > 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 > 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& 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> 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> 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& parameters) { // 将粒子的“启用参数向量”写回项目数据。 // parameters 的维度必须等于用户勾选的参数数量,顺序由 m_enabledParamIndices 决定。 // 这里不直接跑求解器,只负责把 DataManager 调整到该粒子对应的模型状态。 if(parameters.size() != getEnabledParameterCount()) { return; } updateReservoirParameters(parameters); updateWellParameters(parameters); } void nmCalculationAutoFitPSO::updateReservoirParameters(const QVector& 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& parameters) { // 更新目标井上的拟合参数。目前井级可拟合参数主要是: // - skin:写入第一个 perforation; // - wellboreC:写入井筒储集系数。 // 如果目标井不存在或没有射孔数据,这里只记录 debug,不抛异常。 nmDataAnalyzeManager* dataManager = nmDataAnalyzeManager::getCurrentInstance(); // 只更新当前目标井的井参数,避免多井项目中误改其他井。 //QVector wells = dataManager->getWellDataList(); nmDataWellBase* pWell = dataManager->findWellByName(m_targetWellName); if(!pWell) return; //nmDataWellBase* pWell = wells[0]; // 使用第一口井 int paramIndex = 0; for(int i = 0; i < m_parameterSelected.size() && i < 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(pWell); if(pVerticalWell) { QVector wells; wells.append(*pVerticalWell); dataManager->updateVerticalWells(wells); } break; } case NM_WELL_MODEL::Vertical_Fractured_Well: { nmDataVerticalFracturedWell* pVFracturedWell = dynamic_cast(pWell); if(pVFracturedWell) { QVector wells; wells.append(*pVFracturedWell); dataManager->updateVerticalFracturedWells(wells); } break; } case NM_WELL_MODEL::Horizontal_Fractured_Well: { nmDataHorizontalFracturedWell* pHFracturedWell = dynamic_cast(pWell); if(pHFracturedWell) { QVector wells; wells.append(*pHFracturedWell); dataManager->updateHorizontalFracturedWells(wells); } break; } default: break; } } // ==================== 求解器相关方法 ==================== QVector> nmCalculationAutoFitPSO::runSolver() { // 真实求解器入口。目前默认走 DLL 方式,保留 runSolverExe() 是历史/备用路径。 // evaluateFitness() 不关心具体求解器形态,只要求返回可验证的双对数结果。 //return runSolverExe(); return runSolverDll(); } // ==================== 数据处理方法 ==================== QVector nmCalculationAutoFitPSO::interpolateData( const QVector& source, const QVector& targetX) const { // 将一条曲线插值到指定 targetX 网格。 // 双对数误差计算前会分别把目标曲线和模拟曲线插到同一个 commonX 上, // 这样每个时间点才能逐点比较。 QVector result; if(source.isEmpty() || targetX.isEmpty()) { DEBUG_OUT("Warning: Empty data in interpolation"); return result; } // 数据清理和验证合并:跳过 NaN/Inf 和极端大值,尽量保留可用点继续计算。 // 这里不立即把整个粒子判失败,是为了提高求解器偶发异常输出下的鲁棒性。 QVector 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 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(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(qrand()) / RAND_MAX; } void nmCalculationAutoFitPSO::clampToLimits(QVector& 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& 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>& logLogData) const { // 校验双对数曲线结构。约定: // logLogData[0]=time,logLogData[1]=pressure,logLogData[2]=pressure derivative。 // 三列必须长度一致,且至少有足够点数用于插值和误差计算。 if(logLogData.size() < 3) { DEBUG_OUT("LogLog data has less than 3 arrays"); return false; } // 检查数组大小一致性 int size = logLogData[0].size(); if(size == 0) { DEBUG_OUT("Empty LogLog data"); return false; } if(logLogData[1].size() != size || logLogData[2].size() != size) { DEBUG_OUT(QString("LogLog data size mismatch: X=%1, Y1=%2, Y2=%3") .arg(logLogData[0].size()) .arg(logLogData[1].size()) .arg(logLogData[2].size())); return false; } // 检查最小数据点数 if(size < 5) { DEBUG_OUT(QString("Too few LogLog data points: %1").arg(size)); return false; } // 数据有效性检查 for(int i = 0; i < size; ++i) { if(!isFiniteNumber(logLogData[0][i]) || !isFiniteNumber(logLogData[1][i]) || !isFiniteNumber(logLogData[2][i])) { DEBUG_OUT(QString("Invalid LogLog data at index %1").arg(i)); return false; } } return true; } bool 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>& 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& curve1, const QVector& 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 >& target, const QVector >& 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 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 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 alignedTarget1 = interpolateData(targetCurve1, commonX); QVector alignedTarget2 = interpolateData(targetCurve2, commonX); // 插值结果曲线。result 与 target 使用同一 commonX,保证逐点可比。 QVector 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 alignedResult1 = interpolateData(resultCurve1, commonX); QVector 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> 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>(); } ++m_evaluationInProgress; QVector> 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 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> pressureResult = pTargetWell->getResultPressure(); QVector> 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 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> nmCalculationAutoFitPSO::runSolverExe() //{ // DEBUG_OUT("SOLVER EXE START"); // // if(m_evaluationInProgress > 0) { // DEBUG_OUT("EXE Solver already running, skipping"); // return QVector>(); // } // // ++m_evaluationInProgress; // QVector> 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 wells = dataManager->getWellDataList(); // // if(wells.isEmpty()) { // DEBUG_OUT("No wells found in data manager after EXE execution"); // delete exeTask; // --m_evaluationInProgress; // return result; // } // // // 验证结果数据的时间戳或唯一性 // QVector> pressureResult = wells[0]->getResultPressure(); // QVector> 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 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(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(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& 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 > nmCalculationAutoFitPSO::buildSimulatedLogLogData(double progress) const { // 为仿真模式生成一条看起来逐渐贴近目标曲线的 log-log 曲线。 // progress=0 表示偏离目标较多,progress=1 基本等于目标曲线。 QVector > 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 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 nmCalculationAutoFitPSO::calculateSimulatedParams(int iteration) { // 生成随迭代从初始参数平滑靠近目标参数的模拟参数。 // 返回向量仍然是“启用参数向量”,可直接传给 applyParametersToDataManager()。 QVector 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 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); }