2

I am trying to match the simulated result from my model with the experimental result from lab. I have taken the idea from the solution given here: enter link description here

Result before applying GEKKO method:enter image description here

I want the Simulated Curve to be completely match with the Experimental Curve.

Code:

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]))

Output:

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

EDIT: As suggested in the answer I have made changes however, the result I am getting is following: enter image description here

Its still not matching/overlapping the experimental data.

Expected Result:

enter image description here

N_T
  • 33
  • 1
  • 5

1 Answers1

0

Try setting a small positive lower bound for R1 and R2 to avoid divide-by-zero in the equation:

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

The example you provided is missing the data so I included some sample data. This example now converges to a solution.

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

Voltage solution

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()

Switched to IMODE=2 to improve the solution time because there are no differential equations in the model.

John Hedengren
  • 12,068
  • 1
  • 21
  • 25
  • 1
    I have edited my question and taken different initial vlaues of R0,R1,R2,C1 and C2 however, instead of reducing the error it is taking the same values and not optimizing. – N_T Jul 06 '21 at 15:00
  • Does the solver report a successful solution? You may need to add parameter bounds or else give a better initial condition to find the global (not local) optimum. – John Hedengren Jul 09 '21 at 11:46
  • 1
    I think the solver is giving the solution however, it is only taking the 'lb' value of parameter bounds and not really minimizing the error. I have given the best initial values but still not good fit. – N_T Jul 15 '21 at 13:21
  • Try a lower value of the parameter lb, even letting it go negative if needed. – John Hedengren Jul 16 '21 at 17:21