基于集合卡尔曼滤波与贝叶斯优化的蒸汽发生器仿真研究

研读链接：链接

一、背景

1.1 核电发展的宏观背景

论文开篇从能源安全和“双碳”目标切入。中国在2020年提出“碳达峰、碳中和”目标，核能作为一种高效、稳定、低碳的能源形式，在能源结构中占据重要地位。与风能、太阳能相比，核能具有更高的能量密度和更稳定的发电能力，能够弥补可再生能源间歇性发电的不足。

1.2 蒸汽发生器的关键作用

压水堆核电站的核心结构包括核岛和常规岛。蒸汽发生器是一、二回路之间的关键能量交换设备：

一回路：高温高压冷却剂从堆芯吸收裂变能量，流入蒸汽发生器，通过换热管将热量传递给二回路工质；

二回路：工质吸收热量后相变产生蒸汽，推动汽轮机做功，乏汽冷凝后重新回到蒸汽发生器。

据统计，接近25%的压水堆紧急停堆由蒸汽发生器水位控制不良引起。这说明蒸汽发生器的安全、高效运行直接关系到核电站的安全性和经济性。

2.3 数字孪生与参数优化的需求

近年来，核电数字孪生成为研究热点。数字孪生通过将虚拟模型与真实物理对象结合，实现实体精确监测、仿真优化和预测分析。参数优化是数字孪生技术的重要组成部分，但蒸汽发生器模型参数众多、非线性强，传统人工调参效率低下，牛顿法、梯度下降法等又难以获得满意结果。

二、环境条件：研究的基础平台

2.1 研究对象

对象：大亚湾核电站自然循环式U型管蒸汽发生器

结构参数：传热管壁厚1.09 \(mm\)，外径19.05 \(mm\)，管束高度11.5$ m\(，管束数目4474根，换热面积5429\) m^2$

工况水平：30%、50%、70%、100%功率工况

2.2 建模工具

语言：\(Modelica\)（面向对象的物理建模语言）

方法：集总参数法 + 模块化建模

两相流模型：均相流模型

2.3 验证数据

以大亚湾核电站30%、50%、70%、100%工况的稳态运行数据作为验证数据，包括一次侧流量、入口温度、压力，以及二次侧给水流量、给水温度、压力、循环倍率等。

2.4 计算环境

论文未详细说明硬件配置，但从仿真规模来看（\(EnKF\) 100组样本、BO 300次迭代），属于中等计算量，普通工作站即可完成。

三、处理方式：从建模到优化的完整流程

3.1 建模思路：模块化 + 集总参数

蒸汽发生器内部结构复杂，直接建模难度大。论文采用模块化建模思想，将蒸汽发生器划分为多个控制体：

蒸汽室，给水室，下降通道，预热段，沸腾段，上升段，汽水分离器，换热管壁面，一次侧模块

每个控制体基于集总参数法建立质量守恒、能量守恒和压降方程。集总参数法的核心假设是：参数仅是时间的函数，在空间上集中于一点。

3.2 关键假设

论文做了8条简化假设，其中比较重要的有：

1.一、二次侧工质状态参量在同一横截面上相同，忽略径向变化；

2.一次侧工质密度为定值；

3.忽略壁面轴向导热；

4.沸腾段含汽率沿轴向线性分布；

5.汽水分离器分离效率100%，蒸汽室为干饱和蒸汽。

3.3 换热与阻力模型

一次侧： 迪图斯-贝尔特公式\(（Dittus-Boelter）\)

预热段： 单相对流换热，仍用D-B公式

沸腾段： Rohsenow公式或Chen公式

阻力： 单相用范宁公式，两相用均相流模型

3.4 自然循环流量计算

自然循环的驱动压头来自下降通道与加热段工质的密度差：

\[S_{\text{drive}} = \Delta P_f + \Delta P_{\text{loc}} \]

循环流量：

\[G_{\text{cir}} = \sqrt{\frac{S_{\text{drive}}}{\sum \Delta P}} \]

四、算法内容：两大核心方法

4.1 集合卡尔曼滤波（EnKF）

4.1.1 为什么用EnKF？

传统卡尔曼滤波要求系统线性、误差高斯分布，但蒸汽发生器模型是强非线性系统。扩展卡尔曼滤波（EKF）通过泰勒展开线性化，但在强非线性下仍有局限。EnKF用一组样本集合来近似状态分布，避免了显式求解状态转移矩阵，更适合非线性系统。

4.1.2 核心步骤

预测：

