4. Parameter estimation

The module casiopeia.pe contains the classes for parameter estimation. For now, only least squares parameter estimation problems are covered.

4.1. Parameter estimation from single experiments

class casiopeia.pe.LSq(system, time_points, udata=None, qdata=None, ydata=None, pinit=None, xinit=None, wv=None, weps_u=None, discretization_method='collocation', **kwargs)

The class casiopeia.pe.LSq is used to set up least squares parameter estimation problems for systems defined with the casiopeia.system.System class, using a given set of user-provided control data, measurement data and different kinds of weightings.

Raises:

AttributeError, NotImplementedError

Parameters:

• system (casiopeia.system.System) – system considered for parameter estimation, specified using the casiopeia.system.System class
• time_points (numpy.ndarray, casadi.DMatrix, list) – time points \(t_\text{N} \in \mathbb{R}^\text{N}\) used to discretize the continuous time problem. Controls will be applied at the first \(N-1\) time points, while measurements take place at all \(N\) time points.
• udata (numpy.ndarray, casadi.DMatrix) – optional, values for the time-varying controls \(u_\text{N} \in \mathbb{R}^{\text{n}_\text{u} \times \text{N}-1}\) that can change at the switching time points; if no values are given, 0 will be used; note that the the second dimension of \(u_\text{N}\) is \(N-1\) and not \(N\), since there is no control value applied at the last time point
• qdata (numpy.ndarray, casadi.DMatrix) – optional, values for the time-constant controls \(q_\text{N} \in \mathbb{R}^{\text{n}_\text{q}}\); if not values are given, 0 will be used
• ydata (numpy.ndarray, casadi.DMatrix) – values for the measurements at the switching time points \(y_\text{N} \in \mathbb{R}^{\text{n}_\text{y} \times \text{N}}\)
• wv (numpy.ndarray, casadi.DMatrix) – weightings for the measurements \(w_\text{v} \in \mathbb{R}^{\text{n}_\text{y} \times \text{N}}\)
• weps_u (numpy.ndarray, casadi.DMatrix) – weightings for the input errors \(w_{\epsilon_\text{u}} \in \mathbb{R}^{\text{n}_{\epsilon_\text{u}}}\) (only necessary if input errors are used within system)
• pinit (numpy.ndarray, casadi.DMatrix) – optional, initial guess for the values of the parameters that will be estimated \(p_\text{init} \in \mathbb{R}^{\text{n}_\text{p}}\); if no value is given, 0 will be used; note that a poorly or wrongly chosen initial guess can cause the estimation to fail
• xinit (numpy.ndarray, casadi.DMatrix) – optional, initial guess for the values of the states that will be estimated \(x_\text{init} \in \mathbb{R}^{\text{n}_\text{x} \times \text{N}}\); if no value is given, 0 will be used; note that a poorly or wrongly chosen initial guess can cause the estimation to fail
• discretization_method (str) – optional, the method that shall be used for discretization of the continuous time problem w. r. t. the time points given in \(t_\text{N}\); possible values are “collocation” (default) and “multiple_shooting”

Depending on the discretization method specified in discretization_method, the following parameters can be used for further specification:

Parameters:

• collocation_scheme (str) – optional, scheme used for setting up the collocation polynomials, possible values are radau (default) and legendre
• number_of_collocation_points (int) – optional, order of collocation polynomials \(d \in \mathbb{Z}\) (default values is 3)
• integrator (str) – optional, integrator to be used with multiple shooting. See the CasADi documentation for a list of all available integrators. As a default, cvodes is used.
• integrator_options (dict) – optional, options to be passed to the CasADi integrator used with multiple shooting (see the CasADi documentation for a list of all possible options)

The resulting parameter estimation problem has the following form:

\[\begin{split}\begin{aligned} \text{arg}\,\underset{p, x, v, \epsilon_\text{u}}{\text{min}} & & \frac{1}{2} \| R(\cdot) \|_2^2 &\\ \text{subject to:} & & v_\text{k} + y_\text{k} - \phi(x_\text{k}, p; u_\text{k}, q) & = 0 \hspace{1cm} k = 1, \dots, N\\ & & g(x, p, \epsilon_\text{u}; u, q) & = 0 \\ \text{with:} & & \begin{pmatrix} {w_\text{v}}^T & {w_{\epsilon_\text{u}}}^T \end{pmatrix}^{^\mathbb{1}/_\mathbb{2}} \begin{pmatrix} {v} \\ {\epsilon_\text{u}} \end{pmatrix} & = R \\ \end{aligned}\end{split}\]

while \(g(\cdot)\) contains the discretized system dynamics according to the specified discretization method. If the system is non-dynamic, it only contains the user-provided equality constraints.

compute_covariance_matrix()

This function computes the covariance matrix for the estimated parameters from the inverse of the KKT matrix for the parameter estimation problem, which allows for statements on the quality of the values of the estimated parameters [1] [2].

For efficiency, only the inverse of the relevant part of the matrix is computed [3].

The values of the covariance matrix \(\Sigma_{\hat{\text{p}}}\) can afterwards be accessed via the class attribute LSq.covariance_matrix, and the contained standard deviations \(\sigma_{\hat{\text{p}}}\) for the estimated parameters directly via LSq.standard_deviations.

References

[1]	Kostina, Ekaterina and Kostyukova, Olga: Computing Covariance Matrices for Constrained Nonlinear Large Scale Parameter Estimation Problems Using Krylov Subspace Methods, 2012.

[2]	Kostina, Ekaterina and Kriwet, Gregor: Towards Optimum Experimental Design for Partial Differential Equations, SPP 1253 annual conference 2010, slides 12/13.

[3]	Walter, Eric and Prozanto, Luc: Identification of Parametric Models from Experimental Data, Springer, 1997, pages 288/289.

print_estimation_results()

Raises:AttributeError

This function displays the results of the parameter estimation computations. It can not be used before function run_parameter_estimation() has been used. The results displayed by the function contain:

• the values of the estimated parameters \(\hat{p}\) and their corresponding standard deviations \(\sigma_{\hat{\text{p}}},\) (the values of the standard deviations are presented only if the covariance matrix had already been computed),
• the values of the covariance matrix \(\Sigma_{\hat{\text{p}}}\) for the estimated parameters (if it had already been computed), and
• the durations of the estimation and (if already executed) of the covariance matrix computation.

run_parameter_estimation(solver_options={})

Parameters:solver_options (dict) – options to be passed to the IPOPT solver (see the CasADi documentation for a list of all possible options)

This functions will pass the least squares parameter estimation problem to the IPOPT solver. The status of IPOPT printed to the console provides information whether the optimization finished successfully. The estimated parameters \(\hat{p}\) can afterwards be accessed via the class attribute LSq.estimated_parameters.

Note

IPOPT finishing successfully does not necessarily mean that the estimation results for the unknown parameters are useful for your purposes, it just means that IPOPT was able to solve the given optimization problem. You have in any case to verify your results, e. g. by simulation using the casiopeia casiopeia.sim.Simulation class!

4.2. Parameter estimation from multiple experiments

class casiopeia.pe.MultiLSq(pe_setups=[])

The class casiopeia.pe.MultiLSq is used to construct a single least squares parameter estimation problem from multiple least squares parameter estimation problems defined via two or more objects of type casiopeia.pe.LSq.

In this way, the results of multiple independent experimental setups can be used for parameter estimation.

Note

It is assumed that the system description used for setting up the several parameter estimation problems is the same.

Parameters:pe_setups (list) – list of two or more objects of type casiopeia.pe.LSq