Prg_SFunctionEst Class Reference

Estimation of parameters and initial states for a model given as MEX S-function. More...

#include <Prg_SFunctionEst.h>

Inheritance diagram for Prg_SFunctionEst:

Public Member Functions
	Prg_SFunctionEst ()
	constructor
	~Prg_SFunctionEst ()
	destructor
const char *	name ()
	name SFunctionEst
Access methods for program specific members (If prefix: prg_)
int	nex () const
	number of experiments used for estimation
void	set_nex (int val)
	set number of experiments
int	multistage () const
	indicate if problem is treated with one stage per time interval
void	set_multistage (int val)
	set multistage flag
const MATP	ys_ref () const
	reference values for active model outputs (size: KK+1 . ny)
void	set_ys_ref (const MATP val)
	set reference outputs
const MATP	M () const
	measurement matrix M=dy/d(p,x0) (size: (KK+1)*ny . np+nx0)
const MATP	P () const
	precision matrix P=(M^TM)^(-1) (size: np+nx0 . np+nx0)
const MATP	COV () const
	covariance matrix (size: np+nx0 . np+nx0)
Read methods for model specific members (no If prefix).
Note that the prefix mdl_ is omitted in the detailed mathematical problem description.
const VECP	mdl_p () const
	model parameters
const VECP	mdl_p_min () const
	lower bounds for model parameters (note: default is 0)
const VECP	mdl_p_max () const
	upper bounds for model parameters
const IVECP	mdl_p_active () const
	indicate estimated parameters
const VECP	mdl_p_confidence () const
	confidence intervals of estimated parameters
const VECP	mdl_p_nominal () const
	nominal parameter values (for scaling)
const VECP	mdl_x0 () const
	model initial states (read only, note: write initial states of experiments via mdl_x0s)
const VECP	mdl_x0_min () const
	lower bounds for initial states
const VECP	mdl_x0_max () const
	upper bounds for initial states
const IVECP	mdl_x0_active () const
	indicate estimated initial states
const VECP	mdl_x0_confidence () const
	confidence intervals of estimated initial states, averaged over multiple experiments
const VECP	mdl_der_x0_min () const
	minimum for time derivative of x0
const VECP	mdl_der_x0_max () const
	maximum for time derivative of x0
const VECP	mdl_x_nominal () const
	nominal state values (for scaling) More...
const IVECP	mdl_u_order () const
const VECP	mdl_x_min () const
	lower bounds on states at stage boundaries (default: -Inf)
const VECP	mdl_x_max () const
	upper bounds on states at stage boundaries (default: Inf)
const IVECP	mdl_y_active () const
	indicate measured outputs
const VECP	mdl_y_nominal () const
	nominal output values (for scaling)
const MATP	mdl_x0s () const
	initial states for each experiment (size: nex . mdl_nx)
const MATP	mdl_us () const
	model inputs (size: KK+1 . mdl_nu)
const MATP	mdl_xs () const
	model states (size: KK+1 . mdl_nx)
const MATP	mdl_ys () const
	model outputs (read only, size: KK+1 . mdl_ny)
Write methods for model specific members (no If prefix).
void	set_mdl_p (const VECP v)
	set values of model parameter
void	set_mdl_p_min (const VECP v)
	set lower bounds for model parameters
void	set_mdl_p_max (const VECP v)
	set upper bounds for model parameters
void	set_mdl_p_active (const IVECP v)
	set estimated parameters
void	set_mdl_p_nominal (const VECP v)
	set nominal parameters
void	set_mdl_x0_min (const VECP v)
	set lower bounds for initial states
void	set_mdl_x0_max (const VECP v)
	set upper bounds for initial states
void	set_mdl_x0_active (const IVECP v)
	set estimated initial states
void	set_mdl_der_x0_min (const VECP v)
	set minimum for dx0dt
void	set_mdl_der_x0_max (const VECP v)
	set maximum for dx0dt
void	set_mdl_x_nominal (const VECP v)
	set nominal states
void	set_mdl_u_order (const IVECP v)
	set interpolation order for inputs
void	set_mdl_x_min (const VECP v)
	set lower bounds on states
