Infrequent Measurements In Gekko With Extra Simulation Points

I'm solving a more complex dynamic data reconciliation problem that is similar to the example problem posted in the APMonitor dynamic optimization course. Similar to this problem,

Solution 1:

You can combine infrequent measurements with more frequent simulation steps in Pandas with the join() function:

from gekko import GEKKO
import numpy as np
import pandas as pd
# measurements
t_data = [0,0.1,0.2,0.4,0.8,1]
x_data = [2.0,1.6,1.2,0.7,0.3,0.15]
df1 = pd.DataFrame({'time':t_data,'x':x_data})
df1.set_index('time', inplace=True)
# simulation time points
df2 = pd.DataFrame({'time':np.linspace(0,1,51)})
df2.set_index('time', inplace=True)
# merge dataframes
df = df2.join(df1,how='outer')
# get True (1) or False (0) for measurement
df['meas'] = (df['x'].values==df['x'].values).astype(int)
# replace NaN with zeros
df0 = df.fillna(value=0)
print(df.head(15))

This gives the combined set with NaN for the missing data points.

        x  meas
time           
0.002.010.02NaN00.04NaN00.06NaN00.08NaN00.101.610.12NaN00.14NaN00.16NaN00.18NaN00.201.210.22NaN00.24NaN00.26NaN00.28NaN0

You can then create a modified objective function that uses meas to control when to minimize the simulation and measurement only at the defined points.

m = GEKKO(remote=False)
m.time = df.index.values
x = m.Var(value=x_data[0])     # fit to measurement
xd = m.Param(df0['x'].values)
meas = m.Param(df0['meas'].values)
k = m.FV(); k.STATUS = 1# adjustable parameter
m.Equation(x.dt()== -k * x) # differential equation
m.Minimize(meas*(x-xd)**2)
m.options.SOLVER = 1# APOPT solver
m.options.IMODE = 5# dynamic estimation
m.options.NODES = 2# collocation nodes
m.solve(disp=True)   # display solver output
k = k.value[0]

A plot shows the results. I included more data points because NODES=2 isn't very accurate. A switch to NODES=3 helps to improve the accuracy much more than adding more points.

import numpy as np
import matplotlib.pyplot as plt  # plot solution
plt.plot(m.time,x.value,'bo',\
         label='Predicted (k='+str(np.round(k,2))+')')
plt.plot(m.time,df['x'].values,'rs',label='Measured')
# plot exact solution
t = np.linspace(0,1); xe = 2*np.exp(-k*t)
plt.plot(t,xe,'k:',label='Exact Solution')
plt.legend()
plt.xlabel('Time'), plt.ylabel('Value')
plt.ylim([-0.2,2.2])
plt.show()