Least-Squares Parameter Estimation for State-Space Models with State Equality Constraints

TL;DR

This paper develops a method for incorporating state equality constraints into least-squares estimation of state-space models, handling both known and uncertain constraints in time-invariant and time-varying systems.

Contribution

It reformulates state equality constraints as matrix constraints and applies vectorized least squares, extending to uncertain constraints and dynamic systems.

Findings

Effective incorporation of constraints improves estimation accuracy.

Method handles both known and uncertain constraints.

Applicable to time-invariant and time-varying systems.

Abstract

If a dynamic system has active constraints on the state vector and they are known, then taking them into account during modeling is often advantageous. Unfortunately, in the constrained discrete-time state-space estimation, the state equality constraint is defined for a parameter matrix and not on a parameter vector as commonly found in regression problems. To address this problem, firstly, we show how to rewrite the state equality constraints as equality constraints on the state matrices to be estimated. Then, we vectorize the matricial least squares problem defined for modeling state-space systems such that any method from the equality-constrained least squares framework may be employed. Both time-invariant and time-varying cases are considered as well as the case where the state equality constraint is not exactly known.

Equations

$$x_{k+1}$$

$$y_k$$

$$S{x}_k=s,$$

$$y^{T}_{k+1}=\left[\begin{array}[]{cc}y_{k}^{T}&u^{T}_{k}\end{array}\right]\left[\begin{array}[]{c}A^{T}\\ B^{T}\end{array}\right]+e_{k+1}^{T},$$

$$e_{k+1}\triangleq v_{k+1}-A{v}_k+G{w}_k.$$

$$\mathcal{Y}=\mathcal{X}\Theta +\mathcal{E},$$

$$\displaystyle\mathcal{Y}\triangleq\left[\begin{array}[]{c}y_{1}^{T}\\ \vdots\\ y_{N}^{T}\end{array}\right],~{}~{}\mathcal{X}\triangleq\left[\begin{array}[]{cc}y_{0}^{T}&u^{T}_{0}\\ \vdots&\vdots\\ y_{N-1}^{T}&u^{T}_{N-1}\end{array}\right],~{}~{}\mathcal{E}\triangleq\left[\begin{array}[]{c}e_{1}^{T}\\ \vdots\\ e_{N}^{T}\end{array}\right],$$

$$\displaystyle\Theta\triangleq\left[\begin{array}[]{c}A^{T}\\ B^{T}\end{array}\right],$$

$$J_{LS}(\hat{\Theta })\triangleq {\left(\mathcal{Y}-\mathcal{X}\hat{\Theta }\right)}^T\left(\mathcal{Y}-\mathcal{X}\hat{\Theta }\right).$$

$$z_k={\psi }_{k-1}^T\hat{\theta }+{\xi }_k,$$

$$Z=\Psi \hat{\theta }+\Xi .$$

$$D\hat{\theta }=d,$$

$$J_{\text{LS}}(\hat{\theta })\triangleq {\left(Z-\Psi \hat{\theta }\right)}^T\left(Z-\Psi \hat{\theta }\right).$$

$${\hat{\theta }}_{\mathrm{CLS}}={\mathcal{P}}_{\mathcal{N}(D)}{\hat{\theta }}_{\mathrm{LS}}+(I_{n_p}-{\mathcal{P}}_{\mathcal{N}(D)})\bar{d},$$

$${\hat{\theta }}_{\text{LS}}=({\Psi }^T\Psi {)}^{-1}{\Psi }^{T}Z,$$

$${\mathcal{P}}_{\mathcal{N}(D)}\triangleq I_{n_p}-LD,$$

$$L\triangleq ({\Psi }^T\Psi {)}^{-1}{D}^T[D({\Psi }^T\Psi {)}^{-1}{D}^T{]}^{-1},$$

$$\bar{d}\triangleq D^T(D{D}^T{)}^{-1}d$$

$${\hat{\theta }}_{\mathrm{rCLS}}={\left({\Psi }^T\Psi +\mu D^{T}D\right)}^{-1}({\Psi }^{T}Z+\mu D^{T}d),$$

$$K_k$$

$${\hat{\theta }}_{\text{CLS},k}$$

$$P_{\text{CLS},k}$$

$${\hat{\theta }}_{\text{CLS},0}$$

$$P_{\text{CLS},0}$$

$$z_k={\psi }_{k-1}^T{\hat{\theta }}_k+{\xi }_k.$$

$$D_k{\hat{\theta }}_k=d_k.$$

$$K_k$$

$${\hat{\theta }}_{\mathrm{WLS},k}$$

