# MixedLp Estimators Variety for Model Order Reduction in Control Oriented System Identification

Authors: Christophe Corbier; Jean-Claude Carmona
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A new family of MLE typeLp estimators for model order reduction in dynamical systems identification is presented in this paper. A family of Lp distributions proposed in this work combines Lp2 (1<p2<2) and Lp1 (0<p1<1) distributions which are quantified by four parameters. The main purpose is to show that these parameters add degrees of freedom (DOF) in the estimation criterion and reduce the estimated model complexity. Convergence consistency properties of the estimator are analysed and the model order reduction is established. Experimental results are presented and discussed on a real vibration complex dynamical system and pseudo-linear models are considered.


---

## Body

## 1. Introduction

The choice of the norm is a fundamental problem in system identification. For some special cases, the least squares (LS) (p=2), the least absolute deviation (LAD) (p=1), and the Chebyshev estimator (p=∞) methods have been proposed by many authors [1, 2]. It is well known that the L1 norm estimators are efficient when the noise follows the Laplace distribution, the L2 norm estimators are efficient when the noise follows the Gaussian distribution, and the L∞ norm estimators are efficient when the noise follows the uniform distribution. Furthermore, it has been shown that the Lp norm estimator is the MLE if noise is a generalized p-Gaussian [3]. Thus, we can always find a suitable exponent p for the Lp norm method. In the reference paper [4], the authors discuss smooth and sensitive Lp norms for prediction error system identification when disturbances are magnitude-bounded. They show that a necessary condition for Lp norm to be statistically robust with respect to the family P~pω- of the ω--contaminated distribution model with support -β,β for some arbitrary β>0 is that its second derivative does not vanish on the support and must be strictly positive with respect to P~pω-. This latter condition is fundamental in our framework by ensuring the parameter estimate variance convergence rate (see Proposition 4.2 and Corollary 4.3 [4]). However, many drawbacks come from the use of the Lp norm. Firstly, p is the unique DOF of the norm. Secondly, the support is a restricted condition to deal with the outliers [5]. Thirdly, the choice of p may prevent robustness to the large data. Finally, most of the time, L∞ norm is used to the model reduction for robust control [6] with the restricted case where the model parameters number is limited. The Lp norms are widely used in many domains [7, 8]. Even though the formal framework presents some difficulties, previous works show the interest to use these estimators. In system identification, Chen et al. in [9] propose an ARMA robust system identification using a generalized Lp norm estimation algorithm. The authors use this norm to identify a system with a non-Gaussian noise, since the classical LSE is related to the situation when the noise distribution is white-Gaussian. Likewise, these authors in [10] propose the parameter estimation of linear systems with input-output noisy data, supposed to be corrupted by measurement noise unknown distributions.Here, we deal with the problem of theLp estimation in system identification in the presence of atypical data (outliers) in the vibration dynamical system output signal, in order to reduce the model order for the controller design. In such a vibration dynamical system,natural outliers occur in the output signal, increasing considerably the model order in the estimation procedure.To avoid this problem, a new variety of MLE typeLp estimators are considered. We define a mixed Lp distribution with multiple DOF containing 1<p2<2 and 0<p1<1 [11], parametrized by a threshold namedscaling factor, denoted by η~. To ensure both the parameter estimate variance convergence and the model order reduction and apply conditions in [4] related to the second derivative of the Lp norm on the support with respect to P~pω-, we show that there exists a tradeoff between p2, p1, and η~. This reduction is efficient only with η~<1 and p2<2, involving p1<1.Model reduction has been subject to considerable interest in many domains. In [12], a method of the model reduction of the nonlinear complex Ginzburg-Landau equation is presented. In [13] the authors propose dimension reduction and dynamics of a spiking neural network model for decision making under neuromodulation. The model order reduction is a fundamental step, espacially in control design. To estimate a low order model of a system, several possibilities exist. The most obvious one is to directly estimate a lower order model. As known from, for example, Ljung [14], the prediction/output error estimate automatically gives models that are L2 approximations of thetrue system in frequency-weight norm, determined by the input spectrum and noise model. See [14, chapter 8] for more details. A second possibility is to estimate a high order model which is then subjected to model reduction to the desired order. Indeed, some contributions that take into account that the high order model is obtained through an identification experiment when performing model reduction are Porat and Friedlander [15], Porat [16], Söderström et al. [17], Stoica and Söderström [18], Zhu and Backx [19, chapter 7] and Tjärnström and Ljung [20], and Tjärnström [21]. Porat and Friedlander study an ARMA parameter estimation via covariance estimates. Recently, in time domain, [22] presents a relevant paper on a new model order reduction algorithm based on general orthogonal polynomials. The following contributors deal with models having input signals. Söderström et al. [17] focus on model structures that can be embedded in larger structures which are easy to estimate, such as ARX structures. After estimating the high order structure, they reduce the estimate to the low order structure in a weighted nonlinear least-squares sense. The method is calledindirect prediction error method. In [20], the authors focus on L2 model reduction. They study FIR and OE models. The former is estimated directly from data and the latter is computed by reducing a high order model, by L2 model reduction. For OE models, they show that the reduced model has the same variance as the directly estimated one, if the reduced model class used contains the true system. Recently, in [23], the authors illustrate procedures to identify a state-space representation of a lossless or dissipative system from a given noise-free trajectory. The idea is to perform model reduction by obtaining a balanced realization directly from data and truncating it to obtain a reduced order model. However, these different model order reduction methods are performed after estimation phases with high model orders. On the other hand, this reduction is made in favorable conditions, where no outliers are present in the dataset and the prediction errors. In practice, we all know that real conditions lead to model complexity strongly increased. The lack of robustness of these methods damages the estimation and the model order. Different works have been performed, using robust estimators. A first alternative is the use of the L1 norm [24–26], as a possibility to reduce the model complexity to identify a complex dynamical system, for example, an acoustic duct for active noise control (ANC). The authors propose relevant OE (output error) L1 models with complexity dML1 equal to 21. The second one is to use a symmetrix convex Huber’s function [27]. Recently, in [28, 29], the identification of the same plant, using L2-L1 estimators with low values of the scaling factor in Huber’s function, proposes OE-L2-L1 models with complexity dMH equal to 17.Here, we will tackle theLp2-Lp1 distribution approach and, from this, we will propose a Lp parameterized robust estimation criterion (PREC) using a symmetric convex function, offering to the user moreflexibility in terms of an extra DOF, in order to improve the estimation balance robustness/performances. We will show that the number of DOF in this function facilitates the model order reduction. We will see that the reduced order models are provided by the Lp2-Lp1 estimation itself, where classical approaches separate the phases, estimation and then reduction. Indeed, we will demonstrate that the Lp model order thus obtained remains lower than the corresponding L2-L1 estimate. Finally, experimental results on the acoustic duct will confirm these results.We start in Section2 by introducing the mathematical background related to definition and properties of the Lp estimator. We develop the consistency and convergence properties of this estimator. The application to the system identification is presented in Section 3. Section 4 proposes the main results concerning the model order reduction. In Section 5, the improvement of the model order reduction with respect to the number of DOF is presented. Experimental results based on a vibration dynamical system are presented and discussed in Section 6. Conclusions and perspectives are drawn in Section 7.


