下面是一个ODE参数的最小二乘拟合代码。使用了Python的“最小化”和“最小二乘”函数。已经尝试了不同的方法和ODE解算器/步骤(scipy ODE/odeint)。这是一个在MATLAB中很容易解决的问题,但是Python一直返回最初的猜测。我希望您能找到一个编码错误,否则我会对Python优化函数感到失望。Obj表示目标(残差平方和),ode函数(一阶)表示未知参数的方程组。数据集已附加。在
import numpy as np
from scipy.integrate import ode
from scipy.optimize import least_squares
from scipy.optimize import minimize
from scipy.optimize import SR1
import matplotlib.pyplot as plt
import math
Minput=np.loadtxt('C:\\Users\\Ladan\\Documents\\Python Scripts\\Python\\moisturesmoothopt.txt')
Minput=Minput.flatten()
time=np.linspace(0,1800,901)
A=np.zeros(3)
XC,RC,alpha=A
#bnds=([0,0,0],[Minput[0],math.inf,math.inf])
bnds=((0,Minput[0]),(0,math.inf),(0,math.inf))
def firstorder(X,time,A):
if X>=XC:
dX=-RC
if X<XC:
dX=-RC*(X/XC)**alpha
return dX
def obj(A):
X0=Minput[0]
# Xpred=odeint(firstorder,X0,time,args=(A,))
Xpred=ode(firstorder).set_integrator('vode', method='bdf',
order=15).set_initial_value(Minput[0],0).set_f_params(A)
#Xpred=ode(firstorder).set_integrator('lsoda').set_initial_value(Minput[0],0).set_f_params(A)
EPR=Xpred
EPR2=EPR.y.flatten()
ERRone=np.sum(np.power((EPR2-Minput),2))
ERR=ERRone/((901-3)) # residual sum of squares deivided by dof
return ERR
XC=1
RC=0.005
alpha=1.5
A0=[XC,RC,alpha]
Parameters=minimize(obj,A0,method='SLSQP',bounds=bnds,options={'ftol':1e-10,
'maxiter': 1000})
print('parameters',Parameters)
Minput数组的数据在线共享:
虽然我在scipy中使用了较新的ode解算器,但我倾向于使用好的ol'}。虽然取得了一些进展,但我确实收到了一些关于无效值的警告。你需要进一步调查,但这会让你走上正轨
odeint
函数,这个函数稍微老一点,但仍然非常可靠,在许多情况下,由于我不完全理解的原因,它的性能更好。总之,我很大程度上重组了您的分析代码,以同时使用scipy.optimize.least_squares
和{打印输出:
^{pr2}$相关问题 更多 >
编程相关推荐