$$P_{\mathrm{WLS},k}$$

$$L_k$$

$${\hat{\theta }}_{\text{WCLS},k}$$

$$SG$$

$$SA$$

$$SB$$

$$\Theta D_1=D_2,$$

Full text

1 Introduction

In some dynamic systems, dynamics evolve with variables satisfying inequality or equality constraints (Goodwin \BOthers., \APACyear2005). For example, the species concentrations are non-negative in chemical reactions (Massicotte \BOthers., \APACyear1995). Likewise, in the quaternion-based attitude representation, the attitude vector must have unitary norm (Crassidis \BBA Markley, \APACyear2003); and, for ground vehicle tracking problems, the road networks can be viewed as equality constraints on the trajectory (Xu \BOthers., \APACyear2013). The combination of single or multiple cells in biological processes may be represented as a compartment with constant volume (Mohler, \APACyear1974). In addition, compartmental models have applications in classical circuit models, structural models and complex networks, among others, as pointed out in Bernstein \BBA Hyland (\APACyear1993).

In this work, we are specifically concerned with linear state-space dynamic systems satisfying linear equality constraints on the state vector. The scenario we have in mind is the one in which we have dynamical data collected from the dynamic system as well as auxiliary information (written as equality constraints on the state vector) from first principles. Consider the following examples: localization of a land vehicle for which the road map represents a constraint on the trajectory (Xu \BOthers., \APACyear2017); the flight formation of two targets, where the distance between the targets is constant (Xu \BOthers., \APACyear2013); the monitoring of the nitrogen flow in a tropical forest where the amount of nitrogen is constant (Walter \BBA Contreras, \APACyear1999); the experiment of the wet granulation of lactose with deionized water carried out in a ploughshare mixer with constant volume (Lee \BOthers., \APACyear2017); and an interconnected tank system for which prior information on the total volume is available (Hölzel \BBA Bernstein, \APACyear2014).

At this point, one may argue that variable reduction (Hölzel \BBA Bernstein, \APACyear2014; Li, \APACyear2016) may be employed to avoid enforcing the equality constraint on the state vector. However, this approach yields a reduced state vector with a different physical meaning, which is not desirable in many applications. Moreover, if the equality constraint is time varying, keeping a constant state vector parametrization is of interest.

In the last years, the problem of state estimation for both linear and nonlinear equality-constrained dynamic systems has received great attention from the community (Babacan \BOthers., \APACyear2008; Teixeira \BOthers., \APACyear2009; Simon, \APACyear2010; Teixeira \BOthers., \APACyear2008; Xu \BOthers., \APACyear2013; Rengaswamy \BOthers., \APACyear2013; Duan \BBA Li, \APACyear2015; Xu \BOthers., \APACyear2017). The problem of modeling such systems is less often addressed (Xu \BOthers., \APACyear2013; Li, \APACyear2016; Xu \BOthers., \APACyear2017). In the latter works, a two-step modeling procedure is employed. An unconstrained model (auxiliary dynamics) is first obtained and, by projection, such model is fused with the state equality constraint.

Parameter estimation with known equality linear constraints is a solved problem. Auxiliary information such as static function, static gain, and fixed-point location, can be written in the form of linear equality constraints on the parameters of NARX (Nonlinear autoregressive with exogenous inputs) polynomial and RBF (Radial basis function) network models, for instance (Teixeira \BBA Aguirre, \APACyear2011; Aguirre \BOthers., \APACyear2007). If the dynamics are time invariant, one may use the batch equality constrained least squares (Björck, \APACyear1996; Draper \BBA Smith, \APACyear1998). For problems in which auxiliary information is uncertain, the compromise between prediction performance and the equality constraint satisfaction is treated by means a tuning parameter in Teixeira \BBA Aguirre (\APACyear2011). In Arablouei \BBA Dogançay (\APACyear2015) the relaxed solution of the batch equality constrained least squares is addressed to solve the same problem. For a recursive solution, it suffices to use the classical recursive least squares with a proper initialization as shown in Zhou \BOthers. (\APACyear2001); Zhu \BBA Li (\APACyear2007). For convenience, in this work, we present this result using a different perspective in Proposition 3.1. However, for time-varying systems, the equality parameter constraint must be enforced at every time instant in the recursive least squares equations (Alenany \BBA Shang, \APACyear2013). In this regard Vincent \BBA Chaumette (\APACyear2018), exploring connections between Kalman filter and least squares, enforce equality constraint on the Kalman gain (Teixeira \BOthers., \APACyear2008) in order to guarantee that the estimator is unbiased. Conversely, in this manuscript, we enforce constraints on the model matrices in order to guarantee a model whose state vector satisfy an equality constraint.