## 2. Definition and Properties of theLp Estimators

The new family ofLp-distributions proposed in this paper combines Lp2 and Lp1 distributions, quantified by four parameters, p2 (1<p2<2), p1 (0<p1<1), ϕ (ϕ>0), and η~ (0<η~<1). Hence, it can be viewed as a generalization of theredescending mixed distributions, which corresponds to the M-estimates p2=2 and p1=1 (see [27, chapter 4, page 84-85]). A formal definition of a probability density function (pdf) of a random variable is as follows.Definition 1.
Consider a functionfp:R→R+*, X↦fpX such that(1)fpX=fp2X=Ce-Xp2/2ϕp2forX≤η~fp1X=Ce-η~Xp1/ϕp2+η~2/2ϕp2forX>η~,where C=1/2K1+K2 with(2)K1=eη~2/2ϕp2ϕp2/p1p1η~1/p1Γ1p1,η~p1+1ϕp2forX>η~,K2=ϕ21/p2p2Γ1p2-Γ1p2,η~p22ϕp2forX≤η~,where Γa and Γa,X are, respectively, the complete and incomplete Euler’s gamma functions. We can easily verify that, for all X∈R, fpX≥0 and ∫R‍fpXdX=1, which ensure that fp is a probability density. The parameter η~ is defined as a threshold namedscaling factor. Its value depends on the degree of the contamination of the random variables X. Moreover, a low value will ensure a good robustness and a high value a good efficiency of the estimation.The pdffpX belongs to the maximum likelihood estimators (MLE) class. Then, there exists a symmetric function variety ρη~pX such that ρη~pX~-log⁡fpX.Definition 2.
Let one define a new estimator variety as aLp estimator, denoted by θ^Np. Let S,S,P be a probability space and let Xtt=1N be a sequence ofi.i.d.r.v’s with values in S. Let Θ be a Borel subset in Rd. Let ρη~p:S×Θ→R be a symmetric function such that ρη~p•θ:S→R is measurable for each θ∈Θ. The Lp estimator θ^Np is defined by a minimum of the form(3)N-1∑t=1N‍ρη~pXt(θ^Np)=inf⁡θ∈ΘN-1∑t=1N‍ρη~pXt(θ),where(4)ρη~pX=Xp22forX≤η~η~Xp1-η~22forX>η~with 1<p2<2 and 0<p1<1 for all η~>0.In robust statistics, continuity conditions are necessary to ensure good estimations [27, chapter 2, page 24-25]. Therefore, these conditions must be applied for both ρη~p and its first derivative. This latter must be a bounded continuous function. From these considerations, we define the dimensional continuity conditions (DCC) as follows:(5)DCC:ρη~p2±η~=ρη~p1±η~∂∂Xρη~p2X±η~=∂∂Xρη~p1X±η~,where ρη~p2X=Xp2/2 and ρη~p1X=η~Xp1-η~2/2.The dimensionalp1 is given from DCC and, after straightforward calculations, we obtain (6)p1=p21+η~2-p2.Equation (6) allows to show that, for the conditions η~<1 and p2<2, we get p1<1. We now define the gross error model (GEM) in the mixed Lp framework as a generalization of the GEM given in [27, chapter 1, page 12].Definition 3.
LetFp2 be a fixed probability distribution having a twice differentiable density fp2, such that -log⁡fp2 is convex on the convex support of Fp2. Let ω->0 be the level of contamination and let P~pω- be the set of all probability distributions arising from Fp2 through the ω--contamination model(7)P~pω-=Fp∣Fp=1-ω-Fp2+ω-H,H∈M.Here, M is the set of all probability measures on the real line. The corresponding pdf of Fp, contained in P~p, is given by (1), and it puts all contamination outside -η~,η~. A suitable value of η~ should efficiently reduce the propagation of outlier in the prediction errors. This is done by ρη~p1. Accordingly, the Lp1 treatment of the estimated residuals is focused on the outlier occurrence, so the following errors rapidly decrease and they are advantageously treated by ρη~p2.Now, our goal is to provide the relation between the level of contaminationω- and the scaling factor η~. In practice, the robustness is defined by a suitable value of η~ that fixes the ω--contamination model.Theorem 4.
If one adjusts the scaling factorη~ in ρη~p, such that C~=1-ω-C, from the condition ∫-∞∞‍C~e-ρη~pX/σp2dX=1, ω- and η~ are linked by(8)A+B=11-ω-with (9)A=ϕ21/p2p2K1+K2Γ1p2-Γ1p2,η~p22ϕp2B=eη~2/2ϕp2ϕP2/p1K1+K2p1η~1/p1Γ1p1,η~p1+1ϕp2,where Γa,X=∫X∞‍Za-1e-ZdZ, a>0, is incomplete Euler’s gamma function.See AppendixA for the proof. Notice that in the limit case, where (p2,p1)=(2,1) (M-estimates), we obtain, for η=kσ  ϕ=σ, A=1-2Φ-k and B=2φk/k. See [27, chapter 4, page 84].Figure1 shows an example of the level of contamination ω- as a function of η~, p1, parameterized by p2, with ϕ=σ=1. A large value of the scaling factor means that the probability distribution of the prediction errors is weakly disturbed, corresponding to ω-<0.1 (see Figure 1). Conversely, a small value corresponds to a pdf more disturbed. For p2≤1.7, the curves have a minimum, meaning that a particular value of the scaling factor renders the corresponding pdf weakly contaminated.Figure 1
Level of contaminationω- is plotted as a function of η~, p1, parameterized by p2 and ϕ=σ=1.Definition 5.
Letρη~pX be a symmetric function given by(10)ρη~pX=Xp22IX≤η~+η~Xp1-η~22IX>η~,where IY is the unit function with the condition Y. Let us introduce two index subsets νp2 and νp1, respectively, related to the Lp2 and Lp1 contributions, defined by νp2θ=t:εtθ≤η~ and νp1θ=t:εt(θ)>η~. From this, the Lp-PREC to be minimized is given by(11)gNpθ=1N∑t∈νp2(θ)‍εt(θ)p22+η~N∑t∈νp1θ‍εtθp1-η~st2θ2,where εtθ is the estimated modeling error and st(θ) is the sign function such that st(θ)=0 for εt(θ)≤η~, st(θ)=-1 for εt(θ)<-η~, and st(θ)=1 for εt(θ)>η~.In order to prove the consistency and convergence of the robust estimatorθ^Np, we need the following lemma related to the moment of order r of the random variable U.Lemma 6.
For the considered pdffp, the moments of order r(r∈N) of a random variable U are(12)EUr=I1r+I2rwith(13)I1r=ϕr+12r+1/p2p2hhhhhh×Γr+1p2-Γr+1p2,η~p22ϕp2ϕr+12r+1/p2p2×K1+K2-1,I2r=eη~2/2ϕp2ϕp2/p1r+1p1η~r+1/p1hhhhLh×Γr+1p1,η~p1+1ϕp2×K1+K2-1.See AppendixB for the proof.Theorem 7.
LetΘ be a compact subset of Rd and let θ^Np be the estimator obtained as a solution of the problem(14)inf⁡θ∈ΘN-1∑t=1N‍-ln⁡fpεt(θ).Let θ* be the global minimum obtained as a solution of the problem(15)inf⁡θ∈ΘN-1∑t=1N‍-Eln⁡fpεt(θ),asN⟶∞,and then θ^Np is a consistent estimate of θ*; that is,(16)θ^Np⟶θ*,w.p.1.In order to do so, we use the (MLE) consistency result by Newey and McFadden [30] and show that all the assumptions of their theorem hold.Theorem 8.
Suppose thatzii=1,2,… are iid rv’s with pdf fzi∣θ* and(i)
θ≠θ*⇒fzi∣θ≠fzi∣θ*;(ii)
θ*∈Θ, which is compact;(iii)
ln⁡fzi∣θ is continuous at each θ∈Θ w.p.1;(iv)
Esup⁡θ∈Θln⁡fzi∣θ<∞.Then θ^N→probθ*.See AppendixC for the proof.Remark 9.
Before presenting the theorem related to the asymptotic distribution of theLp estimator, let us notice that the standard asymptotic normality results for MLE requires that the PREC be twice continuously differentiable, which is not the case here. There exist, however, asymptotic normality results for nonsmooth functions and we will hereafter use the one proposed by Newey and McFadden [30] and Andrews [31]. The basic insight of their approaches is that the smoothness condition of the PREC, gNp(θ), can be replaced by a smoothness of its limit, which in the standard maximum likelohood case corresponds to the expectation -Eln⁡fpεtθ=gpθ, with the requirement that certain remainder terms are small. Hence, the standard differentiability assumption is replaced by a stochastic differentiability condition, which can then be used to show that the MLE θ^Np is consistent and asymptotically normal. Let us consider the Ψ-function as the derivative with respect to θ of ρη~p denoted by Ψtε;θ. If this function is differentiable in θ, one can establish the asymptotic normality of θ^Np by expanding Nθ^Np-θ* about θ* using element by element mean value expansions. This is the standard way of establishing asymptotic normality of the estimator. In a variety of applications, however, Ψtε;θ is not differentiable in θ, or not even continuous, due to the appearence of a sign function. In such a case, one can still establish asymptotic normality of the estimator provided EΨtε;θ is differentiable in θ. Since the expectation operator is a smoothing operator, EΨtε;θ is often differentiable in θ even though Ψtε;θ is not.
Let us consider the following theorem establishing the asymptotic distribution of the estimator.Theorem 10.
Letε1θ^Np,…,εNθ^Np be iid rv’s from the pdf fp with an unknown parameter θ*, θ*∈B∘ with B a compactness and B∘ interior of B. Then the MLE θ^Np of θ* is consistent and asymptotically normal:(17)Nθ^Np-θ*⟶dN0,Pθ*,where Pθ* is the asymptotic covariance matrix.In order to do so, we use Theorem 7.1 given by Newey and McFadden [30] and show that all the assumptions of their theorem hold.Theorem 11.
SupposegNp(θ^Np)≥sup⁡θ∈DMgNpθ-opN-1, θ^Np→probθ*, and(i)
gp(θ) is maximized on DM at θ*;(ii)
θ* is an interior point of DM;(iii)
gp(θ) is twice differentiable at θ* with nonsingular second derivative g-′′p(θ);(iv)
NE(∂gNpθ/∂θ)θ^Np→dN0,Qθ*;(v)
for anyδN→0, sup⁡θ^Np-θ*≤δNR^Nθ^Np/(1+Nθ^Np-θ*)→prob0 with the remainder(18)R^Nθ=NgNp(θ)-gNp(θ*)-∂gNp(θ)/∂θθ*θ-θ*θ-θ*+-gpθ+gpθ*θ-θ*,and then Nθ^Np-θ*→dN0,Pθ*.The proof is given in AppendixD.


