Chapter 5: SciPy calculation

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import minimize from mpl_toolkits.mplot3d import Axes3D # Rosenbrockfunction（suitable benchmarkproblem） def rosenbrock(x): """ Rosenbrockfunction: f(x, y) = (1-x)^2 + 100(y-x^2)^2 minimum value: f(1, 1) = 0 """ return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2 def rosenbrock_grad(x): """Rosenbrockfunction gradient""" dx = -2 * (1 - x[0]) - 400 * x[0] * (x[1] - x[0]**2) dy = 200 * (x[1] - x[0]**2) return np.array([dx, dy]) print("=" * 60) print("unconstrainedsuitable: Rosenbrockfunction") print("=" * 60) # initial value x0 = np.array([0.0, 0.0]) # method1: Nelder-Mead（functionnot required） result_nm = minimize(rosenbrock, x0, method='Nelder-Mead') print("\n1. Nelder-Mead（functionnot required）") print(f" suitableSolution: x* = {result_nm.x}") print(f" minimum value: f(x*) = {result_nm.fun:.6e}") print(f" iterationiterations: {result_nm.nit}") print(f" functioniterations: {result_nm.nfev}") # method2: BFGS（quasi-Newton method） result_bfgs = minimize(rosenbrock, x0, method='BFGS', jac=rosenbrock_grad) print("\n2. BFGS（quasi-Newton method、gradientuse）") print(f" suitableSolution: x* = {result_bfgs.x}") print(f" minimum value: f(x*) = {result_bfgs.fun:.6e}") print(f" iterationiterations: {result_bfgs.nit}") print(f" functioniterations: {result_bfgs.nfev}") # method3: CG（conjugate gradient method） result_cg = minimize(rosenbrock, x0, method='CG', jac=rosenbrock_grad) print("\n3. CG（conjugate gradient method）") print(f" suitableSolution: x* = {result_cg.x}") print(f" minimum value: f(x*) = {result_cg.fun:.6e}") print(f" iterationiterations: {result_cg.nit}") print(f" functioniterations: {result_cg.nfev}") # fig = plt.figure(figsize=(15, 5)) # : and suitable ax1 = fig.add_subplot(131) x = np.linspace(-2, 2, 200) y = np.linspace(-1, 3, 200) X, Y = np.meshgrid(x, y) Z = (1 - X)**2 + 100 * (Y - X**2)**2 levels = np.logspace(-1, 3, 20) cs = ax1.contour(X, Y, Z, levels=levels, cmap='viridis', alpha=0.6) ax1.clabel(cs, inline=True, fontsize=8) ax1.plot(1, 1, 'r*', markersize=20, label='minimum value (1, 1)') ax1.plot(x0[0], x0[1], 'ko', markersize=10, label='initial value') ax1.plot(result_bfgs.x[0], result_bfgs.x[1], 'bs', markersize=10, label='BFGS') ax1.set_xlabel('x', fontsize=12) ax1.set_ylabel('y', fontsize=12) ax1.set_title('Rosenbrockfunction ', fontsize=14) ax1.legend(fontsize=10) ax1.grid(True, alpha=0.3) # : 3D ax2 = fig.add_subplot(132, projection='3d') ax2.plot_surface(X, Y, np.log10(Z + 1), cmap='viridis', alpha=0.6) ax2.scatter([1], [1], [0], color='red', s=100, label='minimum value') ax2.set_xlabel('x', fontsize=10) ax2.set_ylabel('y', fontsize=10) ax2.set_zlabel('log₁₀(f(x,y) + 1)', fontsize=10) ax2.set_title('3D', fontsize=14) # : method ax3 = fig.add_subplot(133) methods = ['Nelder-Mead', 'BFGS', 'CG'] iterations = [result_nm.nit, result_bfgs.nit, result_cg.nit] func_evals = [result_nm.nfev, result_bfgs.nfev, result_cg.nfev] x_pos = np.arange(len(methods)) width = 0.35 bars1 = ax3.bar(x_pos - width/2, iterations, width, label='iterationiterations', alpha=0.7) bars2 = ax3.bar(x_pos + width/2, func_evals, width, label='functioniterations', alpha=0.7) ax3.set_ylabel('iterations', fontsize=12) ax3.set_title('suitablemethod ', fontsize=14) ax3.set_xticks(x_pos) ax3.set_xticklabels(methods, rotation=15) ax3.legend(fontsize=10) ax3.grid(True, alpha=0.3, axis='y') plt.tight_layout() plt.savefig('optimization_methods_comparison.png', dpi=150, bbox_inches='tight') plt.show() 