void	set_mdl_x_max (const VECP v)
	set upper bounds on states
void	set_mdl_y_active (const IVECP v)
	set measured outputs
void	set_mdl_y_nominal (const VECP v)
	set nominal outputs
void	set_mdl_x0s (const MATP v)
	set initial states
void	set_mdl_us (const MATP v)
	set model inputs
void	set_mdl_xs (const MATP v)
	set model states
Public Member Functions inherited from Prg_SFunction
	Prg_SFunction ()
	constructor
	~Prg_SFunction ()
	destructor
const char *	mdl_name () const
	S-function name.
void	set_mdl_name (const char *str)
	set S-function name
const char *	mdl_path () const
	S-function path, including name, used for dynamic loading. More...
void	set_mdl_path (const char *str)
	set S-function path
const char *	mdl_args () const
	String representation of S-function arguments.
void	set_mdl_args (const char *str)
	set S-function arguments
const VECP	mdl_p () const
	parameters
void	set_mdl_p (const VECP value)
	set parameters
const VECP	mdl_x0 () const
	initial states
void	set_mdl_x0 (const VECP value)
	set initial states
Public Member Functions inherited from Omu_Program
	Omu_Program ()
	constructor
virtual	~Omu_Program ()
	destructor
virtual void	update (int kk, const adoublev &x, const adoublev &u, adoublev &f, adouble &f0, adoublev &c)
	High-level update for the optimization criterion, constraints, and discrete state equations.
virtual void	consistic (int kk, double t, const adoublev &x, const adoublev &u, adoublev &xt)
	High-level consistic (consistent initial conditions) routine for the initialization of continuous-time states from optimization variables x and u. More...

Protected Member Functions
void	write_active_mx_args (VECP p)
	write active parameters that are packed in p to _mx_args
Implementation of predefined methods.
See Also Omu_Program
void	setup_model ()
	load S-function
void	setup_stages (IVECP ks, VECP ts)
	Setup stages and time sampling. More...
void	setup (int k, Omu_VariableVec &x, Omu_VariableVec &u, Omu_VariableVec &c)
	Allocate states x, control parameters u, and constraints c for stage k. More...
void	setup_struct (int k, const Omu_VariableVec &x, const Omu_VariableVec &u, Omu_DependentVec &xt, Omu_DependentVec &F, Omu_DependentVec &f, Omu_Dependent &f0, Omu_DependentVec &c)
void	init_simulation (int k, Omu_VariableVec &x, Omu_VariableVec &u)
	Initialize problem variables using simulation. More...
void	update (int kk, const Omu_StateVec &x, const Omu_Vec &u, const Omu_StateVec &xf, Omu_DependentVec &f, Omu_Dependent &f0, Omu_DependentVec &c)
void	consistic (int kk, double t, const Omu_StateVec &x, const Omu_Vec &u, Omu_DependentVec &xt)
void	continuous (int kk, double t, const Omu_StateVec &x, const Omu_Vec &u, const Omu_StateVec &dx, Omu_DependentVec &F)
Protected Member Functions inherited from Prg_SFunction
void	read_mx_args (VECP p)
	read _mx_args into vector p
void	write_mx_args (VECP p)
	write vector p to _mx_args
bool	setContinuousTask (bool val)
	enable continuous sample time; return true if a continuous sample time exists and has been enabled
bool	setSampleHit (bool val)
	enable hit for all discrete sample times; return true if a discrete sample time exists and has been enabled

Protected Attributes
Omu_VariableVec	_mdl_p
	model parameters (note: default min is 0)
Omu_VariableVec	_mdl_x0
	initial states
Omu_VariableVec	_mdl_x
	state bounds
IVECP	_mdl_p_active
	indicate estimated parameters
VECP	_mdl_p_confidence
	confidence intervals for est parameters
IVECP	_mdl_x0_active
	indicate estimated states
VECP	_mdl_x0_confidence
	confidence intervals for est x0
VECP	_mdl_der_x0_min
	minimum for time derivative of x0
VECP	_mdl_der_x0_max
	maximum for time derivative of x0
IVECP	_mdl_u_order
	interpolation order (default: 1 (linear))
IVECP	_mdl_y_active
	indicate measured outputs