## 3. Application to the System Identification

Having defined the mixedLp framework, our purpose is now to apply this approach to the system identification by introducing the Lp-PREC and its first and second derivatives, to express the asymptotic Lp covariance matrix of the Lp estimator. Before this, let us recall the prediction error context.


### 3.1. Prediction Error Approach

Throughout the paper, we denote the input signal byut, the output signal by yt, and the disturbance signal by et as a random variable sequence with zero mean and variances λ, and N is the total number of measured data. We assume that yt is generated according to(19)yt=G0qut+H0qet,where G0q is a linear time-invariant system, usually referred to as thetrue system and H0q is the noise filter. In the case of OE models, the model set is parameterized by a d-dimensional real-valued parameter vector θ; that is,(20)yt=Gq,θut+Hq,θet=Bq,θFq,θut+et,with Fq,θ=1+∑i=1nFfiq-i, Bq,θ=∑i=1nBbiq-i, and θ=b1…bnBf1…fnFT. Moreover, q-1 is the lag operator such that q-lyt=yt-l, l∈N. The prediction error (also named estimated residuals) is εtθ=yt-y^tθ with y^tθ=φtTθθ the prediction model, φtθ=[ut-1…ut-nB-y^t-1θ…-y^t-nFθ]T the regressor, and ψtθ=∂/∂θy^tθ the gradient of y^tθ [14].