\[x_f^i = A(x_0^i, u)^T \]

\[\bar{x}_f = \frac{1}{N}\sum_{i=1}^{N} x_f^i \]

更新：

\[P = \frac{1}{N-1}\sum_{i=1}^{N}(x_f^i - \bar{x}_f)(x_f^i - \bar{x}_f)^T \]

\[K = PH^T(HPH^T + R)^{-1} \]

\[x_a^i = x_f^i + K(y + w^i - Hx_f^i) \]

关键点：状态向量由汽室压力、沸腾段出口质量含汽率和待同化的修正参数组合\(\beta\)组成：

\[\beta^i = (f_{dc}^i, f_r^i, k_p^i, k_w^i)^T \]

5.2 贝叶斯优化（BO）

5.2.1 为什么用BO？

贝叶斯优化适合目标函数表达式未知、非凸、多峰、评估成本高的优化问题。它通过代理模型拟合目标函数，用采集函数指导下一步采样，只需较少次评估就能接近最优解。

5.2.2 核心组件

  • 代理模型：高斯过程回归（GPR）或树状结构评估器（TPE）
  • 采集函数：期望改进（EI）

高斯过程的核心公式：

\[f(x) \in GP(\mu(x), k(x, x')) \]

\[k(x_i, x_j) = \exp\left(-\frac{1}{2}||x_i - x_j||^2\right) \]

后验分布：

\[P(f_{t+1}|D, x_{t+1}) = N(\mu_t(x_{t+1}), \sigma_t^2(x_{t+1})) \]

期望改进采集函数：

\[EI(x) = (\mu_t(x) - f(x^+) - \xi)\phi(Z) + \sigma_t(x)\phi(Z) \]

六、算法执行逻辑：从参数筛选到优化

6.1 第一步：不确定性与敏感性分析

GLUE（广义似然不确定度估计）

  • 对8个修正系数（\(f_{dc}, f_{sc}, f_{bs}, f_r, k_p, k_w, k_{sc}, k_{bs}\)）进行拉丁超立方抽样，样本数10000；
  • 定义似然函数，筛选符合阈值的参数组合；
  • 计算后验分布和90%不确定度区间。

结果：汽室压力的90%不确定度区间能完全包含仿真对比曲线，说明模型可靠。

Sobol敏感性分析

  • 构造Sobol序列，计算一阶、二阶和总体敏感度指数；
  • 样本数从N=8到N=128，验证收敛性。

结果：

  • 对汽室压力影响最大的参数：\(k_p\)（一次侧换热修正系数）和\(k_w\)（壁面换热修正系数）；
  • 对沸腾段出口质量含汽率影响最大的参数：\(f_r\)（汽水分离段阻力修正系数）和\(f_{dc}\)（下降段阻力修正系数）。

最终筛选出4个关键参数：\(k_p, k_w, f_r, f_{dc}\)。

6.2 第二步：EnKF数据同化

  • 对4个关键参数进行拉丁超立方抽样，样本数100组；
  • 设置观测误差标准差\(\delta = 0.1\)；
  • 最大迭代次数1000，收敛标准为相对误差<0.5%；
  • 分别对30%、50%、70%、100%工况进行同化。

6.3 第三步：贝叶斯优化

  • 迭代次数300次；
  • 分别用GPR和TPE作为代理模型；
  • 目标函数为仿真结果与试验值的相对误差；
  • 同样对4个工况进行优化。

6.4 第四步：水位控制优化

  • 采用单冲量PID控制器；
  • 先用临界振荡法整定PID参数：\(K_u = 3460\)，\(T_u = 2.03\) s，计算得\(K_p = 1557\)，\(T_i = 1.692\) s；
  • 再用贝叶斯优化在临界振荡法结果±50%范围内搜索最优PID参数；
  • 目标函数为超调量与调节时间的乘积。

七、公式解析：关键公式的理解

7.1 蒸汽室压力方程

\[\frac{dP_g}{dt} = \frac{G_g - G_0}{V_{\text{cell}}\frac{\partial \rho_g}{\partial P_g}} \]

理解：汽室压力变化率取决于蒸汽流入流出质量流量差，以及蒸汽密度对压力的偏导数。这个方程将质量守恒转化为压力动态。

7.2 沸腾段能量方程

\[M_{bs}\frac{dh_{bs}}{dt} = G(h_{sl} - h_{bso}) + Q_{bs} \]

理解：沸腾段平均比焓的变化取决于流入流出焓差和换热量。由于假设含汽率线性分布，平均比焓可表示为：

\[h_{bs} = xh_{sg} + (1-x)h_{sl} \]

7.3 两相摩擦压降

\[\Delta P_f = (\Delta P_f)_{lo}\left[1 + x\left(\frac{\rho_f}{\rho_g} - 1\right)\right]\left[1 + x\left(\frac{\mu_l}{\mu_g} - 1\right)\right]^{-0.25} \]

理解：均相流模型中，两相压降通过全液相压降乘以两个修正因子得到，分别考虑密度差和粘度差的影响。

7.4 卡尔曼增益

\[K = PH^T(HPH^T + R)^{-1} \]

理解：卡尔曼增益决定了观测值与预测值的权重。\(P\)是背景误差协方差，\(R\)是观测误差协方差。\(P\)越大（预测越不确定），\(K\)越大，越信任观测；\(R\)越大（观测越不准确），\(K\)越小，越信任预测。

7.5 期望改进采集函数

\[EI(x) = (\mu_t(x) - f(x^+) - \xi)\phi(Z) + \sigma_t(x)\phi(Z) \]

理解：EI平衡了“探索”（高方差区域）和“利用”（高均值区域）。第一项鼓励在均值高于当前最优的区域采样，第二项鼓励在不确定性高的区域采样。

八、得出结果：优化效果显著

8.1 模型验证结果

工况 汽室压力相对误差 质量含汽率相对误差 换热量相对误差
30% 0.97% 0.32% 0.44%
50% 1.69% 3.82% 2.59%
70% 1.29% 4.98% 0.27%
100% 1.62% 3.56% 0.67%

所有相对误差均在5%以内，验证了模型的合理性。

8.2 EnKF优化结果

工况 压力相对误差 干度相对误差
30% 0.55% 0.11%
50% 1.41% 0.40%
70% 0.43% 0.66%
100% 1.03% 0.28%

优化后质量含汽率相对误差不超过1%，压力相对误差基本在1.5%以内。

8.3 贝叶斯优化结果

工况 GPR压力误差 TPE压力误差 GPR干度误差 TPE干度误差
30% 0.03% 0.19% 0.00% 0.00%
50% 0.07% 0.28% 0.07% 0.07%
70% 0.13% 0.23% 0.09% 0.23%
100% 0.06% 0.00% 0.28% 0.40%

优化后相对误差基本在0.5%以内，GPR代理模型略优于TPE。

8.4 方法对比

指标 EnKF 贝叶斯优化
仿真次数 100组 300次
压力误差 <1.5% <0.5%
干度误差 <1% <0.5%
计算成本 较低 较高
优势 计算成本低，适合在线校正 优化精度更高，结合历史信息

8.5 水位控制优化结果

  • 优化前：超调量较大，调节时间较长；
  • 优化后：超调量明显减小，调节时间缩短；
  • 给水流量能较快变化到与蒸汽流量一致，保证水位稳定。

九、代码模拟

# -*- coding: utf-8 -*-
"""
基于集合卡尔曼滤波与贝叶斯优化的蒸汽发生器仿真研究
简化演示代码：模拟完整算法流程
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from SALib.sample import saltelli
from SALib.analyze import sobol
from bayes_opt import BayesianOptimization
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)


# 1. 简化蒸汽发生器模型
def steam_generator_model(params, condition=1.0):
    """
    简化模型：输入8个修正系数，输出汽室压力和沸腾段出口干度
    params: [f_dc, f_sc, f_bs, f_r, k_p, k_w, k_sc, k_bs]
    condition: 工况水平（0.3, 0.5, 0.7, 1.0）
    """
    f_dc, f_sc, f_bs, f_r, k_p, k_w, k_sc, k_bs = params
    
    # 模拟压力：受换热系数和阻力系数影响
    pressure = 7.0 + 0.5*k_p + 0.3*k_w - 0.2*f_dc - 0.1*f_r + 0.05*condition
    # 模拟出口干度：受阻力系数和换热系数影响
    quality = 10.0 + 5.0*f_r + 3.0*f_dc + 2.0*k_p - 1.0*k_w + 0.5*condition
    
    # 加入非线性
    pressure = pressure + 0.1*np.sin(k_p*3) + 0.05*np.cos(f_r*2)
    quality = quality + 0.2*np.sin(f_dc*2) + 0.1*np.cos(k_w*3)
    
    return pressure, quality

# 真实参数（用于生成观测数据）
true_params = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])

# 参数范围（与论文表3.1一致）
param_ranges = np.array([
    [0.3, 1.7],  # f_dc
    [0.3, 1.7],  # f_sc
    [0.3, 1.7],  # f_bs
    [0.3, 1.7],  # f_r
    [0.5, 1.5],  # k_p
    [0.5, 1.5],  # k_w
    [0.5, 1.5],  # k_sc
    [0.5, 1.5]   # k_bs
])

# 生成观测数据（加入噪声）
def generate_observations(condition):
    p_true, q_true = steam_generator_model(true_params, condition)
    p_obs = p_true + np.random.normal(0, 0.02)
    q_obs = q_true + np.random.normal(0, 0.1)
    return p_obs, q_obs

# 2. GLUE 不确定性分析
def glue_analysis(n_samples=2000, condition=0.7):
    """广义似然不确定度估计"""
    print("\n===== GLUE 不确定性分析 =====")
    
    # 拉丁超立方抽样
    from scipy.stats import qmc
    sampler = qmc.LatinHypercube(d=8)
    samples = sampler.random(n=n_samples)
    
    # 缩放到参数范围
    params_samples = qmc.scale(samples, param_ranges[:,0], param_ranges[:,1])
    
    # 计算似然
    p_obs, q_obs = generate_observations(condition)
    likelihoods = []
    
    for p in params_samples:
        p_sim, q_sim = steam_generator_model(p, condition)
        # 压力似然
        L_p = 1 - (p_sim - p_obs)**2 / (0.1**2)
        # 干度似然
        L_q = 1 - (q_sim - q_obs)**2 / (1.0**2)
        L = max(0, L_p * L_q)
        likelihoods.append(L)
    
    likelihoods = np.array(likelihoods)
    
    # 筛选阈值
    threshold = 0.5
    mask = likelihoods > threshold
    if np.sum(mask) < 10:
        threshold = np.percentile(likelihoods, 50)
        mask = likelihoods > threshold
    
    # 后验分布
    posterior_samples = params_samples[mask]
    posterior_weights = likelihoods[mask] / np.sum(likelihoods[mask])
    
    # 计算90%不确定度区间
    p_sim_all = np.array([steam_generator_model(p, condition)[0] for p in posterior_samples])
    q_sim_all = np.array([steam_generator_model(p, condition)[1] for p in posterior_samples])
    
    p_5, p_95 = np.percentile(p_sim_all, [5, 95])
    q_5, q_95 = np.percentile(q_sim_all, [5, 95])
    
    print(f"压力90%不确定度区间: [{p_5:.3f}, {p_95:.3f}] MPa")
    print(f"干度90%不确定度区间: [{q_5:.3f}, {q_95:.3f}] %")
    
    return posterior_samples, posterior_weights

# 3. Sobol 敏感性分析

def sobol_sensitivity(n_samples=128, condition=1.0):
    """Sobol 全局敏感性分析"""
    print("\n===== Sobol 敏感性分析 =====")
    
    problem = {
        'num_vars': 8,
        'names': ['f_dc', 'f_sc', 'f_bs', 'f_r', 'k_p', 'k_w', 'k_sc', 'k_bs'],
        'bounds': param_ranges.tolist()
    }
    
    # 生成样本
    param_values = saltelli.sample(problem, n_samples, calc_second_order=False)
    
    # 计算输出
    Y_pressure = np.array([steam_generator_model(p, condition)[0] for p in param_values])
    Y_quality = np.array([steam_generator_model(p, condition)[1] for p in param_values])
    
    # 分析
    Si_p = sobol.analyze(problem, Y_pressure, print_to_console=False)
    Si_q = sobol.analyze(problem, Y_quality, print_to_console=False)
    
    print("汽室压力总体敏感度:")
    for i, name in enumerate(problem['names']):
        print(f"  {name}: {Si_p['ST'][i]:.4f}")
    
    print("沸腾段出口干度总体敏感度:")
    for i, name in enumerate(problem['names']):
        print(f"  {name}: {Si_q['ST'][i]:.4f}")
    
    # 筛选关键参数（总体敏感度>0.1）
    key_params_p = [problem['names'][i] for i in range(8) if Si_p['ST'][i] > 0.1]
    key_params_q = [problem['names'][i] for i in range(8) if Si_q['ST'][i] > 0.1]
    key_params = list(set(key_params_p + key_params_q))
    print(f"筛选出的关键参数: {key_params}")
    
    return key_params

# 4. 集合卡尔曼滤波 EnKF

def enkf_assimilation(key_params, n_ensemble=100, n_iter=50, condition=0.7):
    """集合卡尔曼滤波参数同化"""
    print("\n===== EnKF 数据同化 =====")
    
    # 关键参数索引
    key_indices = [['f_dc', 'f_sc', 'f_bs', 'f_r', 'k_p', 'k_w', 'k_sc', 'k_bs'].index(k) for k in key_params]
    n_key = len(key_indices)
    
    # 初始化集合
    ensemble = np.zeros((n_ensemble, n_key))
    for i, idx in enumerate(key_indices):
        ensemble[:, i] = np.random.uniform(param_ranges[idx, 0], param_ranges[idx, 1], n_ensemble)
    
    # 观测值
    p_obs, q_obs = generate_observations(condition)
    obs = np.array([p_obs, q_obs])
    R = np.diag([0.02**2, 0.1**2])  # 观测误差协方差
    
    # 迭代
    for it in range(n_iter):
        # 预测
        states = []
        for i in range(n_ensemble):
            full_params = true_params.copy()
            for j, idx in enumerate(key_indices):
                full_params[idx] = ensemble[i, j]
            p_sim, q_sim = steam_generator_model(full_params, condition)
            states.append([p_sim, q_sim])
        states = np.array(states)
        
        # 计算均值和协方差
        mean_state = np.mean(states, axis=0)
        P = np.cov(states.T)
        
        # 卡尔曼增益
        H = np.eye(2)
        K = P @ H.T @ np.linalg.inv(H @ P @ H.T + R)
        
        # 更新
        for i in range(n_ensemble):
            innovation = obs + np.random.multivariate_normal([0,0], R) - states[i]
            update = K @ innovation
            # 更新参数（简化的参数-状态映射）
            ensemble[i] += update[:n_key] * 0.1
        
        # 限制范围
        for i, idx in enumerate(key_indices):
            ensemble[:, i] = np.clip(ensemble[:, i], param_ranges[idx, 0], param_ranges[idx, 1])
        
        if it % 10 == 0:
            mean_params = np.mean(ensemble, axis=0)
            full_params = true_params.copy()
            for j, idx in enumerate(key_indices):
                full_params[idx] = mean_params[j]
            p_sim, q_sim = steam_generator_model(full_params, condition)
            print(f"迭代 {it}: 压力={p_sim:.3f}, 干度={q_sim:.3f}")
    
    # 最终最优参数
    best_params = true_params.copy()
    mean_params = np.mean(ensemble, axis=0)
    for j, idx in enumerate(key_indices):
        best_params[idx] = mean_params[j]
    
    p_best, q_best = steam_generator_model(best_params, condition)
    print(f"EnKF 优化后: 压力={p_best:.3f}, 干度={q_best:.3f}")
    print(f"真实值: 压力={p_obs:.3f}, 干度={q_obs:.3f}")
    
    return best_params

# 5. 贝叶斯优化

def bayesian_optimization(key_params, condition=0.7):
    """贝叶斯优化参数寻优"""
    print("\n===== 贝叶斯优化 =====")
    
    key_indices = [['f_dc', 'f_sc', 'f_bs', 'f_r', 'k_p', 'k_w', 'k_sc', 'k_bs'].index(k) for k in key_params]
    n_key = len(key_indices)
    
    # 观测值
    p_obs, q_obs = generate_observations(condition)
    
    # 目标函数（负误差，最大化）
    def objective(**kwargs):
        full_params = true_params.copy()
        for j, idx in enumerate(key_indices):
            full_params[idx] = kwargs[f'x{j}']
        p_sim, q_sim = steam_generator_model(full_params, condition)
        error = -((p_sim - p_obs)**2 + (q_sim - q_obs)**2)
        return error
    
    # 参数边界
    pbounds = {}
    for j, idx in enumerate(key_indices):
        pbounds[f'x{j}'] = (param_ranges[idx, 0], param_ranges[idx, 1])
    
    # 贝叶斯优化
    optimizer = BayesianOptimization(
        f=objective,
        pbounds=pbounds,
        random_state=42,
        verbose=0
    )
    
    optimizer.maximize(init_points=10, n_iter=50)
    
    # 最优参数
    best_params = true_params.copy()
    for j, idx in enumerate(key_indices):
        best_params[idx] = optimizer.max['params'][f'x{j}']
    
    p_best, q_best = steam_generator_model(best_params, condition)
    print(f"BO 优化后: 压力={p_best:.3f}, 干度={q_best:.3f}")
    print(f"真实值: 压力={p_obs:.3f}, 干度={q_obs:.3f}")
    
    return best_params

# 6. 水位控制 PID + 贝叶斯优化

def pid_control(Kp, Ti, Td, setpoint=13.0, sim_time=100, dt=0.1):
    """模拟PID水位控制"""
    time = np.arange(0, sim_time, dt)
    level = np.zeros_like(time)
    level[0] = 12.5
    integral = 0
    prev_error = 0
    
    for i in range(1, len(time)):
        error = setpoint - level[i-1]
        integral += error * dt
        derivative = (error - prev_error) / dt
        u = Kp * (error + integral / Ti + Td * derivative)
        # 简化水位动态：给水流量影响水位
        level[i] = level[i-1] + 0.01 * u - 0.005 * (level[i-1] - 12.0)
        prev_error = error
    
    # 计算性能指标
    overshoot = max(0, np.max(level) - setpoint) / setpoint * 100
    settling_time = 0
    for i in range(len(time)-1, -1, -1):
        if abs(level[i] - setpoint) > 0.02 * setpoint:
            settling_time = time[i]
            break
    return overshoot, settling_time, level, time

def optimize_pid():
    """贝叶斯优化PID参数"""
    print("\n===== 水位控制 PID 贝叶斯优化 =====")
    
    # 临界振荡法得到的基准参数
    Ku = 3460
    Tu = 2.03
    Kp_base = 0.6 * Ku
    Ti_base = Tu / 2
    Td_base = Tu / 8
    
    def objective(Kp, Ti, Td):
        overshoot, settling_time, _, _ = pid_control(Kp, Ti, Td)
        # 目标：最小化超调量和调节时间的乘积
        return - (overshoot * settling_time)
    
    pbounds = {
        'Kp': (Kp_base * 0.5, Kp_base * 1.5),
        'Ti': (Ti_base * 0.5, Ti_base * 1.5),
        'Td': (Td_base * 0.5, Td_base * 1.5)
    }
    
    optimizer = BayesianOptimization(
        f=objective,
        pbounds=pbounds,
        random_state=42,
        verbose=0
    )
    
    optimizer.maximize(init_points=5, n_iter=30)
    
    best_Kp = optimizer.max['params']['Kp']
    best_Ti = optimizer.max['params']['Ti']
    best_Td = optimizer.max['params']['Td']
    
    print(f"优化前 PID: Kp={Kp_base:.2f}, Ti={Ti_base:.4f}, Td={Td_base:.4f}")
    print(f"优化后 PID: Kp={best_Kp:.2f}, Ti={best_Ti:.4f}, Td={best_Td:.4f}")
    
    # 对比
    over1, set1, level1, time1 = pid_control(Kp_base, Ti_base, Td_base)
    over2, set2, level2, time2 = pid_control(best_Kp, best_Ti, best_Td)
    
    print(f"优化前: 超调量={over1:.2f}%, 调节时间={set1:.2f}s")
    print(f"优化后: 超调量={over2:.2f}%, 调节时间={set2:.2f}s")
    
    # 绘图
    plt.figure(figsize=(10, 5))
    plt.plot(time1, level1, label='优化前')
    plt.plot(time2, level2, label='优化后')
    plt.axhline(y=13.0, color='r', linestyle='--', label='设定值')
    plt.xlabel('时间 (s)')
    plt.ylabel('水位 (m)')
    plt.title('PID 水位控制优化对比')
    plt.legend()
    plt.grid(True)
    plt.savefig('pid_optimization.png')
    plt.show()

# 主程序

if __name__ == "__main__":
    print("="*60)
    print("基于集合卡尔曼滤波与贝叶斯优化的蒸汽发生器仿真研究")
    print("简化演示代码")
    print("="*60)
    
    # 1. GLUE 不确定性分析
    posterior_samples, posterior_weights = glue_analysis(n_samples=2000, condition=0.7)
    
    # 2. Sobol 敏感性分析
    key_params = sobol_sensitivity(n_samples=128, condition=1.0)
    
    # 3. EnKF 数据同化
    best_params_enkf = enkf_assimilation(key_params, n_ensemble=100, n_iter=50, condition=0.7)
    
    # 4. 贝叶斯优化
    best_params_bo = bayesian_optimization(key_params, condition=0.7)
    
    # 5. 水位控制 PID + 贝叶斯优化
    optimize_pid()
    
    print("\n===== 全部流程完成 =====")
posted @ 2026-10-11 19:29  晚叙^_^  阅读(10)  评论(0)    收藏  举报
-->