from scipy.optimize import minimize, LinearConstraint, NonlinearConstraint # function: x^2 + y^2 def objective(x): return x[0]**2 + x[1]**2 def objective_grad(x): return np.array([2*x[0], 2*x[1]]) # : x + y >= 1 def constraint_func(x): return x[0] + x[1] def constraint_grad(x): return np.array([1.0, 1.0]) print("=" * 60) print("constrainedsuitable") print("=" * 60) # problem: min x^2 + y^2 s.t. x + y >= 1 x0 = np.array([0.0, 0.0]) # 1. print("\n1. : 0 <= x, y <= 2") bounds = [(0, 2), (0, 2)] result_bounds = minimize(objective, x0, method='L-BFGS-B', bounds=bounds, jac=objective_grad) print(f" suitableSolution: x* = {result_bounds.x}") print(f" minimum value: f(x*) = {result_bounds.fun:.6f}") # 2. print("\n2. : x + y >= 1") linear_constraint = LinearConstraint([[1, 1]], [1], [np.inf]) result_linear = minimize(objective, x0, method='trust-constr', constraints=linear_constraint, jac=objective_grad) print(f" suitableSolution: x* = {result_linear.x}") print(f" minimum value: f(x*) = {result_linear.fun:.6f}") print(f" : x + y = {np.sum(result_linear.x):.6f}") # 3. print("\n3. : x^2 + y^2 >= 1, x + y = 1") # : x^2 + y^2 >= 1 nonlinear_ineq = NonlinearConstraint( lambda x: x[0]**2 + x[1]**2, 1, np.inf ) # : x + y = 1 nonlinear_eq = NonlinearConstraint( lambda x: x[0] + x[1], 1, 1 ) x0_nonlinear = np.array([0.5, 0.5]) result_nonlinear = minimize(objective, x0_nonlinear, method='trust-constr', constraints=[nonlinear_ineq, nonlinear_eq], jac=objective_grad) print(f" suitableSolution: x* = {result_nonlinear.x}") print(f" minimum value: f(x*) = {result_nonlinear.fun:.6f}") print(f" 1: x^2 + y^2 = {np.sum(result_nonlinear.x**2):.6f}") print(f" 2: x + y = {np.sum(result_nonlinear.x):.6f}") # fig, axes = plt.subplots(1, 3, figsize=(16, 5)) # x = np.linspace(-0.5, 2.5, 200) y = np.linspace(-0.5, 2.5, 200) X, Y = np.meshgrid(x, y) Z = X**2 + Y**2 for idx, (ax, title) in enumerate(zip(axes, [ '1. ', '2. x+y≥1', '3. ' ])): # function cs = ax.contour(X, Y, Z, levels=15, cmap='viridis', alpha=0.6) ax.clabel(cs, inline=True, fontsize=8) # Visualization of if idx == 1: # x + y >= 1 ax.fill_between(x, 1-x, 3, alpha=0.2, color='red', label='') ax.plot(x, 1-x, 'r--', linewidth=2, label='x+y=1') if idx == 2: # x^2 + y^2 >= 1 theta = np.linspace(0, 2*np.pi, 100) ax.plot(np.cos(theta), np.sin(theta), 'r--', linewidth=2, label='x²+y²=1') ax.plot(x, 1-x, 'b--', linewidth=2, label='x+y=1') # suitable if idx == 0: ax.plot(result_bounds.x[0], result_bounds.x[1], 'r*', markersize=20, label='suitable') elif idx == 1: ax.plot(result_linear.x[0], result_linear.x[1], 'r*', markersize=20, label='suitable') else: ax.plot(result_nonlinear.x[0], result_nonlinear.x[1], 'r*', markersize=20, label='suitable') ax.set_xlabel('x', fontsize=12) ax.set_ylabel('y', fontsize=12) ax.set_title(title, fontsize=13) ax.legend(fontsize=10) ax.grid(True, alpha=0.3) ax.set_xlim(-0.5, 2.5) ax.set_ylim(-0.5, 2.5) plt.tight_layout() plt.savefig('constrained_optimization.png', dpi=150, bbox_inches='tight') plt.show() 

