Installation on Ubuntu 16.04¶

The following instructions show the installation on Ubuntu 16.04. If you are planning to install casiopeia on Linux systems different from Ubuntu 16.04, these commands need to be adapted accordingly.

sudo apt-get update
sudo apt-get install python python-numpy python-scipy python-matplotlib python-setuptools --install-recommends

sudo apt-get install python-pip

sudo apt-get install spyder

export PYTHONPATH=$PYTHONPATH:/<path>/<to>/<casadi>/<folder>

echo $PYTHONPATH

>>> from casadi import *
>>> x = MX.sym("x")
>>> print jacobian(sin(x),x)

sudo pip install casiopeia

sudo env PYTHONPATH=$PYTHONPATH pip install casiopeia

sudo pip install casiopeia --upgrade

sudo env PYTHONPATH=$PYTHONPATH pip install casiopeia --upgrade

sudo apt-get update
sudo apt-get install git

git clone git@github.com:adbuerger/casiopeia.git

sudo python setup.py install

sudo env PYTHONPATH=$PYTHONPATH python setup.py install

git pull

sudo python setup.py install

sudo env PYTHONPATH=$PYTHONPATH python setup.py install

Installation on Windows¶

The following instructions have been tested on Windows 7 64 bit.

pip install casiopeia

pip install casiopeia --upgrade

python setup.py install

python setup.py install

Recommendations¶

To speed up computations in casiopeia, it is recommended to install HSL for IPOPT. On how to install the solvers and for further information, see the page Obtaining HSL in the CasADi wiki.

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)[source]¶

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!

Parameter estimation from multiple experiments¶

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

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

Optimum experimental design of single experiments¶

class casiopeia.doe.DoE(system, time_points, uinit=None, umin=None, umax=None, qinit=None, qmin=None, qmax=None, pdata=None, x0=None, xmin=None, xmax=None, wv=None, weps_u=None, discretization_method='collocation', optimality_criterion='A', **kwargs)[source]¶

The class casiopeia.doe.DoE is used to set up Design-of-Experiments-problems for systems defined with the casiopeia.system.System class.

The aim of the experimental design optimization is to identify a set of controls that can be used for the generation of measurement data which allows for a better estimation of the unknown parameters of a system.

To achieve this, an information function on the covariance matrix of the estimated parameters is minimized. The values of the estimated parameters, though they are mostly an initial guess for their values, are not changed during the optimization.

Optimum experimental design and parameter estimation methods can be used interchangeably until a desired accuracy of the parameters has been achieved.

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. umin (numpy.ndarray, casadi.DMatrix) – optional, lower bounds of the time-varying controls \(u_\text{min} \in \mathbb{R}^{\text{n}_\text{u}}\); if not values are given, \(-\infty\) will be used umax (numpy.ndarray, casadi.DMatrix) – optional, upper bounds of the time-vaying controls \(u_\text{max} \in \mathbb{R}^{\text{n}_\text{u}}\); if not values are given, \(\infty\) will be used uinit (numpy.ndarray, casadi.DMatrix) – optional, initial guess for the values of the time-varying controls \(u_\text{N} \in \mathbb{R}^{\text{n}_\text{u} \times \text{N}-1}\) that (might) change at the switching time points; if no values are given, 0 will be used; note that a poorly or wrongly chosen initial guess can cause the optimization to fail, and note that the the second dimension of \(u_N\) is \(N-1\) and not \(N\), since there is no control value applied at the last time point qmin (numpy.ndarray, casadi.DMatrix) – optional, lower bounds of the time-constant controls \(q_\text{min} \in \mathbb{R}^{\text{n}_\text{q}}\); if not values are given, \(-\infty\) will be used qmax (numpy.ndarray, casadi.DMatrix) – optional, upper bounds of the time-constant controls \(q_\text{max} \in \mathbb{R}^{\text{n}_\text{q}}\); if not values are given, \(\infty\) will be used qinit (numpy.ndarray, casadi.DMatrix) – optional, initial guess for the optimal values of the time-constant controls \(q_\text{init} \in \mathbb{R}^{\text{n}_\text{q}}\); if not values are given, 0 will be used; note that a poorly or wrongly chosen initial guess can cause the optimization to fail pdata (numpy.ndarray, casadi.DMatrix) – values of the time-constant parameters \(p \in \mathbb{R}^{\text{n}_\text{p}}\) x0 (numpy.ndarray, casadi.DMatrix, list) – state values \(x_0 \in \mathbb{R}^{\text{n}_\text{x}}\) at the first time point \(t_0\) xmin (numpy.ndarray, casadi.DMatrix) – optional, lower bounds of the states \(x_\text{min} \in \mathbb{R}^{\text{n}_\text{x}}\); if no value is given, \(-\infty\) will be used xmax (numpy.ndarray, casadi.DMatrix) – optional, lower bounds of the states \(x_\text{max} \in \mathbb{R}^{\text{n}_\text{x}}\); if no value is given, \(\infty\) will be used 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) 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” optimality_criterion (str) – optional, the information function \(I_\text{X}(\cdot)\) to be used on the covariance matrix, possible values are A (default) and D, while \[\begin{split}\begin{aligned} I_\text{A}(\Sigma_\text{p}) & = \frac{1}{n_\text{p}} \text{Tr}(\Sigma_\text{p}),\\ I_\text{D}(\Sigma_\text{p}) & = \begin{vmatrix} \Sigma_\text{p} \end{vmatrix} ^{\frac{1}{n_\text{p}}}, \end{aligned}\end{split}\] for further information see e. g. [1]

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)

