在python中使用GEKKO估计电池参数时出错

如何解决在python中使用GEKKO估计电池参数时出错

我正在尝试将模型的模拟结果与实验室的实验结果相匹配。我从这里给出的解决方案中获取了这个想法： enter link description here

应用 GEKKO 方法前的结果：

我希望模拟曲线与实验曲线完全匹配。

m = GEKKO(remote=False)
Vocv = Data.loc[:,'Vocv'].tolist()
Tt=    Data.loc[:,'Tt'].tolist() 
It=    Data.loc[:,'It'].tolist()
Vmeas= Data.loc[:,'Experimental_Voltage(v)']/6000.tolist()
Vocv = m.Param(Vocv); Tt = m.Param(Tt); It = m.Param(It)
##m.time = Tt; time = m.Var(0); m.Equation(time.dt()==1)
    


 R0 = m.FV(lb= 2.448e-07,ub=100); R0.STATUS=1
R1 = m.FV(lb= 3e-07,ub=100);  R1.STATUS=1
R2 = m.FV(lb=3e-07,ub=100); R2.STATUS=1
C1 = m.FV(lb=0.02,ub=1000); C1.STATUS=1
C2 = m.FV(lb=0.02,ub=1000); C2.STATUS=1
ym = m.Param(Vmeas)
yp = m.Var(Vmeas); m.Equation(yp==Vocv+(R0*It) \
                    +(R1*It)*(1-m.exp(-1/(R1*C1))*Tt) \
                    +(R2*It)*(1-m.exp(-1/(R2*C2))*Tt))
m.Minimize((yp-ym)**2)
m.options.IMODE = 2

m.solve(disp=False)

import matplotlib.pyplot as plt
Ex_Time= Data.loc[:,'Experimental_Time(s)']
plt.plot(Ex_Time,Experimental_Voltage/6000)
plt.plot(Tt,yp)
plt.legend([r'$Simulated_Data$',r'$Experimental_Data$'])
plt.ylabel('Voltage')
plt.xlabel('Time')
plt.savefig('sysid.png')
plt.show()

print('R0: ' + str(R0.value[0]))
print('R1: ' + str(R1.value[0]))
print('R2: ' + str(R2.value[0]))
print('C1: ' + str(C1.value[0]))
print('C2: ' + str(C2.value[0]))

输出：

 R0: 2.448e-07
R1: 3e-07
R2: 3e-07
C1: 0.02
C2: 0.02

编辑： 正如答案中所建议的那样，我已经进行了更改，但得到的结果如下：

它仍然不匹配/重叠实验数据。

预期结果：

解决方法

尝试为 R1 和 R2 设置一个小的正下限以避免等式中被零除：

R1 = m.FV(lb=1e-5,ub=100);  R1.STATUS=1
R2 = m.FV(lb=1e-5,ub=100);  R2.STATUS=1

您提供的示例缺少数据，因此我包含了一些示例数据。这个例子现在收敛到一个解决方案。

R0: 0.058215197521
R1: 1.0089645613e-05
R2: 1.0089645613e-05
C1: 42.633186357
C2: 42.633186357

from gekko import GEKKO
m = GEKKO(remote=False)
Vocv = [4.1,4.2,4.05,4.15,4.2]
Tt=    [0,1,2,3,4]
It=    [1.2,1.3,1.4,1.1,1.0]
Vmeas= [4.15,4.25,4.25]
Vocv = m.Param(Vocv); Tt = m.Param(Tt); It = m.Param(It)

R0 = m.FV(lb=0,ub=100);     R0.STATUS=1
R1 = m.FV(lb=1e-5,ub=100);  R2.STATUS=1
C1 = m.FV(lb=0,ub=1000);    C1.STATUS=1
C2 = m.FV(lb=0,ub=1000);    C2.STATUS=1

ym = m.Param(Vmeas)
yp = m.Var(Vmeas); m.Equation(yp==Vocv+R0*It \
                        +R1*It*m.exp(-(Tt/R1)*C1) \
                        +R2*It*m.exp(-(Tt/R2)*C2))
m.Minimize((yp-ym)**2)
m.options.IMODE = 2

m.solve(disp=False)

print('R0: ' + str(R0.value[0]))
print('R1: ' + str(R1.value[0]))
print('R2: ' + str(R2.value[0]))
print('C1: ' + str(C1.value[0]))
print('C2: ' + str(C2.value[0]))

import matplotlib.pyplot as plt
plt.plot(Tt,yp,'b-o')
plt.plot(Tt,ym,'rx')
plt.ylabel('Voltage')
plt.xlabel('Time')
plt.savefig('sysid.png')
plt.show()

切换到 IMODE=2 以缩短求解时间，因为模型中没有微分方程。