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

Question

我正在尝试将我的模型的模拟结果与实验室的实验结果相匹配。我从这里给出的解决方案中得到了这个想法： 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

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

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

预期结果：

Answer 1

尝试为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.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);  R1.STATUS=1
R2 = 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 以缩短求解时间，因为模型中没有微分方程。