If we assume that all state components are directly measured, then least squares methods may be used to estimate the matrices of the linear state-space model. Otherwise, subspace methods must be used (Trnka \BBA Havlena, \APACyear2009; Alenany \BOthers., \APACyear2011; Privara \BOthers., \APACyear2012; Alenany \BBA Shang, \APACyear2013; Wang \BOthers., \APACyear2018) with least squares as a possible internal step. Consider the case of fully measured state vector, for which auxiliary information on the state vector is known in the form of an equality constraint. How to estimate the state matrices of the system such that the free-run prediction of its state vector satisfies the known equality constraint? To address this question (Problem 2.1), one must first be able to mathematically map the equality constraint on the state vector onto an equality constraint on the state space matrices (parameters) to be estimated. In this manuscript, this problem is solved for time-varying linear dynamic systems by adapting a result from Teixeira \BOthers. (\APACyear2009); see Lemmas 4.1 and 4.2. However, the aforementioned equality-constrained least-squares methods cannot be used to enforce such equality constraints because, in this case, the constraint is defined for a parameter matrix and not on a parameter vector as commonly found in regression problems; see Remark 4. To circumvent this problem, we vectorize the matricial least squares problem defined for modeling state-space systems using the vectorization operator and Kronecker product as in Privara \BOthers. (\APACyear2012) such that the existing equality constrained least squares framework may be employed; see Proposition 4.3. The contributions of this manuscript are: (i) to address the two aforementioned problems as a single mathematical problem and to solve such state-space modeling problem with state equality constraints, and (ii) to explore the connections among the papers that address similar problems in the literature (Teixeira \BBA Aguirre, \APACyear2011; Arablouei \BBA Dogançay, \APACyear2015; Zhou \BOthers., \APACyear2001; Zhu \BBA Li, \APACyear2007; Alenany \BBA Shang, \APACyear2013). Here both time-invariant and time-varying cases are considered. Finally, as in Teixeira \BBA Aguirre (\APACyear2011), the case in which the auxiliary information is uncertain is also addressed.

This document is organized as follows. Section 2 formulates the gray-box system identification problem under investigation. In Section 3 we review known equality-constrained parameter estimation methods for both time-invariant and time-varying systems. Section 4 solves the problem formulated in Section 2, presenting the main contributions of this manuscript. In Section 5, numerical examples illustrate the applicability of the proposed approaches. Finally, in Section 6, the concluding remarks are discussed.

Notation is set as follows in this manuscript. $I_{n}$ and $0_{m\times n}$ stand respectively for the $n$-dimensional identity matrix and $m\times n$-dimensional zero matrix. $\otimes$ is the Kronecker product and ${\rm vec}$ is the vectorizer operator.

2 Problem Statement

Consider the linear discrete-time state-space system

[TABLE]

where $x_{k}\in\mathbb{R}^{n}$ is the state vector, $u_{k}\in\mathbb{R}^{p}$ is the input vector, $y_{k}\in\mathbb{R}^{m}$ is the measured output vector, $w_{k}\in\mathbb{R}^{q}$, $q<n$, is the zero-mean process noise with covariance $Q\in\mathbb{R}^{q\times q}$ and $v_{k}\in\mathbb{R}^{m}$ is the zero-mean measurement noise. Note that all the states are assumed to be measured. Assume that the noise terms are mutually uncorrelated. The matrices $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times p}$ and $G\in\mathbb{R}^{n\times q}$ are not assumed to be known. Assume that the system (1) is asymptotically stable. In addition, assume that the state vector satisfies the equality constraint

[TABLE]

where $S\in\mathbb{R}^{n_{r}\times n}$, $s\in\mathbb{R}^{n_{r}\times 1}$ and $1\leq n_{r}\leq n-q$ is the number of constraints. Without loss of generality, we assume that ${\rm rank}(S)=n_{r}$.

The state-space model (1)-(2) can be rewritten as

[TABLE]

where

[TABLE]

Next, assume that sequences of $u_{k}$ and $y_{k}$ are known for $k=0,\ldots,N$, such that we have

[TABLE]

where

[TABLE]

[TABLE]

where $\mathcal{{Y},~{}E}\in\mathbb{R}^{N\times n}$, $\mathcal{{X}}\in\mathbb{R}^{N\times(n+p)}$ and ${\Theta}\in\mathbb{R}^{(n+p)\times n}$.

Define the least squares cost function

[TABLE]

Problem 2.1**.**

