Fit a nonlinear model with python and gekko

Question

I have a dataset of trees. In this dataset, I have the unique number of Plot, the sequence in the order when were take the data "Measurement" and the Height mean in meters and Age mean in years for the trees. Something like this: head of data

next, I define the model to predict the Height using the Age in this way:

Height = B0 * ((1 - exp(-B1 *Age))**B2)

My goal is to determinate the values of B0, B1 & B2 respectively. For this, I use the package gekko to find the parameters of the models with the next code:

num_p = data_gek.Plot.unique()
nmp = 5
number_p = (data_gek.Plot == num_p[nmp])

m = GEKKO()

xm = np.array(data_gek[number_p]['Age'])
x = m.Param(value=xm)

B0 = m.FV(value=38.2) #value=38.2
B0.STATUS = 1

B1 = m.FV(value=0.1) #value=0.1
B1.STATUS = 1

B2 = m.FV(value=2.08) #value=2.08
B2.STATUS = 1

ym =  np.array(data_gek[number_p]['Height'])
z = m.CV(value=ym)

y = m.Var()
m.Equation(y==B0 * ((1 - m.exp(-B1 *x))**B2))
m.Obj(((y-z)/z)**2)

m.options.IMODE = 2
m.options.SOLVER = 1
m.solve(disp=False)

print(B0.value[0],B1.value[0],B2.value[0])
#output 
27.787958561 0.0052435491089 0.21178326158

However, I don't sure that I make in the right way. Is it possible to do this without initial values in parameters? Because I used previous values for B0, B1, and B2 from literature.

If you gonna see my dataset and my process you could access this notebook in Google Colab.

Answer 1

Your script has just one problem. The definition of z needs to be a Param or MV type instead of a CV as z = m.Param(value=ym) because it is an input to your objective function.

You can also use the built-in objective if you define y as a CV instead of a Var. You just need to turn on the feedback status FSTATUS=1 to use an objective function that minimizes the difference between the measurements and model predictions. Here is a modified version of your script.

from gekko import GEKKO
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
url = 'http://apmonitor.com/pdc/uploads/Main/data_2nd_order.txt'
data = pd.read_csv(url)
m  = GEKKO()
xm = np.array(data['time'])
x  = m.Param(value=xm)
B0 = m.FV(1); B1 = m.FV(1); B2 = m.FV(1)
B0.STATUS = 1; B1.STATUS = 1; B2.STATUS = 1
ym =  np.array(data['output (y)'])
y  = m.CV(value=ym)
y.FSTATUS = 1
yi = m.Intermediate(B0 * ((1 - m.exp(-B1 *x))**B2))
m.Equation(y==yi)
m.options.IMODE = 2
m.options.SOLVER = 1
m.solve(disp=True)
print(B0.value[0],B1.value[0],B2.value[0])
plt.plot(xm,ym,'ro')
plt.plot(xm,y.value,'b--')
plt.show()