### 3.2. Limit Hessian andQω-p-Matrix

The main purpose of this section is to express the asymptotic covariance matrix of the estimatorθ^Np for prediction error estimates in case, essentially, a true description of the system is available in the model set. The condition then becomes S∈M, where M is a model structure, a mapping from a parameter space to a set of models. Let us consider the convergence domain defined by Dc=argmin⁡θ∈DMgp(θ)=θ∣θ∈DM,gp(θ)=min⁡θ~∈DMgp(θ~). We put ourselves in the case where the global minimum is equal to the true parameter vector θ0; in this case, Dc=θ*=θ0. To establish this asymptotic covariance matrix, we need the limit Hessian and the Qω-p-matrix. Using Theorem 10, the covariance matrix is given by(21)covθ^NLp~gp-ω-′′θ0-1Qω-pθ0gp-ω-′′θ0-1N=Pθ0Nwith the limit Hessian and the Qω-p-matrix yielded, respectively, by(22)gp-ω-′′θ0=p2p2-12mH,20E-Φνp2,t0+η~p1p1-1mH,10E-Φνp1,t0,Qω-pθ0=p224mH,30E-Φνp2,t0+η~2p12mH,40E-Φνp1,t0,where E-Φνpi,t0=E-ψνpi,t0ψνpi,t0T, i=1,2. See Appendix E for the proof.Remark 12.
It is easy to show that the second derivative ofgNpθ is given by gNp′′θ=1/N∑t=1N‍(∂2/∂θ∂θT)ρη~pεtθ with(23)∂2∂θ∂θTρη~pεtθ=p2p2-12mH,2ψνp2,tθψνp2,tTθ+η~p1p1-1mH,1ψνp1,tθψνp1,tTθ.Since p1<1, the second term on the right-hand side of (23) is negative. Accordingly, there exists a limit value of p1 denoted by p1L such that (23) is positive and then verifies Proposition 4.2 and Corollary 4.3 in [4].


### 3.3.Lp Robust Final Prediction Error