The equality-constrained state-space modeling problem is to obtain the minimizer $\hat{{{\Theta}}}_{\rm CLS}$ of (20) such that the states $x_{k}$ of the corresponding process model (1) satisfy the equality constraint (3).

Recall that this paper does not address the problem of state estimation, although the parameter estimation problem under investigation can be recast as a state estimation problem under proper assumptions.

3 Background on Equality-Constrained

Least Squares

3.1 Time-invariant case

Consider the linear regression model

[TABLE]

where $z_{k}\in\mathbb{R}$ is the measured output, $\psi_{k-1}\in\mathbb{R}^{n_{p}}$ is the known regressor vector, $\xi_{k}\in\mathbb{R}$ is the residue, and $\hat{\theta}\in\mathbb{R}^{n_{p}}$ is the unknown parameter vector to be estimated. Recall that (21) may represent the dynamic model for a linear-in-the-parameters MISO system. Assume that a set of $N$ observations of $z$ and $\psi$ in (21) are available such that we have

[TABLE]

Assume that $\Psi$ has full column rank such that $\Psi$ is left invertible.

Now, assume that the parameters $\hat{\theta}$ must satisfy a set of $n_{r}$ linear equality constraints given by

[TABLE]

where $d\in\mathbb{R}^{n_{r}}$ and $D\in\mathbb{R}^{n_{r}\times n_{p}}$ with ${\rm rank}(D)=n_{p}$. Next, define the least squares cost function

[TABLE]

Then, the minimizer of (24) subject to (23) is given by Björck (\APACyear1996); Draper \BBA Smith (\APACyear1998)

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

is an offset. Therefore, the equality constraint $D\hat{\theta}_{\rm{CLS}}=d$ is exactly satisfied. The estimator (26) is known as the classical least squares (LS) and (25) is known as the equality-constrained least squares (CLS).

Remark 1*.*

Augment the matrices in (22) and (24) by appending a weighted form of the linear constraints (23). Then, the optimal solution $\hat{{\theta}}_{\textrm{LS}}$ (25) is approximated by the relaxed solution (Arablouei \BBA Dogançay, \APACyear2015)

[TABLE]

where $\mu\gg 1$ is the weight associated to the constraints. If one tunes $\mu\to\infty$, then $\hat{{\theta}}_{\rm{rCLS}}\to\hat{{\theta}}_{\rm{CLS}}$. For applications in which the constraints (23) are not precisely known, the relaxed solution (30) is indicated. Indeed, other related least squares approaches may be used to solve this problem; see Teixeira \BBA Aguirre (\APACyear2011, Section 4.2).

In Zhou \BOthers. (\APACyear2001), the recursive counterpart of (25) is investigated. Interestingly, Zhou \BOthers. (\APACyear2001) shows that the CLS and the LS have the same recursive formulas, differing only at the initial values. In Zhu \BBA Li (\APACyear2007) it is shown how to initialize the recursive least square equations in order to guarantee that the corresponding estimates satisfy (23), $\forall k$. For mathematical convenience, Zhu \BBA Li (\APACyear2007) and Zhou \BOthers. (\APACyear2001) derive the recursive equations using Greville formulas, yielding equations in a non-standard format.

Next, for simplicity, we present the recursive counterpart of (25) in a more conventional format. For $k=1,\ldots,N$, we have

[TABLE]

with initial values

[TABLE]

where the projector $\mathcal{P}_{\mathcal{N}(D)}$ (27) guarantees that $\hat{\theta}_{{\textrm{CLS}},0}$ and $P_{{\textrm{CLS}},0}$ are compatible with (23) for any ${\theta}_{0}$ and $P_{0}$. We point out that (31)-(33) correspond to the classical recursive least squares (RLS). The next result proves in a simple way that, if the RLS is properly initialized as in (34)-(35), then its estimates satisfy (23), $\forall k$. As mentioned above, a similar result is presented in Zhu \BBA Li (\APACyear2007) using Greville formulas.

Proposition 3.1**.**

Assume that the initial parameter estimates $\hat{\theta}_{{\textrm{CLS}},0}$ and $P_{{\textrm{CLS}},0}$ are given by (34)-(35) and are used to initialize the RLS equations (31)-(33). Then the estimates given by (31)-(33) satisfy the time-invariant constraint $D\hat{\theta}_{{\rm CLS},k}=d,\forall k$.

Proof.

See Appendix A. ∎

3.2 Time-varying case

Assume now that the parameters may vary with time such that (21) is replaced by

[TABLE]

Also, assume that $\hat{\theta}_{k}$ satisfy the known time-varying constraint

[TABLE]