You do not need to specify initial guesses for the estimated states, since these are obtained with a system simulation using the initial states and the provided initial guesses for the controls.

The resulting optimization problem has the following form:

\[\begin{split}\begin{aligned} \text{arg}\,\underset{u, q, x}{\text{min}} & & I(\Sigma_{\text{p}}(x, u, q; p)) &\\ \text{subject to:} & & g(x, u, q; p) & = 0\\ & & u_\text{min} \leq u_\text{k} & \leq u_\text{max} \hspace{1cm} k = 1, \dots, N-1\\ & & x_\text{min} \leq x_\text{k} & \leq x_\text{max} \hspace{1cm} k = 1, \dots, N\\ & & x_1 \leq x(t_1) & \leq x_1 \end{aligned}\end{split}\]

where \(\Sigma_p = \text{Cov}(p)\) and \(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.

References

[1]	Körkel, Stefan: Numerische Methoden für Optimale Versuchsplanungsprobleme bei nichtlinearen DAE-Modellen, PhD Thesis, Heidelberg university, 2002, pages 74/75.

plot_confidence_ellipsoids(properties='initial')[source]¶

Parameters:	properties (str) – Set whether the experimental properties for the initial setup (“initial”, default), the optimized setup (“optimized”) or for both setups (“all”) shall be plotted. In the later case, both ellipsoids for one pair of parameters will be displayed within one plot.

Plot confidence ellipsoids for all parameter pairs. Since the number of plots is possibly big, all plots will be saved within a folder confidence_ellipsoids_scriptname in you current working directory rather than being displayed directly.

Optimum experimental design of multiple experiments¶

class casiopeia.doe.MultiDoE(doe_setups=[], optimality_criterion='A')[source]¶

The class casiopeia.doe.MultiDoE is used to construct a single experimental design problem from multiple experimental design problems defined via two or more objects of type casiopeia.doe.DoE.

This provides the possibility to design multiple experiments within one single optimization, so that the several experiments can focus on different aspects of the system which in combination then yields more information about the complete system.

Also, this functionality is in particular useful in case an experiment is limited to only small variable bounds, small time horizons, highly depends on the initialization of the system, or any other case when a single experiment might not be enough to capture enough information about a system.

Note

It is assumed that the system description used for setting up the several experimental design problems is the same!

Parameters:	doe_setups (list) – list of two or more objects of type casiopeia.doe.DoE

Parameter estimation for a Lotka-Volterra predator-prey-model¶

The aim of the application lotka_volterra_pe.py is to estimate the unknown parameters of a Lotka-Volterra predator-prey-model for experimentally received measurement data and given standard deviations for the measurements [1]. The predator-prey-model is an ODE of the form \(\dot{x} = f(x,p)\), given by

where \(\alpha = 1.0\) and \(\gamma = 1.0\), the states \(x\) and parameters \(p\) are defined as

and we can measure the full state, i. e.

The values resulting from the parameter estimation are

The results for the system simulation using the estimated parameters in comparison to the measurement data are shown in the figure below.

Figure: Simulation results for the Lotka-Volterra predator-prey-model using the estimated parameters, compared to the given measurement data

Parameter estimation for a pendulum model¶

The aim of the application pendulum_pe.py is to estimate the spring constant \(k\) of a pendulum model for experimentally received measurement data [2]. The pendulum model is an ODE of the form \(\dot{x} = f(x,u,p)\), given by

where \(m = 1.0\), \(L = 3.0\) and \(g = 9.81\), the states \(x\), controls \(u\) and parameters \(p\) are defined as

while the only control \(\psi\) is the initial actuation angle of the pendulum, and therefor stays constant over time. Also, we can measure the full state, i. e.

The value resulting from the parameter estimation is

The results for the system simulation using the estimated parameter in comparison to the measurement data are shown in the figures below.

Figure: Simulation results for the pedulum model using the estimated parameters, compared to the given measurement data

Parameter estimation for a model race car¶

The aim of the application 2d_vehicle_pe.py is to estimate the unknown parameters of a 2D race car model for experimentally received measurement data [3]. The race car and the interpretation of the model states are shown in the figure below [4].

Figure: Depiction of the race car showing the models states

The 2D model of the race car is an ODE of the form \(\dot{x} = f(x,u,p)\), given by

where the states \(x\), controls \(u\) and parameters \(p\) are defined as

and we can measure the full state, i. e.

The values resulting from the parameter estimation are

The results for the system simulation using the estimated parameter in comparison to the measurement data are shown in the figures below.

Figure: Simulation results for the race car model using the estimated parameters, compared to the given measurement data

An evaluation of the covariance matrix for the estimated parameters shows that the standard deviations of \(\hat{C_{1}}\) and \(\hat{C_{2}}\) are relatively small in comparison to their own values, while the standard deviations of the other parameters are relatively big.

This intends that the estimation results for the parameters \(\hat{C_{\text{m}_{1}}}\), \(\hat{C_{\text{m}_{2}}}\), \(\hat{C_{\text{r}_{2}}}\) and \(\hat{C_{\text{r}_{0}}}\) are probably not accurate, and might change substantially for other measurement and control data. Optimum experimental design can be an option to encounter this problem.

Optimum experimental design for a model race car¶

The aim of the application 2d_vehicle_doe_scaled.py is to solve an optimum experimental design problem for the 2D race car model from Parameter estimation for a model race car to obtain control values that allow for a better estimation of the unknown parameters of the model.

Optimized setup¶

Note

Running this optimization takes about 10 min on an Intel(R) Core(TM) i5-4570 3.20GHz CPU.

We use the A-criterion as objective for the experimental design (see Optimum experimental design). The results of the optimization can be analyzed and visualized with the script 2d_vehicle_doe_scaled_validation.py. The figure below shows the optimized control values in comparison to the initially used control values, while the suffix coll indicates that the values were obtained using collocation discretization.

Figure: Optimized control values in comparison to the initially used control values

The figure below shows a comparison of the simulated states values for both initially used and optimized control values, and with this, the effect of the optimization on the route of the race car and it’s velocity during the measurements.

Figure: Comparison of the simulated states values for initial and optimized controls

The quality of the optimized experimental setup in terms of estimated standard deviations of the unknown parameters is evaluated as follows

\[\begin{split}p_\text{O} = \begin{pmatrix} {C_{1, \text{O}}} \\ {C_{2, \text{O}}} \\ {C_{\text{m}_{1},\text{O}}} \\ {C_{\text{m}_{2},\text{O}}} \\ {C_{\text{r}_{2},\text{O}}} \\ {C_{\text{r}_{0},\text{O}}} \end{pmatrix} = \begin{pmatrix} {1.0} \\ {1.0} \\ {1.0} \\ {1.0} \\ {1.0} \\ {1.0} \end{pmatrix} \pm \begin{pmatrix} {1.93054150676} \\ {0.278656552587} \\ {1.96689422255} \\ {1.51815346784} \\ {3.42713773836} \\ {1.88475684297} \end{pmatrix}\end{split}\]

which indicates that the optimized setup is more appropriate for parameter estimation compared to the initial experimental design. Though, the estimated standard deviations are still relatively big in comparison to the scaled parameter values, so it would probably make sense to further refine the experimental design.