import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, Matern, ConstantKernel
from skopt import gp_minimize
from skopt.space import Real
from skopt.utils import use_named_args
from skopt.plots import plot_convergence, plot_objective
from scipy.optimize import minimize
from sklearn.preprocessing import StandardScaler
import warnings
from sklearn.exceptions import ConvergenceWarning
# 设置美观的绘图风格
sns.set(style="whitegrid", palette="muted", font_scale=1.1)
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['axes.labelsize'] = 12
plt.rcParams['xtick.labelsize'] = 10
plt.rcParams['ytick.labelsize'] = 10
# 读取数据
df = pd.read_excel('FEM_data.xlsx')
print("数据预览:")
print(df.head())
# 检查数据分布
plt.figure(figsize=(12, 8))
sns.pairplot(df, diag_kind='kde')
plt.suptitle('Feature and Target Distribution', y=1.02)
plt.tight_layout()
plt.savefig('data_distribution.png', dpi=300)
plt.show()
# 分离特征和目标
X = df[['f_nb', 'E', 'po', 'v_r', 'e_0.2', 'P', 'r']].values
S = df['S'].values.reshape(-1, 1) # 需要最小化
RM = df['RM'].values.reshape(-1, 1) # 需要最大化
# 标准化特征和目标
scaler_X = StandardScaler().fit(X)
scaler_S = StandardScaler().fit(S)
scaler_RM = StandardScaler().fit(RM)
X_scaled = scaler_X.transform(X)
S_scaled = scaler_S.transform(S).ravel()
RM_scaled = scaler_RM.transform(RM).ravel()
# 定义组合目标:S越小越好,RM越大越好 -> 最小化 [S_scaled - RM_scaled]
combined_target = S_scaled - RM_scaled
# 为贝叶斯优化定义参数空间
param_space = [
Real(0.1, 0.9, name='f_nb'), # Continuous variable
Real(120000, 280000, name='E'), # Young's modulus
Real(0.27, 0.34, name='po'), # Poisson's ratio
Real(0.314, 10.314, name='v_r'), # Volume ratio
Real(680, 875, name='e_0.2'), # Yield strength
Real(0, 4000, name='P'), # Load
Real(21, 25, name='r') # Radius
]
# 创建GPR模型 - 进一步扩大核参数边界防止收敛问题
def create_kernel():
return ConstantKernel(1.0, constant_value_bounds=(1e-6, 1e3)) * Matern(
length_scale=1.0, length_scale_bounds=(1e-6, 1e3), nu=2.5
) + ConstantKernel(1.0, constant_value_bounds=(1e-6, 1e3)) * RBF(
length_scale=0.5, length_scale_bounds=(1e-6, 1e3)
)
kernel = create_kernel()
gpr = GaussianProcessRegressor(
kernel=kernel,
alpha=1e-5, # 减小噪声项
n_restarts_optimizer=25, # 增加重启次数
random_state=42,
normalize_y=True # 内部标准化目标值
)
# 训练初始GPR模型
gpr.fit(X_scaled, combined_target)
print(f"初始GPR模型参数: {gpr.kernel_}")
# 定义贝叶斯优化目标函数
@use_named_args(param_space)
def objective(**params):
# 提取参数值
param_values = np.array([params[dim.name] for dim in param_space]).reshape(1, -1)
# 标准化输入
param_scaled = scaler_X.transform(param_values)
# 预测组合目标
mean_pred, std_pred = gpr.predict(param_scaled, return_std=True)
# 获取置信下界 (Lower Confidence Bound) - 平衡探索与开发
lcb = mean_pred - 0.5 * std_pred
return lcb[0] # 返回标量值
# 运行贝叶斯优化 - 使用8个初始点(2的幂)
res_gp = gp_minimize(
func=objective,
dimensions=param_space,
n_calls=50, # 额外评估次数
n_random_starts=8, # 使用2的幂(8个)初始点
initial_point_generator='sobol', # 更好的空间覆盖
random_state=42,
acq_func='EI'
)
print("\n贝叶斯优化结果:")
print(f"最优参数组合: {res_gp.x}")
print(f"最佳目标值: {res_gp.fun}")
# 获取最优参数
optimal_params = np.array(res_gp.x).reshape(1, -1)
optimal_params_scaled = scaler_X.transform(optimal_params)
# 预测最优点的组合目标值
combined_pred = gpr.predict(optimal_params_scaled)
# 分别预测S和RM的标准化值
# 创建新的GPR模型用于S和RM
gpr_S = GaussianProcessRegressor(
kernel=create_kernel(), # 使用相同的核结构
alpha=1e-5,
n_restarts_optimizer=25, # 增加重启次数
random_state=42,
normalize_y=True
)
gpr_S.fit(X_scaled, S_scaled)
gpr_RM = GaussianProcessRegressor(
kernel=create_kernel(), # 使用相同的核结构
alpha=1e-5,
n_restarts_optimizer=25, # 增加重启次数
random_state=42,
normalize_y=True
)
gpr_RM.fit(X_scaled, RM_scaled)
# 预测最优点的S和RM
S_pred_scaled = gpr_S.predict(optimal_params_scaled)
RM_pred_scaled = gpr_RM.predict(optimal_params_scaled)
# 反标准化得到实际值
S_pred = scaler_S.inverse_transform(S_pred_scaled.reshape(-1, 1))
RM_pred = scaler_RM.inverse_transform(RM_pred_scaled.reshape(-1, 1))
print("\n预测最优值:")
print(f"S_min: {S_pred[0][0]:.2f}")
print(f"RM_max: {RM_pred[0][0]:.2f}")
# 可视化优化过程
plt.figure(figsize=(12, 8))
plot_convergence(res_gp)
plt.title('Bayesian Optimization Convergence', fontsize=14)
plt.tight_layout()
plt.savefig('bo_convergence.png', dpi=300)
plt.show()
# 绘制目标函数响应面 - 传递所有维度名称
dimension_names = [dim.name for dim in param_space]
plt.figure(figsize=(16, 12))
ax = plot_objective(res_gp, dimensions=dimension_names)
plt.suptitle('Objective Function Response Surface', fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.96]) # 为标题留出空间
plt.savefig('objective_response.png', dpi=300)
plt.show()
# 分析参数敏感性
sensitivity = []
# 创建包装函数,处理参数传递问题
def wrapped_objective(params):
"""包装函数,处理参数传递问题"""
param_dict = {dim.name: val for dim, val in zip(param_space, params)}
return objective(**param_dict)
# 使用中心差分法计算梯度
delta = 0.05
base_params = np.array(res_gp.x)
for idx, dim in enumerate(param_space):
# 上限参数值
up_params = base_params.copy()
up_params[idx] = min(up_params[idx] + delta * (dim.bounds[1] - dim.bounds[0]), dim.bounds[1])
up_score = -wrapped_objective(up_params)
# 下限参数值
low_params = base_params.copy()
low_params[idx] = max(low_params[idx] - delta * (dim.bounds[1] - dim.bounds[0]), dim.bounds[0])
low_score = -wrapped_objective(low_params)
sensitivity.append({
'parameter': dim.name,
'absolute_sensitivity': up_score - low_score,
'relative_sensitivity': (up_score - low_score) / (abs(up_params[idx] - low_params[idx]) + 1e-10),
'range': dim.bounds
})
# 创建参数敏感性数据框
sensitivity_df = pd.DataFrame(sensitivity).sort_values('absolute_sensitivity', ascending=False)
# 绘制参数敏感性
plt.figure(figsize=(12, 7))
sns.barplot(
y='parameter',
x='absolute_sensitivity',
data=sensitivity_df,
palette='viridis'
)
plt.xlabel('Sensitivity to Objective Function', fontsize=12)
plt.ylabel('Parameter', fontsize=12)
plt.title('Parameter Sensitivity Analysis', fontsize=14)
plt.tight_layout()
plt.savefig('parameter_sensitivity.png', dpi=300)
plt.show()
# 打印完整优化报告
print("\n=================== OPTIMIZATION REPORT ===================")
print(f"最优参数组合 (原始值):")
for dim, val in zip(param_space, res_gp.x):
print(f"- {dim.name}: {val:.4f} (范围: [{dim.bounds[0]}, {dim.bounds[1]}])")
# 计算实际提升百分比
S_mean = df['S'].mean()
RM_mean = df['RM'].mean()
S_improvement = (S_mean - S_pred[0][0]) / S_mean * 100
RM_improvement = (RM_pred[0][0] - RM_mean) / RM_mean * 100
print(f"\n目标提升预测:")
print(f"S最小值 (预测): {S_pred[0][0]:.2f} (降幅: {S_improvement:.1f}%)")
print(f"RM最大值 (预测): {RM_pred[0][0]:.2f} (增幅: {RM_improvement:.1f}%)")
print("\n参数敏感性排序:")
print(sensitivity_df[['parameter', 'absolute_sensitivity']].to_string(index=False))
# 保存优化结果
params_names = [dim.name for dim in param_space]
results_df = pd.DataFrame([res_gp.x], columns=params_names)
results_df['S_pred'] = S_pred[0][0]
results_df['RM_pred'] = RM_pred[0][0]
results_df.to_csv('optimization_results.csv', index=False)
print("\n优化结果已保存到 optimization_results.csv")报错:TypeError: objective() got an unexpected keyword argument 'f_nb'
最新发布