from scipy.interpolate import interp1d, CubicSpline, UnivariateSpline # （） np.random.seed(42) x_data = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) y_true = np.sin(x_data * 0.8) + 0.5 * x_data * 0.1 y_data = y_true + np.random.normal(0, 0.1, len(x_data)) print("=" * 60) print(" ") print("=" * 60) # x_fine = np.linspace(0, 10, 200) y_true_fine = np.sin(x_fine * 0.8) + 0.5 * x_fine * 0.1 # 1. f_linear = interp1d(x_data, y_data, kind='linear') y_linear = f_linear(x_fine) # 2. 3 cs = CubicSpline(x_data, y_data) y_cubic = cs(x_fine) # 3. （） spl_smooth = UnivariateSpline(x_data, y_data, s=0.5) # s: smoothing factor y_smooth = spl_smooth(x_fine) # 4. （） spl_exact = UnivariateSpline(x_data, y_data, s=0) y_exact = spl_exact(x_fine) print("\n:") print(f" : {len(x_data)}") print(f" : {len(x_fine)}") # fig, axes = plt.subplots(2, 2, figsize=(14, 10)) axes = axes.flatten() methods = [ ('', y_linear), ('3', y_cubic), (' (s=0.5)', y_smooth), (' (s=0)', y_exact) ] for idx, (ax, (method_name, y_interp)) in enumerate(zip(axes, methods)): ax.plot(x_fine, y_true_fine, 'g-', linewidth=1, alpha=0.5, label=' function') ax.plot(x_data, y_data, 'ro', markersize=8, label='') ax.plot(x_fine, y_interp, 'b-', linewidth=2, label=method_name) ax.set_xlabel('x', fontsize=12) ax.set_ylabel('y', fontsize=12) ax.set_title(method_name, fontsize=13) ax.legend(fontsize=10) ax.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('spline_interpolation.png', dpi=150, bbox_inches='tight') plt.show() # error print("\n function and error:") for method_name, y_interp in methods: mse = np.mean((y_interp - y_true_fine)**2) print(f" {method_name}: {mse:.6f}") 

# and print("=" * 60) print(": and ") print("=" * 60) # （3function） x_data_poly = np.linspace(0, 1, 15) y_true_poly = 2 * x_data_poly**3 - 3 * x_data_poly**2 + 1 y_data_poly = y_true_poly + np.random.normal(0, 0.1, len(x_data_poly)) x_fine_poly = np.linspace(0, 1, 200) y_true_fine_poly = 2 * x_fine_poly**3 - 3 * x_fine_poly**2 + 1 # degrees = [1, 3, 5, 10] fig, axes = plt.subplots(2, 2, figsize=(14, 10)) axes = axes.flatten() print("\n and error:") for idx, (ax, deg) in enumerate(zip(axes, degrees)): # coeffs = np.polyfit(x_data_poly, y_data_poly, deg) poly = np.poly1d(coeffs) y_fit = poly(x_fine_poly) # error calculation train_error = np.mean((poly(x_data_poly) - y_data_poly)**2) true_error = np.mean((y_fit - y_true_fine_poly)**2) print(f" {deg:2d}: error={train_error:.6f}, error={true_error:.6f}") # ax.plot(x_fine_poly, y_true_fine_poly, 'g-', linewidth=2, alpha=0.5, label=' function') ax.plot(x_data_poly, y_data_poly, 'ro', markersize=8, label='') ax.plot(x_fine_poly, y_fit, 'b-', linewidth=2, label=f'{deg}') ax.set_xlabel('x', fontsize=12) ax.set_ylabel('y', fontsize=12) ax.set_title(f'{deg}（ ）', fontsize=13) ax.legend(fontsize=10) ax.grid(True, alpha=0.3) ax.set_ylim(-1, 2) plt.tight_layout() plt.savefig('polynomial_overfitting.png', dpi=150, bbox_inches='tight') plt.show() print("\n:") print(" 1: （）") print(" 3: suitable（ function and ）") print(" 5-10: （、）") 