VECP	_mdl_p_nominal
	nominal parameter values (for scaling)
VECP	_mdl_x_nominal
	nominal state values (for scaling)
VECP	_mdl_y_nominal
	nominal output values (for scaling)
int	_nx
	number of states for optimizer
int	_np
	number of estimated parameters
int	_nx0
	number of estimated initial states
int	_ny
	number of reference outputs for estimation
int	_multistage
	treat as multistage problem
int	_nex
	number of experiments used for estimation
MATP	_mdl_x0s
	initial states for each experiment
IVECP	_exs
	index identifying experiment in estimation data
MATP	_mdl_us
	given model inputs
MATP	_mdl_xs
	calculated model states
MATP	_mdl_ys
	calculated model outputs
MATP	_ys_ref
	reference values for active outputs
double	_ssr
	sum of squared output residuals
MATP	_M2
	measurement matrix M=dy/d(p,x0)
MATP	_P2
	precision matrix P=(M^TM)^(-1)
MATP	_COV
	covariance matrix COV=
MATP	_dydx
	help variable dy/dx
MATP	_dxdpx0
	help variable dx/d(p,x0)
MATP	_dydpx0
	help variable dy/d(p,x0)
MATP	_dxfdx
	help variable dxf/dx
MATP	_dxfdpx0
	help variable dxf/d(p,x0)
Protected Attributes inherited from Prg_SFunction
char *	_mdl_name
	S-function name.
char *	_mdl_path
	full S-function path
char *	_mdl_args
	S-function parameters.
SimStruct *	_SS
	pointer to SimStruct
mxArray **	_mx_args
	S-function parameters after parsing.
int	_mdl_nargs
	number of S-function arguments
double	_t0_setup_model
	time used for initialization of model
int	_mdl_np
	number of model parameters
int	_mdl_nd
	number of discrete-time states
int	_mdl_nx
	number of model states (incl. discrete)
int	_mdl_nu
	number of model inputs
int	_mdl_ny
	number of model outputs
VECP	_mdl_p
	parameters
VECP	_mdl_x0
	initial states
bool	_mdl_needs_setup
	indicate that setup_model needs to be called as _mdl_name, _mdl_path, or _mdl_args changed

Detailed Description

Estimation of parameters and initial states for a model given as MEX S-function.

The estimation problem is formulated at the discrete time points \(t^{0,l}<t^{1,l}<\ldots<t^{K_l}\) with \(l=1,\ldots,N_{ex}\) a number of experiments. The data of all experiments is concatenated into combined vectors of time points, model inputs, and reference outputs, i.e. \(K_l+1=k_{0,l+1}\). In order to distinguish multiple experiments, it is assumed that \(t^{K_l}>t^{0,l+1}\); for instance each experiment could start with \(t^{0,l}=0\). Parameter values are constrained to be the same for all experiments, whereas initial states are individual for each experiment. Note that this is not a restriction as generally states can be initialized with parameters and as parameters can be modeled as constant states.

In the following all vector operations are defined element wise. The treated estimation problem reads

\[ \begin{array}{l} \displaystyle{} J\ =\ \sum_{l=1}^{N_{ex}} \sum_{k_l=k_{l,0}}^{K_l} \sum_{i=1}^{dim(I)} \left\{ \left[\frac{ys^{k_l}(I)-ys^{k_l}_{ref}}{y_{nominal}(I)}\right]^2 \right\}_i \quad\to\quad \min, \quad I = \mbox{find}(y_{active}), \end{array} \]

subject to the model given with the S-function methods mdlDerivatives \(f\), mdlOutputs \(g\) and mdlUpdate \(h\)

\[ \begin{array}{l} \displaystyle x(t) = \left[\begin{array}{ll} x_d^{kk}, & t\in[t^{kk},t^{kk+1}),\ kk=0,\ldots,KK \\ x_c(t) \end{array}\right], \\[3ex] \displaystyle x_d^{kk} = \frac{h[x_{nominal}\,x(t^-),\ u_{nominal}\,u(t^-)]} {x_{nominal_d}},\quad t=t^{kk},\ kk=1,\ldots,KK,\\[3ex] \displaystyle \frac{dx_c(t)}{dt} = \frac{f[x_{nominal}\,x(t),\ u_{nominal}\,u(t)]} {x_{nominal_c}}, \\[3ex] \displaystyle y(t) = \frac{g[p,\ x_{nominal}\,x(t),\ u(t)]} {y_{nominal}} \end{array} \]

