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.6632270e+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.9029433e+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.9777363e+00 1.09e+01 2.62e+05 -1.7 9.65e+00 - 1.39e-01 5.63e-02h 1
21 2.9781344e+00 1.08e+01 2.63e+05 -1.7 8.00e+00 - 2.30e-03 3.28e-04h 2
22 3.0173368e+00 1.05e+01 2.25e+05 -1.7 6.47e+00 - 4.77e-04 3.45e-02f 2
23 3.0623430e+00 1.01e+01 4.69e+05 -1.7 6.66e+00 - 3.29e-01 3.84e-02h 3
24 3.0978360e+00 9.78e+00 4.85e+05 -1.7 6.78e+00 - 7.45e-02 2.96e-02h 3
25 3.1115805e+00 9.67e+00 6.91e+05 -1.7 6.84e+00 - 3.48e-01 1.13e-02h 4
26 3.2162645e+00 8.90e+00 6.63e+05 -1.7 6.85e+00 - 1.40e-01 8.01e-02h 1
27 3.2173798e+00 8.89e+00 6.66e+05 -1.7 7.20e+00 - 8.13e-03 8.57e-04h 1
28 3.2201923e+00 8.83e+00 6.59e+05 -1.7 3.29e+01 - 5.29e-05 6.29e-03f 1
29 3.2628229e+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.2633787e+00 8.09e+00 5.96e+05 -1.7 2.40e+01 - 1.97e-02 7.38e-04h 1
31 3.2630461e+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.1335637e+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.1335144e+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.1715782e+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.1915039e+00 2.67e+00 2.53e+05 -1.7 6.84e+01 - 9.12e-04 2.90e-03h 1
58 3.2055982e+00 2.63e+00 2.44e+05 -1.7 6.96e+01 - 7.32e-03 1.95e-02f 1
59 3.2058327e+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.2175218e+00 2.58e+00 2.32e+05 -1.7 6.77e+01 - 5.28e-04 2.05e-02f 1
61 3.2177382e+00 2.58e+00 2.33e+05 -1.7 7.99e+01 - 1.92e-03 2.46e-04h 1
62 3.2318552e+00 2.54e+00 2.37e+05 -1.7 7.45e+01 - 4.33e-02 1.91e-02f 1
63 3.2320648e+00 2.54e+00 2.37e+05 -1.7 9.03e+01 - 2.93e-03 2.01e-04h 1
64 3.2606225e+00 2.46e+00 2.13e+05 -1.7 7.38e+01 - 1.52e-04 4.35e-02f 1
65 3.2611843e+00 2.46e+00 2.16e+05 -1.7 1.13e+02 - 1.09e-02 4.56e-04h 1
66 3.2892995e+00 2.39e+00 1.92e+05 -1.7 8.25e+01 - 9.64e-05 4.36e-02f 1
67 3.2898502e+00 2.38e+00 1.93e+05 -1.7 1.22e+02 - 3.87e-03 5.55e-04h 1
68 3.3218379e+00 2.32e+00 1.75e+05 -1.7 1.10e+02 - 6.21e-03 3.84e-02f 1
69 3.3938365e+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.4443210e+00 2.04e+00 6.60e+04 -1.7 9.12e+02 - 2.29e-04 1.43e-03f 1
71 3.4449596e+00 2.04e+00 6.00e+04 -1.7 2.86e+02 - 2.56e-02 7.44e-05h 1
72 4.0535875e+00 1.53e+00 7.81e+04 -1.7 1.36e+02 - 2.11e-01 3.87e-01H 1
73 4.2336486e+00 1.41e+00 9.26e+04 -1.7 6.12e+00 -1.2 5.54e-01 7.73e-02h 1
74 4.3632004e+00 1.33e+00 6.95e+04 -1.7 1.09e+02 - 1.69e-03 5.61e-02h 1
75 4.6869628e+00 1.16e+00 3.63e+05 -1.7 2.98e+01 -1.7 1.00e+00 1.24e-01h 1
76 5.3980665e+00 8.72e-01 2.05e+05 -1.7 4.26e+01 - 1.24e-02 2.43e-01h 1
77 5.5593064e+00 8.25e-01 4.02e+05 -1.7 3.16e+00 -1.3 6.18e-01 5.42e-02h 1
78 8.3227421e+00 2.09e-01 1.64e+05 -1.7 3.10e+01 - 3.25e-03 7.91e-01h 1
79 1.0538910e+01 7.14e-02 3.48e+05 -1.7 6.57e+01 - 5.06e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 1.0017940e+01 7.98e-02 4.35e+05 -1.7 1.65e+02 - 3.27e-02 3.22e-01f 1
81 8.0592267e+00 5.00e-02 4.41e+05 -1.7 5.07e-01 -1.7 3.38e-01 1.00e+00f 1
82 4.4250234e+00 6.34e-01 2.85e+05 -1.7 9.38e+01 - 3.55e-01 1.00e+00F 1
83 3.4756774e+00 7.32e-01 1.01e+05 -1.7 5.79e+01 - 6.46e-01 1.00e+00h 1
84 3.3967064e+00 2.64e-02 4.41e-02 -1.7 8.50e-01 -2.2 1.00e+00 1.00e+00h 1
85 3.2411986e+00 5.26e-02 3.63e+05 -3.8 6.40e-01 -2.7 8.17e-01 1.00e+00h 1
86 3.1731684e+00 6.41e-02 1.10e+04 -3.8 8.14e-01 -3.2 9.70e-01 1.00e+00h 1
87 3.1619840e+00 2.60e-02 3.59e-02 -3.8 2.19e+00 -3.6 1.00e+00 1.00e+00h 1
88 3.1586055e+00 1.48e-02 4.52e-02 -3.8 7.20e+00 -4.1 1.00e+00 1.00e+00h 1
89 3.1412468e+00 2.90e-01 6.33e-02 -3.8 3.34e+01 -4.6 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
90 3.1148750e+00 2.63e-01 6.35e-02 -3.8 2.15e+01 -4.2 1.00e+00 1.00e+00h 1
91 3.0854400e+00 1.60e-01 4.10e-02 -3.8 1.30e+01 -3.7 1.00e+00 1.00e+00h 1
92 3.0613106e+00 6.35e-02 1.63e-02 -3.8 6.81e+00 -3.3 1.00e+00 1.00e+00h 1
93 2.6891173e+00 1.98e+01 5.08e+00 -3.8 1.24e+02 -3.8 3.40e-01 1.00e+00h 1
94 2.7264403e+00 1.60e+01 4.10e+00 -3.8 5.44e+01 - 1.00e+00 2.05e-01h 1
95 2.9276478e+00 4.35e+00 2.41e+00 -3.8 5.13e+01 - 5.35e-01 1.00e+00h 1
96 3.0097281e+00 1.30e+00 1.81e+00 -3.8 3.69e+01 - 2.23e-01 8.66e-01h 1
97 2.1035595e+00 5.62e+01 5.27e+01 -3.8 1.20e+02 -3.4 9.00e-02 1.00e+00h 1
98 2.1744565e+00 4.61e+01 4.33e+01 -3.8 6.67e+01 - 2.33e-01 1.88e-01h 1
99 2.3916968e+00 1.88e+01 1.72e+01 -3.8 5.08e+01 - 1.87e-01 7.34e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
100 2.6368851e+00 4.55e+00 9.76e+00 -3.8 4.13e+01 - 1.00e+00 1.00e+00h 1
101 2.5637768e+00 1.70e+00 4.18e+00 -3.8 2.07e+01 - 3.82e-01 6.92e-01h 1
102 2.2706205e+00 2.98e-01 2.32e+01 -3.8 3.81e+00 - 1.51e-01 9.58e-01h 1
103 1.9720365e+00 3.04e-01 1.71e+01 -3.8 1.67e+01 -2.0 3.82e-01 2.75e-01h 1
104 2.0448003e+00 1.45e-01 3.66e+01 -3.8 1.40e+00 - 2.26e-01 1.00e+00h 1
105 1.7496678e+00 1.67e-01 1.42e+02 -3.8 2.13e+00 - 7.37e-02 1.00e+00h 1
106 1.1518994e+00 1.58e-01 7.06e+00 -3.8 1.28e+00 - 3.02e-01 1.00e+00h 1
107 3.5329621e-01 2.85e-01 3.31e-01 -3.8 1.66e+00 - 4.27e-01 1.00e+00h 1
108 1.4564439e-02 2.18e-01 8.26e-02 -3.8 7.81e-01 - 5.34e-01 1.00e+00h 1
109 3.6173861e-05 8.54e-03 4.02e-03 -3.8 2.63e-01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
110 3.3251411e-05 5.22e-05 1.25e-05 -3.8 2.10e-02 - 1.00e+00 1.00e+00h 1
111 1.1999203e-06 7.06e-05 7.49e-05 -5.7 3.08e-02 - 1.00e+00 1.00e+00h 1
112 1.1912245e-06 1.31e-08 3.12e-08 -5.7 4.68e-04 - 1.00e+00 1.00e+00h 1
113 1.1862889e-06 1.14e-08 1.25e-08 -8.6 3.94e-04 - 1.00e+00 1.00e+00h 1
114 1.1862894e-06 1.20e-14 1.59e-14 -9.0 4.16e-07 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 114
(scaled) (unscaled)
Objective...............: 1.1862894024320862e-06 1.1862894024320862e-06
Dual infeasibility......: 1.5938143090650587e-14 1.5938143090650587e-14
Constraint violation....: 1.1990408665951691e-14 1.1990408665951691e-14
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 9.0909169532454733e-10 9.0909169532454733e-10
Overall NLP error.......: 9.0909169532454733e-10 9.0909169532454733e-10
Number of objective function evaluations = 152
Number of objective gradient evaluations = 115
Number of equality constraint evaluations = 152
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 115
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 114
Total seconds in IPOPT = 0.256
EXIT: Optimal Solution Found.
k1= 5.003509348488082
k2= 0.9999977260457744
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.6632270e+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.9029433e+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.9777363e+00 1.09e+01 2.62e+05 -1.7 9.65e+00 - 1.39e-01 5.63e-02h 1
21 2.9781344e+00 1.08e+01 2.63e+05 -1.7 8.00e+00 - 2.30e-03 3.28e-04h 2
22 3.0173368e+00 1.05e+01 2.25e+05 -1.7 6.47e+00 - 4.77e-04 3.45e-02f 2
23 3.0623430e+00 1.01e+01 4.69e+05 -1.7 6.66e+00 - 3.29e-01 3.84e-02h 3
24 3.0978360e+00 9.78e+00 4.85e+05 -1.7 6.78e+00 - 7.45e-02 2.96e-02h 3
25 3.1115805e+00 9.67e+00 6.91e+05 -1.7 6.84e+00 - 3.48e-01 1.13e-02h 4
26 3.2162645e+00 8.90e+00 6.63e+05 -1.7 6.85e+00 - 1.40e-01 8.01e-02h 1
27 3.2173798e+00 8.89e+00 6.66e+05 -1.7 7.20e+00 - 8.13e-03 8.57e-04h 1
28 3.2201923e+00 8.83e+00 6.59e+05 -1.7 3.29e+01 - 5.29e-05 6.29e-03f 1
29 3.2628229e+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.2633787e+00 8.09e+00 5.96e+05 -1.7 2.40e+01 - 1.97e-02 7.38e-04h 1
31 3.2630461e+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.1335637e+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.1335144e+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.1715782e+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.1915039e+00 2.67e+00 2.53e+05 -1.7 6.84e+01 - 9.12e-04 2.90e-03h 1
58 3.2055982e+00 2.63e+00 2.44e+05 -1.7 6.96e+01 - 7.32e-03 1.95e-02f 1
59 3.2058327e+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.2175218e+00 2.58e+00 2.32e+05 -1.7 6.77e+01 - 5.28e-04 2.05e-02f 1
61 3.2177382e+00 2.58e+00 2.33e+05 -1.7 7.99e+01 - 1.92e-03 2.46e-04h 1
62 3.2318552e+00 2.54e+00 2.37e+05 -1.7 7.45e+01 - 4.33e-02 1.91e-02f 1
63 3.2320648e+00 2.54e+00 2.37e+05 -1.7 9.03e+01 - 2.93e-03 2.01e-04h 1
64 3.2606225e+00 2.46e+00 2.13e+05 -1.7 7.38e+01 - 1.52e-04 4.35e-02f 1
65 3.2611843e+00 2.46e+00 2.16e+05 -1.7 1.13e+02 - 1.09e-02 4.56e-04h 1
66 3.2892995e+00 2.39e+00 1.92e+05 -1.7 8.25e+01 - 9.64e-05 4.36e-02f 1
67 3.2898502e+00 2.38e+00 1.93e+05 -1.7 1.22e+02 - 3.87e-03 5.55e-04h 1
68 3.3218379e+00 2.32e+00 1.75e+05 -1.7 1.10e+02 - 6.21e-03 3.84e-02f 1
69 3.3938365e+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.4443210e+00 2.04e+00 6.60e+04 -1.7 9.12e+02 - 2.29e-04 1.43e-03f 1
71 3.4449596e+00 2.04e+00 6.00e+04 -1.7 2.86e+02 - 2.56e-02 7.44e-05h 1
72 4.0535875e+00 1.53e+00 7.81e+04 -1.7 1.36e+02 - 2.11e-01 3.87e-01H 1
73 4.2336486e+00 1.41e+00 9.26e+04 -1.7 6.12e+00 -1.2 5.54e-01 7.73e-02h 1
74 4.3632004e+00 1.33e+00 6.95e+04 -1.7 1.09e+02 - 1.69e-03 5.61e-02h 1
75 4.6869628e+00 1.16e+00 3.63e+05 -1.7 2.98e+01 -1.7 1.00e+00 1.24e-01h 1
76 5.3980665e+00 8.72e-01 2.05e+05 -1.7 4.26e+01 - 1.24e-02 2.43e-01h 1
77 5.5593064e+00 8.25e-01 4.02e+05 -1.7 3.16e+00 -1.3 6.18e-01 5.42e-02h 1
78 8.3227421e+00 2.09e-01 1.64e+05 -1.7 3.10e+01 - 3.25e-03 7.91e-01h 1
79 1.0538910e+01 7.14e-02 3.48e+05 -1.7 6.57e+01 - 5.06e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 1.0017940e+01 7.98e-02 4.35e+05 -1.7 1.65e+02 - 3.27e-02 3.22e-01f 1
81 8.0592267e+00 5.00e-02 4.41e+05 -1.7 5.07e-01 -1.7 3.38e-01 1.00e+00f 1
82 4.4250234e+00 6.34e-01 2.85e+05 -1.7 9.38e+01 - 3.55e-01 1.00e+00F 1
83 3.4756774e+00 7.32e-01 1.01e+05 -1.7 5.79e+01 - 6.46e-01 1.00e+00h 1
84 3.3967064e+00 2.64e-02 4.41e-02 -1.7 8.50e-01 -2.2 1.00e+00 1.00e+00h 1
85 3.2411986e+00 5.26e-02 3.63e+05 -3.8 6.40e-01 -2.7 8.17e-01 1.00e+00h 1
86 3.1731684e+00 6.41e-02 1.10e+04 -3.8 8.14e-01 -3.2 9.70e-01 1.00e+00h 1
87 3.1619840e+00 2.60e-02 3.59e-02 -3.8 2.19e+00 -3.6 1.00e+00 1.00e+00h 1
88 3.1586055e+00 1.48e-02 4.52e-02 -3.8 7.20e+00 -4.1 1.00e+00 1.00e+00h 1
89 3.1412468e+00 2.90e-01 6.33e-02 -3.8 3.34e+01 -4.6 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
90 3.1148750e+00 2.63e-01 6.35e-02 -3.8 2.15e+01 -4.2 1.00e+00 1.00e+00h 1
91 3.0854400e+00 1.60e-01 4.10e-02 -3.8 1.30e+01 -3.7 1.00e+00 1.00e+00h 1
92 3.0613106e+00 6.35e-02 1.63e-02 -3.8 6.81e+00 -3.3 1.00e+00 1.00e+00h 1
93 2.6891183e+00 1.98e+01 5.08e+00 -3.8 1.24e+02 -3.8 3.40e-01 1.00e+00h 1
94 2.7264413e+00 1.60e+01 4.10e+00 -3.8 5.44e+01 - 1.00e+00 2.05e-01h 1
95 2.9276478e+00 4.35e+00 2.41e+00 -3.8 5.13e+01 - 5.35e-01 1.00e+00h 1
96 3.0097268e+00 1.30e+00 1.81e+00 -3.8 3.69e+01 - 2.23e-01 8.66e-01h 1
97 2.1035070e+00 5.62e+01 5.27e+01 -3.8 1.21e+02 -3.4 9.00e-02 1.00e+00h 1
98 2.1744837e+00 4.61e+01 4.33e+01 -3.8 6.67e+01 - 2.33e-01 1.88e-01h 1
99 2.3915270e+00 1.88e+01 1.72e+01 -3.8 5.08e+01 - 1.87e-01 7.33e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
100 2.6371424e+00 4.56e+00 9.75e+00 -3.8 4.13e+01 - 1.00e+00 1.00e+00h 1
101 2.5640613e+00 1.70e+00 4.16e+00 -3.8 2.07e+01 - 3.83e-01 6.91e-01h 1
102 2.2684719e+00 3.01e-01 2.32e+01 -3.8 3.82e+00 - 1.52e-01 9.58e-01h 1
103 1.9701083e+00 3.03e-01 1.69e+01 -3.8 1.64e+01 -2.0 3.86e-01 2.79e-01h 1
104 2.0392949e+00 1.44e-01 3.63e+01 -3.8 1.40e+00 - 2.27e-01 1.00e+00h 1
105 1.7430278e+00 1.66e-01 1.43e+02 -3.8 2.11e+00 - 7.32e-02 1.00e+00h 1
106 1.1386790e+00 1.64e-01 6.96e+00 -3.8 1.31e+00 - 3.02e-01 1.00e+00h 1
107 3.4192783e-01 2.80e-01 3.20e-01 -3.8 1.66e+00 - 4.27e-01 1.00e+00h 1
108 1.4057099e-02 2.10e-01 8.27e-02 -3.8 7.62e-01 - 5.33e-01 1.00e+00h 1
109 3.5541464e-05 8.02e-03 3.91e-03 -3.8 2.59e-01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
110 3.3254124e-05 4.80e-05 1.15e-05 -3.8 2.01e-02 - 1.00e+00 1.00e+00h 1
111 1.1999249e-06 7.06e-05 7.49e-05 -5.7 3.08e-02 - 1.00e+00 1.00e+00h 1
112 1.1912245e-06 1.32e-08 3.12e-08 -5.7 4.68e-04 - 1.00e+00 1.00e+00h 1
113 1.1862889e-06 1.14e-08 1.25e-08 -8.6 3.94e-04 - 1.00e+00 1.00e+00h 1
114 1.1862894e-06 1.95e-14 1.59e-14 -9.0 4.16e-07 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 114
(scaled) (unscaled)
Objective...............: 1.1862894024318393e-06 1.1862894024318393e-06
Dual infeasibility......: 1.5938170857656184e-14 1.5938170857656184e-14
Constraint violation....: 1.9539925233402755e-14 1.9539925233402755e-14
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 9.0909169532453792e-10 9.0909169532453792e-10
Overall NLP error.......: 9.0909169532453792e-10 9.0909169532453792e-10
Number of objective function evaluations = 152
Number of objective gradient evaluations = 115
Number of equality constraint evaluations = 152
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 115
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 114
Total seconds in IPOPT = 0.253
EXIT: Optimal Solution Found.
k1= 5.003509348488081
k2= 0.9999977260457728