For $k=1,\ldots,N$, the recursive time-varying counterpart of (25) is given by Alenany \BBA Shang (\APACyear2013)

[TABLE]

where $0\ll\lambda\leq 1$ is the forgetting factor. We refer to this method as the recursive weighted constrained LS (RWCLS).

Remark 2*.*

In Alenany \BBA Shang (\APACyear2013), part of the identification data is used to estimate offline the initial parameters $\hat{\theta}_{{\textrm{WLS}},0}$ and $P_{{\textrm{WLS}},0}$ for (38)-(42) by means of the constrained batch algorithm (25). Instead, we suggest to initialize (38)-(42) as in (34)-(35).

4 Equality-Constrained Least Squares for State-Space Modeling

In order to solve Problem 2.1, first it is necessary to map the constraint on the state vector $x_{k}$ given by (3) to a constraint on the parameter vector $\Theta$ given by (19).

The next results address this point by indicating conditions for a state-space model to have a state vector satisfying an equality constraint.

Lemma 4.1**.**

(Teixeira \BOthers., \APACyear2009, Proposition 3.1)* For the system given by (1), assume that*

[TABLE]

Then, for all $k\geq 1$, $Sx_{k}=s$, where $s=Sx_{0}$.

Lemma 4.2**.**

Assume that the state vector of system (1) satisfies the equality constraint Sx_{k}=s,~{}\text{\forall k\geq 0}. Then, the relations (43), (44) (45) are hold.

Proof.

Multiplying (1) by $S$, we obtain $Sx_{k}=SAx_{k-1}+SBu_{k-1}+SGw_{k-1}$. Then $Sx_{k}=s$ implies that $SA=S$, $SB=0_{n_{r}\times p}$ and $SG=0_{n_{r}\times q}$. ∎

Lemma 4.1 gives the conditions (43)-(45) for the dynamic system to satisfy (3), while Lemma 4.2 proves the counterpart. In other words, the previous results provide conditions for process model (1) to be compatible (Li, \APACyear2016) with the state equality constraint (3).

Remark 3*.*

In Teixeira \BOthers. (\APACyear2009, Proposition 3.2), it is proved that if (3) holds, then the system is not controllable in $\mathbb{R}^{n}$ from the process noise $w_{k}$, but it is rather controllable in the subspace defined by (3). Then, we can replace $Gw_{k-1}$ in (1) by $\tilde{w}_{k-1}\triangleq Gw_{k-1}$ with singular noise covariance $\tilde{Q}\triangleq GQG^{T}$ providing (43) is verified as in Lemma 4.1. In so doing, we focus on relations (44)-(45), which are related to the matrices $A$ and $B$ to be estimated.

From (44)-(45), we obtain the equality constraint on the parameter matrix

[TABLE]

where $\Theta$ is given by (19) and

[TABLE]

where $D_{1}\in\mathbb{R}^{n\times n_{r}}$ and ${D}_{2}\in\mathbb{R}^{(n+p)\times n_{r}}$.

Remark 4*.*

Note that the estimator (25) cannot be used to solve Problem 2.1 due to matrix size incompatibility. In (25), the equality constraint (23) is enforced on vector $\hat{\theta}$, whereas in Problem 2.1, the equality constraint (46) is enforced on matrix $\hat{\Theta}$, which is right-multiplied by $D_{1}$.

The next result rewrites the matrix equations (6) and (46) onto vectorized equations like (21) and (23) such that a classical equality-constrained least squares problem is obtained. In so doing, we have a solution for Problem 2.1. Likewise, based on this result, the recursive solution can also be obtained for both time-invariant and time-varying cases.

Proposition 4.3**.**

For the linear regression model (4) with parameter $\hat{\Theta}$ given by (19), the parameter estimate $\hat{\Theta}_{\rm CLS}$ minimizes $J_{\rm LS}(\hat{\Theta})$ given by (20) subject to the equality constraint (46), if and only if $\hat{\theta}_{\rm CLS}$ given by (25) minimizes $J_{\rm LS}(\hat{\theta})$ given by (24) subject to the equality constraint (23) with

[TABLE]

where $\theta\in\mathbb{R}^{(n^{2}+np)}$, $Z\in\mathbb{R}^{Nn}$, $\Psi\in\mathbb{R}^{Nn\times(n^{2}+np)}$, ${{D}}\in\mathbb{R}^{(n_{r}n+n_{r}p)\times(n^{2}+np)}$, and ${{d}}\in\mathbb{R}^{(n_{r}n+n_{r}p)}$.

Proof.

This proof has two parts. First, we rewrite (6) as (21). This is done by using the following relation (Bernstein, \APACyear2005)