Similar to the RFPE in Huber’s framework [29], we now establish a quality measure for a given model. Let m^Np=M~θ^Np be the robust model to be evaluated. It is estimated within the model structure M~. Let us denote J-θ^NpM~;ω-, this quality measure that depends on the actual data record ZN. We classically have(24)J-θ^NpM~;ω-=Eg-ω-pθ^Np,where g-ω-pθ is the limit criterion. Let θ0 be the argument of the minimum of this limit criterion and let us expand g-ω-pθ around θ0, and we obtain(25)g-ω-pθ^Np≈g-ω-pθ0+12θ^Np-θ0Tg′′-ω-pζ^Npθ^Np-θ0.Similarly, since gNp′θ^Np=0,(26)gNpθ^Np≈gNp′θ0-12θ^Np-θ0TgN,ω-′′pζ^Npθ^Np-θ0.From the Hessian uniform convergence theorem and the expectations of (25) and (26), we obtain for the Lp-RFPE(27)J-θ^NpM~;ω-≈EgNpθ^Np+1Ntr⁡Qω-pθ0g′′-ω-pθ0-1.If a single run is taken into account, then we can classically replace the expectation of gNpθ^Np with only one observation available. Therefore,(28)J-θ^NpM~;ω-≈gNpθ^Np+1Ntr⁡Qω-pθ0g′′-ω-pθ0-1.Let us denote Sω-pθ0=Qω-pθ0g′′-ω-pθ0-1 and Tνp1,νp2pθ0=E-Φνp1,t0E-Φνp2,t0-1. We obtain(29)Sω-pθ0=V0I-Xω-0Tνp1,νp2pθ0-Yω-0Tνp1,νp2pθ02,where I is the identity matrix. Now, in order to completely define the expression of the Lp-RFPE, we need the following result due to the Schur complement.Lemma 13.
LetM be a symmetric matrix of the form M=ABBTC, where A, B, and C are real and positive definite matrices. If C is invertible, then the following property holds:(30)(i)M>0iffC>0,A-BC-1BT>0⟹C>A-1.In the mixedLp context, the estimation procedure is mainly a Lp2 estimation since the impulsional effect of the outlier is rapidly treated by the Lp1. Accordingly, cardνp2>cardνp1. Moreover, E-Φνpi,t0=E-ψνpi,t0ψνpi,t0T, i=1,2, are positive definite matrices. Using Lemma 13 with B=Id, A=E-Φνp1,t0-1, and C=E-Φνp2,t0, we have(31)E-Φνp2,t0>E-Φνp1,t0⟹Tνp1,νp2pθ0=α-Id,where d is the model complexity and 0<α-<1. Therefore, Tνp1,νp2pθ02=α-2Id and thetrace of Sω-pθ0 becomes(32)dMpV0-α-Xω-0-α-2Yω-0with dMp being the model complexity. The Lp-RFPE becomes(33)J-θ^NpM~;ω-≈gNpθ^Np+dMpNV0-α-Xω-0-α-2Yω-0.


## 3.1. Prediction Error Approach

Throughout the paper, we denote the input signal byut, the output signal by yt, and the disturbance signal by et as a random variable sequence with zero mean and variances λ, and N is the total number of measured data. We assume that yt is generated according to(19)yt=G0qut+H0qet,where G0q is a linear time-invariant system, usually referred to as thetrue system and H0q is the noise filter. In the case of OE models, the model set is parameterized by a d-dimensional real-valued parameter vector θ; that is,(20)yt=Gq,θut+Hq,θet=Bq,θFq,θut+et,with Fq,θ=1+∑i=1nFfiq-i, Bq,θ=∑i=1nBbiq-i, and θ=b1…bnBf1…fnFT. Moreover, q-1 is the lag operator such that q-lyt=yt-l, l∈N. The prediction error (also named estimated residuals) is εtθ=yt-y^tθ with y^tθ=φtTθθ the prediction model, φtθ=[ut-1…ut-nB-y^t-1θ…-y^t-nFθ]T the regressor, and ψtθ=∂/∂θy^tθ the gradient of y^tθ [14].


## 3.2. Limit Hessian andQω-p-Matrix

The main purpose of this section is to express the asymptotic covariance matrix of the estimatorθ^Np for prediction error estimates in case, essentially, a true description of the system is available in the model set. The condition then becomes S∈M, where M is a model structure, a mapping from a parameter space to a set of models. Let us consider the convergence domain defined by Dc=argmin⁡θ∈DMgp(θ)=θ∣θ∈DM,gp(θ)=min⁡θ~∈DMgp(θ~). We put ourselves in the case where the global minimum is equal to the true parameter vector θ0; in this case, Dc=θ*=θ0. To establish this asymptotic covariance matrix, we need the limit Hessian and the Qω-p-matrix. Using Theorem 10, the covariance matrix is given by(21)covθ^NLp~gp-ω-′′θ0-1Qω-pθ0gp-ω-′′θ0-1N=Pθ0Nwith the limit Hessian and the Qω-p-matrix yielded, respectively, by(22)gp-ω-′′θ0=p2p2-12mH,20E-Φνp2,t0+η~p1p1-1mH,10E-Φνp1,t0,Qω-pθ0=p224mH,30E-Φνp2,t0+η~2p12mH,40E-Φνp1,t0,where E-Φνpi,t0=E-ψνpi,t0ψνpi,t0T, i=1,2. See Appendix E for the proof.Remark 12.
It is easy to show that the second derivative ofgNpθ is given by gNp′′θ=1/N∑t=1N‍(∂2/∂θ∂θT)ρη~pεtθ with(23)∂2∂θ∂θTρη~pεtθ=p2p2-12mH,2ψνp2,tθψνp2,tTθ+η~p1p1-1mH,1ψνp1,tθψνp1,tTθ.Since p1<1, the second term on the right-hand side of (23) is negative. Accordingly, there exists a limit value of p1 denoted by p1L such that (23) is positive and then verifies Proposition 4.2 and Corollary 4.3 in [4].


## 3.3.Lp Robust Final Prediction Error

