Parameter Estimation Example 2: Kinetic Parameter Estimation#
This example demonstrates kinetic parameter estimation for chemical reaction systems.
Mathematical Formulation#
For a chemical reaction system, we typically have:
Material Balance: $\(\frac{dC_i}{dt} = \sum_j \nu_{ij} r_j(C, T, k)\)$
Rate Expression: $\(r_j = k_j f_j(C, T)\)$
Parameter Estimation Objective: $\(\min_{k} \sum_{i,t} w_{i,t} \left(C_{i,\text{exp}}(t) - C_{i,\text{model}}(t)\right)^2\)$
where:
\(C_i\) = concentration of species \(i\)
\(\nu_{ij}\) = stoichiometric coefficient
\(r_j\) = reaction rate
\(k_j\) = kinetic parameters (to be estimated)
\(w_{i,t}\) = weighting factors
Here we work through an example of kinetic parameter estimation.
First we simulate the kinematic behaviour
import pyomo.environ as pyo
from pyomo.dae import ContinuousSet, DerivativeVar, Simulator
import scipy
m = pyo.ConcreteModel()
m.t = ContinuousSet(bounds=(0,1))
m.a = pyo.Var(m.t)
m.b = pyo.Var(m.t)
m.k1 = pyo.Param(initialize=5)
m.k2 = pyo.Param(initialize=1)
m.dadt = DerivativeVar(m.a)
m.dbdt = DerivativeVar(m.b)
m.a[0].fix(1)
m.b[0].fix(0)
def _da(m, t):
return m.dadt[t] == -m.k1*m.a[t]
m.da_con = pyo.Constraint(m.t, rule=_da)
def _db(m, t):
return m.dbdt[t] == m.k1*m.a[t] - m.k2*m.b[t]
m.db_con = pyo.Constraint(m.t, rule=_db)
mysim = Simulator(m, package='scipy')
tsim, profiles = mysim.simulate(integrator='vode', numpoints=100)
import matplotlib.pyplot as plt
varorder = mysim.get_variable_order()
for idx, v in enumerate(varorder):
plt.plot(tsim, profiles[:, idx], label=v)
plt.xlabel('t')
plt.legend(loc='best')
plt.show()
plt.show()
Now let’s estimate the parameters
import pyomo.environ as pyo
from pyomo.dae import ContinuousSet, DerivativeVar
a_conc = {0.1:0.606, 0.2:0.368, 0.3:0.223, 0.4:0.135, 0.5:0.082,
0.6:0.05, 0.7:0.03, 0.8:0.018, 0.9:0.011, 1.0:0.007}
b_conc = {0.1:0.373, 0.2:0.564, 0.3:0.647, 0.4:0.669, 0.5:0.656,
0.6:0.624, 0.7:0.583, 0.8:0.539, 0.9:0.494, 1.0:0.451}
m = pyo.ConcreteModel()
m.meas_time = pyo.Set(initialize=sorted(a_conc.keys()),ordered=True)
m.ameas = pyo.Param(m.meas_time, initialize=a_conc)
m.bmeas = pyo.Param(m.meas_time, initialize=b_conc)
m.time = ContinuousSet(initialize=m.meas_time, bounds=(0,1))
m.a = pyo.Var(m.time, bounds=(0,1))
m.b = pyo.Var(m.time, bounds=(0,1))
m.dadt = DerivativeVar(m.a)
m.dbdt = DerivativeVar(m.b)
m.k1 = pyo.Var()
m.k2 = pyo.Var()
def _a_diffeq(m,t):
return m.dadt[t] == -m.k1*m.a[t]
m.a_diffeq = pyo.Constraint(m.time, rule=_a_diffeq)
def _b_diffeq(m,t):
return m.dbdt[t] == m.k1*m.a[t] - m.k2*m.b[t]
m.b_diffeq = pyo.Constraint(m.time, rule=_b_diffeq)
m.ainit = pyo.Constraint(expr=m.a[0]==1)
m.binit = pyo.Constraint(expr=m.b[0]==0)
def _obj(m):
return sum((m.a[t]-m.ameas[t])**2+(m.b[t]-m.bmeas[t])**2 for t in m.meas_time)
m.obj = pyo.Objective(rule=_obj)
discretizer = pyo.TransformationFactory('dae.collocation')
discretizer.apply_to(m,nfe=10,ncp=3,scheme='LAGRANGE-RADAU')
solver = pyo.SolverFactory('ipopt')
solver.solve(m, tee=True)
print('k1= '+str(pyo.value(m.k1)))
print('k2= '+str(pyo.value(m.k2)))
meas_time = list(m.meas_time)
a_meas = [pyo.value(m.ameas[i]) for i in m.meas_time]
b_meas = [pyo.value(m.bmeas[i]) for i in m.meas_time]
t = list(m.time)
a = [pyo.value(m.a[i]) for i in m.time]
b = [pyo.value(m.b[i]) for i in m.time]
import matplotlib.pyplot as plt
plt.plot(t,a,label='A')
plt.plot(t,b,label='B')
plt.plot(meas_time,a_meas,'o')
plt.plot(meas_time,b_meas,'o')
plt.legend(loc='best')
plt.xlabel('t')
plt.ylabel('concentration')
plt.title('Kinetic Parameter Estimation')
plt.show()
Ipopt 3.14.19:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.1.
Number of nonzeros in equality constraint Jacobian...: 550
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 82
Total number of variables............................: 126
variables with only lower bounds: 0
variables with lower and upper bounds: 62
variables with only upper bounds: 0
Total number of equality constraints.................: 124
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 3.6615861e+00 9.90e-01 5.94e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.6629897e+00 9.88e-01 7.45e-01 -1.7 1.87e+02 - 1.01e-02 2.51e-03h 4
2 3.6603131e+00 9.75e-01 1.30e+00 -1.7 1.72e+02 - 1.69e-02 1.30e-02f 2
3 3.6608441e+00 9.74e-01 1.33e+01 -1.7 2.68e+02 - 5.93e-02 3.38e-04h 6
4 3.5037068e+00 8.92e-01 1.43e+01 -1.7 1.09e+02 - 4.56e-02 8.43e-02f 1
5 3.4924703e+00 8.73e-01 3.43e+02 -1.7 1.49e+02 - 5.89e-01 2.12e-02h 3
6 3.3768345e+00 7.61e-01 6.01e+02 -1.7 8.52e+01 - 7.87e-02 1.29e-01h 2
7 3.3421277e+00 7.11e-01 6.56e+02 -1.7 6.92e+01 - 4.48e-01 6.53e-02H 1
8 3.3131009e+00 6.49e-01 7.61e+02 -1.7 9.22e+01 - 3.90e-02 8.81e-02h 3
9 3.2866091e+00 5.80e-01 7.80e+02 -1.7 1.21e+02 - 2.37e-01 1.06e-01h 3
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 3.0889266e+00 7.63e-01 1.80e+02 -1.7 6.41e+01 - 5.95e-01 9.90e-01H 1
11 3.1019026e+00 2.49e-01 2.25e+02 -1.7 2.44e+01 - 3.31e-01 9.90e-01h 1
12 2.6900564e+00 1.11e+01 2.08e+04 -1.7 8.92e+01 - 3.72e-01 9.90e-01h 1
13 2.5940882e+00 1.04e+01 2.15e+04 -1.7 2.48e+02 - 3.77e-02 3.89e-02h 1
14 2.5828122e+00 1.03e+01 1.94e+04 -1.7 7.45e+01 -2.0 7.81e-03 6.23e-03h 1
15 2.6632269e+00 1.26e+01 3.37e+05 -1.7 1.21e+02 - 9.94e-03 2.69e-01f 1
16 2.7587274e+00 1.19e+01 3.93e+05 -1.7 3.75e+01 -0.7 9.59e-03 6.72e-02h 1
17 2.7634018e+00 1.19e+01 3.32e+05 -1.7 4.07e+01 -1.1 4.56e-02 1.84e-03h 1
18 2.9029434e+00 1.15e+01 2.76e+05 -1.7 1.45e+02 -1.6 5.69e-02 2.07e-02h 1
19 2.9064261e+00 1.15e+01 1.94e+05 -1.7 1.46e+01 - 3.65e-01 1.98e-03h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 2.9777357e+00 1.09e+01 2.62e+05 -1.7 9.65e+00 - 1.39e-01 5.63e-02h 1
21 2.9781338e+00 1.08e+01 2.63e+05 -1.7 8.00e+00 - 2.30e-03 3.28e-04h 2
22 3.0173367e+00 1.05e+01 2.25e+05 -1.7 6.47e+00 - 4.77e-04 3.45e-02f 2
23 3.0623429e+00 1.01e+01 4.69e+05 -1.7 6.66e+00 - 3.29e-01 3.84e-02h 3
24 3.0978359e+00 9.78e+00 4.85e+05 -1.7 6.78e+00 - 7.45e-02 2.96e-02h 3
25 3.1115803e+00 9.67e+00 6.91e+05 -1.7 6.84e+00 - 3.48e-01 1.13e-02h 4
26 3.2162643e+00 8.90e+00 6.63e+05 -1.7 6.85e+00 - 1.40e-01 8.01e-02h 1
27 3.2173796e+00 8.89e+00 6.66e+05 -1.7 7.20e+00 - 8.12e-03 8.57e-04h 1
28 3.2201921e+00 8.83e+00 6.59e+05 -1.7 3.29e+01 - 5.29e-05 6.29e-03f 1
29 3.2628228e+00 8.10e+00 5.89e+05 -1.7 3.17e+01 - 4.70e-02 8.35e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 3.2633786e+00 8.09e+00 5.96e+05 -1.7 2.40e+01 - 1.97e-02 7.38e-04h 1
31 3.2630460e+00 8.08e+00 5.95e+05 -1.7 4.97e+01 - 1.43e-05 1.45e-03f 1
32 3.2712177e+00 7.58e+00 5.52e+05 -1.7 4.12e+01 - 4.76e-02 6.18e-02f 1
33 3.2713519e+00 7.57e+00 5.54e+05 -1.7 3.50e+01 - 7.06e-03 4.92e-04h 1
34 3.1338654e+00 3.99e+00 4.64e+05 -1.7 3.06e+02 - 1.22e-04 9.56e-02h 1
35 3.1335636e+00 3.89e+00 4.55e+05 -1.7 6.09e+01 - 2.86e-02 2.32e-02h 1
36 3.1335174e+00 3.89e+00 4.55e+05 -1.7 6.01e+01 - 4.95e-04 3.26e-04h 1
37 3.1335143e+00 3.89e+00 4.76e+05 -1.7 4.48e+01 - 5.28e-02 6.28e-05h 1
38 3.1241856e+00 3.74e+00 4.43e+05 -1.7 5.23e+01 - 1.60e-04 3.83e-02f 1
39 3.1171846e+00 3.55e+00 4.16e+05 -1.7 4.84e+01 - 3.90e-02 5.14e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 3.1151184e+00 3.38e+00 3.78e+05 -1.7 4.54e+01 - 5.66e-04 4.92e-02h 1
41 3.1151671e+00 3.38e+00 3.79e+05 -1.7 4.23e+01 - 3.53e-03 4.13e-04h 1
42 3.1154077e+00 3.15e+00 3.28e+05 -1.7 4.45e+01 - 2.27e-04 6.71e-02f 1
43 3.1155580e+00 3.15e+00 3.32e+05 -1.7 4.14e+01 - 1.13e-02 4.63e-04h 1
44 3.1167287e+00 3.12e+00 3.26e+05 -1.7 4.37e+01 - 4.19e-05 8.90e-03h 1
45 3.1215623e+00 3.05e+00 3.27e+05 -1.7 4.35e+01 - 4.96e-02 2.54e-02f 1
46 3.1220771e+00 3.04e+00 3.26e+05 -1.7 4.41e+01 - 6.82e-05 1.62e-03h 1
47 3.1288276e+00 2.98e+00 3.15e+05 -1.7 4.42e+01 - 1.02e-02 2.09e-02f 1
48 3.1399751e+00 2.89e+00 2.93e+05 -1.7 4.64e+01 - 1.14e-04 3.23e-02f 1
49 3.1401918e+00 2.89e+00 2.93e+05 -1.7 5.09e+01 - 7.50e-04 3.17e-04h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 3.1487222e+00 2.84e+00 2.81e+05 -1.7 5.06e+01 - 5.51e-04 1.93e-02f 1
51 3.1489213e+00 2.83e+00 2.81e+05 -1.7 5.54e+01 - 7.22e-04 2.81e-04h 1
52 3.1573381e+00 2.80e+00 2.75e+05 -1.7 5.50e+01 - 5.17e-03 1.31e-02f 1
53 3.1707003e+00 2.74e+00 2.61e+05 -1.7 5.78e+01 - 1.58e-04 2.25e-02f 1
54 3.1715780e+00 2.74e+00 2.60e+05 -1.7 6.38e+01 - 1.13e-03 1.29e-03h 1
55 3.1802297e+00 2.71e+00 2.53e+05 -1.7 6.32e+01 - 4.01e-03 1.37e-02f 1
56 3.1895628e+00 2.68e+00 2.54e+05 -1.7 6.82e+01 - 2.38e-02 1.21e-02f 1
57 3.1915046e+00 2.67e+00 2.53e+05 -1.7 6.84e+01 - 9.13e-04 2.91e-03h 1
58 3.2055985e+00 2.63e+00 2.44e+05 -1.7 6.96e+01 - 7.32e-03 1.95e-02f 1
59 3.2058330e+00 2.63e+00 2.44e+05 -1.7 8.39e+01 - 1.43e-03 2.13e-04h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 3.2175221e+00 2.58e+00 2.32e+05 -1.7 6.77e+01 - 5.28e-04 2.05e-02f 1
61 3.2177385e+00 2.58e+00 2.33e+05 -1.7 7.99e+01 - 1.92e-03 2.46e-04h 1
62 3.2318555e+00 2.54e+00 2.37e+05 -1.7 7.45e+01 - 4.33e-02 1.91e-02f 1
63 3.2320651e+00 2.54e+00 2.37e+05 -1.7 9.03e+01 - 2.93e-03 2.01e-04h 1
64 3.2606227e+00 2.46e+00 2.13e+05 -1.7 7.38e+01 - 1.52e-04 4.35e-02f 1
65 3.2611845e+00 2.46e+00 2.16e+05 -1.7 1.13e+02 - 1.09e-02 4.56e-04h 1
66 3.2892997e+00 2.39e+00 1.92e+05 -1.7 8.25e+01 - 9.64e-05 4.36e-02f 1
67 3.2898504e+00 2.38e+00 1.93e+05 -1.7 1.22e+02 - 3.87e-03 5.55e-04h 1
68 3.3218381e+00 2.32e+00 1.75e+05 -1.7 1.10e+02 - 6.21e-03 3.84e-02f 1
69 3.3938561e+00 2.04e+00 6.58e+04 -1.7 7.37e+01 - 4.32e-04 4.11e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 3.4442615e+00 2.04e+00 6.60e+04 -1.7 9.17e+02 - 2.28e-04 1.42e-03f 1
71 3.4448998e+00 2.04e+00 5.99e+04 -1.7 2.84e+02 - 2.60e-02 7.50e-05h 1
72 4.0533179e+00 1.53e+00 7.78e+04 -1.7 1.36e+02 - 2.12e-01 3.87e-01H 1
73 4.2336999e+00 1.41e+00 9.34e+04 -1.7 6.13e+00 -1.2 5.56e-01 7.75e-02h 1
74 4.3639078e+00 1.33e+00 7.01e+04 -1.7 1.08e+02 - 1.70e-03 5.64e-02h 1
75 4.6886245e+00 1.16e+00 3.63e+05 -1.7 2.98e+01 -1.7 1.00e+00 1.25e-01h 1
76 5.4005466e+00 8.72e-01 2.05e+05 -1.7 4.26e+01 - 1.24e-02 2.44e-01h 1
77 5.5621033e+00 8.24e-01 4.01e+05 -1.7 3.15e+00 -1.3 6.16e-01 5.43e-02h 1
78 8.3100244e+00 2.07e-01 1.62e+05 -1.7 3.18e+01 - 3.07e-03 7.88e-01h 1
79 1.0537702e+01 7.16e-02 3.28e+05 -1.7 6.58e+01 - 5.32e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 1.0009575e+01 7.87e-02 4.37e+05 -1.7 1.34e+02 - 3.74e-02 3.88e-01f 1
81 8.1592884e+00 4.42e-02 3.27e+05 -1.7 3.61e-01 -1.7 4.94e-01 1.00e+00f 1
82 6.6779695e+00 9.04e-02 2.18e+05 -1.7 1.09e+02 - 3.35e-01 2.47e-01f 3
83 3.7582674e+00 5.48e-01 4.23e+04 -1.7 5.19e+00 - 8.06e-01 9.88e-01f 1
84 3.6375713e+00 6.83e-02 5.70e-02 -1.7 1.00e+00 -2.2 1.00e+00 1.00e+00h 1
85 3.3189973e+00 7.33e-02 3.02e+05 -2.5 7.78e-01 -2.7 8.24e-01 1.00e+00h 1
86 3.1884913e+00 7.15e-02 3.00e-02 -2.5 2.35e+00 -3.2 1.00e+00 1.00e+00h 1
87 3.1596974e+00 2.95e-02 5.85e-02 -2.5 7.67e+00 -3.6 1.00e+00 1.00e+00h 1
88 3.1283114e+00 4.86e-01 1.49e-01 -2.5 3.64e+01 -4.1 1.00e+00 1.00e+00h 1
89 3.0508764e+00 8.07e-01 2.40e-01 -2.5 2.92e+01 -3.7 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
90 2.6664966e+00 7.19e+00 2.06e+00 -2.5 5.46e+01 -3.3 5.12e-01 1.00e+00h 1
91 2.7832045e+00 6.95e+00 1.81e+00 -2.5 2.60e+01 - 1.00e+00 1.00e+00h 1
92 1.2716500e+00 4.73e+00 2.41e+00 -2.5 2.87e+01 - 3.50e-02 4.38e-01h 1
93 1.4372239e+00 3.27e+00 1.07e+01 -2.5 9.74e+00 -1.0 2.47e-01 1.00e+00h 1
94 5.9591171e-01 2.08e-01 2.84e+00 -2.5 3.00e+00 -0.6 5.48e-01 1.00e+00h 1
95 4.4091074e-01 4.01e-02 4.23e-01 -2.5 6.37e-01 -0.2 8.30e-01 1.00e+00h 1
96 3.6146725e-01 2.22e+00 4.41e-01 -2.5 3.20e+00 - 1.00e+00 1.00e+00F 1
97 4.1014737e-02 2.17e+00 2.73e-01 -2.5 6.74e+00 - 1.00e+00 1.00e+00h 1
98 3.2376617e-03 3.07e-01 6.16e-02 -2.5 2.56e+00 - 1.00e+00 1.00e+00h 1
99 8.3326962e-03 3.85e-03 2.24e-03 -2.5 2.62e-01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
100 2.9806291e-04 1.09e-02 6.49e+03 -3.8 3.66e-01 - 9.76e-01 1.00e+00h 1
101 3.1011459e-05 4.45e-04 7.37e-04 -3.8 8.37e-02 - 1.00e+00 1.00e+00h 1
102 1.2002431e-06 6.73e-05 7.43e-05 -5.7 3.05e-02 - 1.00e+00 1.00e+00h 1
103 1.1912248e-06 1.29e-08 2.97e-08 -5.7 4.62e-04 - 1.00e+00 1.00e+00h 1
104 1.1862889e-06 1.14e-08 1.25e-08 -8.6 3.94e-04 - 1.00e+00 1.00e+00h 1
105 1.1862894e-06 1.20e-14 1.59e-14 -9.0 4.16e-07 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 105
(scaled) (unscaled)
Objective...............: 1.1862894024317846e-06 1.1862894024317846e-06
Dual infeasibility......: 1.5938567068277406e-14 1.5938567068277406e-14
Constraint violation....: 1.1990408665951691e-14 1.1990408665951691e-14
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 9.0909169534166345e-10 9.0909169534166345e-10
Overall NLP error.......: 9.0909169534166345e-10 9.0909169534166345e-10
Number of objective function evaluations = 147
Number of objective gradient evaluations = 106
Number of equality constraint evaluations = 147
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 106
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 105
Total seconds in IPOPT = 0.242
EXIT: Optimal Solution Found.
k1= 5.003509348488083
k2= 0.9999977260457749
import pyomo.environ as pyo
from pyomo.dae import ContinuousSet, DerivativeVar
a_conc = {0.1:0.606, 0.2:0.368, 0.3:0.223, 0.4:0.135, 0.5:0.082,
0.6:0.05, 0.7:0.03, 0.8:0.018, 0.9:0.011, 1.0:0.007}
b_conc = {0.1:0.373, 0.2:0.564, 0.3:0.647, 0.4:0.669, 0.5:0.656,
0.6:0.624, 0.7:0.583, 0.8:0.539, 0.9:0.494, 1.0:0.451}
m = pyo.ConcreteModel()
m.meas_time = pyo.Set(initialize=sorted(a_conc.keys()),ordered=True)
m.ameas = pyo.Param(m.meas_time, initialize=a_conc)
m.bmeas = pyo.Param(m.meas_time, initialize=b_conc)
m.time = ContinuousSet(initialize=m.meas_time, bounds=(0,1))
m.a = pyo.Var(m.time, bounds=(0,1))
m.b = pyo.Var(m.time, bounds=(0,1))
m.dadt = DerivativeVar(m.a)
m.dbdt = DerivativeVar(m.b)
m.k1 = pyo.Var()
m.k2 = pyo.Var()
def _a_diffeq(m,t):
return m.dadt[t] == -m.k1*m.a[t]
m.a_diffeq = pyo.Constraint(m.time, rule=_a_diffeq)
def _b_diffeq(m,t):
return m.dbdt[t] == m.k1*m.a[t] - m.k2*m.b[t]
m.b_diffeq = pyo.Constraint(m.time, rule=_b_diffeq)
m.ainit = pyo.Constraint(expr=m.a[0]==1)
m.binit = pyo.Constraint(expr=m.b[0]==0)
def _obj(m):
return sum((m.a[t]-m.ameas[t])**2+(m.b[t]-m.bmeas[t])**2 for t in m.meas_time)
m.obj = pyo.Objective(rule=_obj)
discretizer = pyo.TransformationFactory('dae.collocation')
discretizer.apply_to(m,nfe=10,ncp=3,scheme='LAGRANGE-RADAU')
solver = pyo.SolverFactory('ipopt')
solver.solve(m, tee=True)
print('k1= '+str(pyo.value(m.k1)))
print('k2= '+str(pyo.value(m.k2)))
meas_time = list(m.meas_time)
a_meas = [pyo.value(m.ameas[i]) for i in m.meas_time]
b_meas = [pyo.value(m.bmeas[i]) for i in m.meas_time]
t = list(m.time)
a = [pyo.value(m.a[i]) for i in m.time]
b = [pyo.value(m.b[i]) for i in m.time]
import matplotlib.pyplot as plt
plt.plot(t,a,label='A')
plt.plot(t,b,label='B')
plt.plot(meas_time,a_meas,'o')
plt.plot(meas_time,b_meas,'o')
plt.legend(loc='best')
plt.xlabel('t')
plt.ylabel('concentration')
plt.title('Kinetic Parameter Estimation')
plt.show()
Ipopt 3.14.19:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.1.
Number of nonzeros in equality constraint Jacobian...: 550
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 82
Total number of variables............................: 126
variables with only lower bounds: 0
variables with lower and upper bounds: 62
variables with only upper bounds: 0
Total number of equality constraints.................: 124
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 3.6615861e+00 9.90e-01 5.94e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.6629897e+00 9.88e-01 7.45e-01 -1.7 1.87e+02 - 1.01e-02 2.51e-03h 4
2 3.6603131e+00 9.75e-01 1.30e+00 -1.7 1.72e+02 - 1.69e-02 1.30e-02f 2
3 3.6608441e+00 9.74e-01 1.33e+01 -1.7 2.68e+02 - 5.93e-02 3.38e-04h 6
4 3.5037068e+00 8.92e-01 1.43e+01 -1.7 1.09e+02 - 4.56e-02 8.43e-02f 1
5 3.4924703e+00 8.73e-01 3.43e+02 -1.7 1.49e+02 - 5.89e-01 2.12e-02h 3
6 3.3768345e+00 7.61e-01 6.01e+02 -1.7 8.52e+01 - 7.87e-02 1.29e-01h 2
7 3.3421277e+00 7.11e-01 6.56e+02 -1.7 6.92e+01 - 4.48e-01 6.53e-02H 1
8 3.3131009e+00 6.49e-01 7.61e+02 -1.7 9.22e+01 - 3.90e-02 8.81e-02h 3
9 3.2866091e+00 5.80e-01 7.80e+02 -1.7 1.21e+02 - 2.37e-01 1.06e-01h 3
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 3.0889266e+00 7.63e-01 1.80e+02 -1.7 6.41e+01 - 5.95e-01 9.90e-01H 1
11 3.1019026e+00 2.49e-01 2.25e+02 -1.7 2.44e+01 - 3.31e-01 9.90e-01h 1
12 2.6900564e+00 1.11e+01 2.08e+04 -1.7 8.92e+01 - 3.72e-01 9.90e-01h 1
13 2.5940882e+00 1.04e+01 2.15e+04 -1.7 2.48e+02 - 3.77e-02 3.89e-02h 1
14 2.5828122e+00 1.03e+01 1.94e+04 -1.7 7.45e+01 -2.0 7.81e-03 6.23e-03h 1
15 2.6632269e+00 1.26e+01 3.37e+05 -1.7 1.21e+02 - 9.94e-03 2.69e-01f 1
16 2.7587274e+00 1.19e+01 3.93e+05 -1.7 3.75e+01 -0.7 9.59e-03 6.72e-02h 1
17 2.7634018e+00 1.19e+01 3.32e+05 -1.7 4.07e+01 -1.1 4.56e-02 1.84e-03h 1
18 2.9029434e+00 1.15e+01 2.76e+05 -1.7 1.45e+02 -1.6 5.69e-02 2.07e-02h 1
19 2.9064261e+00 1.15e+01 1.94e+05 -1.7 1.46e+01 - 3.65e-01 1.98e-03h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 2.9777357e+00 1.09e+01 2.62e+05 -1.7 9.65e+00 - 1.39e-01 5.63e-02h 1
21 2.9781338e+00 1.08e+01 2.63e+05 -1.7 8.00e+00 - 2.30e-03 3.28e-04h 2
22 3.0173367e+00 1.05e+01 2.25e+05 -1.7 6.47e+00 - 4.77e-04 3.45e-02f 2
23 3.0623429e+00 1.01e+01 4.69e+05 -1.7 6.66e+00 - 3.29e-01 3.84e-02h 3
24 3.0978359e+00 9.78e+00 4.85e+05 -1.7 6.78e+00 - 7.45e-02 2.96e-02h 3
25 3.1115803e+00 9.67e+00 6.91e+05 -1.7 6.84e+00 - 3.48e-01 1.13e-02h 4
26 3.2162643e+00 8.90e+00 6.63e+05 -1.7 6.85e+00 - 1.40e-01 8.01e-02h 1
27 3.2173796e+00 8.89e+00 6.66e+05 -1.7 7.20e+00 - 8.12e-03 8.57e-04h 1
28 3.2201921e+00 8.83e+00 6.59e+05 -1.7 3.29e+01 - 5.29e-05 6.29e-03f 1
29 3.2628228e+00 8.10e+00 5.89e+05 -1.7 3.17e+01 - 4.70e-02 8.35e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 3.2633786e+00 8.09e+00 5.96e+05 -1.7 2.40e+01 - 1.97e-02 7.38e-04h 1
31 3.2630460e+00 8.08e+00 5.95e+05 -1.7 4.97e+01 - 1.43e-05 1.45e-03f 1
32 3.2712177e+00 7.58e+00 5.52e+05 -1.7 4.12e+01 - 4.76e-02 6.18e-02f 1
33 3.2713519e+00 7.57e+00 5.54e+05 -1.7 3.50e+01 - 7.06e-03 4.92e-04h 1
34 3.1338654e+00 3.99e+00 4.64e+05 -1.7 3.06e+02 - 1.22e-04 9.56e-02h 1
35 3.1335636e+00 3.89e+00 4.55e+05 -1.7 6.09e+01 - 2.86e-02 2.32e-02h 1
36 3.1335174e+00 3.89e+00 4.55e+05 -1.7 6.01e+01 - 4.95e-04 3.26e-04h 1
37 3.1335143e+00 3.89e+00 4.76e+05 -1.7 4.48e+01 - 5.28e-02 6.28e-05h 1
38 3.1241856e+00 3.74e+00 4.43e+05 -1.7 5.23e+01 - 1.60e-04 3.83e-02f 1
39 3.1171846e+00 3.55e+00 4.16e+05 -1.7 4.84e+01 - 3.90e-02 5.14e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 3.1151184e+00 3.38e+00 3.78e+05 -1.7 4.54e+01 - 5.66e-04 4.92e-02h 1
41 3.1151671e+00 3.38e+00 3.79e+05 -1.7 4.23e+01 - 3.53e-03 4.13e-04h 1
42 3.1154077e+00 3.15e+00 3.28e+05 -1.7 4.45e+01 - 2.27e-04 6.71e-02f 1
43 3.1155580e+00 3.15e+00 3.32e+05 -1.7 4.14e+01 - 1.13e-02 4.63e-04h 1
44 3.1167287e+00 3.12e+00 3.26e+05 -1.7 4.37e+01 - 4.19e-05 8.90e-03h 1
45 3.1215623e+00 3.05e+00 3.27e+05 -1.7 4.35e+01 - 4.96e-02 2.54e-02f 1
46 3.1220771e+00 3.04e+00 3.26e+05 -1.7 4.41e+01 - 6.82e-05 1.62e-03h 1
47 3.1288276e+00 2.98e+00 3.15e+05 -1.7 4.42e+01 - 1.02e-02 2.09e-02f 1
48 3.1399751e+00 2.89e+00 2.93e+05 -1.7 4.64e+01 - 1.14e-04 3.23e-02f 1
49 3.1401918e+00 2.89e+00 2.93e+05 -1.7 5.09e+01 - 7.50e-04 3.17e-04h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 3.1487222e+00 2.84e+00 2.81e+05 -1.7 5.06e+01 - 5.51e-04 1.93e-02f 1
51 3.1489213e+00 2.83e+00 2.81e+05 -1.7 5.54e+01 - 7.22e-04 2.81e-04h 1
52 3.1573381e+00 2.80e+00 2.75e+05 -1.7 5.50e+01 - 5.17e-03 1.31e-02f 1
53 3.1707003e+00 2.74e+00 2.61e+05 -1.7 5.78e+01 - 1.58e-04 2.25e-02f 1
54 3.1715780e+00 2.74e+00 2.60e+05 -1.7 6.38e+01 - 1.13e-03 1.29e-03h 1
55 3.1802297e+00 2.71e+00 2.53e+05 -1.7 6.32e+01 - 4.01e-03 1.37e-02f 1
56 3.1895628e+00 2.68e+00 2.54e+05 -1.7 6.82e+01 - 2.38e-02 1.21e-02f 1
57 3.1915046e+00 2.67e+00 2.53e+05 -1.7 6.84e+01 - 9.13e-04 2.91e-03h 1
58 3.2055985e+00 2.63e+00 2.44e+05 -1.7 6.96e+01 - 7.32e-03 1.95e-02f 1
59 3.2058330e+00 2.63e+00 2.44e+05 -1.7 8.39e+01 - 1.43e-03 2.13e-04h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 3.2175221e+00 2.58e+00 2.32e+05 -1.7 6.77e+01 - 5.28e-04 2.05e-02f 1
61 3.2177385e+00 2.58e+00 2.33e+05 -1.7 7.99e+01 - 1.92e-03 2.46e-04h 1
62 3.2318555e+00 2.54e+00 2.37e+05 -1.7 7.45e+01 - 4.33e-02 1.91e-02f 1
63 3.2320651e+00 2.54e+00 2.37e+05 -1.7 9.03e+01 - 2.93e-03 2.01e-04h 1
64 3.2606227e+00 2.46e+00 2.13e+05 -1.7 7.38e+01 - 1.52e-04 4.35e-02f 1
65 3.2611845e+00 2.46e+00 2.16e+05 -1.7 1.13e+02 - 1.09e-02 4.56e-04h 1
66 3.2892997e+00 2.39e+00 1.92e+05 -1.7 8.25e+01 - 9.64e-05 4.36e-02f 1
67 3.2898504e+00 2.38e+00 1.93e+05 -1.7 1.22e+02 - 3.87e-03 5.55e-04h 1
68 3.3218381e+00 2.32e+00 1.75e+05 -1.7 1.10e+02 - 6.21e-03 3.84e-02f 1
69 3.3938561e+00 2.04e+00 6.58e+04 -1.7 7.37e+01 - 4.32e-04 4.11e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 3.4442615e+00 2.04e+00 6.60e+04 -1.7 9.17e+02 - 2.28e-04 1.42e-03f 1
71 3.4448998e+00 2.04e+00 5.99e+04 -1.7 2.84e+02 - 2.60e-02 7.50e-05h 1
72 4.0533179e+00 1.53e+00 7.78e+04 -1.7 1.36e+02 - 2.12e-01 3.87e-01H 1
73 4.2336999e+00 1.41e+00 9.34e+04 -1.7 6.13e+00 -1.2 5.56e-01 7.75e-02h 1
74 4.3639078e+00 1.33e+00 7.01e+04 -1.7 1.08e+02 - 1.70e-03 5.64e-02h 1
75 4.6886245e+00 1.16e+00 3.63e+05 -1.7 2.98e+01 -1.7 1.00e+00 1.25e-01h 1
76 5.4005466e+00 8.72e-01 2.05e+05 -1.7 4.26e+01 - 1.24e-02 2.44e-01h 1
77 5.5621033e+00 8.24e-01 4.01e+05 -1.7 3.15e+00 -1.3 6.16e-01 5.43e-02h 1
78 8.3100244e+00 2.07e-01 1.62e+05 -1.7 3.18e+01 - 3.07e-03 7.88e-01h 1
79 1.0537702e+01 7.16e-02 3.28e+05 -1.7 6.58e+01 - 5.32e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 1.0009575e+01 7.87e-02 4.37e+05 -1.7 1.34e+02 - 3.74e-02 3.88e-01f 1
81 8.1592885e+00 4.42e-02 3.27e+05 -1.7 3.61e-01 -1.7 4.94e-01 1.00e+00f 1
82 6.6779695e+00 9.04e-02 2.18e+05 -1.7 1.09e+02 - 3.35e-01 2.47e-01f 3
83 3.7582675e+00 5.48e-01 4.23e+04 -1.7 5.19e+00 - 8.06e-01 9.88e-01f 1
84 3.6375713e+00 6.83e-02 5.70e-02 -1.7 1.00e+00 -2.2 1.00e+00 1.00e+00h 1
85 3.3189973e+00 7.33e-02 3.02e+05 -2.5 7.78e-01 -2.7 8.24e-01 1.00e+00h 1
86 3.1884913e+00 7.15e-02 3.00e-02 -2.5 2.35e+00 -3.2 1.00e+00 1.00e+00h 1
87 3.1596974e+00 2.95e-02 5.85e-02 -2.5 7.67e+00 -3.6 1.00e+00 1.00e+00h 1
88 3.1283114e+00 4.86e-01 1.49e-01 -2.5 3.64e+01 -4.1 1.00e+00 1.00e+00h 1
89 3.0508764e+00 8.07e-01 2.40e-01 -2.5 2.92e+01 -3.7 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
90 2.6664966e+00 7.19e+00 2.06e+00 -2.5 5.46e+01 -3.3 5.12e-01 1.00e+00h 1
91 2.7832045e+00 6.95e+00 1.81e+00 -2.5 2.60e+01 - 1.00e+00 1.00e+00h 1
92 1.2716511e+00 4.73e+00 2.41e+00 -2.5 2.87e+01 - 3.50e-02 4.38e-01h 1
93 1.4372148e+00 3.27e+00 1.07e+01 -2.5 9.74e+00 -1.0 2.47e-01 1.00e+00h 1
94 5.9592382e-01 2.07e-01 2.84e+00 -2.5 3.00e+00 -0.6 5.48e-01 1.00e+00h 1
95 4.4091992e-01 4.01e-02 4.23e-01 -2.5 6.37e-01 -0.2 8.30e-01 1.00e+00h 1
96 3.6167742e-01 2.22e+00 4.42e-01 -2.5 3.20e+00 - 1.00e+00 1.00e+00F 1
97 4.1155476e-02 2.17e+00 2.73e-01 -2.5 6.74e+00 - 1.00e+00 1.00e+00h 1
98 3.2286718e-03 3.08e-01 6.17e-02 -2.5 2.56e+00 - 1.00e+00 1.00e+00h 1
99 8.3331542e-03 3.87e-03 2.25e-03 -2.5 2.62e-01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
100 2.9800146e-04 1.09e-02 6.50e+03 -3.8 3.66e-01 - 9.76e-01 1.00e+00h 1
101 3.1011523e-05 4.45e-04 7.37e-04 -3.8 8.37e-02 - 1.00e+00 1.00e+00h 1
102 1.2002429e-06 6.73e-05 7.43e-05 -5.7 3.05e-02 - 1.00e+00 1.00e+00h 1
103 1.1912248e-06 1.29e-08 2.97e-08 -5.7 4.62e-04 - 1.00e+00 1.00e+00h 1
104 1.1862889e-06 1.14e-08 1.25e-08 -8.6 3.94e-04 - 1.00e+00 1.00e+00h 1
105 1.1862894e-06 1.24e-14 1.59e-14 -9.0 4.16e-07 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 105
(scaled) (unscaled)
Objective...............: 1.1862894024318386e-06 1.1862894024318386e-06
Dual infeasibility......: 1.5938567041807626e-14 1.5938567041807626e-14
Constraint violation....: 1.2434497875801753e-14 1.2434497875801753e-14
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 9.0909169534166076e-10 9.0909169534166076e-10
Overall NLP error.......: 9.0909169534166076e-10 9.0909169534166076e-10
Number of objective function evaluations = 147
Number of objective gradient evaluations = 106
Number of equality constraint evaluations = 147
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 106
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 105
Total seconds in IPOPT = 0.254
EXIT: Optimal Solution Found.
k1= 5.003509348488082
k2= 0.9999977260457752