[TABLE]

where $\mathcal{M}$, $\mathcal{N}$ and $\mathcal{O}$ are real matrices of appropriate size.

The deterministic part of the model (4) can be rewritten as

[TABLE]

which is the form of (21).

Second, we rewrite (46) as (23). For a data set of size $N$, we have

[TABLE]

which is an equation of the form (21), that is, $Z=\Psi\theta$, with $Z$ and $\Psi$ set as in (56)-(61).

Thus, the cost function $J_{\rm LS}(\hat{\theta})$ given by (24) is equivalent to the cost function $J_{\rm LS}(\hat{\Theta})$ given by (20).

Likewise, the equality constraint (46) can be rewritten as

[TABLE]

which is an equation of the form (23), with $D$ and $d$ set as in (66). ∎

Corollary 4.4**.**

Given the vectorizing definitions (56)-(66), the recursive equality-constrained parameter estimates $\hat{\Theta}_{{\rm CLS},k}$ are given by (31)-(35).

Proof.

Given the result of Proposition 4.3 and that the recursive equations (31)-(35) are equivalent to the equality constrained least squares (25) (see Proposition 3.1), we prove this result. ∎

Remark 5*.*

For the time-varying counterpart of (1)-(3), using definitions similar to (56)-(66), the recursive equality-constrained parameter estimates $\hat{\Theta}_{{\rm CLS},k}$ can be obtained from (38)-(42).

Note that (1)-(2) characterizes an output-error model. So, the next result proves that the LS estimator is biased for such type of model.

Proposition 4.5**.**

For the state-space model (1)-(2), the LS estimator given by (26) and (56)-(61) is biased, that is, $\rm{E}[\hat{\theta}]-\theta\neq 0$.

Proof.

See Appendix B. ∎

Thus, in this work, we use algorithms based on the extended LS (Ljung, \APACyear1987); however, for brevity, we omit the term “extended”. Other unbiased estimators could be used instead.

5 Simulated Results

5.1 Compartmental system: time-invariant case

Consider the linear discrete-time compartmental model (Teixeira \BOthers., \APACyear2009) represented by (1)-(2) involving mass exchange among compartments whose matrices are given by

[TABLE]

with state vector $x_{k}$ $\in$ $\mathbb{R}^{3}$ composed by the amount of mass in each compartment, initial condition $x_{0}$ $=$ $[1~{}~{}1~{}~{}1]^{T}$ , and process noise and observation noise covariance matrices $\tilde{Q}$ $=$ $\sigma_{w}^{2}GG^{T}$, where G=\left[\begin{array}[]{cc}0.05&-0.03\\ -0.02&0.01\\ -0.03&0.02\end{array}\right], and $R_{k}=\sigma_{v}^{2}I_{2\times 2}$.

One realization of simulated identification data for this system is shown in Fig. 1 for $\sigma_{w}$ $=$ $1.0$ and $\sigma_{v}$ $=$ $0.1$. Note that conditions of Lemmas 4.1 and 4.2 hold for (69) such that the trajectory of $x_{k}$ $\in$ $\mathbb{R}^{3}$ lies on the plane (3), whose parameters are assumed to be known and are given by

[TABLE]

that is, mass conservation is verified. The validation data is simulated with different initial condition $x_{0}$ $=$ $[2~{}~{}1~{}~{}0]^{T}$.

We investigate a 1000-run Monte Carlo simulation testing identification data with different noise realizations for $\sigma_{w}$ $=$ $1.0$ and $\sigma_{v}$ $=$ $0.1$. The LS given by (26) for the state-space system, for which (6)-(20) are defined, is used to yield estimates for the matrix $A$ as discussed in Section 2. Likewise, as indicated by Proposition 4.3, CLS was also employed.

Fig. 2 shows the results regarding the Monte Carlo validation of the obtained models (mean values with two standard-devation confidence interval). In order to quantify the fit between the simulation of the system and the identified models, we use the root-mean-square error for each $i$th state component, $i=1,\ldots,n$,

[TABLE]

where $N$ is the length of the measured data and $m_{r}=1,\ldots,M$, where $M$ is the number of realizations. Table 5.1 shows the mean $\overline{\textrm{RMSE}}_{n}$ and standard deviation $\sigma_{\rm{RMSE}_{n}}$ of the RMSE for each state. Note that the performance of the model obtained with CLS is better than the model estimated by LS. That is, the auxiliary information about mass conservation was useful.

In addition, we consider the case where the auxiliary information is uncertain. Suppose that the uncertain state equality constraint is assumed to be given by (3) with

[TABLE]

