利用灵敏性分析确定模型参数的敏感性

   import numpy as np

   def calculate_combined_difference(Y, Y_b):
      # 位置差异 (Position difference)
      delta_p = np.mean(Y) - np.mean(Y_b)
      
      # 尺度差异 (Scale difference)
      delta_s = (np.max(Y) - np.min(Y)) - (np.max(Y_b) - np.min(Y_b))
      
      # 形状差异 (Shape difference - RMSE)
      delta_r = np.sqrt(np.mean((np.array(Y) - np.array(Y_b))**2)) # 确保Y, Y_b为numpy数组
      
      # 综合差异 (Combined difference)
      Delta_i = delta_p + delta_s + delta_r
      return Delta_i

   param_definitions = {
      "Negative electrode thickness [m]": {
         "default": default_param["Negative electrode thickness [m]"],
         "bounds_factor": 0.2,  # 变动 +/- 20%
      },
      "Positive electrode thickness [m]": {
         "default": default_param["Positive electrode thickness [m]"],
         "bounds_factor": 0.2,
      },
      "Separator thickness [m]": {
         "default": default_param["Separator thickness [m]"],
         "bounds_factor": 0.2,
      },
      "Negative particle radius [m]": {
         "default": default_param["Negative particle radius [m]"],
         "bounds_factor": 0.2,
      },
      "Positive particle radius [m]": {
         "default": default_param["Positive particle radius [m]"],
         "bounds_factor": 0.2,
      },
   }

   problem = {
      "num_vars": len(param_definitions),
      "names": list(param_definitions.keys()),
      "bounds": [],
   }

   for name, props in param_definitions.items():
      lower_bound = props["default"] * (1 - props["bounds_factor"])
      upper_bound = props["default"] * (1 + props["bounds_factor"])
      problem["bounds"].append([lower_bound, upper_bound])
      print(
         f"Parameter: {name}, Default: {props['default']:.2e}, Bounds: [{lower_bound:.2e}, {upper_bound:.2e}]"
      )

   from SALib.sample import saltelli
   # import pybamm # 假设pybamm已安装
   # import numpy as np # 假设numpy已安装

   # 生成参数样本
   # N 通常取2的幂次，例如1024。样本总数将是 N * (D + 2) 或 N * (2D + 2)
   # D 是参数数量 (problem['num_vars'])
   param_values = saltelli.sample(problem, 1024) 

   # 运行一次基线仿真
   model_baseline = pybamm.lithium_ion.SPMe()  
   param_baseline = pybamm.ParameterValues("Chen2020")  
   experiment_baseline = pybamm.Experiment(["Discharge at 1C until 2.5 V"])  
   sim_baseline = pybamm.Simulation(  
      model_baseline, experiment=experiment_baseline, parameter_values=param_baseline  
   )  
   sol_baseline = sim_baseline.solve()  
   t_eval_baseline = sol_baseline["Time [s]"].entries  
   voltage_baseline = sol_baseline["Terminal voltage [V]"].entries

   # 为了确保所有仿真都在相同的时间点上进行比较，定义一个固定的时间向量 t_eval_common
   # 使用基线仿真的最大时间，并插值到 200 个点
   max_time_baseline = t_eval_baseline[-1]  
   t_eval_common = np.linspace(0, max_time_baseline, 200) # 200 个评估点

   # 对基线解在 t_eval_common 上进行插值
   voltage_baseline_interp = np.interp(t_eval_common, t_eval_baseline, voltage_baseline)

   print(  
      "Baseline simulation complete and common evaluation time points (t_eval_common) prepared."  
   )  
   print(f"Baseline simulation ran for {max_time_baseline:.2f} seconds.")  
   print(  
      f"Common evaluation time points range from {t_eval_common[0]:.2f}s to {t_eval_common[-1]:.2f}s with {len(t_eval_common)} points."  
   )

   def evaluate_model_delta_i(parameter_sample):  
      """  
      运行 PyBaMM SPMe 模型并计算与基线电压曲线的综合差异指标 Delta_i.  
      Delta_i = delta_p + delta_s + delta_r  
      其中:  
         delta_p: 位置差异 (mean_run - mean_baseline)  
         delta_s: 尺度差异 ((max_run - min_run) - (max_baseline - min_baseline))  
         delta_r: 形状差异 (RMSE(run, baseline))  
      参数:  
         parameter_sample (np.array): SALib 生成的参数样本。  
      返回:  
         float: 电压曲线的综合差异指标 Delta_i。  
      """  
      model_run = pybamm.lithium_ion.SPMe()  
      param_run = pybamm.ParameterValues("Chen2020") # 每次都从新的默认参数开始

      current_params = {}
      for i, name in enumerate(problem["names"]):
         current_params[name] = parameter_sample[i]

      try:
         param_run.update(current_params)
         sim_run = pybamm.Simulation(
               model_run, experiment=experiment_baseline, parameter_values=param_run
         )
         sol_run = sim_run.solve(t_eval=t_eval_common)
         voltage_run = sol_run["Terminal voltage [V]"].entries

         if len(voltage_run) < len(t_eval_common):
               padding_value = (
                  voltage_run[-1]
                  if len(voltage_run) > 0
                  else param_run["Lower voltage cut-off [V]"]
               )
               voltage_run_padded = np.full_like(t_eval_common, padding_value)
               voltage_run_padded[: len(voltage_run)] = voltage_run
               voltage_run = voltage_run_padded

         # 计算 delta_p (位置差异)
         delta_p = np.mean(voltage_run) - np.mean(voltage_baseline_interp)

         # 计算 delta_s (尺度差异)
         range_run = np.max(voltage_run) - np.min(voltage_run)
         range_baseline = np.max(voltage_baseline_interp) - np.min(voltage_baseline_interp)
         delta_s = range_run - range_baseline

         # 计算 delta_r (形状差异 - RMSE)
         delta_r = np.sqrt(np.mean((voltage_run - voltage_baseline_interp) ** 2))

         # 计算综合差异指标 Delta_i
         # 注意：直接相加可能会因为各项尺度不同导致某一项主导。可以考虑加权或归一化，但按当前要求直接相加。
         delta_i = delta_p + delta_s + delta_r

         # print(f"Run with {current_params}, d_p={delta_p:.3f}, d_s={delta_s:.3f}, d_r={delta_r:.3f}, Delta_i={delta_i:.3f}")
         return delta_i

      except Exception as e:
         # print(f"Error during simulation with params {current_params}: {e}")
         return 1e6  # 一个较大的惩罚值

   from SALib.analyze import sobol

   # 执行 Sobol 分析
   # Y_outputs 应为上一步仿真得到的实际输出结果数组
   # 注意处理Y_outputs中的NaN值（如果仿真中可能出现）
   # Y_outputs_cleaned = Y_outputs[~np.isnan(Y_outputs)]
   # param_values_cleaned = param_values[~np.isnan(Y_outputs)] 
   # 如果清理了Y_outputs，对应的param_values也需要清理，但这会破坏Saltelli样本结构
   # 更好的做法是在仿真失败时赋一个合理的惩罚值，而不是NaN，以保持样本完整性
   Si = sobol.analyze(problem, Y_outputs, calc_second_order=True, print_to_console=False)

   # 提取一阶和总效应指数
   S1_indices = Si['S1']
   ST_indices = Si['ST']
   
   print("一阶 Sobol 指数 (S1):")
   for name, s1_val in zip(problem['names'], S1_indices):
      print(f"- {name}: {s1_val:.4f}")

   print("\\n总 Sobol 指数 (ST):")
   for name, st_val in zip(problem['names'], ST_indices):
      print(f"- {name}: {st_val:.4f}")

   # 如果 calc_second_order=True, 还可以获取二阶指数 Si['S2'] 和对应的参数对 Si['S2_conf']