Dynamic parameter estimation of 4 parameters with 4 ODEs

user_19611701 0
Viewed 37

I would like to perform dynamic parameter estimation of 4 parameters in a model as shown in example 4 of this link: https://apmonitor.com/do/index.php/Main/DynamicEstimation. The difference I have with my application is there are 4 differential equations which collectively contain the 4 parameters. There is experimental data for only 1 of the equations. I have adjusted the code given in the example to work for this application, but it is returning 0 for each parameter.

How would I get this to work with 4 equations instead of 1? Here is the code so far:

from gekko import GEKKO
import matplotlib.pyplot as plt  # plot solution
import numpy as np

t_data = [1.22, 5.15, 23.67, 51.17, 74.58, 97.83, 118.97, 143.33, 166.73, 192.08, 222.83, 245.08, 266.67, 286.18, 309.90]
m4_data = [3634.8, 7035.9, 7797.8, 9351.4, 10041.0, 10674.6, 11339.5, 11115.7, 11225.1, 11465.4, 11383.2, 11456.5, 11506.8, 11683.6, 11588.2]


m = GEKKO()
m.time = t_data
m4 = m.CV(value=m4_data); m4.FSTATUS = 1  # fit to measurement
m1,m2,m3 = m.Array(m.Var,3,value=3)
a,b,c,d = m.Array(m.FV,4)
t = 50 
v = 7 
s = 10
#DiffEqs
#m1' = -aS*(b-m3/v)
#m2' = -cs*(b-m3/v)
#m3' = (-d*t/v)*m3 -m1'-m2'
#m4' = d*t*m3/v

m.Equations([m1.dt() == -a*s*(b-m3/v),m2.dt() == -c*s*(b-m3/v),m3.dt() ==(-d*t/v)*m3-m1.dt()-m2.dt() , m4.dt() == d*t*m3/v])

m.options.IMODE = 5
m.options.NODES = 5  
m.solve()   
print(a.value[0],b.value[0],c.value[0],d.value[0])
gekko python
1 Answers

Turn on the STATUS for the FV to switch:

0 = user-specified value (fixed constant)
1 = degree of freedom that the solver can adjust to minimize the objective

Here is one way to turn on the STATUS for each. Another way is to use a.STATUS=1, b.STATUS=1, c.STATUS=1, d.STATUS=1

p = m.Array(m.FV,4); a,b,c,d=p
for pi in p:
    pi.STATUS=1

Here is the complete code:

from gekko import GEKKO
import matplotlib.pyplot as plt  # plot solution
import numpy as np
t_data = [1.22, 5.15, 23.67, 51.17, 74.58,
          97.83, 118.97, 143.33, 166.73,
          192.08, 222.83, 245.08, 266.67,
          286.18, 309.90]
m4_data = [3634.8, 7035.9, 7797.8, 9351.4,
           10041.0, 10674.6, 11339.5, 11115.7,
           11225.1, 11465.4, 11383.2, 11456.5,
           11506.8, 11683.6, 11588.2]
m = GEKKO()
m.time = t_data
m4 = m.CV(value=m4_data); m4.FSTATUS = 1  # fit to measurement
m1,m2,m3 = m.Array(m.Var,3,value=3)
p = m.Array(m.FV,4); a,b,c,d=p
for pi in p:
    pi.STATUS=1
t = 50; v = 7; s = 10
m.Equations([m1.dt() == -a*s*(b-m3/v),
             m2.dt() == -c*s*(b-m3/v),
             m3.dt() ==(-d*t/v)*m3-m1.dt()-m2.dt(),
             m4.dt() == d*t*m3/v])
m.options.IMODE = 5
m.options.NODES = 5  
m.solve()   
print(a.value[0],b.value[0],c.value[0],d.value[0])

The solution is:

305601.02084 0.00010747038506 -305687.94961 17.323603876
user_23669410
Related