The rCLS given by (30) is used to estimate the state-space model with the uncertain auxiliary information (72). As in Arablouei \BBA Dogançay (\APACyear2015), we tuned the parameter $\mu$ in order to obtain models with good prediction performance. The results are shown in Fig. 2 and Table 5.1 for $\mu_{1}=5\times 10^{3}$ ($\rm{rCLS}_{1}$) and $\mu_{2}=5\times 10^{4}$ ($\rm{rCLS}_{2}$). Note that the results yielded by ($\rm{rCLS}_{1}$) are better than those from LS. Moreover, results from $\rm{rCLS}_{2}$ almost coincide to those from CLS. Then, the appropriate use of uncertain prior information may improve the quality of the estimated model, as discussed in Teixeira \BBA Aguirre (\APACyear2011).

We also test the recursive solution to this problem as indicated by Corollary 4.4. RLS is properly initialized as in (34)-(35), with $\theta_{0}$ $=$ $0_{9\times 1}$ and $P_{0}$ $=$ $10^{3}I_{9\times 9}$, yielding RCLS and is compared with the batch LS and CLS estimates in Fig. 2 and Table 5.1. Note that, when compared to LS, the use of auxiliary information (70) in the initialization of RCLS improves the performance of the estimated model.

5.2 Compartmental system: time-varying case

We now consider a time-varying compartmental system. As in a reconfigurable system, we consider the case in which the linear dynamics switches among three different modes. For instance, this may be the case for a multi-tank system with reconfigurable valves. The first mode is simulated with $A_{1}$ as in (69). The second and third modes are described by the matrices

[TABLE]

[TABLE]

The matrices $B$ and $C$ are defined as in (69) for all modes.

A typical realization of the simulated identification data is shown in Fig. 3 for $\sigma_{w}$ $=$ $10$ and $\sigma_{v}$ $=$ $1$ and $x_{0}$ $=$ $[20~{}~{}20~{}~{}10]^{T}$. The mass conservation is verified for all operating points and the assumedly known parameters of (3) are given by

[TABLE]

Note that the conditions of Lemmas 4.1 and 4.2 are verified for all modes. A new initial condition is set as $x_{0}$ $=$ $[15~{}~{}10~{}~{}25]^{T}$ to simulate the validation data.

We generate $1000$ Monte Carlo simulations with different noise realizations for $\sigma_{w}$ $=$ $10$ and $\sigma_{v}$ $=$ $1$ in order to obtain the identification data. For each running simulation we employ both RWLS given by (38)-(40) and RWCLS given by (38)-(42) to estimate the time-varying model. These recursive estimators are randomly initialized with an arbitrary initial condition ${\rm vec}(\Theta_{0})\in\mathbb{R}^{{9\times 1}}$ given by a normal distribution with $\sigma_{\Theta}=1$ and $P_{0}=10^{4}I_{9\times 9}$. Recall that in Alenany \BBA Shang (\APACyear2013) part of the identification data is used to estimate de initial conditions by means of the batch algorithm (25); see Remark 2. Here, we use the result given by (34)-(35) to more conveniently proceed the identification procedure using the RWCLS. The forgetting factor $\lambda$ is set to $0.95$.

The Monte Carlo validation results for RWLS and RWCLS are shown in Fig. 4. We would like to draw attention to the variance of the estimated models. For the three different modes, we verify that the performance of the model obtained with RWCLS is better than the model estimated by RWLS. That is, the prior information about mass conservation improved the quality of the estimated model.

5.3 Practical application: forest ecosystem

Consider the compartmental model of nitrogen flow in a tropical forest studied in Walter \BBA Contreras (\APACyear1999) and summarized in Figure 5, where $k_{ij}$ accounts for the flow rates between compartments, associated with the mass leaving the $i$th compartment and arriving at the $j$th compartment. Note that each compartment allows interaction from its neighbour compartments in both directions and that the parameters are given for a continuous-time model. Choosing the sampling-period as $T_{s}=0.1$ years, the discrete-time matrices are given by

[TABLE]