Similar to the RFPE in Huber’s framework [29], we now establish a quality measure for a given model. Let m^Np=M~θ^Np be the robust model to be evaluated. It is estimated within the model structure M~. Let us denote J-θ^NpM~;ω-, this quality measure that depends on the actual data record ZN. We classically have(24)J-θ^NpM~;ω-=Eg-ω-pθ^Np,where g-ω-pθ is the limit criterion. Let θ0 be the argument of the minimum of this limit criterion and let us expand g-ω-pθ around θ0, and we obtain(25)g-ω-pθ^Np≈g-ω-pθ0+12θ^Np-θ0Tg′′-ω-pζ^Npθ^Np-θ0.Similarly, since gNp′θ^Np=0,(26)gNpθ^Np≈gNp′θ0-12θ^Np-θ0TgN,ω-′′pζ^Npθ^Np-θ0.From the Hessian uniform convergence theorem and the expectations of (25) and (26), we obtain for the Lp-RFPE(27)J-θ^NpM~;ω-≈EgNpθ^Np+1Ntr⁡Qω-pθ0g′′-ω-pθ0-1.If a single run is taken into account, then we can classically replace the expectation of gNpθ^Np with only one observation available. Therefore,(28)J-θ^NpM~;ω-≈gNpθ^Np+1Ntr⁡Qω-pθ0g′′-ω-pθ0-1.Let us denote Sω-pθ0=Qω-pθ0g′′-ω-pθ0-1 and Tνp1,νp2pθ0=E-Φνp1,t0E-Φνp2,t0-1. We obtain(29)Sω-pθ0=V0I-Xω-0Tνp1,νp2pθ0-Yω-0Tνp1,νp2pθ02,where I is the identity matrix. Now, in order to completely define the expression of the Lp-RFPE, we need the following result due to the Schur complement.Lemma 13.
LetM be a symmetric matrix of the form M=ABBTC, where A, B, and C are real and positive definite matrices. If C is invertible, then the following property holds:(30)(i)M>0iffC>0,A-BC-1BT>0⟹C>A-1.In the mixedLp context, the estimation procedure is mainly a Lp2 estimation since the impulsional effect of the outlier is rapidly treated by the Lp1. Accordingly, cardνp2>cardνp1. Moreover, E-Φνpi,t0=E-ψνpi,t0ψνpi,t0T, i=1,2, are positive definite matrices. Using Lemma 13 with B=Id, A=E-Φνp1,t0-1, and C=E-Φνp2,t0, we have(31)E-Φνp2,t0>E-Φνp1,t0⟹Tνp1,νp2pθ0=α-Id,where d is the model complexity and 0<α-<1. Therefore, Tνp1,νp2pθ02=α-2Id and thetrace of Sω-pθ0 becomes(32)dMpV0-α-Xω-0-α-2Yω-0with dMp being the model complexity. The Lp-RFPE becomes(33)J-θ^NpM~;ω-≈gNpθ^Np+dMpNV0-α-Xω-0-α-2Yω-0.


## 4. Model Order Reduction

Let two model classesM and M~ be given, respectively, related to Huber’s and Lp estimators. To avoid some drawbacks with the notation of the parameter vectors, model classes M and M~ are parameterized by θ and θ~, respectively. Recall that the PREC in Huber’s approach is defined as follows:(34)WNθ=1N∑t∈ν2θ‍εt2θ2+ηN∑t∈ν1θ‍εtθ-ηst2θ2.See [29, chapter 4, page 50]. Let us consider the following conditions C1:(1)
η<1 and ω<0.4;(2)
η~<1 and ω-<0.4;(3)
1.8≤p2<2 with p1<1;(4)
ϕ<2.Theorem 14.
Assume that the conditions C1 are verified. LetΩi=sup⁡t,θ,θ~εtθ,εt(θ~) be the innovation outlier; then(35)sup⁡θ~gNpθ~<sup⁡θWNθ.See the proof in AppendixF. Let θ^NH be Huber’s estimator. Using J-θ^NHM;ω in [29, chapter 7, page 125] and Lemma 13, since α<1 and ω<0.4, then (36)sup⁡θ^NHJ-θ^NHM;ω=1+Ωi+2λ^ωdMHN.From (33) and conditions C1, since Xω-0<1/2 and Yω-0<1/2, we obtain(37)sup⁡θ^NLpJ-θ^NLpM~;ω-=1+Ωip1+dMpN1+V0.From the curves of V0 and using the conditions C1, then(38)V0≈2p225p2p2-12η~+1.Let us denote Rd=2λ^ω/(1+V0), where ω=gη and η~,p2 are the DOF, respectively, related to Huber’s and Lp symmetric functions.Theorem 15.
Ifsup⁡θ^NpJ-θ^NpM~;ω-<sup⁡θ^NHJ-θ^NHM;ω and Rd<1, then(39)dMp<2λ^ω1+V0dMH⟹dMp<dMH.The proof is given in AppendixG. This tuning offers to the user the necessary DOF to estimate and validate the reduced order models, as described in the following section.


## 5. Effect of DOF in the Symmetric Function to the Model Order Reduction

