Scipy优化函数求MSE最小值时参数数组未更新问题求助
2026-8-11
问题描述
使用Scipy寻找函数最小值优化粘弹性模型，目标是找到使均方误差（MSE）最小的数组参数g和tau，但优化仅迭代1次就终止，完全未修改初始猜测值（g=[4,4,4,4,4]，tau=[0.8,0.8,0.8,0.8,0.8]）。
原代码
import numpy as np import pandas as pd import scipy.optimize as spo ## Insert the hihest numner of Prony Parameters for optimization.. N_Prony = 5 #########Intial_guess of g_i and tau_i################ #### Define Array for the guessing ######## Prony_0 = np.ones(N_Prony*2) ## Five for the relaxtion time (taui) and five for g_i # Initial guess of g_i and Tau_i gini = 4 Taui = 0.8 Prony_0[0:N_Prony] = Prony_0[0:N_Prony]*gini Prony_0[N_Prony:2*N_Prony] = Prony_0[N_Prony:2*N_Prony]*Taui g = Prony_0[0:N_Prony] ## The first part of the array is g tau = Prony_0[N_Prony:2*N_Prony] # The second part of the array is Tau df1 = pd.read_excel(r'C:\Users\Mahmoud Khadijeh\Desktop\DSR Application\Testdata_Einf.xlsx') ## Read the data from Excel file w = df1.iloc[:,1] ## Read the frequency from the Excel file E_INF = df1.iloc[4,2]; NU = df1.iloc[0,5] ## Read E_INF & Poission's Ratio from the EXCEL FILE G_INF = (E_INF)/2*(1+NU) # Calculate G_INF from G0 = G_INF/(1-sum(g)) # Calculate G0 from G_INF TANW_MEAS = (df1.iloc[:,3])/(df1.iloc[:,2]) # Degree of Viscoelasticity list1 = [] # This list is to store G' from the loop in an array {For each Frequnecy} list2 = [] # This list is to store G'' from the loop in an array {For each Frequnecy} ## Calculation.. for K in range(len(w)): #Second part of Equation 5 --> G' for L in range(N_Prony): GPrime_1 = G_INF + G0*((g[L]*((tau[L])**2)*(w[K])**2)/(1+((tau[L]**2))*w[K]**2)*N_Prony) list1.append(GPrime_1) print(list1[0]) df1["G'"] = list1 df1["E'"] = np.dot(2*(1+NU),list1) #Convert G' to E' and add it to the table print('df1 = ', df1) for J in range(len(w)): #Second part of Equation 5 --> G' for i in range(N_Prony): GPrime_2 = G0*((g[i]*((tau[i]))*w[J])/(1+((tau[i]**2))*w[J]**2)*N_Prony) list2.append(GPrime_2) print(list2[0]) df1["G''"] = list2 df1["E''"] = np.dot(2*(1+NU),list2) #Convert G'' to E'' and add it to the table #### Initial Guess array x0 = np.array([g, tau]) def objective(SE): global new_df g = SE[0] # Variable 1 that we have to optimize tau = SE[1] # Variable 2 that we have to optimize print('g::', g) print('tau::', tau) new_df = pd.DataFrame() new_df["E'_meas"] = df1.iloc[:,2] new_df["E''_meas"] =df1.iloc[:,3] #list1 = new_df["E'_meas"] #list2 = new_df["E''_meas"] new_df["E'_cal"] = list1 # Where list1 is E' new_df["E''_cal"] = list2 # Where list2 is E'' new_df["Tan(d)_meas"] = TANW_MEAS new_df["Tan(d)_cal"] = new_df["E''_meas"]/new_df["E'_meas"] MSE = np.square(np.subtract(new_df["E'_meas"],new_df["E'_cal"])).mean() #minimize = (((new_df["E'_meas"] - new_df["E'_cal"])**2)/np.std(new_df["E'_meas"])) + \ # (((new_df["E''_meas"] - new_df["E''_cal"])**2)/np.std(new_df["E''_meas"])) + \ # (((new_df["Tan(d)_meas"] - new_df["Tan(d)_cal"])**2)/np.std(new_df["Tan(d)_cal"])) return MSE sol = spo.minimize(objective, x0, method='SLSQP', options={'disp': True}) print(sol)
运行输出
Optimization terminated successfully (Exit mode 0) Current function value: 3714530.31378857 Iterations: 1 Function evaluations: 11 Gradient evaluations: 1 fun: 3714530.31378857 jac: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) message: 'Optimization terminated successfully' nfev: 11 nit: 1 njev: 1 status: 0 success: True x: array([4. , 4. , 4. , 4. , 4. , 0.8, 0.8, 0.8, 0.8, 0.8])
问题根源与修正方案
核心问题
- 目标函数未关联优化变量：
objective函数直接使用全局的list1和list2作为计算值，这两个数组是用初始g/tau计算的，和优化输入参数SE无关。导致MSE对所有优化变量的梯度为0，Scipy认为当前点就是最小值，直接终止迭代。 - 初始猜测数组形状错误：
x0 = np.array([g, tau])生成(2,5)的二维数组，但scipy.optimize.minimize要求输入必须是一维数组。 - Prony项计算逻辑错误：双层循环计算G'和G''时，每次迭代都覆盖结果变量，最终只保留最后一个Prony项的贡献，没有累加所有项的和，公式实现完全错误。
- G0计算未随迭代更新：G0依赖g的和，应该在每次优化迭代时用当前的g值重新计算，而非初始值。
- 公式括号错误：原代码中
G_INF = (E_INF)/2*(1+NU)的括号位置错误，正确公式应为G_INF = E_INF / (2*(1+NU))。
修正后的代码
import numpy as np import pandas as pd import scipy.optimize as spo N_Prony = 5 # 读取数据 df1 = pd.read_excel(r'C:\Users\Mahmoud Khadijeh\Desktop\DSR Application\Testdata_Einf.xlsx') w = df1.iloc[:,1].values E_INF = df1.iloc[4,2] NU = df1.iloc[0,5] G_INF = E_INF / (2 * (1 + NU)) # 修正公式括号位置 E_prime_meas = df1.iloc[:,2].values E_double_prime_meas = df1.iloc[:,3].values # 初始猜测：转为一维数组 gini = 4 Taui = 0.8 x0 = np.concatenate([np.full(N_Prony, gini), np.full(N_Prony, Taui)]) def calculate_viscoelastic_params(g, tau, w, G_INF, NU): """根据当前g和tau计算E'和E''""" sum_g = np.sum(g) if np.isclose(1 - sum_g, 0): raise ValueError("sum(g)接近1，G0会无穷大") G0 = G_INF / (1 - sum_g) # 向量化计算所有频率的G'和G''，替代循环 tau_sq = tau ** 2 w_sq = w[:, np.newaxis] ** 2 denom = 1 + tau_sq * w_sq # G' = G_INF + G0 * sum(g_i * tau_i² w² / (1 + tau_i² w²)) G_prime = G_INF + G0 * np.sum(g * tau_sq * w_sq / denom, axis=1) # G'' = G0 * sum(g_i * tau_i w / (1 + tau_i² w²)) G_double_prime = G0 * np.sum(g * tau * w[:, np.newaxis] / denom, axis=1) # 转换为E'和E'' E_prime = 2 * (1 + NU) * G_prime E_double_prime = 2 * (1 + NU) * G_double_prime return E_prime, E_double_prime def objective(SE): """目标函数：计算E'的MSE""" g = SE[:N_Prony] tau = SE[N_Prony:] # 用当前优化参数计算预测值 E_prime_cal, E_double_prime_cal = calculate_viscoelastic_params(g, tau, w, G_INF, NU) # 计算MSE mse = np.mean((E_prime_meas - E_prime_cal) ** 2) return mse # 添加约束：sum(g) < 1（避免G0无穷大） constraints = ({'type': 'ineq', 'fun': lambda x: 1 - np.sum(x[:N_Prony]) - 1e-6}) # 加小偏移避免数值问题 # 执行优化 sol = spo.minimize(objective, x0, method='SLSQP', constraints=constraints, options={'disp': True}) print("优化结果：") print(f"最优g值: {sol.x[:N_Prony]}") print(f"最优tau值: {sol.x[N_Prony:]}") print(f"最小MSE: {sol.fun}")
关键修改说明
- 目标函数关联优化变量：每次迭代时用输入的
SE参数重新计算E'和E''，确保MSE随参数变化。 - 一维初始猜测数组：用
np.concatenate将g和tau合并为10元素的一维数组，符合Scipy要求。 - 向量化计算替代循环：用NumPy向量化操作实现Prony项的累加，既修正了原循环的逻辑错误，又提升计算效率。
- 添加约束条件：限制
sum(g) < 1，避免G0计算出现无穷大的情况。 - 修正公式括号：调整G_INF的计算式，保证粘弹性参数转换的正确性。
内容的提问来源于stack exchange，提问作者Mhmoud Khadija
相关产品推荐
相关产品推荐
