pyodefit - A state and parameter estimation algorithm

pyodefit is a implementation of the variable and parameter estimation algorithm described in Ref. [1]. With this optimization based approach, a cost function is minimized to estimate the parameters and model variables of a mathematical model \( \mathrm{d}\mathbf{y}(t) / \mathrm{d}t = \mathbf{F}(t,\mathbf{y}(t),\mathbf{p}) \), given as a set of ordinary differential equations, from given observed or measured time series. The goal is to find a solution for the model variables and parameters, so that on the one hand the model variables describe the data as good as possible and on the other hand the model equations are fulfilled as exactly as possible for the estimated variables and parameters. The used cost function (see [1]) has the form \[ C(\mathbf{x}) = \sum_i H_i(\mathbf{x})^2 \; , \] where \(\mathbf{x}\) is the input vector and \( \mathbf{H}(\mathbf{x}) \) is a vector valued function. The cost function penalizes a deviation between the time series and a user defined measurement function (which describes the data) on the one hand and the error in approximating the model equations \( \mathbf{u}(t) = \mathrm{\Delta}\mathbf{y}(t) / \mathrm{\Delta}t - \mathbf{F}(t,\mathbf{y}(t),\mathbf{p}) \) on the other hand (\(\mathrm{\Delta}\mathbf{y}(t) / \mathrm{\Delta}t \) is the finite difference approximation of the time derivative). Two more terms in the cost function are used to obtain a smooth solution for the model variables (employing Hermite interpolation) and to restrict variables and parameters to predefined bounds. Only the model equations (ODEs), a measurement function and the data have to be providedto the estimation algorithm.

To minimize \( C(\mathbf{x}) \) sparseLM [2] (a sparse implementation of the Levenberg-Marquardt algorithm) is used here. To call sparseLM from Python the Python wrapper pysparseLM is needed. Furthermore, sparseLM requires the sparse Jacobian matrix of \( \mathbf{H}(\mathbf{x}) \) with respect to  \( \mathbf{x} \) and its sparsity pattern (indices of non-zero elements). To compute the sparsity pattern and the values of the corresponding elements of the Jacobian matrix the automatic differentiation tool ADOL-C is used. The ADOL-C library is called from Python by the wrapper pyadolc.

References

  1. J. Schumann-Bischoff and U. Parlitz, "State and parameter estimation using unconstrained optimization", Phys. Rev. E 84, 056214 (2011)
  2. M. I. A. Lourakis, "Sparse Non-linear Least Squares Optimization for Geometric Vision", European Conference on Computer Vision, 2, pp. 43-56 (2010)