with state vector $x_{k}$ $\in$ $\mathbb{R}^{5}$ composed by the amount of nitrogen in each compartment, initial condition $x_{0}\!\!=\!\![\!-\!3.5~{}\!-\!2.52~{}~{}0~{}~{}520~{}~{}26.5]^{T}\!\!+\!\!\mathcal{K}[3.82~{}~{}316~{}~{}1~{}~{}576~{}~{}41\!]^{T}$, where $\mathcal{K}$ depends on the total amount of nitrogen in the system in the beginning and is set as $\mathcal{K}=1.5$, the process noise and observation noise covariance matrices $\tilde{Q}$ $=$ $\sigma_{w}^{2}GG^{T}$, where G=\left[\begin{array}[]{cccc}0.1220&0.1634&0.0249&-0.0383\\ -0.0420&-0.0048&-0.1430&0.0235\\ 0.1640&-0.0317&-0.0057&0.0571\\ -0.1871&-0.0877&0.1697&-0.0098\\ -0.0569&-0.0392&-0.0459&-0.0325\end{array}\right], and $R_{k}=\sigma_{v}^{2}I_{4\times 4}$. The conditions of Lemmas 4.1 and 4.2 hold for (78) such that the trajectory of $x_{k}$ $\in$ $\mathbb{R}^{5}$ lies on the plane (3), whose parameters are assumed to be known and are given by

[TABLE]

that is, mass conservation is verified.

We investigate a 1000-run Monte Carlo simulation testing identification data with $u_{k}=1+\sigma_{u}w_{k}^{u}$ and different noise realizations for $\sigma_{u}$ $=$ $0.1$, $\sigma_{w}$ $=$ $1.0$ and $\sigma_{v}$ $=$ $1.0$ (not shown for brevity), where $w^{u}$ is a zero mean white noise to ensure the persistence of excitation of the input. For each running simulation we employ both LS given by (26) and CLS given by Proposition 4.3. The validation data is simulated with different initial condition $x_{0}\!\!=\!\![72.2~{}~{}381.5~{}~{}101.5~{}~{}1264.0~{}~{}128.0]^{T}$. Figure 6 shows the results regarding the RMSE of the Monte Carlo validation of the obtained models. For all the states the use of auxiliary information in the CLS improves the performance of the estimated models. Specifically, the CLS estimation of the $4$th and $5$th state components always outperform the LS estimation. For in the CLS estimator, $\overline{\textrm{RMSE}}_{4}$ is around $15$ times smaller than the same index for LS. For the $5$th state, CLS improves $\overline{\textrm{RMSE}}_{5}$ by a factor of 3.

6 Concluding Remarks

We address the problem of modeling state-space dynamic systems for which the state vector satisfies an exactly known or an uncertain equality constraint (prior information). We assume that all state components are measured such that least-square methods can be used to estimate the state-space matrices. Both batch and recursive algorithms are considered. By means of the latter, the time-varying case is also addressed.

First, we show how to map the known equality constraint on the state vector on an equality constraint on the parameter matrix to be estimated by the least-square based method. Then, we show how to rewrite the corresponding least squares problem into a vectorized form such that existing equality-constrained least squares methods may be used.

In addition to obtaining state-space matrices that yield an equality-constrained model on the state vector, we observe that the usage of both exactly known and uncertain prior information improves the prediction quality of the model compared to the case in which the equality constraint is not enforced. Such results are consistent with those from Teixeira \BBA Aguirre (\APACyear2011). The algorithms here investigated are also of interest for gray-box subspace identification methods that employ least squares as an internal step; see Trnka \BBA Havlena (\APACyear2009); Alenany \BOthers. (\APACyear2011); Privara \BOthers. (\APACyear2012); Alenany \BBA Shang (\APACyear2013).

Appendix A

Next, we present the proof of Proposition 3.1.

Proof.

Given that (27) is a projector onto the null space of $D$, initialize the RLS equations given by (31)-(33) with the initial values (34)-(35) such that

[TABLE]

For $k=1$, multiplying (31)-(33) by $D$ we have

[TABLE]

By symmetry, from (84)-(85), we verify that $DK_{2}=0$, $D\hat{\theta}_{\textrm{CLS},2}=d$ and $DP_{\textrm{CLS},2}=0$.

Likewise, at time $k+1$, we have $DK_{k+1}=0$, $D\hat{\theta}_{\textrm{CLS},k+1}=d$ and $DP_{\textrm{CLS},k+1}=0$.

Thus, by induction, we have $DK_{k}=0$, $D\hat{\theta}_{\textrm{CLS},k}=d$ and $DP_{\textrm{CLS},k}=0$, $\forall$ $k$, completing the proof. ∎

Appendix B

In the following, we show that the classical LS estimator is biased for the estimation of the matrices of the state-space model (1)-(2).

Proof.

Replace $y_{k}$ in (4), write this result as the linear regression model (21) and apply the vectorization operador (68), then we obtain

[TABLE]

where

[TABLE]

The terms $v_{k}$ underlined are identical. So, the regressor $y_{k}$ between the braces is correlated with $e_{k+1}$. In this case, the LS estimator is biased. ∎