Parameter Estimation 2#
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.13.2:
******************************************************************************
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 http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
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.5828123e+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.7587273e+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.9029431e+00 1.15e+01 2.76e+05 -1.7 1.45e+02 -1.6 5.69e-02 2.07e-02h 1
19 2.9064258e+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.9777374e+00 1.09e+01 2.62e+05 -1.7 9.65e+00 - 1.39e-01 5.63e-02h 1
21 2.9781355e+00 1.08e+01 2.63e+05 -1.7 8.00e+00 - 2.30e-03 3.28e-04h 2
22 3.0173366e+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.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.2162648e+00 8.90e+00 6.63e+05 -1.7 6.85e+00 - 1.40e-01 8.01e-02h 1
27 3.2173800e+00 8.89e+00 6.66e+05 -1.7 7.20e+00 - 8.13e-03 8.57e-04h 1
28 3.2201926e+00 8.83e+00 6.59e+05 -1.7 3.29e+01 - 5.29e-05 6.29e-03f 1
29 3.2628230e+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.2633788e+00 8.09e+00 5.96e+05 -1.7 2.40e+01 - 1.97e-02 7.38e-04h 1
31 3.2630462e+00 8.08e+00 5.95e+05 -1.7 4.97e+01 - 1.43e-05 1.45e-03f 1
32 3.2712178e+00 7.58e+00 5.52e+05 -1.7 4.12e+01 - 4.76e-02 6.18e-02f 1
33 3.2713520e+00 7.57e+00 5.54e+05 -1.7 3.50e+01 - 7.06e-03 4.92e-04h 1
34 3.1338655e+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.1335175e+00 3.89e+00 4.55e+05 -1.7 6.01e+01 - 4.95e-04 3.26e-04h 1
37 3.1335145e+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.1171847e+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.1151185e+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.1220772e+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.1401919e+00 2.89e+00 2.93e+05 -1.7 5.09e+01 - 7.51e-04 3.17e-04h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 3.1487223e+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.1715785e+00 2.74e+00 2.60e+05 -1.7 6.38e+01 - 1.13e-03 1.29e-03h 1
55 3.1802298e+00 2.71e+00 2.53e+05 -1.7 6.32e+01 - 4.01e-03 1.37e-02f 1
56 3.1895629e+00 2.68e+00 2.54e+05 -1.7 6.82e+01 - 2.38e-02 1.21e-02f 1
57 3.1915021e+00 2.67e+00 2.53e+05 -1.7 6.84e+01 - 9.11e-04 2.90e-03h 1
58 3.2055975e+00 2.63e+00 2.44e+05 -1.7 6.96e+01 - 7.31e-03 1.95e-02f 1
59 3.2058320e+00 2.63e+00 2.44e+05 -1.7 8.39e+01 - 1.44e-03 2.13e-04h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 3.2175211e+00 2.58e+00 2.32e+05 -1.7 6.77e+01 - 5.29e-04 2.05e-02f 1
61 3.2177375e+00 2.58e+00 2.33e+05 -1.7 7.99e+01 - 1.92e-03 2.46e-04h 1
62 3.2318545e+00 2.54e+00 2.37e+05 -1.7 7.45e+01 - 4.33e-02 1.91e-02f 1
63 3.2320641e+00 2.54e+00 2.37e+05 -1.7 9.03e+01 - 2.93e-03 2.01e-04h 1
64 3.2606218e+00 2.46e+00 2.13e+05 -1.7 7.38e+01 - 1.52e-04 4.35e-02f 1
65 3.2611836e+00 2.46e+00 2.16e+05 -1.7 1.13e+02 - 1.09e-02 4.56e-04h 1
66 3.2892990e+00 2.39e+00 1.92e+05 -1.7 8.25e+01 - 9.64e-05 4.37e-02f 1
67 3.2898497e+00 2.38e+00 1.93e+05 -1.7 1.22e+02 - 3.87e-03 5.55e-04h 1
68 3.3218374e+00 2.32e+00 1.75e+05 -1.7 1.10e+02 - 6.21e-03 3.84e-02f 1
69 3.3937825e+00 2.04e+00 6.58e+04 -1.7 7.37e+01 - 4.33e-04 4.11e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 3.4444856e+00 2.04e+00 6.60e+04 -1.7 8.98e+02 - 2.32e-04 1.46e-03f 1
71 3.4451251e+00 2.04e+00 6.02e+04 -1.7 2.89e+02 - 2.45e-02 7.28e-05h 1
72 4.0543379e+00 1.52e+00 7.87e+04 -1.7 1.36e+02 - 2.09e-01 3.87e-01H 1
73 4.2335621e+00 1.41e+00 9.07e+04 -1.7 6.09e+00 -1.2 5.50e-01 7.69e-02h 1
74 4.3614927e+00 1.33e+00 6.80e+04 -1.7 1.10e+02 - 1.67e-03 5.54e-02h 1
75 4.6828370e+00 1.16e+00 3.62e+05 -1.7 2.97e+01 -1.7 1.00e+00 1.24e-01h 1
76 5.3918616e+00 8.74e-01 2.05e+05 -1.7 4.27e+01 - 1.24e-02 2.43e-01h 1
77 5.5522789e+00 8.27e-01 4.03e+05 -1.7 3.18e+00 -1.3 6.23e-01 5.39e-02h 1
78 8.3533766e+00 2.13e-01 1.69e+05 -1.7 2.91e+01 - 3.77e-03 7.99e-01h 1
79 1.0541745e+01 7.07e-02 3.89e+05 -1.7 6.55e+01 - 4.53e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 1.0046735e+01 7.91e-02 4.40e+05 -1.7 2.78e+02 - 1.92e-02 1.94e-01f 1
81 7.8972966e+00 6.09e-02 5.36e+05 -1.7 7.10e-01 -1.7 2.43e-01 1.00e+00f 1
82 4.0836282e+00 9.72e-01 3.90e+05 -1.7 1.15e+02 - 2.73e-01 1.00e+00F 1
83 3.5507361e+00 9.75e-02 5.61e-01 -1.7 7.00e-01 -2.2 1.00e+00 1.00e+00h 1
84 3.2930543e+00 9.28e-02 3.53e+05 -2.5 8.16e-01 -2.7 7.94e-01 1.00e+00h 1
85 3.2146757e+00 4.93e-02 6.95e-03 -2.5 1.02e+00 -3.2 1.00e+00 1.00e+00h 1
86 3.2013289e+00 5.60e-03 1.74e-02 -2.5 3.05e+00 -3.6 1.00e+00 1.00e+00h 1
87 3.1802213e+00 2.52e-02 2.14e+04 -3.8 3.64e+00 -4.1 9.20e-01 1.00e+00h 1
88 3.1754618e+00 2.37e-02 1.42e-01 -3.8 1.18e+01 -4.6 1.00e+00 1.00e+00h 1
89 3.1588649e+00 5.23e-01 1.00e-01 -3.8 6.04e+01 -5.1 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.1154685e+00 1.23e+00 2.80e-01 -3.8 5.74e+01 -4.6 1.00e+00 1.00e+00h 1
91 3.0483853e+00 1.84e-01 4.87e-02 -3.8 7.01e+00 -3.3 1.00e+00 1.00e+00h 1
92 3.0219307e+00 3.81e-02 9.21e-03 -3.8 3.88e+00 -2.9 1.00e+00 1.00e+00h 1
93 2.8828318e+00 1.48e+00 3.67e-01 -3.8 2.80e+01 -3.4 1.00e+00 1.00e+00h 1
94 2.4687798e+00 5.58e-01 1.35e-01 -3.8 4.98e+00 -2.0 1.00e+00 1.00e+00h 1
95 1.7126984e+00 1.36e+00 3.10e-01 -3.8 7.42e+00 -1.6 3.62e-01 1.00e+00h 1
96 6.1388536e-01 6.05e+00 1.14e+00 -3.8 1.34e+01 -1.2 3.01e-02 5.19e-01f 1
97 6.1413847e-01 6.05e+00 1.53e+00 -3.8 1.95e+01 - 2.30e-01 6.08e-05h 1
98 5.1440213e-01 4.01e+00 1.17e+00 -3.8 6.04e+00 - 4.69e-02 3.37e-01h 1
99 5.2496160e-01 3.65e+00 1.79e+01 -3.8 3.85e+00 - 5.56e-03 8.81e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
100 5.2516703e-01 3.65e+00 2.70e+01 -3.8 6.32e+00 - 1.94e-02 3.91e-04h 1
101 1.7963225e-01 7.77e-01 1.17e+01 -3.8 3.17e+00 - 4.96e-03 1.00e+00f 1
102 2.0189004e-01 1.47e-01 1.66e+02 -3.8 1.58e+00 - 1.45e-01 1.00e+00h 1
103 1.1250787e-01 6.17e-02 1.32e+01 -3.8 1.13e+00 - 1.69e-01 1.00e+00h 1
104 2.4033538e-02 1.52e-01 1.45e+00 -3.8 1.72e+00 - 2.74e-01 1.00e+00h 1
105 2.8637398e-03 1.76e-02 1.19e-01 -3.8 4.12e-01 - 4.96e-01 1.00e+00h 1
106 2.0192356e-05 7.62e-03 3.57e-03 -3.8 3.18e-01 - 1.00e+00 1.00e+00h 1
107 3.3274318e-05 4.77e-06 2.27e-06 -3.8 5.72e-03 - 1.00e+00 1.00e+00h 1
108 1.1999865e-06 7.04e-05 7.50e-05 -5.7 3.08e-02 - 1.00e+00 1.00e+00h 1
109 1.1912245e-06 1.32e-08 3.13e-08 -5.7 4.68e-04 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
110 1.1862889e-06 1.14e-08 1.25e-08 -8.6 3.94e-04 - 1.00e+00 1.00e+00h 1
111 1.1862894e-06 1.27e-14 1.59e-14 -9.0 4.16e-07 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 111
(scaled) (unscaled)
Objective...............: 1.1862894024323384e-06 1.1862894024323384e-06
Dual infeasibility......: 1.5938203946741843e-14 1.5938203946741843e-14
Constraint violation....: 1.2656542480726785e-14 1.2656542480726785e-14
Complementarity.........: 9.0909169532399198e-10 9.0909169532399198e-10
Overall NLP error.......: 9.0909169532399198e-10 9.0909169532399198e-10
Number of objective function evaluations = 150
Number of objective gradient evaluations = 112
Number of equality constraint evaluations = 150
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 112
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 111
Total CPU secs in IPOPT (w/o function evaluations) = 0.097
Total CPU secs in NLP function evaluations = 0.007
EXIT: Optimal Solution Found.
k1= 5.003509348488082
k2= 0.9999977260457736
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.13.2:
******************************************************************************
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 http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
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.5828123e+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.7587273e+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.9029431e+00 1.15e+01 2.76e+05 -1.7 1.45e+02 -1.6 5.69e-02 2.07e-02h 1
19 2.9064258e+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.9777374e+00 1.09e+01 2.62e+05 -1.7 9.65e+00 - 1.39e-01 5.63e-02h 1
21 2.9781355e+00 1.08e+01 2.63e+05 -1.7 8.00e+00 - 2.30e-03 3.28e-04h 2
22 3.0173366e+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.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.2162648e+00 8.90e+00 6.63e+05 -1.7 6.85e+00 - 1.40e-01 8.01e-02h 1
27 3.2173800e+00 8.89e+00 6.66e+05 -1.7 7.20e+00 - 8.13e-03 8.57e-04h 1
28 3.2201926e+00 8.83e+00 6.59e+05 -1.7 3.29e+01 - 5.29e-05 6.29e-03f 1
29 3.2628230e+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.2633788e+00 8.09e+00 5.96e+05 -1.7 2.40e+01 - 1.97e-02 7.38e-04h 1
31 3.2630462e+00 8.08e+00 5.95e+05 -1.7 4.97e+01 - 1.43e-05 1.45e-03f 1
32 3.2712178e+00 7.58e+00 5.52e+05 -1.7 4.12e+01 - 4.76e-02 6.18e-02f 1
33 3.2713520e+00 7.57e+00 5.54e+05 -1.7 3.50e+01 - 7.06e-03 4.92e-04h 1
34 3.1338655e+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.1335175e+00 3.89e+00 4.55e+05 -1.7 6.01e+01 - 4.95e-04 3.26e-04h 1
37 3.1335145e+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.1171847e+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.1151185e+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.1220772e+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.1401919e+00 2.89e+00 2.93e+05 -1.7 5.09e+01 - 7.51e-04 3.17e-04h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 3.1487223e+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.1715785e+00 2.74e+00 2.60e+05 -1.7 6.38e+01 - 1.13e-03 1.29e-03h 1
55 3.1802298e+00 2.71e+00 2.53e+05 -1.7 6.32e+01 - 4.01e-03 1.37e-02f 1
56 3.1895629e+00 2.68e+00 2.54e+05 -1.7 6.82e+01 - 2.38e-02 1.21e-02f 1
57 3.1915021e+00 2.67e+00 2.53e+05 -1.7 6.84e+01 - 9.11e-04 2.90e-03h 1
58 3.2055975e+00 2.63e+00 2.44e+05 -1.7 6.96e+01 - 7.31e-03 1.95e-02f 1
59 3.2058320e+00 2.63e+00 2.44e+05 -1.7 8.39e+01 - 1.44e-03 2.13e-04h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 3.2175211e+00 2.58e+00 2.32e+05 -1.7 6.77e+01 - 5.29e-04 2.05e-02f 1
61 3.2177375e+00 2.58e+00 2.33e+05 -1.7 7.99e+01 - 1.92e-03 2.46e-04h 1
62 3.2318545e+00 2.54e+00 2.37e+05 -1.7 7.45e+01 - 4.33e-02 1.91e-02f 1
63 3.2320641e+00 2.54e+00 2.37e+05 -1.7 9.03e+01 - 2.93e-03 2.01e-04h 1
64 3.2606218e+00 2.46e+00 2.13e+05 -1.7 7.38e+01 - 1.52e-04 4.35e-02f 1
65 3.2611836e+00 2.46e+00 2.16e+05 -1.7 1.13e+02 - 1.09e-02 4.56e-04h 1
66 3.2892990e+00 2.39e+00 1.92e+05 -1.7 8.25e+01 - 9.64e-05 4.37e-02f 1
67 3.2898497e+00 2.38e+00 1.93e+05 -1.7 1.22e+02 - 3.87e-03 5.55e-04h 1
68 3.3218374e+00 2.32e+00 1.75e+05 -1.7 1.10e+02 - 6.21e-03 3.84e-02f 1
69 3.3937825e+00 2.04e+00 6.58e+04 -1.7 7.37e+01 - 4.33e-04 4.11e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 3.4444856e+00 2.04e+00 6.60e+04 -1.7 8.98e+02 - 2.32e-04 1.46e-03f 1
71 3.4451251e+00 2.04e+00 6.02e+04 -1.7 2.89e+02 - 2.45e-02 7.28e-05h 1
72 4.0543379e+00 1.52e+00 7.87e+04 -1.7 1.36e+02 - 2.09e-01 3.87e-01H 1
73 4.2335621e+00 1.41e+00 9.07e+04 -1.7 6.09e+00 -1.2 5.50e-01 7.69e-02h 1
74 4.3614927e+00 1.33e+00 6.80e+04 -1.7 1.10e+02 - 1.67e-03 5.54e-02h 1
75 4.6828370e+00 1.16e+00 3.62e+05 -1.7 2.97e+01 -1.7 1.00e+00 1.24e-01h 1
76 5.3918616e+00 8.74e-01 2.05e+05 -1.7 4.27e+01 - 1.24e-02 2.43e-01h 1
77 5.5522789e+00 8.27e-01 4.03e+05 -1.7 3.18e+00 -1.3 6.23e-01 5.39e-02h 1
78 8.3533766e+00 2.13e-01 1.69e+05 -1.7 2.91e+01 - 3.77e-03 7.99e-01h 1
79 1.0541745e+01 7.07e-02 3.89e+05 -1.7 6.55e+01 - 4.53e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 1.0046735e+01 7.91e-02 4.40e+05 -1.7 2.78e+02 - 1.92e-02 1.94e-01f 1
81 7.8972966e+00 6.09e-02 5.36e+05 -1.7 7.10e-01 -1.7 2.43e-01 1.00e+00f 1
82 4.0836282e+00 9.72e-01 3.90e+05 -1.7 1.15e+02 - 2.73e-01 1.00e+00F 1
83 3.5507361e+00 9.75e-02 5.61e-01 -1.7 7.00e-01 -2.2 1.00e+00 1.00e+00h 1
84 3.2930543e+00 9.28e-02 3.53e+05 -2.5 8.16e-01 -2.7 7.94e-01 1.00e+00h 1
85 3.2146757e+00 4.93e-02 6.95e-03 -2.5 1.02e+00 -3.2 1.00e+00 1.00e+00h 1
86 3.2013289e+00 5.60e-03 1.74e-02 -2.5 3.05e+00 -3.6 1.00e+00 1.00e+00h 1
87 3.1802213e+00 2.52e-02 2.14e+04 -3.8 3.64e+00 -4.1 9.20e-01 1.00e+00h 1
88 3.1754618e+00 2.37e-02 1.42e-01 -3.8 1.18e+01 -4.6 1.00e+00 1.00e+00h 1
89 3.1588649e+00 5.23e-01 1.00e-01 -3.8 6.04e+01 -5.1 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.1154685e+00 1.23e+00 2.80e-01 -3.8 5.74e+01 -4.6 1.00e+00 1.00e+00h 1
91 3.0483853e+00 1.84e-01 4.87e-02 -3.8 7.01e+00 -3.3 1.00e+00 1.00e+00h 1
92 3.0219307e+00 3.81e-02 9.21e-03 -3.8 3.88e+00 -2.9 1.00e+00 1.00e+00h 1
93 2.8828318e+00 1.48e+00 3.67e-01 -3.8 2.80e+01 -3.4 1.00e+00 1.00e+00h 1
94 2.4687798e+00 5.58e-01 1.35e-01 -3.8 4.98e+00 -2.0 1.00e+00 1.00e+00h 1
95 1.7126984e+00 1.36e+00 3.10e-01 -3.8 7.42e+00 -1.6 3.62e-01 1.00e+00h 1
96 6.1388536e-01 6.05e+00 1.14e+00 -3.8 1.34e+01 -1.2 3.01e-02 5.19e-01f 1
97 6.1413847e-01 6.05e+00 1.53e+00 -3.8 1.95e+01 - 2.30e-01 6.08e-05h 1
98 5.1440213e-01 4.01e+00 1.17e+00 -3.8 6.04e+00 - 4.69e-02 3.37e-01h 1
99 5.2496160e-01 3.65e+00 1.79e+01 -3.8 3.85e+00 - 5.56e-03 8.81e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
100 5.2516703e-01 3.65e+00 2.70e+01 -3.8 6.32e+00 - 1.94e-02 3.91e-04h 1
101 1.7963225e-01 7.77e-01 1.17e+01 -3.8 3.17e+00 - 4.96e-03 1.00e+00f 1
102 2.0189004e-01 1.47e-01 1.66e+02 -3.8 1.58e+00 - 1.45e-01 1.00e+00h 1
103 1.1250787e-01 6.17e-02 1.32e+01 -3.8 1.13e+00 - 1.69e-01 1.00e+00h 1
104 2.4033538e-02 1.52e-01 1.45e+00 -3.8 1.72e+00 - 2.74e-01 1.00e+00h 1
105 2.8637398e-03 1.76e-02 1.19e-01 -3.8 4.12e-01 - 4.96e-01 1.00e+00h 1
106 2.0192356e-05 7.62e-03 3.57e-03 -3.8 3.18e-01 - 1.00e+00 1.00e+00h 1
107 3.3274318e-05 4.77e-06 2.27e-06 -3.8 5.72e-03 - 1.00e+00 1.00e+00h 1
108 1.1999865e-06 7.04e-05 7.50e-05 -5.7 3.08e-02 - 1.00e+00 1.00e+00h 1
109 1.1912245e-06 1.32e-08 3.13e-08 -5.7 4.68e-04 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
110 1.1862889e-06 1.14e-08 1.25e-08 -8.6 3.94e-04 - 1.00e+00 1.00e+00h 1
111 1.1862894e-06 1.27e-14 1.59e-14 -9.0 4.16e-07 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 111
(scaled) (unscaled)
Objective...............: 1.1862894024323384e-06 1.1862894024323384e-06
Dual infeasibility......: 1.5938203946741843e-14 1.5938203946741843e-14
Constraint violation....: 1.2656542480726785e-14 1.2656542480726785e-14
Complementarity.........: 9.0909169532399198e-10 9.0909169532399198e-10
Overall NLP error.......: 9.0909169532399198e-10 9.0909169532399198e-10
Number of objective function evaluations = 150
Number of objective gradient evaluations = 112
Number of equality constraint evaluations = 150
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 112
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 111
Total CPU secs in IPOPT (w/o function evaluations) = 0.091
Total CPU secs in NLP function evaluations = 0.006
EXIT: Optimal Solution Found.
k1= 5.003509348488082
k2= 0.9999977260457736