英文:
Non-Linear data fitting for kinetic data
问题
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
# 步骤 1：读取实验数据
data1 = pd.read_excel("本地文件地址")
# 提取偏压（A、B、C）和反应速率的数据
pp = data1.iloc[:49, 2:5].values
TT = data1.iloc[:49, 16].values
rr = data1.iloc[:49, 17].values
F0 = data1.iloc[:49, 5:10].values
W = data1.iloc[:49, 15].values
FF = data1.iloc[:49, 10:15].values
data = np.hstack((pp, F0, FF, W.reshape(-1, 1), TT.reshape(-1, 1)))
#------------------------------------------------------------------------
class Kinetic:
"""Langmuir-Hinshelwood Two-Site Kinetic Model Solver"""
R = 8.314 # J/mol.K
Ea = 45000 # J/mol
theta_g = [0.9, 0.1] # 1
w = np.linspace(0, 1, 100)
@staticmethod
def k3(flow, data, K1, K2, k4, K5, k0):
"""Arrhenius estimation for k3"""
pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = data
return k0 * np.exp(-Kinetic.Ea / (Kinetic.R * T))
@staticmethod
def s(flow, data, K1, K2, k4, K5, k0):
"""Partial kinetic rate"""
pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = data
# 计算反应器每个位置的气体压力
pA_t = flow[0] / np.sum(flow)
pB_t = flow[1] / np.sum(flow)
pC_t = flow[2] / np.sum(flow)
k3 = Kinetic.k3(flow, data, K1, K2, k4, K5, k0)
return (K1 * K2 * k3 / k4) * pA_t * pB_t
@staticmethod
def equation(x, y, C, s):
return 1/(1 + C + np.sqrt(s*y/x)) - x
@staticmethod
def system(theta, flow, data, K1, K2, k4, K5, k0):
"""等温覆盖系统"""
pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = data
t0, t1 = theta
s = Kinetic.s(flow, data, K1, K2, k4, K5, k0)
pA_t = flow[0] / np.sum(flow)
pB_t = flow[1] / np.sum(flow)
pC_t = flow[2] / np.sum(flow)
C0 = K2 * pB_t
C1 = K1 * pA_t + pC_t / K5
return np.array([
Kinetic.equation(t0, t1, C0, s),
Kinetic.equation(t1, t0, C1, s),
])
@staticmethod
def _theta(flow, data, K1, K2, k4, K5, k0):
"""单等温覆盖"""
solution = fsolve(
Kinetic.system,
Kinetic.theta_g,
args=(flow, data, K1, K2, k4, K5, k0),
full_output=False
)
if not(all((solution >= 0) & (solution <= 1.0))):
raise ValueError("Theta 不受限于域: %s" % solution)
return solution
@staticmethod
def theta(flow, data, k1, k2, k4, k5, k0):
"""等温覆盖"""
return np.apply_along_axis(Kinetic._theta, 0, flow, data, k1, k2, k4, k5, k0)
@staticmethod
def r1(flow, w, data, K1, K2, k4, K5, k0):
"""全局动力学速率"""
pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = data
pA_t = flow[0] / np.sum(flow)
pB_t = flow[1] / np.sum(flow)
pC_t = flow[2] / np.sum(flow)
k3 = Kinetic.k3(flow, data, K1, K2, k4, K5, k0)
t0, t1 = Kinetic.theta(flow, data, K1, K2, k4, K5, k0)
v=np.array([-1,-1, 1, 0, 0])
return ((K1 * K2 * k3) * pA_t * pB_t * t0 * t1) * v
@staticmethod
def dFdm(data, K1, K2, k4, K5, k0):
"""全局动力学速率"""
residual =[]
sq=[]
for i in data:
pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = i
F0 = np.array([F0_a, F0_b, F0_c, F0_ar, F0_h])
FF = np.array([FF_a, FF_b, FF_c, FF_ar, FF_h]) # 最终实验流量
flow = odeint(Kinetic.r1, F0, Kinetic.w, args = (i, K1, K2, k4, K5, k0))
flow_pred = flow[-1,2] # 预测的最终模型流量
res = (flow_pred - FF_c)*1000
residual.append(res)
return residual
bounds=((0, 0, 0, 0, 0), (np.inf,
<details>
<summary>英文:</summary>
![Langmuir Hinshelwood Model Image][1]
[Data](https://docs.google.com/spreadsheets/d/1-M1DpuI4XNfSJktVG_ibzy-RzSAeUSIj/edit?usp=sharing&ouid=116419220787995800483&rtpof=true&sd=true)
*Update: I modified the code after the suggestions of @jlandercy. Further, the code now solves a differential equation shown below to optimize the flow of my product gas. Sample data is updated.
[![PFR differential equation optimized][2]][2]
I am trying to fit a two-site Langmuir Hinshelwood model which I derived, (basically, rate as a function of pressures of reactants and products) to my experimental data. The model converges, but changing the initial parameters can make a stark difference. The problem is those parameters matter for my analysis. Since I am comparing different models, the model vs experiment % fit depends on these parameters (K1, K2 ...).
I wanted to know if there is any other method to check that a particular set of parameters are better than the others (obviously a particular combination will give a slightly better-fit percentage, but that, in my opinion, may not the deciding criteria). The MSE landscape calculation, as suggested by @jlandercy, has been useful in understanding the convexity relationship between different parameters, however, I could not extract more information from them to move ahead and fix my K parameters.
P.S.: The other models are way simpler, in terms of mathematics, than this model. So, if I this is well understood, I should be able to compare the models with decent parameter accuracy.
Thanks in advance for the help!
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
Step 1: Read experimental data
data1 = pd.read_excel("Local address of file")
Extract the partial pressures (A, B, C) and the reaction rate
pp = data1.iloc[:49, 2:5].values
TT = data1.iloc[:49, 16].values
rr = data1.iloc[:49, 17].values
F0 = data1.iloc[:49, 5:10].values
W = data1.iloc[:49, 15].values
FF = data1.iloc[:49, 10:15].values
data = np.hstack((pp, F0, FF, W.reshape(-1,1), TT.reshape(-1,1)))
#-----------------------------------------------------------------------------
class Kinetic:
"""Langmuir-Hinshelwood Two-Site Kinetic Model Solver"""
R = 8.314 # J/mol.K
Ea = 45000 # J/mol
theta_g = [0.9, 0.1] # 1
w = np.linspace(0,1,100)
@staticmethod
def k3(flow, data, K1, K2, k4, K5, k0):
"""Arrhenius estimation for k3"""
pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = data
return k0 * np.exp(-Kinetic.Ea / (Kinetic.R * T))
@staticmethod
def s(flow, data, K1, K2, k4, K5, k0):
"""Partial kinetic rate"""
pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = data
#Calculating gas pressures at every position of PFR reactor
pA_t = flow[0] / np.sum(flow)
pB_t = flow[1] / np.sum(flow)
pC_t = flow[2] / np.sum(flow)
k3 = Kinetic.k3(flow, data, K1, K2, k4, K5, k0)
return (K1 * K2 * k3 / k4) * pA_t * pB_t
@staticmethod
def equation(x, y, C, s):
return 1/(1 + C + np.sqrt(s*y/x)) - x
@staticmethod
def system(theta, flow, data, K1, K2, k4, K5, k0):
"""Isothermal coverages system"""
pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = data
t0, t1 = theta
s = Kinetic.s(flow, data, K1, K2, k4, K5, k0)
pA_t = flow[0] / np.sum(flow)
pB_t = flow[1] / np.sum(flow)
pC_t = flow[2] / np.sum(flow)
C0 = K2 * pB_t
C1 = K1 * pA_t + pC_t / K5
return np.array([
Kinetic.equation(t0, t1, C0, s),
Kinetic.equation(t1, t0, C1, s),
])
@staticmethod
def _theta(flow, data, K1, K2, k4, K5, k0):
"""Single Isothermal coverages"""
solution = fsolve(
Kinetic.system,
Kinetic.theta_g,
args=(flow, data, K1, K2, k4, K5, k0),
full_output=False
)
if not(all((solution >= 0) & (solution <= 1.0))):
raise ValueError("Theta are not constrained to domain: %s" % solution)
return solution
@staticmethod
def theta(flow, data, k1, k2, k4, k5, k0):
"""Isothermal coverages"""
return np.apply_along_axis(Kinetic._theta, 0, flow, data, k1, k2, k4, k5, k0)
@staticmethod
def r1(flow, w, data, K1, K2, k4, K5, k0):
"""Global kinetic rate"""
pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = data
pA_t = flow[0] / np.sum(flow)
pB_t = flow[1] / np.sum(flow)
pC_t = flow[2] / np.sum(flow)
#print(pA)
k3 = Kinetic.k3(flow, data, K1, K2, k4, K5, k0)
t0, t1 = Kinetic.theta(flow, data, K1, K2, k4, K5, k0)
v=np.array([-1,-1, 1, 0, 0])
#v = v[:, np.newaxis] # Convert v to a 2D array with shape (5, 1)
return ((K1 * K2 * k3) * pA_t * pB_t * t0 * t1) * v
@staticmethod
def dFdm(data, K1, K2, k4, K5, k0):
"""Global kinetic rate"""
residual =[]
sq=[]
for i in data:
pA, pB, pC, F0_a, F0_b, F0_c, F0_ar, F0_h, FF_a, FF_b, FF_c, FF_ar, FF_h, mg, T = i
F0 = np.array([F0_a, F0_b, F0_c, F0_ar, F0_h])
FF = np.array([FF_a, FF_b, FF_c, FF_ar, FF_h]) # Final Experimental flows
flow = odeint(Kinetic.r1, F0, Kinetic.w, args = (i, K1, K2, k4, K5, k0))
flow_pred = flow[-1,2] # Predicted final model flows
res = (flow_pred - FF_c)*1000
residual.append(res)
return residual
bounds=((0, 0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf, np.inf))
initial_guesses = [8, 8, 50, 16, 6]
result = least_squares(lambda k: Kinetic.dFdm(data, *k), x0=initial_guesses, bounds=bounds, ftol=1e-14)
K1_fit, K2_fit, k4_fit, K5_fit, k0_fit = result.x
print(K1_fit, K2_fit, k4_fit, K5_fit, k0_fit)
[1]: https://i.stack.imgur.com/kULtU.jpg
[2]: https://i.stack.imgur.com/KypHf.jpg
</details>
# 答案1
**得分**: 2
以下是翻译好的部分：
### TL;DR
无法通过您的11个数据点确定它们是否适合所提出的模型。模型收敛，但参数错误较大，对初始猜测非常敏感。因此，在这一点上不要信任模型的输出。
但是，可以显示，至少使用50个跨越变量空间的数据点，模型会以合理的精度收敛到预期的参数，并且在数据上具有合理的误差项，过程是稳定的。
这意味着您要么需要以特定的精度收集更多数据点，要么您的动力学可能不合理地遵循所提出的模型。
### Procedure
让我们将全局任务分解为更小的任务。要解决动力学问题，我们需要解决以下问题：
- 具有域约束的非线性方程组求解（部分覆盖度位于0和1之间）：
- 确保解的存在性和唯一性（参见算盘）
- 对初始参数进行合理的猜测（收敛和解的选择）
- 具有域约束的非线性最小二乘（NLLS）求解（化学常数是正定的）：
- 确保均方误差（MSE）地形合理，以便达到全局最小值（在变量空间上具有足够的数据点和可接受的误差）
- 使过程对初始参数的变化具有鲁棒性
- 参数归一化和解耦（如果需要）
在这两个任务中，我们需要确保收敛到全局最小值或所需的解，而且需要保持数值稳定性。
`scipy`包是解决此类优化问题的好工具，可使用以下函数：
- [`optimize.fsolve`][1] 用于非线性方程组；
- [`optimize.curve_fit`][2] 用于NLLS拟合。
#### Solver
下面是一个详细处理此问题的类示例：
```python
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
class Kinetic:
"""Langmuir-Hinshelwood Two-Site Kinetic Model Solver"""
R = 8.314 # J/mol.K
Ea = 37000 # J/mol
T = 473.15 # K
theta_ = [0.9, 0.9] # 1
@staticmethod
def k3(p, k1, k2, k4, k5, k0):
"""Arrhenius estimation for k3"""
return k0*np.exp(-Kinetic.Ea/(Kinetic.R*Kinetic.T))
# 其他方法的翻译省略
# ...
请注意，此类使用TRF求解器，而不是默认的LM求解器，以确保强制执行参数约束。对于您的数据集，使用LM会导致某些初始猜测产生负参数。
Experimental dataset
如果我们使用提议的初始猜测检查模型与您的实验数据集：
import pandas as pd
df = pd.read_excel("data.xlsx")
p = df.filter(regex="p").values
r = df["rate"].values
Kinetic.solve(p, r, 5, 1e5, 50, 10, 100)
我们得到以下调整结果：
(array([1.39075854e-01, 2.15668993e+03, 1.16513471e+00, 2.69648513e-01,
2.64057386e+00]),
array([[ 1.48622028e+03, -6.20830186e+06, 2.86587881e+05,
1.03973854e+02, -2.89149484e+04],
[-6.20830186e+06, 1.73740772e+12, -3.36688691e+09,
9.40415116e+06, 1.17499454e+08],
[ 2.86587881e+05, -3.36688691e+09, 5.84376597e+07,
5.50069835e+03, -5.57235249e+06],
[ 1.03973854e+02, 9.40415116e+06, 5.50069835e+03,
8.40397226e+01, -2.04455154e+03],
[-2.89149484e+04, 1.17499454e+08, -5.57235249e+06,
-2.04455154e+03, 5.62563731e+05]]),
{'nfev': 10,
'fvec': array([ 1.63900711e-06, 1.66564630e-06, 1.55523211e-06, 2.83708013e-06,
2.40778598e-06, 2.85927736e-06, 3.37890387e-08, -1.19231145e-06,
-4.03211541e-07, -2.37714019e-06, -3.98891395e-06])},
'`gtol` termination condition is satisfied.',
1)
调整收敛，但绝对不令人满意：误差比参数大数个数量级。
此外，更改初始猜测会返回完全不同的参数，这意味着MSE地形没有正确塑造，无法确保稳定收敛。
这两种症状通常都是缺乏数据和/或在回归变量上具有太高误差的证据。但事实上，也可能是错误的模型选择和错误的实施，因此需要分析过程的敏感性。
Synthetic dataset
让我们构建一个跨越变量空
英文:
TL;DR
It's not possible to tell with your 11 points if they fits properly the proposed model. The model converges but with high errors on parameters and high sensitivity w.r.t. initial guess. So at this point don't trust the output of the model.
But it is possible to show that with at least 50 points spanning the variable space the model converges towards expected parameters with decent precision and the process is stable with reasonable error term on data.
It means that you either need to collect more points with ad hoc precision or your kinetic may not reasonably follow the proposed model.
Procedure
Let's break down the global task into smaller tasks. To solve the kinetic we need to address:
- Non linear system of equations solving with domain constraints (partial coverages lie between 0 and 1):
- Ensure solution existence and uniqueness (see abacus)
- Educated guess for initial parameters (convergence and solution selection)
- Non linear least squares (NLLS) solving with domain constraints (chemical constants are positive definite):
- Ensure reasonable MSE landscape to make global minimum reachable (enough points with bearable errors over variable space)
- Make procedure robust against initial parameters variation
- Parameters normalization and decoupling (if necessary)
In both task we need to ensure convergence towards global minimum or desired solution in addition of numerical stability.
The scipy package is the good companion to solve this class of optimization problems using:
optimize.fsolvefor non linear equations system;optimize.curve_fitfon NLLS adjustment.
Solver
Below a class handling this problem in details:
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
class Kinetic:
"""Langmuir-Hinshelwood Two-Site Kinetic Model Solver"""
R = 8.314 # J/mol.K
Ea = 37000 # J/mol
T = 473.15 # K
theta_ = [0.9, 0.9] # 1
@staticmethod
def k3(p, k1, k2, k4, k5, k0):
"""Arrhenius estimation for k3"""
return k0*np.exp(-Kinetic.Ea/(Kinetic.R*Kinetic.T))
@staticmethod
def s(p, k1, k2, k4, k5, k0):
"""Partial kinetic rate"""
pA, pB, pC = p
k3 = Kinetic.k3(p, k1, k2, k4, k5, k0)
return (k1*k2*k3/k4)*pA*pB
@staticmethod
def equation(x, y, C, s):
return 1/(1 + C + np.sqrt(s*y/x)) - x
@staticmethod
def isopleth(x, y, s):
return (x + np.sqrt(s*x*y) - 1)/x
@staticmethod
def system(theta, p, k1, k2, k4, k5, k0):
"""Isothermal coverages system"""
pA, pB, pC = p
t0, t1 = theta
s = Kinetic.s(p, k1, k2, k4, k5, k0)
C0 = k2*pB
C1 = k1*pA + pC/k5
return np.array([
Kinetic.equation(t0, t1, C0, s),
Kinetic.equation(t1, t0, C1, s),
])
@staticmethod
def _theta(p, k1, k2, k4, k5, k0):
"""Single Isothermal coverages"""
solution = optimize.fsolve(
Kinetic.system,
Kinetic.theta_,
args=(p, k1, k2, k4, k5, k0),
full_output=False
)
if not(all((solution >= 0.) & (solution <= 1.0))):
raise ValueError("Theta are not constrained to domain: %s" % solution)
return solution
@staticmethod
def theta(p, k1, k2, k4, k5, k0):
"""Isothermal coverages"""
return np.apply_along_axis(Kinetic._theta, 0, p, k1, k2, k4, k5, k0)
@staticmethod
def r(p, k1, k2, k4, k5, k0):
"""Global kinetic rate"""
pA, pB, pC = p
k3 = Kinetic.k3(p, k1, k2, k4, k5, k0)
t0, t1 = Kinetic.theta(p, k1, k2, k4, k5, k0)
return (k1*k2*k3)*pA*pB*t0*t1
@staticmethod
def solve(p, r, *k):
"""Global kinetic constants adjustment"""
return optimize.curve_fit(
Kinetic.r, p.T, r, p0=k,
bounds=((0, 0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf, np.inf)),
method="trf",
gtol=1e-10,
full_output=True
)
@staticmethod
def abaque(s, C=None):
lin = np.linspace(0.0001, 1.1, 200)
X, Y = np.meshgrid(lin, lin)
if C is None:
C = np.arange(0, 26, 1)
Cx = Kinetic.isopleth(X, Y, s)
Cy = Kinetic.isopleth(Y, X, s)
fig, axe = plt.subplots()
xlabels = axe.contour(X, Y, Cx, C, cmap="jet")
axe.clabel(xlabels, xlabels.levels, inline=True, fontsize=7)
ylabels = axe.contour(X, Y, Cy, C, cmap="jet")
axe.clabel(ylabels, ylabels.levels, inline=True, fontsize=7)
axe.set_title("Isothemal Coverages System\nConstant Isopleths ($s={}$)".format(s))
axe.set_xlabel(r"Partial Coverage, $\theta_0$ [-]")
axe.set_ylabel(r"Partial Coverage, $\theta_1$ [-]")
axe.set_aspect("equal")
axe.set_xlim([0, 1.1])
axe.set_ylim([0, 1.1])
axe.grid()
return axe
Notice it uses the TRF solver instead of the default LM solver to ensure parameters constraints are enforced. With your dataset, using LM return negative parameters for some initial guesses.
Experimental dataset
If we check the model against your experimental dataset with proposed initial guess:
import pandas as pd
df = pd.read_excel("data.xlsx")
p = df.filter(regex="p").values
r = df["rate"].values
Kinetic.solve(p, r, 5, 1e5, 50, 10, 100)
We get the following adjustment:
(array([1.39075854e-01, 2.15668993e+03, 1.16513471e+00, 2.69648513e-01,
2.64057386e+00]),
array([[ 1.48622028e+03, -6.20830186e+06, 2.86587881e+05,
1.03973854e+02, -2.89149484e+04],
[-6.20830186e+06, 1.73740772e+12, -3.36688691e+09,
9.40415116e+06, 1.17499454e+08],
[ 2.86587881e+05, -3.36688691e+09, 5.84376597e+07,
5.50069835e+03, -5.57235249e+06],
[ 1.03973854e+02, 9.40415116e+06, 5.50069835e+03,
8.40397226e+01, -2.04455154e+03],
[-2.89149484e+04, 1.17499454e+08, -5.57235249e+06,
-2.04455154e+03, 5.62563731e+05]]),
{'nfev': 10,
'fvec': array([ 1.63900711e-06, 1.66564630e-06, 1.55523211e-06, 2.83708013e-06,
2.40778598e-06, 2.85927736e-06, 3.37890387e-08, -1.19231145e-06,
-4.03211541e-07, -2.37714019e-06, -3.98891395e-06])},
'`gtol` termination condition is satisfied.',
1)
The adjustment converges but it is definitely not satisfying: errors are several decade higher than parameters.
In addition, changing the initial guess returns drastically different parameters meaning the MSE landscape is not properly shaped to ensure stable convergence.
Both symptoms are generally an evidence against the lack of data and/or too high errors on regressed variable. But indeed it can be also bad model selection and bad implementation, hence the need to analyze procedure sensitivity.
Synthetic dataset
Let's build a synthetic dataset spanning the variable space with normal noise on it:
plin = np.linspace(0.1, 0.3, 5)
PA, PB, PC = np.meshgrid(plin, plin, plin)
P = np.vstack([
PA.flatten(), PB.flatten(), PC.flatten()
]).T
R = Kinetic.r(P.T, 5, 1e5, 50, 10, 100)
R += np.random.randn(R.shape[0])*1e-7
Following calls decently converge towards the expected parameters with fair errors:
Kinetic.solve(P, R, 1, 1, 1, 1, 1)
Kinetic.solve(P, R, 1e2, 1e2, 1e2, 1e2, 1e2)
Kinetic.solve(P, R, 1e4, 1e4, 1e4, 1e4, 1e4)
(array([4.99738701e+00, 1.02579705e+05, 4.74129338e+01, 9.98744998e+00,
1.00075578e+02]),
array([[ 8.22367638e-06, 8.36365990e-02, 8.07407481e-03,
-5.91768686e-07, -2.40408987e-04],
[ 8.36365990e-02, 8.97492687e+07, 8.22981373e+01,
-3.59165834e-04, -7.40000336e+00],
[ 8.07407481e-03, 8.22981373e+01, 7.94766011e+00,
-1.77033582e-04, -2.36480247e-01],
[-5.91768686e-07, -3.59165834e-04, -1.77033582e-04,
4.05361075e-05, 5.41574698e-06],
[-2.40408987e-04, -7.40000336e+00, -2.36480247e-01,
5.41574698e-06, 7.03833667e-03]]),
{'nfev': 175,
'fvec': array([ 1.63383888e-07, -4.98390311e-08, -3.59565336e-08, 6.62435527e-08,
-7.48972275e-08, -9.51285315e-08, 8.58923659e-10, -9.53560561e-08,
1.92097574e-07, 1.72359976e-08, 1.06487677e-07, -8.48179470e-08,
...
7.36713698e-08, -4.52283602e-08, 1.21064513e-07, -1.40410586e-08,
-8.25508331e-08])},
'`gtol` termination condition is satisfied.',
1)
Parameters are close to expected values and errors on parameters are low, leading at least 3 significant digits on each parameters.
MSE Landscape
When investigating such kind of problem it is important to sketch the MSE Landscape to check the convergence towards global optimum.
Because your problem has a high dimensionality, it requires to project it a lot. I leave the toolbox to do it here:
partial = functools.partial(Kinetic.r, P.T)
vectorized = np.vectorize(partial, otypes=[np.ndarray])
def error(y):
def wrapped(yhat):
return np.log10(np.sum((y-yhat)**2)/y.size)
return wrapped
N = 100
k1, k2, k4, k5, k0 = 5, 1e5, 50, 10, 100
k1l = np.logspace(-2, 2, N, base=10)
k2l = np.logspace(3, 7, N, base=10)
k4l = np.logspace(0, 4, N, base=10)
k5l = np.logspace(0, 4, N, base=10)
k0l = np.logspace(0, 4, N, base=10)
K1, K5 = np.meshgrid(k1l, k5l)
Rk = vectorized(K1, k2, k4, K5, k0)
mse = np.vectorize(error(R))
MSEk = mse(Rk)
And I'll take snapshot for some projections. For instance we see the non linear coupling between k4 and k0:
And between k1 and k0:
Other projections are more regular:
By inspecting the MSE Landscape, most of the projections are sufficiently smooth, steep and convex to ensure stable convergence once we have enough points with acceptable errors.
But the coupling with parameters k0 and k4 introduce some extra complexity as it shapes canyons in the MSE Landscape that may slower or disturb the convergence. Anyway it seems an educated guess can prevent the algorithm to pass by those canyons.
Conclusions
Solving this kinetic requires two critical steps:
- Stable system solving with constraint (calculus + educated guess)
- Stable NLLS fit solving with constraint (TRF)
We can be confident, if your experiment follow the proposed model and experimental dataset contains enough data points with acceptable errors we should be able to regress kinetic constants from it using the above procedure.
Miscellaneous
Details of the procedure are shown in the note below:
System Abacus
To solve the non linear system, we can investigate abacus to visually assess the existence of the solution within the domain:
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论