How to increase predicted data when estimating parameters of set of complex differential equations by fitting experimental data?

user_14544264 0
Viewed 66

I am trying to estimate the parameters of my set of ODEs in my program, always minimizing the error between my experimental data and predicted data.

The problem is, I can obtain a good prediction and a very good fit, but I can only estimate the same number of points as my experimental data, which returns a very strange output.

Can you please provide me more information in how can I obtain a more accurate set of predicted points?

Code can be found below.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
xm = np.array([0,1,2,3,4,5])
ym = np.array([2.0,1.5,np.nan,2.2,3.0,5.0])
m = GEKKO(remote=False)
m.time = xm
a = m.FV(lb=0.1,ub=2.0)
a.STATUS=1
y = m.CV(value=ym,name='y',fixed_initial=False)
y.FSTATUS=1
m.Equation(y.dt()==a*y)
m.options.IMODE = 5
m.options.SOLVER = 1
m.solve(disp=True)
print('Optimized, a = ' + str(a.value[0]))
plt.figure(figsize=(6,2))
plt.plot(xm,ym,'bo',label='Meas')
plt.plot(xm,y.value,'r-',label='Pred')
plt.ylabel('y')
plt.ylim([0,6])
plt.legend()
plt.show()

If I replace variable m.time to obtain more data predicted with :
m.time = np.linspace(0,5,30) I get the error: raise Exception('Data arrays must have the same length, and match time discretization in dynamic problems') Exception: Data arrays must have the same length, and match time discretization in dynamic problems

gekko python
1 Answers

There are two options (Methods 1 and 2) that I've shown below. You can either plot the interpolating nodes to give you more resolution or create a new model for simulation.

Prediction Options with Gekko

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
xm = np.array([0,1,2,3,4,5])
ym = np.array([2.0,1.5,np.nan,2.2,3.0,5.0])
m = GEKKO(remote=False)
m.time = xm
a = m.FV(lb=0.1,ub=2.0)
a.STATUS=1
y = m.CV(value=ym,name='y',fixed_initial=False)
y.FSTATUS=1
m.Equation(y.dt()==a*y)
m.options.IMODE = 5
m.options.SOLVER = 1
m.options.CSV_WRITE = 2  # For Method 1
m.options.NODES = 3      # For Method 1 (options 3-6)
m.solve(disp=True)
print('Optimized, a = ' + str(a.value[0]))

# Method 1: Plot interpolating nodes
import json
with open(m.path+'//results_all.json') as f:
    results = json.load(f)

# Method 2: Re-simulate with more points
sim = GEKKO(remote=False)
ap = a.value[0]
xp = np.linspace(0,7); sim.time=xp
yp = sim.Var(y.value[0])
sim.Equation(yp.dt()==ap*yp)
sim.options.NODES = 3
sim.options.IMODE=4; sim.solve()

plt.figure(figsize=(6,2))
plt.plot(xm,ym,'bo',label='Meas')
plt.plot(xm,y.value,'gs-.',label='Pred Original')
plt.plot(results['time'],results['y'],'kx-',\
         MarkerSize=10,label='Pred Method 1')
plt.plot(xp,yp,'r.--',label='Pred Method 2')
plt.ylabel('y')
plt.ylim([0,10])
plt.legend()
plt.show()

A third option is to reset the .value of the original model but that can be tedious. Instead, you can also create both the estimation and simulation models in a loop as is done with an example of Moving Horizon Estimation and Model Predictive Control that use the same model but transfer parameters between them:

MHE and MPC

# use remote=True for MacOS
mhe = GEKKO(name='tclab-mhe',remote=False)
mpc = GEKKO(name='tclab-mpc',remote=False)

# create 2 models (MHE and MPC) in one loop
for m in [mhe,mpc]:
    # Parameters with bounds
    m.K1 = m.FV(value=0.607,lb=0.1,ub=1.0)
    m.K2 = m.FV(value=0.293,lb=0.1,ub=1.0)
    m.K3 = m.FV(value=0.24,lb=0.1,ub=1.0)
    m.tau12 = m.FV(value=192,lb=100,ub=200)
    m.tau3 = m.FV(value=15,lb=10,ub=20)
    m.Ta = m.Param(value=23.0) # degC

    m.Q1 = m.MV(value=0,lb=0,ub=100,name='q1')
    m.Q2 = m.MV(value=0,lb=0,ub=100,name='q2')

    # Heater temperatures
    m.TH1 = m.SV(value=T1m[0])
    m.TH2 = m.SV(value=T2m[0])
    # Sensor temperatures
    m.TC1 = m.CV(value=T1m[0],name='tc1')
    m.TC2 = m.CV(value=T2m[0],name='tc2')

    # Temperature difference between two heaters
    m.DT = m.Intermediate(m.TH2-m.TH1)

    # Equations
    m.Equation(m.tau12*m.TH1.dt()+(m.TH1-m.Ta)==m.K1*m.Q1+m.K3*m.DT)
    m.Equation(m.tau12*m.TH2.dt()+(m.TH2-m.Ta)==m.K2*m.Q2-m.K3*m.DT)
    m.Equation(m.tau3*m.TC1.dt()+m.TC1==m.TH1)
    m.Equation(m.tau3*m.TC2.dt()+m.TC2==m.TH2)
user_23669410
Related