Here, we focus on a new vision on the link between the model reduction objective and its link to the estimation chosen, in particular the number of DOF of this symmetric function. LetDf be the number of DOF. We present in Table 1 different symmetric function related to the MLE context and the corresponding Df with their associated tuning parameters.Table 1
DOF in the symmetric function.
Estimator
Function
D
f
DOF
LSE
X
2
2
∑
454
541
D
F
21454
f
d
0
0
LSAD
|
X
|
0
0
Huber
{
X
2
2
for
|
X
|
≤
η
η
|
X
|
-
η
2
2
for
|
X
|
>
η
1
η
Mixed-Lp
{
|
X
|
p
2
2
for
|
X
|
≤
η
~
η
~
|
X
|
p
1
-
η
~
2
2
for
|
X
|
>
η
~
2
η
~, p2First, the reason for which the LSE and LSAD estimators do not allow the model order reduction seems to be due to the number of DOF; here,Df=0. From the comparison between Akaike’s FPE and the L1-FPE [24], we obtain dML1<λ^NdML2. The condition dML1<dML2 is ensured only by the value of the variance λ^N of the L2 estimation. Sometimes, this one is obtained with difficulty, especially in the presence of some outliers in the prediction errors. Since Df=0, the model complexity reduction remains difficult. From our point of view, the introduction of DOF in the estimation symmetric function facilitates the model complexity reduction. For example, the DOF in Huber’s symmetric function may ensure both the robustness and the model order reduction. Therefore, comparing (36) and the L1-FPE, we have dMH<(1/2λ^(ω))dML1. Consequently, the model order reduction may be obtained by tuning η in Huber’s function to have the condition 1/2λ^(ω)<1. By extention, with three DOF, the condition Rd<1 in (G.2) may be obtained by tuning ω (through η), η~, and p2. Thus, theflexibility of the Lp symmetric function ensures the model order reduction. Table 2 summarizes these main results related to the effect of the norm’s DOF on the model order reduction.Table 2
Effect of the DOF in the model order reduction.
Dimensional
D
f
Condition
DOF
d
M
L
1
<
d
M
L
2
0
λ
^
N
<
1
0
d
M
H
<
d
M
L
1
1
1
2
λ
^
(
ω
)
<
1,  ω=g(η)
η
d
M
p
<
d
M
H
3
R
d
<
1
η, η~, p2


## 6. Experimental Results

### 6.1. Presentation

This section presents different results related to the identification of a real system, namely, an acoustic duct plant, in order to estimate and validate reduced order models, usingLp estimators. We compare the Lp estimation results with those given by Huber’s [28, 29], LSAD, LSE [24], and L∞ estimators (θ^N∞=min⁡θmax⁡tεtθ). As shown in Figure 3, we have to find discrete-time dynamical models with reduced orders of the plant including the secondary loudspeaker, its amplifier, the acoustic secondary path, the measurement microphone, and its amplifier/antialiasing filter. For more details, the reader could refer to [24]. The identification has been realized from a PRBS (pseudrandom binary sequence) as excitation input signal, with a length L=210-1 and level ±3V, approximating a white noise upon the frequency range 0;1000Hz. The sampling period is Ts=500 μs and the number of data is N=1024. The pure delay d of the duct is given as an integer of the sampling period. Here, d=7. Figure 2 shows the output signal given by the microphone during the recording phase of the data set ZN of the process. An estimations campaign has been conducted to get the value of p1L ensuring Proposition 4.2 and Corollary 4.3 in [4]. p1 was varied over the range [0.5;1[ with an incremental step of 0.01. The convergence is ensured for p1≥0.88. We deduce that p1L=0.88. In the sequel, for the Lp estimation context, 1.8≤p2≤2, η~<1, and 0.88≤p1<1, and we will denote all the models by (α,β)-(nF,nB).Figure 2
Output signal given by the microphone during the recording phase of the data setZN. It clearly shows some outliers with distinguished levels.Figure 3
Experimental setup.


### 6.2. Reduced Order Models

Table3 presents the main results from L∞, LSE, LSAD, M, and Lp estimators. The model order reduction is clearly shown and we can notice that the complexity reduction between LSE and the low dimension of the new mixed Lp estimator is 40%. In Table 4, the parameters of the (2,1)-(9,8) and (1.8,0.939)-(7,7) models are given. Figures 4 and 5 show the frequency responses of the (2,1)-(9,8) and (1.8,0.939)-(7,7) models, respectively, versus the spectral estimation of the process. Even though the Lp model presents a lower complexity than Huber’s model, its frequency response especially in low frequencies particularly interesting in controller design is good. Remark that, for L∞ estimator, the fit is low and the model order increases considerably with dM∞=42.Table 3
Main results fromL∞, LSE, LSAD, M, and Lp estimators for model order reduction.
Estimator
Model
Fit (%)
k
η
ω
λ
^
(
ω
)
η
~
ω
-
R
d
d
M
L
∞
(31, 11)
12
—
—
—
—
—
—
—
42
LSE
(2, 0)-(14, 9)
48
—
—
—
—
—
—
—
23
LSAD
(0, 1)-(14, 7)
72
—
—
—
—
—
—
—
21
M
(2, 1)-(9, 8)
69.45
0.75
0.6159
0.22
0.2819
—
—
—
17
L
p
(1.9, 0.9924)-(8, 7)
67.95
—
—
—
—
0.71
0.265
0.185
15
(1.8, 0.939)-(7, 7)
65.52
—
—
—
—
0.64
0.351
0.182
14Table 4
Parameters from the Huber estimation(2,1)-(9,8) (bnH, fnH) and from the Lp estimation (1.8, 0.939)-(7, 7) (bnp, fnp).
n
1
2
3
4
5
6
7
8
9
b
n
H
-
0.0856
-
0.1642
0.1062
0.2241
0.1460
0.0440
-
0.1929
-
0.0707
—
f
n
H
0.4241
-
0.5992
-
0.4307
-
0.0207
0.2362
0.3232
-
0.1615
0.0462
0.0744
b
n
p
-
0.0927
-
0.1415
0.1829
0.1553
0.0447
0.0412
-
0.1857
—
—
f
n
p
-
0.0699
-
0.7319
0.0578
0.0387
-
0.0021
0.3249
-
0.2332
—
—Figure 4
Frequency response of the(2,1)-(9,8) model (solid line) versus spectral estimation of the process (dashed line).Figure 5
Frequency response of the(1.8,0.939)-(7,7) model (solid line) versus spectral estimation of the process (dashed line).


## 6.1. Presentation

This section presents different results related to the identification of a real system, namely, an acoustic duct plant, in order to estimate and validate reduced order models, usingLp estimators. We compare the Lp estimation results with those given by Huber’s [28, 29], LSAD, LSE [24], and L∞ estimators (θ^N∞=min⁡θmax⁡tεtθ). As shown in Figure 3, we have to find discrete-time dynamical models with reduced orders of the plant including the secondary loudspeaker, its amplifier, the acoustic secondary path, the measurement microphone, and its amplifier/antialiasing filter. For more details, the reader could refer to [24]. The identification has been realized from a PRBS (pseudrandom binary sequence) as excitation input signal, with a length L=210-1 and level ±3V, approximating a white noise upon the frequency range 0;1000Hz. The sampling period is Ts=500 μs and the number of data is N=1024. The pure delay d of the duct is given as an integer of the sampling period. Here, d=7. Figure 2 shows the output signal given by the microphone during the recording phase of the data set ZN of the process. An estimations campaign has been conducted to get the value of p1L ensuring Proposition 4.2 and Corollary 4.3 in [4]. p1 was varied over the range [0.5;1[ with an incremental step of 0.01. The convergence is ensured for p1≥0.88. We deduce that p1L=0.88. In the sequel, for the Lp estimation context, 1.8≤p2≤2, η~<1, and 0.88≤p1<1, and we will denote all the models by (α,β)-(nF,nB).Figure 2
Output signal given by the microphone during the recording phase of the data setZN. It clearly shows some outliers with distinguished levels.Figure 3
Experimental setup.


## 6.2. Reduced Order Models

Table3 presents the main results from L∞, LSE, LSAD, M, and Lp estimators. The model order reduction is clearly shown and we can notice that the complexity reduction between LSE and the low dimension of the new mixed Lp estimator is 40%. In Table 4, the parameters of the (2,1)-(9,8) and (1.8,0.939)-(7,7) models are given. Figures 4 and 5 show the frequency responses of the (2,1)-(9,8) and (1.8,0.939)-(7,7) models, respectively, versus the spectral estimation of the process. Even though the Lp model presents a lower complexity than Huber’s model, its frequency response especially in low frequencies particularly interesting in controller design is good. Remark that, for L∞ estimator, the fit is low and the model order increases considerably with dM∞=42.Table 3
Main results fromL∞, LSE, LSAD, M, and Lp estimators for model order reduction.
Estimator
Model
Fit (%)
k
η
ω
λ
^
(
ω
)
η
~
ω
-
R
d
d
M
L
∞
(31, 11)
12
—
—
—
—
—
—
—
42
LSE
(2, 0)-(14, 9)
48
—
—
—
—
—
—
—
23
LSAD
(0, 1)-(14, 7)
72
—
—
—
—
—
—
—
21
M
(2, 1)-(9, 8)
69.45
0.75
0.6159
0.22
0.2819
—
—
—
17
L
p
(1.9, 0.9924)-(8, 7)
67.95
—
—
—
—
0.71
0.265
0.185
15
(1.8, 0.939)-(7, 7)
65.52
—
—
—
—
0.64
0.351
0.182
14Table 4
Parameters from the Huber estimation(2,1)-(9,8) (bnH, fnH) and from the Lp estimation (1.8, 0.939)-(7, 7) (bnp, fnp).
n
1
2
3
4
5
6
7
8
9
b
n
H
-
0.0856
-
0.1642
0.1062
0.2241
0.1460
0.0440
-
0.1929
-
0.0707
—
f
n
H
0.4241
-
0.5992
-
0.4307
-
0.0207
0.2362
0.3232
-
0.1615
0.0462
0.0744
b
n
p
-
0.0927
-
0.1415
0.1829
0.1553
0.0447
0.0412
-
0.1857
—
—
f
n
p
-
0.0699
-
0.7319
0.0578
0.0387
-
0.0021
0.3249
-
0.2332
—
—Figure 4
Frequency response of the(2,1)-(9,8) model (solid line) versus spectral estimation of the process (dashed line).Figure 5
Frequency response of the(1.8,0.939)-(7,7) model (solid line) versus spectral estimation of the process (dashed line).


## 7. Conclusion

In this paper, a new robust estimator based on a mixedLp distribution has been proposed in order to reduce the estimated model complexity. The main consistency and convergence properties of the robust estimator have been established and it has been shown that the tuning of DOF in the symmetric function improves the model order reduction. In effect, the additional DOF in the symmetric function of the estimation criterion present the advantage to tune its parameters. This flexibility offers to the user the choice to estimate and reduce the estimated model complexity in the same procedure. We presented experimental results from a vibration real process and we showed the effective reductions. Many aspects of these studies are opened to further researches. First, from a formal point of view, we should study the generalization of multiple DOF in a mixed symmetric function by investigating the link between the number of DOF and the model reduction. Then, it would be relevant to apply our Lp framework to a process with a level of contamination of the data more significant.



---

*Source: 349070-2015-03-19.xml*