from scipy.fft import fft, fftfreq, ifft from scipy.signal import find_peaks # （multiple + ） fs = 1000 # (Hz) T = 1.0 # () N = int(fs * T) # t = np.linspace(0, T, N, endpoint=False) # : 50Hz + 120Hz + 250Hz + freq1, freq2, freq3 = 50, 120, 250 signal = (np.sin(2 * np.pi * freq1 * t) + 0.5 * np.sin(2 * np.pi * freq2 * t) + 0.3 * np.sin(2 * np.pi * freq3 * t)) noise = 0.2 * np.random.randn(N) signal_noisy = signal + noise print("=" * 60) print(" analysis") print("=" * 60) # FFT yf = fft(signal_noisy) xf = fftfreq(N, 1/fs)[:N//2] # power = 2.0/N * np.abs(yf[:N//2]) # peaks, properties = find_peaks(power, height=0.1) peak_freqs = xf[peaks] peak_powers = power[peaks] print(f"\n: {fs} Hz") print(f": {N}") print(f": {fs/N:.2f} Hz") print("\n was :") for freq, pwr in zip(peak_freqs, peak_powers): print(f" {freq:.1f} Hz (: {pwr:.3f})") # fig, axes = plt.subplots(3, 1, figsize=(14, 10)) # （ ） axes[0].plot(t[:500], signal[:500], 'b-', linewidth=1, alpha=0.7, label=' ') axes[0].set_xlabel(' [s]', fontsize=12) axes[0].set_ylabel('', fontsize=12) axes[0].set_title(' （）', fontsize=13) axes[0].legend(fontsize=10) axes[0].grid(True, alpha=0.3) # （） axes[1].plot(t[:500], signal_noisy[:500], 'r-', linewidth=1, alpha=0.7, label='') axes[1].set_xlabel(' [s]', fontsize=12) axes[1].set_ylabel('', fontsize=12) axes[1].set_title('', fontsize=13) axes[1].legend(fontsize=10) axes[1].grid(True, alpha=0.3) # （） axes[2].plot(xf, power, 'g-', linewidth=2, label='') axes[2].plot(peak_freqs, peak_powers, 'r^', markersize=10, label='') for freq, pwr in zip(peak_freqs, peak_powers): axes[2].annotate(f'{freq:.0f} Hz', xy=(freq, pwr), xytext=(freq+10, pwr+0.1), fontsize=10, arrowprops=dict(arrowstyle='->', color='red')) axes[2].set_xlabel(' [Hz]', fontsize=12) axes[2].set_ylabel('', fontsize=12) axes[2].set_title('FFT ', fontsize=13) axes[2].set_xlim(0, 400) axes[2].legend(fontsize=10) axes[2].grid(True, alpha=0.3) plt.tight_layout() plt.savefig('fft_frequency_analysis.png', dpi=150, bbox_inches='tight') plt.show() # （ ） print("\n" + "=" * 60) print("FFT （）") print("=" * 60) cutoff_freq = 300 # (Hz) yf_filtered = yf.copy() yf_filtered[np.abs(xf) > cutoff_freq] = 0 # # FFT signal_filtered = np.real(ifft(yf_filtered)) plt.figure(figsize=(14, 5)) plt.subplot(1, 2, 1) plt.plot(t[:500], signal_noisy[:500], 'r-', alpha=0.5, linewidth=1, label='') plt.plot(t[:500], signal[:500], 'b-', linewidth=2, label=' ') plt.xlabel(' [s]', fontsize=12) plt.ylabel('', fontsize=12) plt.title(' vs ', fontsize=13) plt.legend(fontsize=10) plt.grid(True, alpha=0.3) plt.subplot(1, 2, 2) plt.plot(t[:500], signal_filtered[:500], 'g-', linewidth=2, label='') plt.plot(t[:500], signal[:500], 'b--', alpha=0.5, linewidth=1, label=' ') plt.xlabel(' [s]', fontsize=12) plt.ylabel('', fontsize=12) plt.title(f'FFT（: {cutoff_freq} Hz）', fontsize=13) plt.legend(fontsize=10) plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('fft_noise_filtering.png', dpi=150, bbox_inches='tight') plt.show() # mse_noisy = np.mean((signal_noisy - signal)**2) mse_filtered = np.mean((signal_filtered - signal)**2) print(f"\nerror（MSE）:") print(f" : {mse_noisy:.6f}") print(f" : {mse_filtered:.6f}") print(f" : {mse_noisy / mse_filtered:.2f}times") 

