如何解决参数估计和最小二乘最小化 - Python 和 Gekko
我正在尝试估计 4 个参数（kdoc1、kdoc2、S1、S2），使用 GEKKO 求解代数常方程组，最小化观察值和预测值之间残差的平方和。我在测量数据和预测数据之间没有得到很好的调整，你对我在代码中做错了什么有什么建议吗？
此外，我还尝试添加另一组测量数据。我希望模型使用两个测量数据集来估计参数，我需要使用两个数据集最小化残差之和。你知道如何在代码中做到这一点吗？
我已经添加了额外的点来帮助求解器使用建议 python - Infrequent measurements in Gekko with extra simulation points - Stack Overflow
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import math as math
import pandas as pd
#TOC=5.93 mg/L
#M2 C1 5mg/L
#time: [0,0.08333,0.5,1,4,22.61667]; C2 measured: [0,9.33e-5,8.55e-5,7.77e-5,7.00e-5]
#M4 C1 20mg/l
#time: [0,2.92e-4,2.72e-4,2.14e-4]
#measurements
t_data = [0.08333,22.61667]
x_data = [9.33e-5,7.00e-5]
#x_data2 = [2.92e-4,2.14e-4]
df1 = pd.DataFrame({'time':t_data,'x':x_data})
df1.set_index('time',inplace=True)
# simulation time points
df2 = pd.DataFrame({'time':np.linspace(0,25,51)}) #(0,51)
df2.set_index('time',inplace=True)
# merge dataframes
df = df2.join(df1,how='outer')
# get True (1) or False (0) for measurement
df['meas'] = (df['x'].values==df['x'].values).astype(int)
# replace NaN with zeros
df0 = df.fillna(value=0)
#Estimator Model
m = GEKKO()
m.time = df.index.values
# adjustable parameters
kdoc1 = m.FV() #m-1h-1
kdoc2 = m.FV() #m-1h-1
S1 = m.FV()
S2 = m.FV()
#State variables
# M (mol/L)
C0 = m.Var(value=8.00e-5)
C1 = m.Var(value=1.36e-3)
C2 = m.Var(value=x_data[0])
C2m = m.Param(df0['x'].values)
meas = m.Param(df0['meas'].values)
m.Minimize(meas*(C2-C2m)**2)
C3 = m.Var(value=0)
C4 = m.Var(value=3.16e-8)
C5 = m.Var(value=3.83e-7)
C6 = m.Var(value=0)
C7 = m.Var(value=0)
C8 = m.Var(value=0)
C9 = m.Var(value=0)
C10 = m.Var(value=0)
C11 = m.Var(value=0)
C12 = m.Var(value=3.21e-3)
TOC = m.Var(value=5.93)#mg-C/L
cC1 = m.Var(value=0)
cC2 = m.Var(value=0)
#temperature ºC
T = 20
#Rate constants
k1 = 2040000000*math.exp(-1887/(273.15+T))*3600
k2 = 138000000*math.exp(-8800/(273.15+T))*3600
k3 = 300000*math.exp(-2010/(273.15+T))*3600
k4 = 0.00000065*3600
k6 = 26700*3600
k7 = 167*3600
k8 = 27700*3600
k9 = 8300*3600
k10 = 0
k5 = 37800000000*math.exp(-2169/(273.15+T))*C4+2.52e25*math.exp(-16860/(273.15+T))*C10+0.87*math.exp(-503/(273.15+T))*C9
r1 = k1 * C0 * C1
r2 = k2 * C2
r3 = k3 * C0 * C2
r4 = k4 * C3
r5 = k5 * C2 * C2
r6 = k6 * C3 * C1* C4
r7 = k7 * C3 * C5
r8 = k8 * C6 * C3
r9 = k9 * C6 * C2
r10 = k10 * C2 * C3
r11 = kdoc1*S1*TOC*C2/12000
r12 = kdoc2*S2*TOC*C0/12000
t = m.Param(value=m.time)
m.Equation(C0.dt()== -r1 + r2 - r3 + r4 + r8 - r12)
m.Equation(C1.dt()== -r1 + r2 + r5 - r6 + r11)
m.Equation(C2.dt()== r1 - r2 - r3 + r4 - r5 + r6 - r9 - r10 - r11)
m.Equation(C3.dt()== r3 - r4 + r5 - r6 - r7 - r8 - r10)
m.Equation(C4.dt()== 0)
m.Equation(C5 == 1e-14/C4)
m.Equation(C6.dt()== r7 - r8 - r9)
m.Equation(C7 == (3.16e-8*C0)/C4)
m.Equation(C1 == (5.01e-10*C8)/C4)
m.Equation(C11 == (5.01e-11*C10)/C4)
m.Equation(C9 == (5.01e-7*C9)/C4)
m.Equation(C10 == C12 - 2*C11 - C5 + C4)
m.Equation(C12.dt()== 0)
m.Equation(TOC.dt()== 0)
m.Equation(cC1 == 17000*C1)
m.Equation(cC2 == 51500*C2)
#Application options
m.options.soLVER = 1 #APOPT solver
m.options.IMODE = 5 #Dynamic Simultaneous - estimation
m.options.EV_TYPE = 2 #absolute error
m.options.NODES = 5 #collocation nodes (2,5)
if True:
kdoc1.STATUS=1
kdoc2.STATUS=1
S1.STATUS=1
S2.STATUS=1
#Solve
m.solve(disp=False)
#Results
data={'Time (h)':t,'C0 (M)':C0.value,'C1 (M)':C1.value,'C2 (M)':C2.value,'C3 (M)':C3.value,'C4 (M)':C4.value,'C5 (M)':C5.value,'C6 (M)':C6.value,'C7 (M)':C7.value,'C8 (M)':C8.value,'C9 (M)':C9.value,'C10 (M)':C10.value,'C11 (M)':C11.value,'c12 (M)':C12.value}
dfr = pd.DataFrame(data,columns=['Time (h)','C0 (M)','C1 (M)','C2 (M)','C3 (M)','C4 (M)','C5 (M)','C6 (M)','C7 (M)','C8 (M)','C9 (M)','C10 (M)','C11 (M)','C12 (M)'])
#dfr.to_csv(r'C:\dfestimation.csv',index = False)
print('Solution')
print('kdoc1 = ' + str(kdoc1.value[0]))
print('kdoc2 = ' + str(kdoc2.value[0]))
print('S1 = ' + str(S1.value[0]))
print('S2 = ' + str(S2.value[0]))
plt.figure(1,figsize=(12,8))
plt.subplot(3,1)
plt.plot(m.time,C2.value,'bo',label='Predicted')
plt.plot(m.time,df['x'].values,'rs',label='Meas')
plt.plot(m.time,label ='C2')
plt.legend(loc='best')
plt.subplot(3,2)
plt.plot(m.time,cC2.value,label ='C2')
plt.plot(m.time,cC1.value,label ='C1')
plt.legend(loc='best')
plt.subplot(3,3)
plt.plot(m.time,C0.value,label ='C0')
plt.plot(m.time,C1.value,label ='C1')
plt.plot(m.time,C3.value,label ='C3')
plt.plot(m.time,C4.value,label ='C4')
plt.plot(m.time,C5.value,label ='C5')
plt.plot(m.time,C6.value,label ='C6')
plt.plot(m.time,C7.value,label ='C7-')
plt.plot(m.time,C8.value,label ='C8')
plt.plot(m.time,C9.value,label ='C9')
plt.plot(m.time,C10.value,label ='C10')
plt.plot(m.time,C11.value,label ='C11')
plt.plot(m.time,C12.value,label ='C12')
plt.xlabel('time (h)')
plt.ylabel('Concentration (mol/L)')
plt.legend(loc='best')
解决方法
以下是更接近的新结果。
我向参数添加了约束，并为它们提供了非零起始值 (1)。
kdoc1 = m.FV(1,lb=0.01,ub=10) #m-1h-1
kdoc2 = m.FV(1,ub=10) #m-1h-1
S1 = m.FV(1,ub=10)
S2 = m.FV(1,ub=10)
计算 C0 和 C1 的初始条件，因为它们在 C2 中产生了巨大的瞬时跳跃。如果这些初始浓度已知，那么您可能需要切换回固定的初始条件。
m.free_initial(C0)
m.free_initial(C1)
缩放目标函数，因为 1e-4**2 已经是一个很小的数字，可能已经满足客观收敛标准。原问题的目标函数值很小。我让每个平方值大约等于 1 作为起点。
m.Minimize(meas*(1e4*(C2-C2m))**2)
这提供了一个合理的解决方案，您可以开始使用它来获得您需要的回归结果。
kdoc1 = 3.7210285146
kdoc2 = 2.5743961631
S1 = 3.7210285146
S2 = 2.5743961631
参数值不会在约束条件下结束，因此这意味着该解决方案是最佳的。即使参数约束最后没有激活，约束确实有助于指导求解器在适当的搜索区域。
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import math as math
import pandas as pd
#measurements
t_data = [0.08333,0.5,1,4,22.61667]
x_data = [9.33e-5,8.55e-5,7.77e-5,7.00e-5]
#x_data2 = [2.92e-4,2.72e-4,2.92e-4,2.14e-4]
df1 = pd.DataFrame({'time':t_data,'x':x_data})
df1.set_index('time',inplace=True)
# simulation time points
df2 = pd.DataFrame({'time':np.linspace(0,25,51)})
df2.set_index('time',inplace=True)
# merge dataframes
df = df2.join(df1,how='outer')
# get True (1) or False (0) for measurement
df['meas'] = (df['x'].values==df['x'].values).astype(int)
# replace NaN with zeros
df0 = df.fillna(value=0)
#Estimator Model
m = GEKKO()
m.time = df.index.values
# adjustable parameters
kdoc1 = m.FV(1,ub=10)
#State variables
# M (mol/L)
C0 = m.Var(value=8.00e-5,lb=0,ub=1e-3)
m.free_initial(C0)
C1 = m.Var(value=1.36e-3,ub=2e-3)
m.free_initial(C1)
C2 = m.Var(value=x_data[0])
C2m = m.Param(df0['x'].values,lb=0)
meas = m.Param(df0['meas'].values)
m.Minimize(meas*(1e4*(C2-C2m))**2)
C3 = m.Var(value=0)
C4 = m.Var(value=3.16e-8)
C5 = m.Var(value=3.83e-7)
C6 = m.Var(value=0)
C7 = m.Var(value=0)
C8 = m.Var(value=0)
C9 = m.Var(value=0)
C10 = m.Var(value=0)
C11 = m.Var(value=0)
C12 = m.Var(value=3.21e-3)
TOC = m.Var(value=5.93)#mg-C/L
cC1 = m.Var(value=0)
cC2 = m.Var(value=0)
#temperature ºC
T = 20
#Rate constants
k1 = 2040000000*math.exp(-1887/(273.15+T))*3600
k2 = 138000000*math.exp(-8800/(273.15+T))*3600
k3 = 300000*math.exp(-2010/(273.15+T))*3600
k4 = 0.00000065*3600
k6 = 26700*3600
k7 = 167*3600
k8 = 27700*3600
k9 = 8300*3600
k10 = 0
k5 = 37800000000*math.exp(-2169/(273.15+T))*C4\
+2.52e25*math.exp(-16860/(273.15+T))*C10\
+0.87*math.exp(-503/(273.15+T))*C9
r1 = k1 * C0 * C1
r2 = k2 * C2
r3 = k3 * C0 * C2
r4 = k4 * C3
r5 = k5 * C2 * C2
r6 = k6 * C3 * C1* C4
r7 = k7 * C3 * C5
r8 = k8 * C6 * C3
r9 = k9 * C6 * C2
r10 = k10 * C2 * C3
r11 = kdoc1*S1*TOC*C2/12000
r12 = kdoc2*S2*TOC*C0/12000
t = m.Param(value=m.time)
m.Equation(C0.dt()== -r1 + r2 - r3 + r4 + r8 - r12)
m.Equation(C1.dt()== -r1 + r2 + r5 - r6 + r11)
m.Equation(C2.dt()== r1 - r2 - r3 + r4 - r5 + r6 - r9 - r10 - r11)
m.Equation(C3.dt()== r3 - r4 + r5 - r6 - r7 - r8 - r10)
m.Equation(C4.dt()== 0)
m.Equation(C5 == 1e-14/C4)
m.Equation(C6.dt()== r7 - r8 - r9)
m.Equation(C7 == (3.16e-8*C0)/C4)
m.Equation(C1 == (5.01e-10*C8)/C4)
m.Equation(C11 == (5.01e-11*C10)/C4)
m.Equation(C9 == (5.01e-7*C9)/C4)
m.Equation(C10 == C12 - 2*C11 - C5 + C4)
m.Equation(C12.dt()== 0)
m.Equation(TOC.dt()== 0)
m.Equation(cC1 == 17000*C1)
m.Equation(cC2 == 51500*C2)
#Application options
m.options.SOLVER = 3 #APOPT solver
m.options.IMODE = 5 #Dynamic Simultaneous - estimation
m.options.EV_TYPE = 2 #absolute error
m.options.NODES = 3 #collocation nodes (2,5)
if True:
kdoc1.STATUS=1
kdoc2.STATUS=1
S1.STATUS=1
S2.STATUS=1
#Solve
m.solve(disp=True)
print('Solution')
print('kdoc1 = ' + str(kdoc1.value[0]))
print('kdoc2 = ' + str(kdoc2.value[0]))
print('S1 = ' + str(S1.value[0]))
print('S2 = ' + str(S2.value[0]))
plt.figure(1,figsize=(8,5))
plt.subplot(2,1)
plt.plot(m.time,C2.value,'bo',label='Predicted')
plt.plot(m.time,df['x'].values,'rs',label='Meas')
plt.plot(m.time,label ='C2')
plt.legend(loc='best')
plt.subplot(2,2)
plt.plot(m.time,cC2.value,label ='C2')
plt.plot(m.time,cC1.value,label ='C1')
plt.legend(loc='best')
plt.xlabel('time (h)')
plt.figure(2,figsize=(12,8))
plt.subplot(4,3,C0.value,label ='C0')
plt.legend()
plt.subplot(4,C1.value,label ='C1')
plt.legend()
plt.subplot(4,3)
plt.plot(m.time,C3.value,label ='C3')
plt.legend()
plt.subplot(4,4)
plt.plot(m.time,C4.value,label ='C4')
plt.legend()
plt.subplot(4,5)
plt.plot(m.time,C5.value,label ='C5')
plt.legend()
plt.subplot(4,6)
plt.plot(m.time,C6.value,label ='C6')
plt.legend()
plt.subplot(4,7)
plt.plot(m.time,C7.value,label ='C7-')
plt.legend()
plt.subplot(4,8)
plt.plot(m.time,C8.value,label ='C8')
plt.legend()
plt.subplot(4,9)
plt.plot(m.time,C9.value,label ='C9')
plt.legend()
plt.xlabel('time (h)')
plt.subplot(4,10)
plt.plot(m.time,C10.value,label ='C10')
plt.legend()
plt.xlabel('time (h)')
plt.subplot(4,11)
plt.plot(m.time,C11.value,label ='C11')
plt.legend()
plt.xlabel('time (h)')
plt.subplot(4,12)
plt.plot(m.time,C12.value,label ='C12')
plt.legend()
plt.xlabel('time (h)')
plt.show()
@John Hedengren，感谢您的支持。我按照您的想法修改了我的代码，并按照您在此处的解释添加了另一组数据：
How to set up GEKKO for parameter estimation from multiple independent sets of data?
但是，我没有找到一个好的解决方案，特别是对于测量 4。您对如何更好地调整它有什么想法吗？感谢您的时间和帮助。
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import math as math
import pandas as pd
#data set 1
t_data1 = [0.08333,22.61667]
x_data1 = [7.77e-5,7e-5,6.22e-5,3.89e-5]
x_data1mgl = [4,3.6,3.2,2]
#data set 2
t_data4 = [0.08333,22.61667]
x_data4 = [2.92e-4,2.14e-4]
x_data4mgl = [15,14,15,11]
#combine data and time in the same dataframe
data1 = pd.DataFrame({'time':t_data1,'x1':x_data1})
data4 = pd.DataFrame({'time':t_data4,'x4':x_data4})
data1.set_index('time',inplace=True)
data4.set_index('time',inplace=True)
# merge dataframes
data = data1.join(data4,how='outer')
# simulation time points
dftp = pd.DataFrame({'time':np.linspace(0,51)})
dftp.set_index('time',inplace=True)
# merge dataframes
df = data.join(dftp,how='outer')
# get True (1) or False (0) for measurement
#df['meas'] = (df['x'].values==df['x'].values).astype(int)
z1 = (df['x1']==df['x1']).astype(int)
z4 = (df['x4']==df['x4']).astype(int)
# replace NaN with zeros
df0 = df.fillna(value=0)
#Estimator Model
m = GEKKO(remote=True)
#m = GEKKO(remote=False) # remote=False to produce local folder with results
m.time = df.index.values
# measurements - Gekko arrays to store measured data
xm = m.Array(m.Param,2)
xm[0].value = df0['x1'].values
xm[1].value = df0['x4'].values
cl=m.Array(m.Var,2)
cl[0].value= 4e-5
cl[1].value= 1.33e-4
TOC=m.Array(m.Var,2)
TOC[0].value= 11.85
TOC[1].value= 11.85
C0=m.Array(m.Var,2)
C1=m.Array(m.Var,2)
cC1=m.Array(m.Var,2)
C2=m.Array(m.Var,2)
C2m=m.Array(m.Var,2)
cC2=m.Array(m.Var,2)
C3=m.Array(m.Var,2)
C4=m.Array(m.Var,2)
C5=m.Array(m.Var,2)
C6=m.Array(m.Var,2)
C7=m.Array(m.Var,2)
C8=m.Array(m.Var,2)
C9=m.Array(m.Var,2)
C10=m.Array(m.Var,2)
C11=m.Array(m.Var,2)
C12=m.Array(m.Var,2)
r1=m.Array(m.Var,2)
r2=m.Array(m.Var,2)
r3=m.Array(m.Var,2)
r4=m.Array(m.Var,2)
r5=m.Array(m.Var,2)
r6=m.Array(m.Var,2)
r7=m.Array(m.Var,2)
r8=m.Array(m.Var,2)
r9=m.Array(m.Var,2)
r10=m.Array(m.Var,2)
r11=m.Array(m.Var,2)
r12=m.Array(m.Var,2)
k5=m.Array(m.Var,2)
#Define GEKKO variables that determine if time point contains data to be used in regression
#index for objective (0=not measured,1=measured)
zm = m.Array(m.Param,2)
zm[0].value=z1
zm[1].value=z4
# fit to measurement
x=m.Array(m.Var,2)
x[0].value=x_data1[0]
x[1].value=x_data4[0]
#adjustable parameters
kdoc1 = 0 #m-1h-1
kdoc2 = m.FV(1,ub=10) #m-1h-1
S1 = 0
S2 = m.FV(1,ub=10)
#Define Variables
for i in range(2):
#State variables
# M (mol/L)
#C0 = m.Var(cl[i])
C0[i] = m.Var(value=cl[i],ub=1e-3)
C1[i] = m.Var(value=6.9e-6,ub=2e-3)
C2[i] = m.Var(value=x[i],ub=1e-3)
C2m[i] = m.Param(xm[i],lb=0)
meas = m.Param(zm[i])
m.Minimize(meas*((C2[i]-C2m[i]))**2)
C3[i] = m.Var(value=0)
C4[i] = m.Var(value=3.16e-8)
C5[i] = m.Var(value=3.83e-7)
C6[i] = m.Var(value=0)
C7[i]= m.Var(value=0)
C8[i] = m.Var(value=0)
C9[i] = m.Var(value=0)
C10[i] = m.Var(value=0)
C11[i] = m.Var(value=0)
C12[i] = m.Var(value=3.21e-3)
TOC[i] = m.Var(value=5.93)#mg-C/L
cC1[i] = m.Var(value=0)
cC2[i] = m.Var(value=0)
#temperature ºC
T = 20
#Rate constants
k1 = 2040000000*math.exp(-1887/(273.15+T))*3600
k2 = 138000000*math.exp(-8800/(273.15+T))*3600
k3 = 300000*math.exp(-2010/(273.15+T))*3600
k4 = 0.00000065*3600
k6 = 26700*3600
k7 = 167*3600
k8 = 27700*3600
k9 = 8300*3600
k10 = 0
k5[i] = 37800000000*math.exp(-2169/(273.15+T))*C4[i]\
+2.52e25*math.exp(-16860/(273.15+T))*C10[i]\
+0.87*math.exp(-503/(273.15+T))*C9[i]
r1[i] = k1 * C0[i] * C1[i]
r2[i] = k2 * C2[i]
r3[i] = k3 * C0[i] * C2[i]
r4[i] = k4 * C3[i]
r5[i] = k5[i] * C2[i] * C2[i]
r6[i] = k6 * C3[i] * C1[i]* C4[i]
r7[i] = k7 * C3[i] * C5[i]
r8[i] = k8 * C6[i] * C3[i]
r9[i] = k9 * C6[i] * C2[i]
r10[i] = k10 * C2[i] * C3[i]
r11[i] = kdoc1*S1*TOC[i]*C2[i]/12000
r12[i] = kdoc2*S2*TOC[i]*C0[i]/12000
t = m.Param(value=m.time)
m.Equation(C0[i].dt()== -r1[i] + r2[i] - r3[i] + r4[i] + r8[i] - r12[i])
m.Equation(C1[i].dt()== -r1[i] + r2[i] + r5[i] - r6[i] + r11[i])
m.Equation(C2[i].dt()== r1[i] - r2[i] - r3[i] + r4[i] - r5[i] + r6[i] - r9[i] - r10[i] - r11[i])
m.Equation(C3[i].dt()== r3[i] - r4[i] + r5[i] - r6[i] - r7[i] - r8[i] - r10[i])
m.Equation(C4[i].dt()== 0)
m.Equation(C5[i] == 1e-14/C4[i])
m.Equation(C6[i].dt()== r7[i] - r8[i] - r9[i])
m.Equation(C7[i] == (3.16e-8*C0[i])/C4[i])
m.Equation(C1[i] == (5.01e-10*C8[i])/C4[i])
m.Equation(C11[i] == (5.01e-11*C10[i])/C4[i])
m.Equation(C9[i] == (5.01e-7*C9[i])/C4[i])
m.Equation(C10[i] == C12[i] - 2*C11[i] - C5[i] + C4[i])
m.Equation(C12[i].dt()== 0)
m.Equation(TOC[i].dt()== 0)
m.Equation(cC1[i] == 17000*C1[i])
m.Equation(cC2[i] == 51500*C2[i])
m.Equation(S1+S2<1)
#Application options
m.options.SOLVER = 3 #IPOPT solver
m.options.IMODE = 5 #Dynamic Simultaneous - estimation = MHE
m.options.EV_TYPE = 2 #absolute error
m.options.NODES = 3 #collocation nodes (2,5)
if True:
#kdoc1.STATUS=1
kdoc2.STATUS=1
#S1.STATUS=1
S2.STATUS=1
#Solve
#m.open_folder() # open folder if remote=False to see infeasibilities.txt
m.solve(disp=True)
print('Final SSE Objective: ' + str(m.options.objfcnval))
print('Solution')
#print('kdoc1 = ' + str(kdoc1.value[0]))
print('kdoc2 = ' + str(kdoc2.value[0]))
#print('S1 = ' + str(S1.value[0]))
print('S2 = ' + str(S2.value[0]))
plt.figure(1,C2[0],'b.--',label='Predicted 1')
plt.plot(m.time,C2[1],'r.--',label='Predicted 4')
plt.plot(t_data1,x_data1,'bx',label='Meas 1')
plt.plot(t_data4,x_data4,'rx',label='Meas 4')
plt.legend(loc='best')
plt.subplot(2,cC2[0].value,'b',label ='C2_M1')
plt.plot(m.time,cC2[1].value,'r',label ='C2_M4')
plt.plot(t_data1,x_data1mgl,x_data4mgl,label='Meas 4')
#plt.plot(m.time,C0[i].value,C1[i].value,C3[i].value,C4[i].value,C5[i].value,C6[i].value,C7[i].value,label ='C7')
plt.legend()
plt.subplot(4,C8[i].value,C9[i].value,C10[i].value,C11[i].value,C12[i].value,label ='C12')
plt.legend()
plt.xlabel('time (h)')
plt.show()
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