with piecewise constant or linear interpolation of the inputs

\[ \begin{array}{ll} u_i(t)\ = & \left\{\begin{array}{ll} us_i^{k_l}, & i \in \mbox{find}(u_{order} = 0) \\[1ex] \displaystyle \frac{t^{k_l+1}-t}{t^{k_l+1}-t^{k_l}}\ us^{k_l} + \frac{t-t^{k_l}}{t^{k_l+1}-t^{k_l}}\ us^{k_l+1}, & \mbox{else} \end{array}\right. \\[5ex] & \quad t\in[t^{k_l},t^{k_l+1}), \quad k_l=k_{l,0},\ldots,K_l-1, \quad l=1,\ldots,N_{ex}, \end{array} \]

and subject to the constraints

\[ \begin{array}{rcccll} \displaystyle \left\{ \frac{p_{min}}{p_{nominal}} \right. &\le& \displaystyle \frac{p}{p_{nominal}} &\le& \left. \displaystyle \frac{p_{max}}{p_{nominal}} \right\}_i, \quad & i \in \mbox{find}(p_{active}), \\[3ex] \displaystyle \left\{ \frac{x^0_{min}}{x_{nominal}} \right. &\le& x(t^{0,l}) &\le& \displaystyle \left. \frac{x^0_{max}}{x_{nominal}} \right\}_i, \quad & i \in \mbox{find}(x^0_{active}), \quad l=1,\ldots,N_{ex}, \\[3ex] \displaystyle \left\{ \frac{der\_x^0_{min}}{x_{nominal}} \right. &\le& \dot{x}(t^{0,l}) &\le& \displaystyle \left. \frac{der\_x^0_{max}}{x_{nominal}} \right\}_i, \quad & i \in \mbox{find}(x^0_{active}), \quad l=1,\ldots,N_{ex}, \\[3ex] && \{ x(t^{0,l}) &=& \displaystyle \frac{x0s^l}{x_{nominal}} \}_i, \quad & i \notin \mbox{find}(x^0_{active}), \quad l=1,\ldots,N_{ex}, \\[3ex] \displaystyle \frac{x_{min}}{x_{nominal}} &\le& x(t^{k_l}) &\le& \displaystyle \frac{x_{max}}{x_{nominal}}, \quad & k_l=k_{l,0},\ldots,K_l, \quad l=1,\ldots,N_{ex}. \end{array} \]

The problem is treated as multistage problem with \(K_{N_{ex}}\) stages and one sample period per stage per default (multistage=1 or multistage=2). Discrete-time state variables with unknown initial value are introduced for the estimated parameters p, in order to preserve a sparse structure. Additional \(K_{N_{ex}}\) junction conditions (equality constraints) are introduced for the estimated parameters and \(K_{N_{ex}}-N_{ex}+1\) junction conditions are introduced for the state variables x. Alternatively the problem can be treated without stages hiding model states from the optimizer (multistage=0).

The states of all stages are either initialized with the initial states \(x0s^l\)

\[ \begin{array}{rcll} x_{initial}(t^{k_l}) &=& \displaystyle \frac{x0s^l}{x_{nominal}}, & k_l=k_{l,0},\ldots,K_l, \quad l=1,\ldots,N_{ex} \end{array} \]

or with individually given initial guesses \(xs\)

\[ \begin{array}{rcll} x_{initial}(t^{k_l}) &=& \displaystyle \frac{xs^{k_l}}{x_{nominal}}, & k_l=k_{l,0},\ldots,K_l, \quad l=1,\ldots,N_{ex}. \end{array} \]

All but the initial states of the experiments may be initialized with the results of an initial-value simulation for given initial states and inputs (multistage=0 or multistage=1). With multistage=2, individual initial states are used for each stage, resulting in the multiple shooting method. This might be useful if no sensible initial guesses can be given for unknown parameters, so that a simulation fails.

states and inputs. If the problem is treated as multistage problem, then not performing the simulation results in the multiple shooting method. This might be useful if no sensible initial guesses can be given for unknown parameters, so that a simulation fails.

Model inputs, states and outputs can be accessed through

\[ \begin{array}{l} us^{k_l} = u_{nominal}\,u(t^{k_l}), \\[1ex] xs^{k_l} = x_{nominal}\,x(t^{k_l}), \\[1ex] ys^{k_l} = y_{nominal}\,y(t^{k_l}), \quad k_l=k_{l,0},\ldots,K_l, \quad l=1,\ldots,N_{ex}. \end{array} \]

In addition to the estimated parameters and initial states, the measurement matrix M is calculated in each optimization iteration. It is defined as

\[ M = \left[\begin{array}{cc} M_p^0, & M_{x^0}^0 \\[2ex] M_p^1, & M_{x^0}^1 \\[2ex] \vdots, & \vdots \\[2ex] M_p^{K_{N_{ex}}}, & M_{x^0}^{K_{N_{ex}}} \end{array}\right] \]

with the sub-matrices

\[ \begin{array}{r} \displaystyle M_p^{k,l} = \frac{d\,\displaystyle\frac{ys^{k,l}} {y_{nominal}}(I_y)} {d\,\displaystyle\frac{p}{p_{nominal}}(I_p)}, \quad M_{x^0}^{k,l} = \frac{d\,\displaystyle\frac{ys^{k,l}}{y_{nominal}}(I_y)} {d\,\displaystyle\frac{x0s^l}{x_{nominal}}(I_{x^0})}, \quad k_l=k_{l,0},\ldots,K_l, \quad l=1,\ldots,N_{ex}, \\[7ex] I_y = \mbox{find}(y_{active}),\ I_p = \mbox{find}(p_{active}), \ I_{x^0} = \mbox{find}(x^0_{active}). \end{array} \]

Note that the sub-matrices for estimated initial states of all experiments are stored in compressed form in one column, even though individual initial states are being estimated for each experiment.

Some simple quality measures are directly calculated in order to support a first interpretation of estimation results. These are the precision matrix

\[ P = (M^TM)^{-1}, \]

and the covariance matrix

\[ COV = s_R^2 P, \qquad s_R^2=\frac{J}{n}, \qquad n=\mbox{dim}(I_y)(K+1) - \mbox{dim}(I_p) - N_{ex}\mbox{dim}(I_{x^0}) - 1, \]

assuming that the residual variance \(s_R^2\) is equal to the measurement variance. Last but not least confidence intervals for estimated parameters and initial states \(px^0\) are obtained applying the t-Test:

\[ px^0_{confidence}\ =\ |px^0_{estimated} - px^0_{actual}|\ = \ px_{nomainal} t_{\alpha/2,n}\sqrt{\mbox{diag(COV)}}, \qquad \alpha=0.05. \]

Member Function Documentation

void Prg_SFunctionEst::init_simulation ( int  k, Omu_VariableVec &  x, Omu_VariableVec &  u  )

protectedvirtual

Initialize problem variables using simulation.

A problem specification may set state variables x and/or control parameters u prior to the evaluation of stage k. The default version causes an initial value simulation based on initial values given in setup.

Reimplemented from Omu_Program.

const VECP Prg_SFunctionEst::mdl_x_nominal ( ) const

inline

nominal state values (for scaling)

interpolation order (0 (constant) or 1 (linear), default: 1)

References _mdl_x_nominal.

void Prg_SFunctionEst::setup ( int  k, Omu_VariableVec &  x, Omu_VariableVec &  u, Omu_VariableVec &  c  )

protectedvirtual

Allocate states x, control parameters u, and constraints c for stage k.

Additionally setup bounds and initial values.

Implements Omu_Program.

void Prg_SFunctionEst::setup_stages ( IVECP  ks, VECP  ts  )

protectedvirtual

Setup stages and time sampling.

The default implementation sets up an optimization problem without stages.

Reimplemented from Omu_Program.

---

The documentation for this class was generated from the following file:

• Prg_SFunctionEst.h