from scipy.integrate import solve_ivp from scipy.optimize import differential_evolution # problem: suitable # : 、 and def cooling_ode(t, T, k): """ Newton : dT/dt = -k(T - T_ambient) T: [K] k: （） """ T_ambient = 300 # (K) return -k * (T - T_ambient) def hardness_model(T_min, cooling_rate): """ （） and """ # : 500K Ms = 500 if T_min > Ms: return 200 # （） else: # （） hardness = 400 + 100 * np.tanh(cooling_rate / 50) return hardness def residual_stress_model(cooling_rate): """ """ return 50 * cooling_rate**0.8 def evaluate_heat_treatment(k, target_hardness=450): """ function Parameters: ----------- k : float （） target_hardness : float Returns: -------- cost : float function（ do ） """ # T0 = 1000 # K（） # t_span = (0, 100) # 0-100 sol = solve_ivp(cooling_ode, t_span, [T0], args=(k,), dense_output=True, max_step=0.5) # T_min = np.min(sol.y[0]) # （） cooling_rate = np.mean(np.abs(np.diff(sol.y[0]) / np.diff(sol.t))) # and calculation hardness = hardness_model(T_min, cooling_rate) stress = residual_stress_model(cooling_rate) # function: + cost = (hardness - target_hardness)**2 + 0.1 * stress**2 return cost, hardness, stress, sol print("=" * 60) print(" suitable") print("=" * 60) # suitable（suitable: Differential Evolution） print("\nsuitable...") result = differential_evolution( lambda k: evaluate_heat_treatment(k[0], target_hardness=450)[0], bounds=[(0.01, 0.5)], # maxiter=50, seed=42 ) k_optimal = result.x[0] cost_optimal, hardness_optimal, stress_optimal, sol_optimal = \ evaluate_heat_treatment(k_optimal) print(f"\nsuitable:") print(f" suitable: k* = {k_optimal:.4f}") print(f" : {hardness_optimal:.2f} HV (: 450 HV)") print(f" : {stress_optimal:.2f} MPa") print(f" function: {cost_optimal:.4f}") # : suitable k_slow = 0.05 k_fast = 0.3 cost_slow, hardness_slow, stress_slow, sol_slow = evaluate_heat_treatment(k_slow) cost_fast, hardness_fast, stress_fast, sol_fast = evaluate_heat_treatment(k_fast) print("\n:") print(f" slow (k={k_slow}):") print(f" : {hardness_slow:.2f} HV, : {stress_slow:.2f} MPa") print(f" (k={k_fast}):") print(f" : {hardness_fast:.2f} HV, : {stress_fast:.2f} MPa") # fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # ax1 = axes[0, 0] t_plot = np.linspace(0, 100, 500) ax1.plot(t_plot, sol_slow.sol(t_plot)[0], 'b-', linewidth=2, label=f'slow (k={k_slow})') ax1.plot(t_plot, sol_optimal.sol(t_plot)[0], 'g-', linewidth=2, label=f'suitable (k={k_optimal:.3f})') ax1.plot(t_plot, sol_fast.sol(t_plot)[0], 'r-', linewidth=2, label=f' (k={k_fast})') ax1.axhline(y=500, color='k', linestyle='--', alpha=0.5, label=' Ms') ax1.set_xlabel(' [s]', fontsize=12) ax1.set_ylabel(' [K]', fontsize=12) ax1.set_title(' ', fontsize=13) ax1.legend(fontsize=10) ax1.grid(True, alpha=0.3) # and ax2 = axes[0, 1] k_range = np.linspace(0.01, 0.5, 50) hardness_range = [] stress_range = [] for k in k_range: _, h, s, _ = evaluate_heat_treatment(k) hardness_range.append(h) stress_range.append(s) ax2.plot(stress_range, hardness_range, 'b-', linewidth=2) ax2.plot(stress_optimal, hardness_optimal, 'r*', markersize=20, label='suitable') ax2.axhline(y=450, color='g', linestyle='--', alpha=0.5, label='') ax2.set_xlabel(' [MPa]', fontsize=12) ax2.set_ylabel(' [HV]', fontsize=12) ax2.set_title(' vs （）', fontsize=13) ax2.legend(fontsize=10) ax2.grid(True, alpha=0.3) # and ax3 = axes[1, 0] ax3.plot(k_range, hardness_range, 'b-', linewidth=2) ax3.plot(k_optimal, hardness_optimal, 'r*', markersize=20, label='suitable') ax3.axhline(y=450, color='g', linestyle='--', alpha=0.5, label='') ax3.set_xlabel(' k', fontsize=12) ax3.set_ylabel(' [HV]', fontsize=12) ax3.set_title(' and ', fontsize=13) ax3.legend(fontsize=10) ax3.grid(True, alpha=0.3) # and function ax4 = axes[1, 1] cost_range = [] for k in k_range: c, _, _, _ = evaluate_heat_treatment(k) cost_range.append(c) ax4.plot(k_range, cost_range, 'b-', linewidth=2) ax4.plot(k_optimal, cost_optimal, 'r*', markersize=20, label='minimum value') ax4.set_xlabel(' k', fontsize=12) ax4.set_ylabel('function', fontsize=12) ax4.set_title('function ', fontsize=13) ax4.legend(fontsize=10) ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('heat_treatment_optimization.png', dpi=150, bbox_inches='tight') plt.show() 

