Design-Optimization-Workflow/Function_BO.py at main · Jinge91/Design-Optimization-Workflow · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def runBO(input, ):
    import numpy as np
    import matplotlib.pyplot as plt
    from skopt.plots import plot_gaussian_process
    from skopt.learning import GaussianProcessRegressor
    from sklearn.gaussian_process.kernels import Matern
    import pandas as pd


    width = data['width']
    # print("\nChannel width:")
    # print(width)

    flowrate = data.iloc[:,1]
    # print("\nFlow rate")
    # print(flowrate)

    En = data.iloc[:,2]

    # Create the 2D input array (X) and 1D output array (Y)
    # .values converts the DataFrame selection into a NumPy array ------------------------------
    input_columns = ['width', 'flow rate']
    output_column = 'entropy'
    X_in = data[input_columns].values
    Y_out = data[output_column].values
    Y_out = np.log(Y_out)
    # print(X_in)
    # print(Y_out)

    # Fit GP model ---------------------------------------------------------------------
    kernel = Matern(nu=2.5)
    gpr_model = GaussianProcessRegressor(kernel=kernel,
                                        n_restarts_optimizer=10,
                                        random_state=0)

    En_gpmodel= gpr_model.fit(X_in,Y_out)

    # creat meshgrid for design space-------------------------------------------------
    W_range = np.linspace(0.2, 2, 200)
    F_range = np.linspace(0.0001, 1, 200)
    W,F = np.meshgrid(W_range,F_range)
    X_plot = np.vstack((W.ravel(), F.ravel())).T

    # get predicted mean and std -----------------------------------------------------
    predicted_mean, predicted_std = En_gpmodel.predict(X_plot, return_std=True)

    print(f"the average of std is{np.mean(predicted_std)}")


    # calculate standard deviation bound ---------------
    visualization_factor = 1
    upper_bound = predicted_mean + predicted_std * visualization_factor
    lower_bound = predicted_mean - predicted_std * visualization_factor

    # --- 2. 创建一个3D图像 ---
    fig = plt.figure(figsize=(12, 9))
    ax = fig.add_subplot(111, projection='3d')

    # --- 3. 绘制核心图层 ---

    # a) 绘制预测均值曲面 (颜色鲜艳，不透明)
    #    cmap='viridis' 是一个常用的漂亮色系
    ax.plot_surface(W, F, predicted_mean.reshape(F.shape), cmap='viridis',
                    edgecolor='none', label='Predicted Mean')

    # b) 绘制置信区间曲面 (灰色，半透明)
    #    alpha=0.3 控制透明度，让我们可以透过它看到均值曲面
    ax.plot_surface(W, F, upper_bound.reshape(F.shape), color='gray',
                    alpha=0.5, edgecolor='none', rstride=10, cstride=10)
    ax.plot_surface(W, F, lower_bound.reshape(F.shape), color='gray',
                    alpha=0.5, edgecolor='none', rstride=10, cstride=10)

    # c) 绘制原始样本点 (红色散点，突出显示)
    #    X_in[:, 0] 是第一个输入特征 (width)
    #    X_in[:, 1] 是第二个输入特征 (flow rate)
    #    Y_out 是对应的输出 (entropy)
    ax.scatter(X_in[:, 0], X_in[:, 1], Y_out, color='red', s=50, edgecolor='black', depthshade=True, label='Sample Points')


    # --- 4. 设置图表信息 ---
    ax.set_title('GP Model of entropy', fontsize=16)
    ax.set_xlabel('Width', fontsize=12, labelpad=10)
    ax.set_ylabel('Flow Rate', fontsize=12, labelpad=10)
    ax.set_zlabel('Entropy', fontsize=12, labelpad=10)
    ax.legend()

    plt.show()

    from skopt.acquisition import gaussian_ei

    # --------------- EI 函数计算 ------------------

    # 1. 找到当前已观测到的最优值 (最小值)
    #    EI 函数需要这个值作为基准
    y_opt = np.min(Y_out)

    # 2. 计算网格上每个点的预期改进值
    #    这就是您要的采集函数指令
    ei_values = gaussian_ei(X=X_plot, model=En_gpmodel, y_opt=y_opt)

    print(f"\n已成功为网格上的 {len(ei_values)} 个点计算了EI值。")


    # --- 3. 可视化预期改进函数曲面 ---
    # 将一维的EI值数组重塑为二维网格
    Z_ei = ei_values.reshape(W.shape)

    fig = plt.figure(figsize=(10, 7))
    ax = fig.add_subplot(111, projection='3d')

    surf = ax.plot_surface(W, F, Z_ei, cmap='plasma', edgecolor='none')
    ax.set_title('Expected Improvement (EI) Acquisition Surface', fontsize=16)
    ax.set_xlabel('Width', fontsize=12)
    ax.set_ylabel('Flow Rate', fontsize=12)
    ax.set_zlabel('EI Value', fontsize=12)
    fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)

    # 找到EI值最大的点，并在图上标记出来
    max_ei_index = np.argmax(ei_values)
    best_next_point = X_plot[max_ei_index]
    ax.scatter(best_next_point[0], best_next_point[1], np.max(ei_values),
            color='red', s=100, edgecolor='black', depthshade=True,
            label=f'Max EI (Next Point)')

    ax.legend()
    plt.show()

    next_width = best_next_point[0]
    next_flow = best_next_point[1]
    print(f"next point to sample are channel width = {best_next_point[0]:.4f} mm, flow rate = {best_next_point[1]:.4f} kg/s")

    return next_width, next_flow

if __name__ =="__main__":
    result_data = runBO()