# suitableproblem # : A → B → C # : B 、 def reactor_model(params): """ （CSTR） Parameters: ----------- params : array [T, tau, C_A0] T: [K] tau: [min] C_A0: A [mol/L] Returns: -------- C_B : float B [mol/L] cost : float [$/h] """ T, tau, C_A0 = params # （Arrhenius） k1_ref = 0.5 # A → B @ 350K k2_ref = 0.2 # B → C @ 350K Ea1 = 50000 # [J/mol] Ea2 = 60000 R = 8.314 # k1 = k1_ref * np.exp(-Ea1/R * (1/T - 1/350)) k2 = k2_ref * np.exp(-Ea2/R * (1/T - 1/350)) # （ ） # dC_A/dt = 0 = (C_A0 - C_A)/tau - k1*C_A # dC_B/dt = 0 = -C_B/tau + k1*C_A - k2*C_B C_A = C_A0 / (1 + k1 * tau) C_B = k1 * tau * C_A / (1 + k2 * tau) # # （ ）+ （ ） heating_cost = 0.01 * (T - 300)**2 # $/h reactor_cost = 5 * tau # $/h feedstock_cost = 10 * C_A0 # $/h total_cost = heating_cost + reactor_cost + feedstock_cost return C_B, total_cost def multi_objective_cost(params, weight_yield=0.7, weight_cost=0.3): """ suitable function 1: B （-C_B ） 2: """ C_B, cost = reactor_model(params) # （） C_B_normalized = C_B / 2.0 # 2 mol/L and cost_normalized = cost / 100 # 100 $/h and return -weight_yield * C_B_normalized + weight_cost * cost_normalized print("=" * 60) print(" suitable") print("=" * 60) # # : 300K ≤ T ≤ 400K（） # : 1 ≤ tau ≤ 20 min（） # : 0.5 ≤ C_A0 ≤ 5.0 mol/L（） bounds = [ (300, 400), # T [K] (1, 20), # tau [min] (0.5, 5.0) # C_A0 [mol/L] ] # suitable（suitable） print("\nsuitable（: 70%, : 30%）...") result = differential_evolution( multi_objective_cost, bounds, args=(0.7, 0.3), maxiter=100, seed=42 ) T_opt, tau_opt, C_A0_opt = result.x C_B_opt, cost_opt = reactor_model(result.x) print(f"\nsuitable:") print(f" : T* = {T_opt:.2f} K") print(f" : τ* = {tau_opt:.2f} min") print(f" : C_A0* = {C_A0_opt:.2f} mol/L") print(f" B: {C_B_opt:.4f} mol/L") print(f" : {cost_opt:.2f} $/h") # calculation（ and ） print("\n" + "=" * 60) print(" calculation（ vs ）") print("=" * 60) weights = np.linspace(0, 1, 11) pareto_C_B = [] pareto_cost = [] pareto_T = [] for w_yield in weights: w_cost = 1 - w_yield result_pareto = differential_evolution( multi_objective_cost, bounds, args=(w_yield, w_cost), maxiter=50, seed=42 ) C_B, cost = reactor_model(result_pareto.x) pareto_C_B.append(C_B) pareto_cost.append(cost) pareto_T.append(result_pareto.x[0]) # fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # ax1 = axes[0, 0] ax1.plot(pareto_cost, pareto_C_B, 'bo-', linewidth=2, markersize=8) ax1.plot(cost_opt, C_B_opt, 'r*', markersize=20, label='') ax1.set_xlabel(' [$/h]', fontsize=12) ax1.set_ylabel('B [mol/L]', fontsize=12) ax1.set_title('（ vs ）', fontsize=13) ax1.legend(fontsize=10) ax1.grid(True, alpha=0.3) # ax2 = axes[0, 1] T_range = np.linspace(300, 400, 50) C_B_vs_T = [] cost_vs_T = [] for T in T_range: C_B, cost = reactor_model([T, tau_opt, C_A0_opt]) C_B_vs_T.append(C_B) cost_vs_T.append(cost) ax2_twin = ax2.twinx() ax2.plot(T_range, C_B_vs_T, 'b-', linewidth=2, label='B') ax2_twin.plot(T_range, cost_vs_T, 'r-', linewidth=2, label='') ax2.axvline(x=T_opt, color='g', linestyle='--', alpha=0.7, label='suitable') ax2.set_xlabel(' [K]', fontsize=12) ax2.set_ylabel('B [mol/L]', fontsize=12, color='b') ax2_twin.set_ylabel(' [$/h]', fontsize=12, color='r') ax2.set_title(' ', fontsize=13) ax2.tick_params(axis='y', labelcolor='b') ax2_twin.tick_params(axis='y', labelcolor='r') ax2.grid(True, alpha=0.3) # ax3 = axes[1, 0] tau_range = np.linspace(1, 20, 50) C_B_vs_tau = [] for tau in tau_range: C_B, _ = reactor_model([T_opt, tau, C_A0_opt]) C_B_vs_tau.append(C_B) ax3.plot(tau_range, C_B_vs_tau, 'b-', linewidth=2) ax3.axvline(x=tau_opt, color='r', linestyle='--', linewidth=2, label='suitable') ax3.set_xlabel(' [min]', fontsize=12) ax3.set_ylabel('B [mol/L]', fontsize=12) ax3.set_title(' ', fontsize=13) ax3.legend(fontsize=10) ax3.grid(True, alpha=0.3) # ax4 = axes[1, 1] C_A0_range = np.linspace(0.5, 5.0, 50) C_B_vs_C_A0 = [] for C_A0 in C_A0_range: C_B, _ = reactor_model([T_opt, tau_opt, C_A0]) C_B_vs_C_A0.append(C_B) ax4.plot(C_A0_range, C_B_vs_C_A0, 'b-', linewidth=2) ax4.axvline(x=C_A0_opt, color='r', linestyle='--', linewidth=2, label='suitable') ax4.set_xlabel(' C_A0 [mol/L]', fontsize=12) ax4.set_ylabel('B [mol/L]', fontsize=12) ax4.set_title(' ', fontsize=13) ax4.legend(fontsize=10) ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('reactor_optimization.png', dpi=150, bbox_inches='tight') plt.show() print("\n:") print(" - and 、") print(" - do and do 、 ") print(" - than 、 and ")