## A unified framework for solving a general class of conditional and robust set-membership estimation problems /diamondmath

∗ /diamondmath

Vito Cerone , Member, IEEE, Jean-Bernard Lasserre , Dario Piga , § Member, IEEE, Diego Regruto , ∗ Member, IEEE

## Abstract

In this paper we present a unified framework for solving a general class of problems arising in the context of set-membership estimation/identification theory. More precisely, the paper aims at providing an original approach for the computation of optimal conditional and robust projection estimates in a nonlinear estimation setting where the operator relating the data and the parameter to be estimated is assumed to be a generic multivariate polynomial function and the uncertainties affecting the data are assumed to belong to semialgebraic sets. By noticing that the computation of both the conditional and the robust projection optimal estimators requires the solution to min-max optimization problems that share the same structure, we propose a unified two-stage approach based on semidefinite-relaxation techniques for solving such estimation problems. The key idea of the proposed procedure is to recognize that the optimal functional of the inner optimization problems can be approximated to any desired precision by a multivariate polynomial function by suitably exploiting recently proposed results in the field of parametric optimization. Two simulation examples are reported to show the effectiveness of the proposed approach.

/diamondmath Accepted for publication in IEEE Transactions on Automatic Control

∗ V. Cerone and D. Regruto are with Dipartimento di Automatica e Informatica, Politecnico di Torino, 10129 Torino, Italy vito.cerone@polito.it, diego.regruto@polito.it

/diamondmath J. B. Lasserre is with LAAS-CNRS and Institute of Mathematics, University of Toulouse, 7 Avenue du Colonel Roche 31 077 Toulouse Cedex 4, France lasserre@laas.fr

§ D. Piga is with the Control Systems Group, Department of Electrical Engineering, Eindhoven University of Technology, P.O. Box 513, 5600 MB, Eindhoven, The Netherlands. D.Piga@tue.nl

## I. INTRODUCTION

Estimation theory can roughly be defined as a branch of mathematics dealing with the problem of inferring the values of some unknown variables, usually called parameters , from a set of empirical data related to the unknown parameters through a given, possibly uncertain, mathematical relation. Experimental data are usually obtained by means of measurement procedures that are known to be affected by uncertainty. Most of the results available in the estimation theory literature are based on a statistical description of the uncertainty affecting the data.

A worthwhile alternative to the stochastic description, is the so-called bounded-error or setmembership characterization where measurement errors are assumed to be unknown but bounded (UBB), i.e., the measurement uncertainties are assumed to belong to a given bounded set. Such a description seems to be more suitable in those cases where either a priori statistical information is not available or the errors are better characterized in a deterministic way (e.g., systematic and class errors in measurement equipments, rounding and truncation errors in digital devices). Based on the UBB description of the uncertainty, a new paradigm called bounded-error or setmembership estimation has been proposed starting with the seminal work of Schweppe [1]. In the last three decades, set-membership estimation theory has been the subject of extensive research efforts which led to a number of relevant results with emphasis on the application of the set-membership paradigm in the context of system identification. The interested reader is referred to the book [2], the survey papers [3], [4] and the reference therein for a thorough review of the fundamental principles of the theory. Set-membership estimation algorithms can roughly be divided in two main categories: (i) set-valued estimators (see, e.g., [5]-[12] and the references therein), aimed at deriving either exact or approximate descriptions of the so-called feasible parameter set , i.e., the set of all possible parameter values consistent with the collected experimental data and a set of a-priori assumptions; (ii) pointwise estimators (see, e.g., [13]-[18] and the references therein), that return a single element of the parameter space according to a given selection criteria.

In this paper we focus on the latter category and, in particular, on two classes of pointwise estimation algorithms called conditional estimators and projection estimators respectively. In a nutshell, a set-membership estimation algorithm is called a conditional estimator when the sought estimate is constrained to belong to a given set (see, e.g., [15], [17]-[19]), while it is called a

projection estimator (see, e.g., [14]-[16]) when the parameter estimate is sought by minimizing a certain norm of the so-called regression error. To the best of the authors' knowledge, most of the results presented in the bounded-error literature about conditional and/or projection estimation are derived under a set of simplifying hypotheses, including the assumptions that: (i) the operator relating the parameter and the experimental data is linear and is not affected by uncertainty, (ii) the error affecting the measured data belongs to simple-shaped convex sets (e.g. boxes, ellipsoids) and (iii) the parameter estimate to be computed is looked for in the entire parameter space or at most in a linearly parameterized subset of the parameter space.

In this work, by recognizing that the problems of computing the conditional and projection estimates require the solution to min-max optimization problems that share essentially the same structure, a unified approach is proposed to approximate to any desired precision the optimal (either conditional or projection) estimate by assuming that (i) the operator relating the parameter and the experimental data is a generic nonlinear polynomial function possibly dependent on a set of uncertain variables assumed to belong to a given semialgebraic set, (ii) the error affecting the measured data belongs to a semialgebraic set and (iii) the parameter estimate to be computed is sought in a semialgebraic subset of the parameter space. It is worth noticing that in full generality, solving nonconvex min-max optimization problems is a real challenge for which no general methodology is available. An exception is a certain class of robust versions of some convex optimization problems when the uncertainty set has some special form. In this case, computationally tractable robust counterparts of these convex problems may exist. See for instance [20], [21] and the references therein.

The paper is organized as follows. The addressed estimation problem is formulated in Section II, where the proposed unified framework is also presented. A two-stage approach based on semidefinite-relaxation techniques for the solution of the considered class of estimation problems is then presented in Section III. The effectiveness of the proposed approach is demonstrated by means of two simulation examples in Section IV. Concluding remarks end the paper.

## II. PROBLEM FORMULATION

In this paper we consider a class of parametric nonlinear set-membership estimation problems where a given nonlinear operator F maps the parameter θ ∈ P θ ⊆ R /lscript to be estimated into the

output vector w ∈ R N as follows

$$\mathbf w = \mathcal { F } ( \theta, \varepsilon _$$

where ε F is an uncertain variable. The set P θ takes into account possible prior information on the parameter θ to be estimated. In this work P θ is assumed to be a semialgebraic set of the form where k , z z = 1 , . . . , m are multivariate polynomials in the /lscript components of the vector θ . Output measurements y ∈ R N are assumed to be corrupted by bounded noise as follows

$$\mathcal { P } _ { \theta } = \left \{ \theta \in$$

$$y = g ( \mathbf w, \varepsilon _ { y } )$$

where g is a polynomial function in the variable w and ε y . The uncertain variables ε F and ε y are assumed to belong to the following semialgebraic set

$$\begin{array} {$$

with h , i i = 1 , . . . , r being multivariate polynomials in the q components of the vector ε F and the N components of the vector ε y . In this work we restrict our attention to the case where the nonlinear operator F is a multivariate polynomial function of variables θ and ε F .

In the set-membership estimation framework, all the values of θ that are consistent with the assumed model structure described in (1), collected measurements y (3) and bounds on the uncertainty variables (4) are considered as feasible solutions to the estimation problem. The set D θ of all such values is called the feasible parameter set (FPS) and can be defined as the projection onto the parameter space R /lscript of the following set D :

$$\mathcal { D } = \left \{ ( \theta, \varepsilon _$$

In this paper we provide a unified approach to address some minmax estimation problems arising in set-membership identification. In particular we will refer throughout the paper to the general formulation of the considered identification problems presented in Section II-A.

## A. General min-max formulation of the considered class of SM estimation problems

The contribution of the paper is to provide an approach to solve the following general nonlinear set-membership estimation problems:

$$\text{P1} \colon \theta _ { \text{rob} } = \arg \min _ { \theta \in \mathcal { M } } \max _ { \alpha \in \mathcal { S } _ { \alpha, \theta } } J ( \theta, \alpha )$$

where M can be either D θ or P θ or any other possible subset of R /lscript described by a set of polynomial inequalities, α ∈ R T ,

$$\mathcal { S } _ { \alpha, \theta } = \left \{ \alpha \in \mathbb { R } ^ { T } \colon d _ { \mu } ( \alpha, \theta ) \geq 0, \mu = 1, \dots, M \right \} \\ \text{gebraic set, and } d _ { \mu }, \, \mu = 1, \dots, M \text{ are multivariate polynomials in the components}$$

is a semialgebraic set, and d , µ µ = 1 , . . . , M are multivariate polynomials in the components of the vectors α and θ . In the rest of the paper we will refer to (6) as robust SM estimation problem P1 . To the best of the authors' knowledge this is the first attempt towards the solution of a robust SM estimation problem in such a general form.

It is worth noting that computation of the global optimal solution θ rob of problem (6) is a difficult and challenging problem since (6) is an NP-hard robust nonconvex optimization problem. As already mentioned, in full generality there is no methodology to solve (6) except for robust versions of some convex optimization problems when the uncertainty set has some special form. Indeed, such problems have computationally tractable robust counterparts as described, for example, in [20], [21].

In sections II-B and II-C reported below, we show that the two classes of set-membership estimation problems considered in the paper ( conditional central estimation and robust projection estimation ) can be interpreted as two specific instances of problem P1 .

As will be discussed in details in the following, the approach proposed in this paper to solve problem P1 relies on the results presented in [22] (see Theorem 1 of this paper) that were derived under the following assumption:

Assumption 1: For each fixed value of θ = θ the set S α,θ is nonempty.

In sections II-B and II-C we show that Assumption 1 is always satisfied for the specific classes of estimation problems considered in the paper.

B. Conditional central estimation for set-membership nonlinear errors-in-variables parametric identification

Consider a single-input single-output (SISO) nonlinear dynamic system which transforms the noise-free input x t into the noise-free output w t according to

$$w _ { t } = f ( \theta, w _ { t - 1 }, w _ { t - 2 }, \dots w _ { t - n a }, x _ { t - 1 }, x _ { t - 2 }, \dots x _ { t - n b } )$$

where θ is the parameter vector to be estimated and f is assumed to be a multivariate polynomial function. Both input and output data sequences are corrupted by additive noise, ξ t and η t respectively, i.e.

$$u _ { t } = x _ { t } + \xi _ { t }, \ \ y _ { t } = w _ { t } + \eta _ { t }.$$

The noise samples ξ t and η t are bounded by given ∆ ξ t and ∆ η t respectively, that is:

$$| \, \xi _ { t } \, | \leq \Delta \xi _ { t }, \ \ | \, \eta _ { t } \, | \leq \Delta \eta _ { t }.$$

The nonlinear errors-in-variables (NEIV) model structure described by (8)-(10) can be written in the form (1) by setting:

$$\theta = \left [ a _ { 1 } \ \dots \ \ a _ { n a } \ \ b _ { 0 } \ \ b _ { 1 } \ \dots \ \ b _ { n b } \right ] ^ { T },$$

$$\mathbf w = [ w _ { 1 } \ w _ { 2 } \dots w _ { N } ],$$

$$\mathbf y = [ y _ { 1 } \ y _ { 2 } \dots y _ { N } ],$$

$$\varepsilon _ { F } = \left [ \varepsilon _ { F, 1 } \ \varepsilon _ { F, 2 } \ \dots \ \varepsilon _ { F, N } \right ] ^ { T }$$

$$\varepsilon _ { y, t } = [ \eta _ { 1 } \ \eta _ { t - 2 } \dots \eta _ { N } ],$$

$$\mathcal { F } ( \theta, \varepsilon _ { F } ) = \left [ \mathcal { F } _ { 1 } ( \theta, \varepsilon _ { F, 1 } ) \ \mathcal { F } _ { 2 } ( \theta, \varepsilon _ { F, 2 } ) \ \dots \ \mathcal { F } _ { N } ( \theta, \varepsilon _ { F, N } ) \right ] ^ { \mathbb { T } }$$

where, for all t = 1 2 , , . . . , N ,

$$\mathcal { F } _ { t } ( \theta, \varepsilon _ { F } ) = f ( \theta, y _ { t - 1 } - \eta _ { t - 1 }, \dots y _ { t - n a } - \eta _ { t - n a }, u _ { t - 1 } - \xi _ { t - 1 }, \dots u _ { t - n b } - \xi _ { t - n b } ). \quad ( 1 7 )$$

and

$$\varepsilon _ { F, t } & = [ \eta _ { t - 1 } \ \eta _ { t - 2 } \dots \eta _ { t - n a } \ \xi _ { t } \ \xi _ { t - 1 } \dots \xi _ { t - n b } ].$$

The NEIV model structure considered in (8) is quite general and comprises many important nonlinear model classes usually considered in system identification including, among the other,

Hammerstein, Wiener and Lur'e models (see, e.g., [23]-[25]), linear-parameter-varying (LPV) models with polynomial dependence on the scheduling variables [26], [27], polynomial nonlinear autoregressive (NARX) and nonlinear output error (NOE) models.

In the NEIV bounded error problem considered here, the set S ε F ,ε y is simply described by the linear inequalities in (10), while D in (5) is the set of all parameter values and noise samples consistent with the collected experimental data, the model structure in (8) and (9), and noise bounds in (10).

As far as the class of linear time-invariant dynamic systems is considered, the nonlinear operator F ( θ, ε F ) simplifies to:

$$\mathcal { F } ( \theta, \varepsilon _ { F } ) = F ( \varepsilon _ { F } ) ^ { \intercal } \theta$$

where

$$F ( \varepsilon _ { F } ) = [ F _ { 1 } \ F _ { 2 } \ \dots \ F _ { N } ] ^ { \mathbb { T } }$$

and, for all t = 1 2 , , . . . , N ,

$$F _ { t } = & \left [ - y _ { t - 1 } + \eta _ { t - 1 } \right ] - y _ { t - 2 } + \eta _ { t - 2 } \dots - y _ { t - n a } + \eta _ { t - n a } \\ & u _ { t } - \xi _ { t } \left ( u _ { t - 1 } - \xi _ { t - 1 } \dots u _ { t - n b } - \xi _ { t - n b } \right ].$$

As is well known, the linear EIV identification set-up (see [28] for details) is quite general in the sense that many other common linear identification problems can be written in this framework. In fact, the problem of identifying an output error (OE) model is obtained by setting ξ t = 0 , the case of finite-impulse-response (FIR) models is obtained for n a = 0 , while the structure in (8) and (9) turns out to be an equation error (EE) model when ξ t = 0 and η t = ∑ na i =1 a η i t -i .

It is worth noting that in the general NEIV problems D is a nonconvex semialgebraic set since the constraints y -F ( θ, ε F ) = ε y in (5) are polynomial functions of θ and ε F and, moreover, the same property holds true in the simplified linear-time-invariant case, where constraints y -F ( θ, ε F ) = ε y in (5) are bilinear in θ and ε F due to (19).

Although D is the set of all nonlinear dynamic models with structure (8) that are consistent with experimental data and measurement error bounds, neither the feasible parameter set nor the

tight outerbounding box derived in [11], [12] can be straightforwardly exploited for controller design or system behavior simulation. Thus, in many applications, the problem of selecting a single model among the feasible ones arises. One of the most common choices in the SM literature is to look for the value of the parameter θ that minimizes the worst case /lscript p estimation error computed over the entire feasible set, i.e.

$$\theta _ { c } \doteq \arg \min _ { \theta \in \mathbb { R } ^ { \ell } \, ( \theta _ { \nu }, \varepsilon _ { F }, \varepsilon _ { y } ) \in \mathcal { D } _ { \nu } } \| \theta _ { \nu } - \theta \| _ { p }.$$

where

$$\mathcal { D } _ { \nu } = & \left \{ ( \theta _ { \nu }, \varepsilon _ { F }, \varepsilon _ { y } ) \in \mathcal { P } _ { \theta _ { \nu } } \times \mathbb { R } ^ { q } \times \mathbb { R } ^ { N } \, \colon \\ & y = g ( \mathcal { F } ( \theta _ { \nu }, \varepsilon _ { F } ), \varepsilon _ { y } ), \ ( \varepsilon _ { F }, \varepsilon _ { y } ) \in \mathcal { S } _ { \varepsilon _ { F }, \varepsilon _ { y } } \right \}, \\ \text{s the $\ell_{p}$-norm of a vector.} \quad \text{$\cdot$} \quad \cdot \quad \cdot \quad \cdot \quad \frac { \infty } { \infty } \cdot \cdot \cdot \cdot \quad \frac { \infty } { \infty } \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \colon$$

and ‖ · ‖ p is the /lscript p -norm of a vector.

The estimate θ c computed by solving (22) is the so-called /lscript p -Chebyshev center of D , also called central estimate in the SM literature.

Remark 1: In the case p = ∞ , the central estimate is the center of the minimum-volume-box outerbounding D θ and can be computed by exploiting the convex relaxation approach proposed in [12].

Although the central estimate provides the minimum of the worst-case estimation error, it may show some undesirable features in the case of EIV identification or, more generally, when the set D is nonconvex. More precisely, in those cases, the Chebyshev center θ c is neither guaranteed to belong to the set D θ nor to the set P θ and, as a consequence, the identified LTI system could result inconsistent either with the experimental data or with some of the a-priori physical information on the parameter θ . In order to avoid such drawbacks, it is most desirable to force the computed parameter estimate to belong to a given set M by modifying the optimization problem (22) as follows

$$\theta _ { c } ^ { \mathcal { M } } \doteq \arg \min _ { \theta \in \mathcal { M } } \max _ { ( \theta _ { \nu }, \varepsilon _ { F }, \varepsilon _ { y } ) \in \mathcal { D } _ { \nu } } \| \theta _ { \nu } - \theta \| _ { p }$$

where M ⊂ R /lscript is assumed to be a semialgebraic set described by polynomial inequalities. Such a set is: (a) D θ if our aim is to constrain the computed estimate to belong to the feasible parameter set, (b) P θ if it is required to guarantee that the identified system satisfies the set of

available a-priori information, (c) D ∩P θ θ , or (d) any other semialgebraic set if, more generally, we want to force the identified system to belong to a particular model class. Problem (24) falls into the class of Conditional set-membership estimation problems and, in particular, θ M c is referred to as the conditional Chebyshev center of the feasible parameter set D with respect to the model class M . The problem of conditional central estimation is still a challenging problem in the field of set-membership identification/information based-complexity and a number of papers have appeared in the literature in the last decades on the subject (we refer the reader to the paper [18] and the references therein for a thorough review). In particular, conditional central algorithms have been proposed to effectively address the problems of reduced order modeling [15], set-membership state smoothing and filtering [29] and worst-case identification [18]. For such problems, computationally efficient and/or closed-form solutions to the problem of conditional central estimation have been derived assuming that: (i) F is a linear operator in both θ and ε F , (ii) ε F = 0 , (iii) M is a linear manifold, and (iv) S ε y is a simple-shaped convex set (usually a box, an ellipsoid or a polytope). Unfortunately, such assumptions are not satisfied in many relevant identification problems including, for example, the EIV problem considered in this section. As an additional motivating example leading to the class of estimation problems defined in (24), we mention the problem of identifying input-output linear systems that are apriori known to be bounded-input bounded-output (BIBO) stable. In this case, we are interested in computing the optimal estimate of the system parameter, in the Chebyshev center sense, over the set P stab of all the parameter values that guarantee BIBO stability of the system. Since, as shown in [30], the set P stab is semialgebraic and described by polynomial inequalities, such a problem naturally leads to a conditional estimation problem of the general form (24) where the set M = P stab .

Remark 2: It is worth remarking that in problem (24) S α,θ coincides with D ν where α = ( θ , ε ν F , ε y ) . Therefore, since D ν does not depend on θ , Assumption 1 is satisfied as long as D ν is a nonempty set, that is a common assumption in Set-membership identification, often satisfied in practice unless the identification problem is not well posed (e.g. the considered a-priori assumption on the system to be identified are completely wrong).

## C. Robust conditional projection estimation

Another class of estimator of particular interest in the set-membership/information-based complexity (IBC) framework is given by the so-called projection algorithms (see, e.g., [14], [16] and the references therein) where the parameter estimate is computed by solving an optimization problem of the following form:

$$\theta _ { p } = \arg \min _ { \theta \in \mathcal { M } } \| e _ { R } ( \varepsilon _ { F }, \theta ) \| _ { p }$$

where e R = y -F ( θ, ε F ) is the so-called regression error (see, e.g., [31], [32]) and the set M either coincides with R /lscript or is a subset of R /lscript . In the latter case the obtained estimator is called a conditional projection algorithm (see, e.g., [15]). The optimization criterion ‖ e R ‖ p is widely adopted in the identification literature (see, e.g., [31], [33], [34] and the references therein) and, in particular, the popular least-square estimator (see, e.g., [31], [32]) is obtained by setting p = 2 , M = R /lscript and assuming that ε F = 0 where 0 is the null element of the space R /lscript . Projection estimators and their optimality properties have been extensively investigated in the SM framework and a number of interesting results have been derived (see, e.g., [14], [16], [19]). However, to the best of the authors' knowledge, most of such results have been obtained under the assumptions that: (i) F is a linear operator in both θ and ε F , (ii) M is a linear manifold, (iii) S ε y is a simple-shaped convex set (usually a box, an ellipsoid or a polytope) and, most important, (iv) assuming that the operator F is not affected by uncertainty, i.e., ε F = 0 . In this work, we consider the following generalization of (25)

$$\theta _ { p } ^ { r } = \arg \min _ { \theta \in \mathcal { M } } \max _ { ( \varepsilon _ { F }, \varepsilon _ { y } ) \in \mathcal { S } _ { \varepsilon _ { F }, \varepsilon _ { y } } } \| e _ { R } ( \theta, \varepsilon _ { F } ) \| _ { p }$$

where, in order to take care of the effects of the uncertainty affecting the problem, we look for the parameter estimate θ r p that minimizes the worst case /lscript p regression error computed over the entire uncertainty set S ε F ,ε y . As far as the set M is concerned, we only assume that M be a subset of R /lscript described by polynomial inequalities. In such a way, the user is allowed, for example, to constrain the optimal estimate to belong to the feasible parameter set ( M = D θ ), to guarantee that θ r p satisfies the available a-priori physical information ( M = P θ ) or, more generally, to force the estimated model to belong to a specific, possibly reduced-order, model class (see, e.g., [15], [19]). In the rest of the paper we refer to θ r p as the robust p-norm projection estimate or RPE for short. It is worth noticing that the problem of computing the optimal projection estimate for the

/negationslash case p = 2 (least squares estimate) in the presence of an uncertainty ε F = 0 under the restricting assumptions that F is linear operator and the set S ε F ,ε y is convex, has been widely studied also outside the context of set-membership estimation and a number of different approaches have been proposed (see, e.g., paper [35] and the references therein for a thorough review of the available methods and results).

Remark 3: It is worth remarking that in problem (26) S α,θ coincides with S ε F ,ε y where α = ( ε F , ε y ) . Therefore, since S ε F ,ε y does not depend on θ , Assumption 1 is satisfied for all θ ∈ R /lscript as long as S ε F ,ε y is a nonempty set, that is a common assumption in Set-membership identification, often satisfied in practice unless the considered a-priori assumption on the measurement errors are completely wrong.

## III. A SEMIDEFINITE RELAXATION APPROACH

In this section a two-stage approach is proposed to approximate to any desired precision the global optimal solution to the general SM robust identification problem (6). The proposed approach is based on the following basic observations:

- (i) problem (6) is a two-players non-cooperative game (see, e.g., [36]) where M and S α,θ are the action sets of the first and the second player respectively;
- (ii) from the point of view of the second player, P1 is a parametric optimization problem (see [22] and the references therein) in the sense that the optimal value of the inner maximization problem in (6) is a function of the decision of player 1, i.e. the value of the parameter θ ;
- (iii) the optimal value function J of the parametric inner maximization problem is given by

and it is a function of parameter θ only.

$$s a \text{ union or use asion or player } 1, i. e. \text{ use value or use parameter } \sigma, \\ e \text{ function } \widetilde { J } \text{ of the parametric inner maximization problem is given by } \\ \widetilde { J } ( \theta ) \doteq J ( \theta, \alpha ^ { * } ( \theta ) ) = \max _ { \alpha \in \mathcal { S } _ { \alpha, \theta } } J ( \theta, \alpha ) \quad \text{ (27)} \\ \text{on of parameter } \theta \text{ only.} \\ \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \colon$$

Thanks to observations (i)-(iii) above, once ˜ J θ ( ) is known, problem P1 simplifies to the following optimization problem:

Unfortunately, a general methodology to derive an exact closed-form expression for function ˜ J θ ( ) is not available and, therefore, a two-stage procedure is proposed here to approximate the global optimal solution of problem P1 . In the first stage, a polynomial function ˜ J ∗ τ ( θ ) of degree August 5, 2014 DRAFT

$$P 2 \colon \theta _ { _ { \text{rob} } } = \arg \min _ { \theta \in \mathcal { M } } \widetilde { J } ( \theta ). \text{ \quad \ \ } ( 2 8 ) \text{ \quad \ } \text{
\text{methodology to derive an exact closed-form expression for function}$$

2 τ , upper approximating (in a strong sense) the optimal value function ˜ J θ ( ) of the parametric inner maximization problem is computed. Then, in the second stage, problem P2 is replaced with the following polynomial optimization problem:

$$\mathbf P 3 \colon \theta _ { \tau } = \arg \min _ { \theta \in \mathcal { M } } \widetilde { J } _ { \tau } ^ { * } ( \theta ). \\ \intertext { a t i o n of the function } \widetilde { J } ( \theta )$$

A. Polynomial approximation of the function ˜ J θ ( ) Polynomial approximation of the function ˜ J θ ( ) is performed here by exploiting the methodology for parametric polynomial optimization proposed in [22]. However, in order to apply the results presented in [22], we first need to compute a set R ⊂ θ R /lscript outerbounding M , whose shape is simple enough to allow one to easily compute all the moments of a Borel probability measure ϕ with uniform distribution on R θ . In this work we exploit the SDP-relaxation based procedure proposed in [12] to compute the minimum-volume axis-aligned box containing the set M . Once the box R θ = [ θ, θ ] ⊂ R /lscript , i.e.,

$$\mathcal { R } _ { \theta } \, = \, \{ \theta \in \mathbb { R } ^ { \ell } \ \colon \, \phi _ { k } ( \theta ) \geq 0, \, \phi _ { k } ( \theta ) \doteq ( \overline { \theta _ { k } } - \theta _ { k } ) ( \theta _ { k } - \underline { \theta } _ { k } ), \ \ k = 1, \dots, \ell \}$$

has been computed, we can formulate the following optimization problem where we look for the upper polynomial approximation ˜ J τ ( θ ) of the optimal value function ˜ J θ ( ) such that the integral ∫ R θ ˜ J τ ( θ dϕ θ ) ( ) is minimized:

$$\Lambda \text{ is minimized} & \quad \begin{matrix} \min _ { \tilde { J } _ { \tau } ( \theta ) } \int _ { \mathcal { R } _ { \theta } } \tilde { J } _ { \tau } ( \theta ) d \varphi \\ \text{ s.t. } \tilde { J } _ { \tau } ( \theta ) \geq J ( \theta, \alpha ) & \forall ( \theta, \alpha ) \in \mathcal { R } _ { \theta } \times \mathcal { S } _ { \alpha, \theta }. \end{matrix}$$

By noticing that (i) the objective function can be written as a linear combination of the moments of the uniform distribution measure supported on R θ and (ii) the inequality constraint can be approximately replaced by a SOS constraint, the following semidefinite relaxed problem

is obtained ( [22]):

where for each β = [ β 1 . . . β /lscript ] ∈ N /lscript and θ = [ θ 1 . . . θ /lscript ] the notation θ β stands for the monomial θ β 1 1 θ β 2 2 . . . θ β /lscript /lscript , N /lscript 2 τ = { β ∈ N /lscript : ∑ j β j ≤ 2 τ } , Σ[ θ, α ] is the set of SOS polynomials in the variables θ and α , while γ β are the moments of the Borel probability measure ϕ with uniform distribution on R θ , defined as (see, e.g., [37]):

$$\begin{cases} \sin _ { \lambda _ { \beta }, \sigma _ { \mu }, \psi _ { k } } \sum _ { \beta \in \mathbb { N } _ { 2 r } ^ { 2 } } \lambda _ { \beta } \gamma _ { \beta } } \\ \sin _ { \beta \in \mathbb { N } _ { 2 i } ^ { 2 } } \sum _ { \beta \in \mathbb { N } _ { 2 r } ^ { 2 } } \lambda _ { \beta } \theta ^ { \beta } - J ( \theta, \alpha ) = \sigma _ { 0 } ( \theta, \alpha ) + \sum _ { \mu = 1 } ^ { M } \sigma _ { \mu } ( \theta, \alpha ) d _ { \mu } ( \alpha, \theta ) + \sum _ { k = 1 } ^ { \ell } \psi _ { k } ( \theta, \alpha ) \, \phi _ { k } ( \theta ) \\ \sigma _ { \mu } \subset \Sigma [ \theta, \alpha ], \, \mu = 1, \dots, M \\ \psi _ { k } \subset \Sigma [ \theta, \alpha ], \, k = 1, \dots, \ell \\ \deg ( \sigma _ { 0 } ) \leq 2 \tau ; \, \deg ( \sigma _ { \mu } d _ { \mu } ) \leq 2 \tau, \, \mu = 1, \dots, M, \, \deg ( \psi _ { k } \phi _ { k } ) \leq 2 \tau, \, k = 1, \dots, \ell, \\ \intertext { where for each } \beta = [ \beta _ { 1 } \dots \beta _ { k } ] \in \mathbb { N } ^ { \ell } \, \text{and } \theta = [ \theta _ { 1 } \dots \theta _ { k } ] \, \text{the notation } \theta ^ { \beta } \, \text{ stands for the monomial} \\ \theta _ { 1 } ^ { \beta _ { 1 } } \theta _ { 2 } ^ { \beta _ { 2 } } \dots \theta _ { k } ^ { \beta _ { k } }, \, \mathbb { N } _ { 2 r } ^ { \ell } = \{ \beta \, \in \mathbb { N } ^ { \ell } \, \colon \sum _ { j } \beta _ { j } \leq 2 \tau \}, \, \Sigma [ \theta, \alpha ] \, \text{is the set of SOS polynomials in the} \end{cases}$$

$$\text{as (see, e.g., } \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \ iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota\iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iia \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \jmath \,. \text{as (see, e.g., } \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iia \iota \iota \iota \jmath \,.$$

Lemma 1: If R ×S θ α,θ contains an open set, then the semidefinite program (31) has an optimal solution ( λ , σ ∗ β ∗ µ , ψ ∗ k ) , µ = 0 , . . . , M , k = 1 , . . . , /lscript , provided that tau is sufficiently large.

A detailed proof of Lemma 1 is postponed to the Appendix.

Next, by applying the results presented in [22] about parametric polynomial optimization, it is possible to show that the optimal solution of (31) enjoys the important property stated in the following theorem:

Theorem 1: ([22]) Let ( λ , σ ∗ β ∗ µ , ψ ∗ k ) , µ = 0 , . . . , M , k = 1 , . . . , /lscript , be an optimal solution of problem (31) for a given degree τ and let us define the polynomial ˜ J ∗ τ ( θ ) = ∑ β ∈ N /lscript 2 τ λ θ ∗ β β . Then, ˜ J ∗ τ ( θ ) converges to the optimal value function ˜ J θ ( ) for the L 1 ( R θ , ϕ ) -norm as τ goes to infinity, i.e.:

For the proof of Theorem 1 we refer the reader to the paper [22]. We also have the following property.

$$\int _ { \mathcal { R } _ { \theta } } | \widetilde { J } _ { \tau } ^ { * } ( \theta ) - \widetilde { J } ( \theta ) | \, d \varphi ( \theta ) & \to 0. \\ \text{orem 1 we refer the reader to the paper } [ 2 ]. \, \text{We also have the following}$$

↦

Proposition 1: The optimal value function θ → ˜ J θ ( ) is upper semicontinuous (u.s.c.) on R θ .

Proof: Let ( θ n ) ⊂ R θ be a sequence such that θ n → θ as n →∞ , and

$$\lim \sup _ { z \to \theta } \widetilde { J } ( z ) \, = \, \lim _ { n \to \infty } \widetilde { J } ( \theta _ { n } ). \\ \ell ^ { * } ( \theta _ { n } ) \, \text{ be an arbitrary } \maximizer \, \text{in} \,. \\ \cdot \, \hat { \sim } \, \cdot \, \cdot \, \cdot \, \cdot \,.$$

Next, for each n ∈ N , let α ∗ ( θ n ) be an arbitrary maximizer in S α,θ for the max problem in (27). By compactness of R θ and S α,θ , there is a subsequence denoted ( n /lscript ) , /lscript ∈ N , and a point ( α, θ ) ∈ S α,θ ×R θ such that ( α ∗ ( θ n /lscript ) , θ n /lscript ) → ( α, θ ) as /lscript →∞ . Consequently, using continuity of J , lim sup z → θ ˜ J z ( ) = lim n →∞ ˜ J θ ( n ) = lim /lscript →∞ J θ ( n /lscript , α ∗ ( θ n /lscript )) = J θ, α ( ) ≤ max α ∈S α,θ J θ, α ( ) = ˜ J θ , ( ) which proves that ˜ J is u.s.c. Thanks to Theorem 1 we are in the position of proving the following result, which shows that the solution to problem P3 converges to the solution of P2 (and hence problem P1 ) as τ goes to infinity.

Theorem 2: Let ˜ J ∗ τ ( θ ) , τ ∈ N , be the polynomial defined in Theorem 1. Consider the polynomial optimization problem P3 in (29) with optimal value denoted by J ∗ τ , and let θ ∗ τ ∈ M be an optimal solution of P3 . Let ̂ J τ = min k ≤ τ J ∗ k = ˜ J ∗ k τ ( ) ( θ ∗ k τ ( ) ) for some k τ ( ) ∈ [1 , . . . , τ ] . Then:

Moreover, if ˜ J θ ( ) is continuous on M and θ rob in (28) is unique, then θ ∗ k τ ( ) → θ rob as τ →∞ . If θ rob is not unique then any accumulation point of the sequence ( θ ∗ k τ ( ) ) , τ ∈ N , is a global minimizer of problem min { ˜ ( J θ ) : θ ∈ M} .

$$\text{an} \colon \\ \lim _ { \tau \to \infty } \left ( \min _ { k \leq \tau } J _ { k } ^ { * } \right ) = \lim _ { \tau \to \infty } \widehat { J } _ { \tau } = \min _ { \theta \in \mathcal { M } } \max _ { \alpha \in \mathcal { S } _ { \alpha, \theta } } J ( \theta, \alpha ) = \colon J ^ { * }. \\ \text{over}, \, \text{if } \widetilde { J } ( \theta ) \text{ is continuous on } \mathcal { M } \text{ and } \theta _ { \text{rob} } \text{ in } (28) \text{ is unique, then } \theta _ { k ( \tau ) } ^ { * } \to \theta _ { \text{rob} } \text{ as } \tau \to \infty. \\ \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot$$

Proof: Observe that being ˜ J ∗ τ ( θ ) continuous on R θ (hence on M ), it has a global minimizer θ ∗ τ ∈ M , for every τ . From Theorem 1, ˜ J ∗ τ ( θ ) L 1 ( R θ ,ϕ ) → ˜ J θ ( ) (i.e., convergence in the L 1 ( R θ , ϕ ) -norm). Hence, by [38, Theorem 2.5.3], there exists a subsequence ( τ /lscript ) such that ˜ J ∗ τ /lscript ( θ ) → ˜ J θ ( ) , ϕ -almost uniformly on R θ .

/negationslash

Next, by Proposition 1, the optimal value mapping ˜ J is u.s.c. on R θ (hence on M ). With /epsilon1 &gt; 0 fixed, arbitrary, let B ( /epsilon1 ) = . { θ ∈ M : ˜ J θ ( ) &lt; J ∗ + /epsilon1 } and let κ . = ϕ ( B ( /epsilon1 )) . As ˜ J is u.s.c., B ( /epsilon1 ) is nonempty, open, and therefore κ &gt; 0 . As ˜ J ∗ τ /lscript ( θ ) → ˜ J θ ( ) , ϕ -almost uniformly on R θ , there exists a Borel set A κ ∈ B R ( θ ) such that ϕ A ( κ ) &lt; κ and ˜ J ∗ τ /lscript ( θ ) → ˜ J θ ( ) , uniformly on R \ θ A κ . Hence, as ∆ = ( . R \ θ A κ ) ∩ B ( /epsilon1 ) = ∅ , one has

$$\lim _ { \ell \to \infty } \widetilde { J } _ { \tau _ { \ell } } ^ { * } ( \theta ) \, = \, \widetilde { J } ( \theta ) \, \leq \, J ^ { * } + \epsilon, \quad \forall \, \theta \in \Delta,$$

and so, as J ∗ τ ≤ ˜ J ∗ τ ( θ ) on ∆ , one obtains lim /lscript →∞ J ∗ τ /lscript ≤ J ∗ + /epsilon1 . As /epsilon1 &gt; 0 was arbitrary, one finally gets lim /lscript →∞ J ∗ τ /lscript = J ∗ . On the other hand, by monotonicity of the sequence ( ̂ J τ ) , and so (34) holds.

$$J ^ { * } & \leq \lim _ { \tau \to \infty } \widehat { J } _ { \tau } \, = \, \lim _ { \ell \to \infty } \widehat { J } _ { \tau _ { \ell } } \, \leq \, \lim _ { \ell \to \infty } J _ { \tau _ { \ell } } ^ { * } \, = \, J ^ { * }, \\ \text{oids}. & \quad \sim$$

Next, let θ ∗ τ ∈ M be a global minimizer of ˜ J ∗ τ ( θ ) on M . As M is compact, there exists θ ∈ M and a subsequence τ /lscript such that θ ∗ k τ ( /lscript ) → θ as /lscript →∞ . In addition, from ˜ J ∗ τ ( θ ) ≥ ˜ J θ ( ) for every τ and every θ ∈ R θ ,

$$J ^ { * } \, & \leq \, \widetilde { J } ( \theta _ { k ( \tau _ { \ell } ) } ^ { * } ) \, \leq \, \widetilde { J } _ { k ( \tau _ { \ell } ) } ^ { * } ( \theta _ { k ( \tau _ { \ell } ) } ^ { * } ) \, = \, \widehat { J } _ { \tau _ { \ell } }. \\ \text{ting } \ell \to \infty \text{ yields the desired result } J ^ { * } = \widetilde { J } ( \overline { \theta } ). \, S \text{}$$

˜ Remark 4: Even though the sequence ˜ J ∗ τ → ˜ J for the L 1 -norm, the sequence ˜ J ∗ τ , τ ∈ N , is not necessarily monotone (meaning ˜ ( J ∗ τ θ ) ≥ ˜ J ∗ τ +1 ( θ ) for all τ ∈ N , θ ∈ R θ , does not necessarily holds). For this reason, it would be useful to know bounds on the distance between ˜ J ∗ τ and ˜ J for each fixed value of τ . Unfortunately, computation of such bounds is a difficult open problem that requires further investigation.

So using (35) and letting /lscript → ∞ yields the desired result J ∗ = ˜ J θ ( ) . So if the minimizer of ˜ J on M is unique, one has θ = θ rob , and as the converging subsequence ( τ /lscript ) was arbitrary, the desired result follows. If the minimizer is not unique then every accumulation point θ is a global minimizer since J ∗ = J θ ( ) , as just shown above.

## B. Solution to problem P3 via SDP relaxation

Once a polynomial approximation ˜ J ∗ τ ( θ ) of the optimal value function ˜ J θ ( ) of the inner maximization problem in (6) has been computed as discussed in III-A, we are in the position of solving problem P3 which is a multivariate polynomial optimization problem in the variable θ on the compact semi-algebraic set M . By applying the moments-based relaxation approach proposed in [39], a hierarchy of SDP relaxations ( Q t ) , t ∈ N , can be constructed with the following properties:

- -The resulting sequence (inf Q t ) , t ∈ N , of optimal values is monotone non decreasing and converges to the optimal value J ∗ τ of problem P3 .
- -If the global minimizer θ τ ∈ M of P3 is unique then the vector of 'first-order' moments of an optimal solution y of Q t converges to θ τ as t →∞ .

For more details on SDP relaxations for generalized moment problems, the interested reader is referred, e.g., to [37]. In fact, in view of recent results in [40] and [41], the convergence is finite provided that the problem satisfies a set of mild conditions (see [41] and the references therein for details), that is, generically , the optimal value J ∗ τ is attained at a particular relaxation in the hierarchy, i.e., J ∗ τ = inf Q t for some t . Finally, another recent result by Nie [42] ensures that generically , eventually some rank test is passed at some step t in the hierarchy, which permits to detect finite convergence at step t , and extract global minimizers (which generically are finitely many). The reader is referred to the papers [41], [42] for a discussion on the precise technical meaning of the word generically in this context.

Remark 5: [Exploiting sparsity] At a first sight, the applicability of the relaxation-based procedure proposed in this paper seems to be limited in practice to small-size identification problems, due to large dimensions of the SDP problems involved in the two stages of the proposed approach. That is certainly true for problem P1 in its general form (6), where the functional J θ, α ( ) is a generic multivariate polynomial and the sets M and S α,θ are generic semialgebraic sets. However, in view of the discussion and results reported in works [12], [23], it is possible to show that a number of identification problems arising from real-word applications enjoy a peculiar sparsity structure, called correlative sparsity in the framework of large-scale optimization (see, e.g., [43], [44]), which can be exploited to significantly reduce the computational complexity and the size of the involved SDP optimization problems either by means of the approach proposed in [43], [45] or by means of the ad-hoc procedure presented in [11]. More specifically, it is possible to show that a number of set-membership identification problems leads to semialgebraic optimization where the constraints and the functional satisfy the so-called running intersection property (see [45]), a condition that guarantees convergence of the solution of the relaxed problem to the global optimum of the polynomial problem also when the correlative sparsity pattern is used to derive semidefinite relaxations of reduced complexity (see, e.g., [43], [44]). Analysis of the correlative sparsity structure of problem P1 cannot be performed in general, since it requires to precisely specify the mathematical structure of the sets M and S α,θ . At the same time, providing a general discussion on the subject of sparsity exploitation in the context of SDP relaxation for polynomial problems is far beyond the scope of the paper, and the interested reader is referred to papers [43]-[45]. However, we will try to provide here a sketch of the main

ideas, covering the subject mostly at the level of intuition.

Let { 1 , . . . n } be the union ⋃ p k =1 I k of subsets I k ⊂ { 1 , . . . n } . A polynomial optimization problem is said to enjoy a correlative type of sparsity structure if: (i) each polynomial involved in the description of the set of constraints is only concerned with variables { X i : i ∈ I k } for some k ; (ii) the functional to be optimized J can be written as the sum J = J 1 + . . . + J p such that each J k only involves variables { X i : i ∈ I k } . Furthermore, the problem satisfies the running intersection property if the following condition is fulfilled:

$$I _ { k + 1 } \cap \bigcup _ { j = 1 } ^ { k } I _ { j } \subseteq I _ { s }, \text{for some } s \leq k \\ \text{detected either } h v \text{ insertion or } h v \text{ exnolitino}$$

The subsets { I k } can be detected either by inspection or by exploiting the systematic approach proposed in [43] and implemented in the software package [46]. If the problem enjoys a correlative sparsity structure, this can be used to derive SDP relaxations of lower complexity, as described in [43], [45]. Essentially, the intuitive idea underlying the approaches proposed in [43], [45] is the following: if the constraints and the objective function can be properly decomposed in subsets/subfunctionals depending only on a small subset of variables, then 'sparse' SDP relaxations can be constructed. This means that the involved SOS polynomials depend, each one, only on a small subset of variables of the original polynomial optimization problem. The fact that the linear EIV identification and the nonlinear Hammerstein identification problems enjoy a correlative sparsity structure satisfying the running intersection property, has been proved in previous papers [12], [23]. The same arguments/reasoning can be used to show that the polynomial approximation/optimization problems obtained by applying the approach proposed in this paper to the problem of conditional central estimation problem (24), enjoy the correlative sparsity structure and satisfy the running intersection property for many different choices of the set M including, e.g., the case M = D . This is also true for a number of problems in the class of robust conditional projection estimators including the nonlinear nonconvex robust least squares problem considered in Example 2 of Section IV.

Remark 6: It is worth-remarking that, by exploiting recent results presented in [47], the twostage relaxation-based procedure proposed in this section can be extended to a more general class of problems where the function J θ, α ( ) in problem P1 is a non-polynomial semialgebraic function.

with ε y ( ) t s.t.

$$\varepsilon _ { y } ( t ) & \leq 0 \quad \text{ if } \mathbf y ( t ) = 1, \\ \varepsilon _ { y } ( t ) & \geq - 1 \quad \text{ if } \mathbf y ( t ) = 0.$$

Based on eqs. (39), the uncertainty set S ε y can be written in terms of nonnegative inequality constraints as with

$$h _ { t } ( \varepsilon _ { y } ) \doteq \begin{cases} \ - \varepsilon _ { y } ( t ) & \text{ if } y ( t ) = 1, \\ \varepsilon _ { y } ( t ) + 1 & \text{ if } y ( t ) = 0. \end{cases} \\ \text{ DRAFT}$$

$$\mathcal { S } _ { \varepsilon _ { y } } = \left \{ \varepsilon _ { y } \in \mathbb { R } ^ { N } \colon h _ { t } ( \varepsilon _ { y } ) \geq 0, \ t = 1, \dots, N \right \}, \\ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

## IV. SIMULATION EXAMPLES

The capabilities of the presented approach are shown in this section by means of two simulation examples.

## Example 1

The first illustrative example comes from the problem addressed in [48] on the identification of ARX models based on quantized measurements. Consider the system analyzed in [48], i.e.,

$$\mathbf w ( t ) = \theta _ { 1 } ^ { o } \mathbf w ( t - 1 ) + \theta _ { 2 } ^ { o } u ( t ) + d ( t ) = 0. 6 \mathbf w ( t - 1 ) + 0. 6 u ( t ) + d ( t ),$$

where u t ( ) and w ( ) t are the input and output signals at time t , respectively, and d t ( ) is an unknown additive disturbance which is assumed to belong to the interval [ -0 1 . , 0 1] . . The system is simulated using a white input signal u t ( ) uniformly distributed within [ -2 5 . , 2 5] . and a disturbance d t ( ) with uniform distribution in the interval [ -0 1 . , 0 1] . . The output w ( ) t is measured by a binary sensor with threshold C = 1 , i.e., where y ( ) t is the output of the binary sensor. Indeed, the system output w ( ) t is not accessible and only its measurement y ( ) t is available. The estimate of the parameters θ of the system in (36) is computed based on a collection of N = 200 input/output measurements. Note that y ( ) t can be written in the form of (3) as follows:

$$\text{sensor with threshold $C=1$, i.e.,} y ( t ) = \begin{cases} \ 1 & \text{if $w(t) \geq 1$} \\ 0 & \text{otherwise} \end{cases} ( \text{$37$} )$$

$$y ( t ) = \mathbf w ( t ) + \varepsilon _ { y } ( t ),$$

Substitution of eq. (38) into (36) leads to the following relation between input and noise-corrupted output y ( ) t :

$$\mathbf y ( t ) = \theta _ { 1 } \left ( \mathbf y ( t - 1 ) - \varepsilon _ { y } ( t - 1 ) \right ) + \theta _ { 2 } u ( t ) + d ( t ) + \varepsilon _ { y } ( t ).$$

The FPS D θ for the considered system is thus defined as the projection over the parameter space of the set D defined by (42), (40) and the a-priori assumption on the disturbance d t ( ) , i.e.

$$\text{or} \, \text{ne set } \mathcal { V } \, \text{einen} \, \text{by} \, ( 4 2 ), ( 4 0 ) \, \text{and} \, \text{the} \, a \text{-prior} \, \text{assumption on} \, \text{the} \, a \text{/}, \, i. e. \\ \mathcal { D } = \left \{ ( \theta, d, \varepsilon _ { y } ) \in \mathbb { R } ^ { 2 + 2 N } \, \colon \, \text{y} ( t ) = \theta _ { 1 } \, ( \text{y} ( t - 1 ) - \varepsilon _ { y } ( t - 1 ) ) + \theta _ { 2 } u ( t ) + d ( t ) + \varepsilon _ { y } ( t ), \\ h _ { t } ( \varepsilon _ { y } ) \geq 0, \ \ - 0. 1 \leq d ( t ) \leq 0. 1, \ \ t = 1, \dots, N \right \}. \\ \dots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \end{$$

Note that D is described by polynomial constraints because of the product between the unknown parameter θ 1 and the noise ε y ( t -1) in the equality constraint appearing in (43). In this example we will compute the /lscript 2 -norm conditional Chebyshev center θ D c of the FPS D with respect to D θ itself, i.e.,

$$\theta _ { c } ^ { \mathcal { D } } \doteq \arg \min _ { \theta \in \mathcal { D } _ { \theta } } \max _ { ( \theta _ { \nu }, d, \varepsilon _ { y } ) \in \mathcal { D } } \| \theta _ { \nu } - \theta \| _ { 2 } ^ { 2 }.$$

In order to compute a solution to problem (44) through the procedure discussed in the paper, an outer-bounding box R θ of the FPS D θ is first evaluated by means of the approach proposed in [12] for bounding the parameters of linear systems in the bounded-error EIV framework. The computed outer-bounding box R θ is reported in Fig. 1, together with the true FPS D θ . Then, a polynomial ˜ J τ ( θ ) ∗ of degree 2 τ (with τ = 2 ) upper approximating the function is computed by solving the SDP problem (31). It is worth remarking that problem (31) enjoys a particular structured sparsity which is used to reduce the computational complexity in constructing the SOS polynomials in (31). In fact, the objective function ‖ θ ν - ‖ θ 2 2 in (45) only depends on the model parameters θ ν , while each constraint defining D in (43) only depends on a small subset of variables, namely, the model parameters θ ν , the disturbance d t ( ) and the noise samples ε y ( t -1) and ε y ( ) t . A correlative sparsity structure satisfying conditions in Remark 5 can be easily detected through a procedure similar to the one discussed in [12] in the context of set-membership EIV identification. The obtained 4-degree polynomial J τ ( θ ) ∗ given by

$$\widetilde { J } ( \theta ) = \max _ { ( \theta _ { \nu }, d, \varepsilon _ { y } ) \in \mathcal { D } } \| \theta _ { \nu } - \theta \| _ { 2 } ^ { 2 }, & ( 4 5 ) \\ \intertext { g the SDP problem (31). It is worth remarking that problem (31) enjoys }. \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot$$

$$\text{can be easily detected through a procedure similar to the one discussed in [12] in the context} \\ \text{of set-membership EIV identification. The obtained 4-degree polynomial } \mathbb { J } _ { \tau } ( \theta ) ^ { * } \text{ given by} \\ \mathbb { J } _ { \tau } ( \theta ) ^ { * } = & 0. 9 3 9 - 0. 7 9 5 \theta _ { 1 } - 0. 0 3 7 \theta _ { 2 } - 1. 0 3 9 \theta _ { 1 } ^ { 2 } - 3. 7 3 1 \theta _ { 1 } \theta _ { 2 } + 4. 3 1 5 \theta _ { 2 } ^ { 2 } + \\ & 2. 6 1 7 \theta _ { 1 } ^ { 3 } + 0. 7 9 0 \theta _ { 1 } ^ { 2 } \theta _ { 2 } + 2. 4 7 3 \theta _ { 1 } \theta _ { 2 } ^ { 2 } - 4. 5 6 7 \theta _ { 2 } ^ { 3 } + \\ & - 0. 9 6 8 \theta _ { 1 } ^ { 4 } + 1. 7 4 0 \theta _ { 1 } ^ { 3 } \theta _ { 2 } - 7. 2 6 7 \theta _ { 1 } ^ { 2 } \theta _ { 2 } ^ { 2 } + 7. 9 1 7 \theta _ { 1 } \theta _ { 2 } ^ { 3 } - 1. 1 6 6 \theta _ { 2 } ^ { 4 },$$

Fig. 1. Exact Feasible Parameter Set D θ (grey region), outer-bounding box R θ (region inside the box), (unconditional) Chebyshev center ( /squaresolid ), exact conditional Chebyshev center θ D c ( × ), approximation of the conditional Chebyshev center computed with the proposed two-stage approach (O).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0532v1-figure-e8cf2a76-8263-4e53-b987-a2a2d0155de1.jpg)

is plotted in Fig. 2, together with the true function ˜ J θ ( ) in (45), which in turn has been obtained by gridding.

Fig. 2. True function ˜ J θ ( ) in (45) (gray) and computed polynomial approximation ˜ J ∗ τ ( θ ) (black).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0532v1-figure-9e2e635d-ce0e-4e36-a1dc-0e4faa51d4c2.jpg)

The (unconditional) Chebyshev center of the FPS is computed by minimizing ˜ J τ ( θ ) ∗ over the outer-bounding box R θ , while an approximation of the conditional Chebyshev center is computed by solving problem P3 via the SDP relaxation approach discussed in Section III-B. Both the unconditional and the conditional Chebyshev center are reported in Fig. 1, which shows that the unconditional one does not belong to the FPS, while the computed approximation of the conditional Chebyshev center does. In the same figure, the exact conditional Chebyshev center θ D θ c , that is, the minimum of the true function ˜ J θ ( ) over the exact FPS D θ is also reported showing that the proposed relaxation approach is able to provide a good approximation of the global optimal solution to problem (44). The CPU time taken to compute the conditional Chebyshev center θ D θ c is about 1320 seconds on a 2.40-GHz Intel Pentium IV with 3 GB of RAM. More specifically, the time required to compute the polynomial approximation ˜ J θ ( ) (i.e. time required to compute the solution to problem (31)) is about 450 seconds, while the second step (solution to minimization problem P3 with order of relaxation 2 ) takes about 870 seconds. The maximum amount of memory used by Matlab during the computation was about 891 MB. The solver SeDuMi has been used to solve the SDP problem (31) and the SDP problems relaxing P3 .

## Example 2

In this example, the method is applied to the problem of robust estimation of a non-linear-inthe-parameter static model when both the input and the output measurements are corrupted by bounded noise.

The multi-input-single-output (MISO) data-generating system is given by

$$\mathbf w ( t ) = \theta _ { 1 } ^ { 0 } x _ { 1 } ( t ) + \theta _ { 1 } ^ { 0 } \theta _ { 2 } ^ { 0 } x _ { 2 } ( t ) + \theta _ { 3 } ^ { 0 } x _ { 3 } ( t ) + ( \theta _ { 1 } ^ { 0 } ) ^ { 2 } x _ { 4 } ( t ) + \theta _ { 4 } ^ { 0 } \theta _ { 5 } ^ { 0 } x _ { 5 } ( t ) + ( \theta _ { 5 } ^ { 0 } ) ^ { 2 } x _ { 6 } ( t ) + \theta _ { 4 } ^ { 0 } \theta _ { 6 } ^ { 0 } x _ { 7 } ( t ) =$$

$$\langle t \rangle + ( \theta _ { 1 } ^ { o } ) ^ { 2 } \, x _ { 4 } ( t ) + \theta _ { 4 } ^ { o } \theta _ { 5 } ^ { o } x _ { 5 } ( t ) + ( \theta _ { 5 } ^ { o } ) ^ { 2 } \, x _ { 6 } ( t ) + \theta _ { 4 } ^ { o } \theta _ { 6 } ^ { o } x _ { 7 } ( t ) =$$

$$= 1 x _ { 1 } ( t ) + ( 1 \cdot 0. 6 ) x _ { 2 } ( t ) - 0. 5 x _ { 3 } ( t ) + 1 x _ { 4 } ( t ) + 0. 3 \cdot 0. 8 x _ { 5 } ( t ) + ( 0. 8 ) ^ { 2 } x _ { 6 } ( t ) - 0. 3 \cdot 0. 5 x _ { 7 } ( t ),$$

where x i ( ) t , with i = 1 , . . . , 7 , is the i -th noise-free input and w ( ) t is the noise-free output at time t . The inputs x i ( ) t are i.i.d. random processes uniformly distributed in the interval [ -1 , 1] with length N = 400 . Both the inputs x i ( ) t and the output w ( ) t are corrupted by additive

uncertainties ξ i ( ) t and η t ( ) , respectively, i.e.,

$$u _ { i } ( t ) = & x _ { i } ( t ) + \xi _ { i } ( t ), \quad i = 1, \dots, 7,$$

$$\mathbf y ( t ) = & \mathbf w ( t ) + \eta ( t ),$$

where ξ i ( ) t and η t ( ) are white-noise processes uniformly distributed in the intervals [ -∆ ξ , i ∆ ] = ξ i [ -0 2 . , 0 2] . (for all i = 1 , . . . , 7 ) and [ -∆ η, ∆ ] = [ η -0 25 . , 0 25] . , respectively. The signal-tonoise ratio on the inputs, SNR x i , and on output, SNR w , defined as

$$\mathfrak { L } _ { t = 1 } \quad / \quad t = 1 \quad \mathfrak { J } \\ \text{SNR} _ { \mathfrak { w } } = 1 0 \log \left \{ \sum _ { t = 1 } ^ { N } \mathfrak { w } ^ { 2 } ( t ) \Big / \sum _ { t = 1 } ^ { N } \eta ^ { 2 } ( t ) \right \}, \\ \text{all } i = 1, \dots, 7 \text{ and } 1 6 \, \text{db, respectively. Let us denote with } \theta = [ \theta _ { 1 }, \, \theta _ { 2 }, \, \theta _ { 3 }, \, \theta _ { 4 }, \, \theta _ { 5 }, \, \theta _ { 6 } ]$$

$$\text{the inputs, $\mathbb{NK}_{x_{i}}$, and on output, $\mathbb{NK}_{w}$, dennead as} \\ \text{SNR} _ { x _ { i } } = 1 0 \log \left \{ \sum _ { t = 1 } ^ { N } x _ { i } ^ { 2 } ( t ) \Big / \sum _ { t = 1 } ^ { N } \xi _ { i } ^ { 2 } ( t ) \right \}, \\ \left ( \underline { \ N } \right ) / \underline { \ N } \right )$$

are 13 db (for all i = 1 , . . . , 7 ) and 16 db, respectively. Let us denote with θ = [ θ , 1 θ , 2 θ , 3 θ , 4 θ , 5 θ 6 ] the parameters of the model to be estimated. The FPS D θ is then given by the projection over the parameter space of the following set:

$$& \text{the parameter space of the following set:} \\ & \mathcal { D } = \left \{ ( \theta, \xi, \eta ) \in \mathbb { R } ^ { 3 + 5 N } \, \colon \, y ( t ) = \theta _ { 1 } \left ( x _ { 1 } \right ) - \xi _ { 1 } ( t ) \right ) + \theta _ { 1 } \theta _ { 2 } \left ( x _ { 2 } ( t ) - \xi _ { 2 } ( t ) \right ) + \\ & \quad + \theta _ { 3 } \left ( x _ { 3 } ( t ) - \xi _ { 3 } ( t ) \right ) + \theta _ { 1 } ^ { 2 } \left ( x _ { 4 } ( t ) - \xi _ { 4 } ( t ) \right ) + \\ & \quad + \theta _ { 4 } \theta _ { 5 } \left ( x _ { 5 } ( t ) - \xi _ { 5 } ( t ) \right ) + \left ( \theta _ { 5 } \right ) ^ { 2 } \left ( x _ { 6 } ( t ) - \xi _ { 6 } ( t ) \right ) + \theta _ { 4 } \theta _ { 6 } \left ( x _ { 7 } ( t ) - \xi _ { 7 } ( t ) \right ) + \\ & \quad + \eta ( t ), \\ & \quad \left | \eta ( t ) \right | \leq \Delta \eta, \ \left | \xi _ { i } ( t ) \right | \leq \Delta \xi _ { i }, \ \ t = 1, \dots, N, \ \ i = 1, \dots, 7 \right \}.$$

$$( 5 3 )$$

Now, let ˆ( y t, θ ) be the output of the model to be estimated, given by:

$$\hat { \mathbf y } ( t, \theta ) = \theta _ { 1 } u _ { 1 } ( t ) + \theta _ { 1 } \theta _ { 2 } u _ { 2 } ( t ) + \theta _ { 3 } u _ { 3 } ( t ) + \theta _ { 1 } ^ { 2 } u _ { 4 } ( t ) + \theta _ { 4 } \theta _ { 5 } u _ { 5 } ( t ) + \left ( \theta _ { 5 } \right ) ^ { 2 } u _ { 6 } ( t ) + \theta _ { 4 } \theta _ { 6 } u _ { 7 } ( t ).$$

In this example, we compute the parameter estimate θ ∗ = [ θ , ∗ 1 θ , ∗ 2 θ , ∗ 3 θ , ∗ 4 θ , ∗ 5 θ ∗ 6 ] that minimizes the worst-case /lscript 2 -loss function V ( θ, ξ ) , defined as

$$\mathcal { V } ( \theta, \xi ) = \sum _ { t = 1 } ^ { N } \left ( \mathbf y ( t ) - \hat { \mathbf y } ( t, \theta ) \right ) ^ { 2 },$$

over all possible realizations of the input uncertainties ξ i ( ) t in the interval [ -∆ ξ , i ∆ ] ξ i under the constraint that the identified parameters belong to the FPS. The considered estimation problem can be formulated as the following min-max optimization problem:

$$\hat { \theta } ^ { * } = \arg \min _ { \theta \in \mathcal { D } } \max _ { \xi \in \mathcal { S } _ { \xi } } \mathcal { V } ( \theta, \xi ),$$

where S ξ is defined as

$$\mathcal { S } _ { \xi } = \left \{ \xi \, \colon | \xi _ { i } ( t ) | \leq \Delta \xi _ { i }, \text{ \ for all } \, i = 1, \dots, 7, \text{ \ } t = 1, \dots, N \right \}.$$

It is worth noting that problem (56) is: (i) a nonlinear nonconvex least squares problems, due to the nonlinear-in-parameter structure of the system to be estimated; (ii) a robust nonlinear least-square problem, due to the presence of uncertainty in all the explanatory variables; (iii) a nonconvex constrained least square problem, since the optimal estimate is looked for over the feasible parameter set D θ . Therefore, problem (56) is a challenging estimation problem for which, to the best of the authors' knowledge, no solution has been previously proposed in the literature.

Here, the solution to Problem (56) is computed by applying the two-stage relaxation based method presented in the paper, which leads to the following estimate of the model parameters:

$$\hat { \theta } ^ { * } = \left [ \hat { \theta } ^ { * } _ { 1 }, \, \hat { \theta } ^ { * } _ { 2 }, \, \hat { \theta } ^ { * } _ { 3 }, \, \hat { \theta } ^ { * } _ { 4 }, \, \hat { \theta } ^ { * } _ { 5 }, \, \hat { \theta } ^ { * } _ { 6 } \right ] = [ 0. 9 8, \, 0. 5 7, \, - 0. 5 4, \, 0. 3 6, \, 0. 7 9, \, - 0. 5 9 ] \,.$$

It is worth remarking that, in order to apply the proposed method, an outer-bounding box of the feasible set D has been computed by suitable modifications of the algorithm proposed in [11]. Furthermore, as in Example 1, problem (56) enjoys a particular sparsity structure which has been exploited to reduce the computational load in solving (56). In fact, the objective function V ( θ, ξ ) is given by the sum of N terms ( y ( ) t -ˆ( y t, θ )) 2 , each one involving only the model parameters θ and the noise samples ξ i ( ) t (with i = 1 , . . . , 7 ) as unknown variables. Furthermore, each constraint defining S ξ in (57) only depends on the noise variable ξ i ( ) t . Similarly, each constraint defining the set D in (53) only involves a small subset of optimization variables, namely: the model parameters θ , the input noise samples ξ i ( ) t (with i = 1 , . . . , 7 ) and the output noise sample η t ( ) . By stacking the variables θ, ξ i ( ) t , η ( ) t in the vector

$$X = [ \theta _ { 1 }, \, \dots, \, \theta _ { 6 }, \, \xi _ { 1 } ( 1 ), \dots, \, \xi _ { 7 } ( 1 ), \dots, \, \xi _ { 1 } ( N ), \dots, \, \xi _ { 7 } ( N ), \, \eta ( 1 ), \, \dots, \, \eta ( N ) ] \,.$$

The index sets I t (with t = 1 , . . . , N ) introduced in Remark 5 and satisfying the running intersection property can be defined as

$$I _ { t } = \{ 1, \dots, 6, 6 + 7 ( t - 1 ) + 1, \dots, 6 + 7 ( t - 1 ) + 7, 6 + 7 N + t \} \,. \quad \quad ( 6 0 )$$

In this way, each constraint defining D in (53) is only concerned with variables { X i : i ∈ I t } , that is θ , ξ i ( ) t (with i = 1 , . . . , 7 ) and η t ( ) .

The performance of the estimated model is tested on a validation set with N val = 100 input/output measurements. The noise-free output w ( ) t and the estimated output signal ˆ( y t, θ ˆ ) ∗ are plotted in Fig. 3, while the difference between w ( ) t and ˆ( y t, θ ˆ ) ∗ is depicted in Fig. 4 showing a good agreement between the two signals. The CPU time taken to compute the parameter estimate ˆ θ ∗ is about 7 hours. More specifically, the time required to compute the solution to problem (31) with τ = 2 is about 2 5 . hours, while the second step (solution to minimization problem P3 with order of relaxation 2 ) takes about 4 5 . hours. The maximum amount of memory used by Matlab during the computation was about 1.9 GB. Based on the authors' experience, although sparsity is exploited, the identification problem considered in this example becomes computationally intractable (in commercial workstations and using general purpose SDP solvers like SeDuMi) when models with more that 7 parameters are considered.

Fig. 3. Noise-free output signal w ( ) t (thick line) and estimated output ˆ( y t, θ ˆ ) ∗ (thin line).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0532v1-figure-45e639c3-ab7e-4823-9a85-2e3760225fcb.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0532v1-figure-fbc9a127-b070-4c74-86aa-ef9519855008.jpg)

Fig. 4. Estimate output error w ( ) t - ˆ( y t, θ ˆ ) ∗ .

## V. CONCLUSIONS

In this work we have presented a two-stage approach, based on suitable convex semidefinite relaxations, for approximating the global optimal solution to a general class of min-max constrained semialgebraic optimization problems arising in the framework of set-membership estimation theory. We have shown that the proposed methodology can be profitably applied to the problem of computing both conditional central and robust projection estimators in a nonlinear setting where the operator relating the data and the parameter to be estimated is assumed to be a generic multivariate polynomial function and the uncertainties affecting the data are assumed to belong to semialgebraic sets. The key idea of the approach is to first compute a convergent polynomial approximation of the optimal value function of the inner maximization problem. Once such an approximation has been computed, the outer minimization problem reduces to a standard polynomial optimization problem solved by constructing a convergent hierarchy of semidefinite relaxations. Two simulation examples have been reported to show the effectiveness of the proposed approach. In particular, in the first example we have demonstrated that the presented two-stage algorithm provides good approximation of the global optimum of

the considered min-max estimation problems, while in the second example we have shown that the proposed methodology can be applied to compute the solution to a challenging nonconvex constrained robust least square estimation problem.

## APPENDIX

## Proof of Lemma 1.

To prove that the semidefinite program (31) has an optimal solution, we prove that Slater's condition holds for its dual, which is the semidefinite program:

$$\text{as for us} \, \text{and}, \, \text{when} \, is \, \text{use semuquemeine program} \colon \\ \max _ { z } \, & L _ { z } ( J ( \theta, \alpha ) ) \\ & M _ { \tau } ( z ) \succeq 0 \\ & M _ { \tau - r _ { \mu } } ( d _ { \mu } \, z ) \succeq 0, \quad \mu = 1, \dots, M \quad \text{(61)} \\ & M _ { \tau - 1 } ( \phi _ { k } \, z ) \succeq 0, \quad k = 1, \dots, \ell \\ & L _ { z } ( \theta ^ { \beta } ) \, = \, \gamma _ { \beta }, \quad \forall \beta \in \mathbb { N } _ { 2 \tau } ^ { \ell }$$

where

- · z = ( z κ ) , κ ∈ N /lscript + T 2 τ , is a sequence indexed in the canonical basis of monomials ( θ β α ν ) , of R [ θ, α ] 2 τ (the vector space of polynomials of degree at most 2 τ ).
- · M z τ ( ) is the moment matrix of order τ , associated with the sequence z .
- · M τ -r µ ( d µ z ) is the localizing matrix of order τ -r µ , associated with the sequence z and the polynomial d µ ∈ R [ θ, α ] (and where r µ . = (deg /ceilingleft d µ ) / 2 /ceilingright ).
- · L z : R [ θ, α ] → R is the so-called Riesz functional:

↦

For more details on moment and localizing matrices, and how they are used in polynomial optimization, the interested reader is referred to [37]. Now let O be an open set contained in R ×S θ α,θ , with projection O 1 on R θ . Let ϕ be the Borel probability measure uniformly distributed on R θ with moments ( γ β ) . Let ψ be the Borel probability measure on R ×S θ α,θ defined by:

$$\bullet \ L _ { z } \colon \mathbb { R } [ \theta, \alpha ] & \to \mathbb { R } \text{ is the so-called Riesz functional} \colon \\ p \left ( = \sum _ { ( \beta, \nu ) \in \mathbb { N } ^ { \ell + T } } p _ { \beta \nu } \theta ^ { \beta } \alpha ^ { \nu } \right ) \quad \mapsto \ L _ { z } ( p ) = \sum _ { ( \beta, \nu ) \in \mathbb { N } ^ { \ell + T } } p _ { \beta \nu } \, z _ { \beta, \nu }, \quad p \in \mathbb { R } [ \theta, \alpha ]. \\ \text{For more details on moment and localizing matrices, and how they are used in polynomial}$$

$$\varkappa _ { \theta } \, \text{with} \, \text{moms} \, ( \gamma _ { \beta } ) \, \text{Le} \, \psi \, \text{ve} \, \text{one} \, \text{don} \, \text{provision} \, \text{measure on} \, \varkappa _ { \theta } \times \varkappa _ { \alpha, \theta } \, \text{meuse} \, \text{t} \\ \psi ( A \times B ) \, \doteq \, \int _ { A } \mathbf Q ( B \, | \, \theta ) \, \varphi ( d \theta ), \quad \, B \in \mathcal { B } ( \mathcal { S } _ { \alpha, \theta } ), \, A \in \mathcal { B } ( \mathcal { R } _ { \theta } ),$$

where Q ( ·|· ) is a stochastic kernel on R × S θ α,θ such that Q ( · | θ ) is the probability measure uniformly distributed on S α,θ if θ ∈ O 1 , and Q ( ·| θ ) is any probability measure on S α,θ if θ ∈ R θ \ O 1 . Then let z be the sequence of moments associated with ψ , i.e.,

$$\text{$\text{$item_zero$ use sequence on moments associated with $\psi$, $i.e.,} \\ z _ { \beta, \nu } \doteq \int \theta ^ { \beta } \, \alpha ^ { \nu } \, d \psi ( \theta, \alpha ), \quad \forall \, ( \beta, \nu ) \in \mathbb { N } _ { 2 \tau } ^ { \ell + T }. \\ \mathbb { O M } \quad ( d \tau ) \succ \mathbb { O M } _ { 2 \nu } = 1 \quad M \text{ and } M \quad. ( d \theta, \tau ) \succ \mathbb { O M } _ { 2 \nu } \, k = \cdot$$

/negationslash

Then M z τ ( ) /follows 0 , M τ -r µ ( d µ z ) /follows 0 , µ = 1 , . . . , M , and M τ -1 ( φ k z ) /follows 0 , k = 1 , . . . , /lscript as well. Indeed suppose for instance that M z τ ( ) p = 0 for some vector p = 0 . Let ˜ p ∈ R [ θ, α ] τ be the polynomial with coefficient vector p .

$$\text{comment vector $p$}. \\ 0 \, = \, \langle p, \mathbf M _ { \tau } ( \mathbf z ) \, p \rangle \, = \, \int \tilde { p } ( \theta, \alpha ) ^ { 2 } \, d \psi ( \theta, \alpha ) = 0, \\ \text{that $\tilde{$n$ vanishes on the whole onen set $\O$ in contradiction with $t$} }$$

/negationslash but this implies that ˜ p vanishes on the whole open set O , in contradiction with p = 0 . Similarly, let µ ∈ { 1 , . . . , M } be arbitrary, and suppose that M τ -r µ ( z ) p = 0 for some vector p = 0 , and let ˜ p ∈ R [ θ, α ] τ -r µ be the polynomial with coefficient vector p .

/negationslash

$$\Im [ \sigma, \alpha ] _ { \tau - r _ { \mu } } \, \text{ use polynomial with even even } p. \\ 0 \, = \, \langle p, \mathbf M _ { \tau - r _ { \mu } } ( \mathbf z ) \, p \rangle \, = \, \int \tilde { p } ( \theta, \alpha ) ^ { 2 } \, d _ { \mu } ( \theta, \alpha ) \, d \psi ( \theta, \alpha ) = 0, \\ \text{minies that } \tilde { n } \, \text{ vanishes on the whole open set } \mathcal { O } \, \text{ in contradiction with } n \neq 0 \text{ } \mathcal { A }$$

/negationslash but this implies that ˜ p vanishes on the whole open set O , in contradiction with p = 0 . A similar argument shows that M τ -1 ( φ k z ) /follows 0 for every k = 1 , . . . , /lscript . Moreover, from the definitions of ψ and ϕ ,

$$\nu \, \text{and} \, \varphi, \\ L _ { z } ( \theta ^ { \beta } ) \, = \, \int _ { \mathcal { R } _ { \theta } \times \mathcal { S } _ { \alpha, \theta } } \theta ^ { \beta } \, d \psi ( \theta, \alpha ) \, = \, \int _ { \mathcal { R } _ { \theta } } \theta ^ { \beta } \, \varphi ( d \theta ) \, = \, \gamma _ { \beta }, \quad \beta \in \mathbb { N } _ { 2 \tau } ^ { \ell }, \\ \text{and} \, \alpha \, \pi \, \text{is oddiclable for } \mathcal { G } 1 \, \text{that for } \mathcal { G } 1 \, \text{and} \, \nu \, \text{is small} \, \text{now}$$

and so z is admissible for (61). Therefore Slater's condition holds for (61), and by a well-known result of convex optimization, the dual of (61) (i.e. (31)) has a optimal solution if its value is finite.

But the value of the primal semidefinite program (31) is bounded below by ∫ R θ J θ, α dϕ θ ( ) ( ) . Moreover, there exists M &gt; 0 such that M -J θ, α ( ) &gt; 0 on R ×S θ α,θ . Therefore by Putinar's Positivstellensatz, there exists some integer τ 0 such that

$$M - J ( \theta, \alpha ) = \sigma _ { 0 } ( \theta, \alpha ) + \sum _ { \mu = 1 } ^ { M } \sigma _ { \mu } ( \theta, \alpha ) d _ { \mu } ( \alpha, \theta ) + \sum _ { k = 1 } ^ { \ell } \psi _ { k } ( \theta, \alpha ) \, \phi _ { k } ( \theta ) \\ \text{there} \, \text{de} \sigma ( \sigma _ { \cap } ) \, < \, 2 \pi _ { \cap } \, \text{ de} \sigma ( \sigma _ { \dots } ) \, < \, 2 \pi _ { \cap } \, \text{ $u$ = 1}$} \quad M \text{ and } \text{de} \sigma ( \sigma _ { \dots } ) \, < \, 2 \pi _ { \cap } \, k = 1$$

where deg( σ 0 ) ≤ 2 τ 0 , deg( σ d µ µ ) ≤ 2 τ 0 , µ = 1 , . . . , M , and deg( ψ φ k k ) ≤ 2 τ 0 , k = 1 , . . . , /lscript . Hence the optimal value of (31) is finite whenever τ ≥ τ 0 and so (31) has an optimal solution.

/square

## REFERENCES

- [1] F. Schweppe, 'Recursive state estimation: unknown but bounded error and system input,' IEEE Trans. Aut. Control , vol. AC-13, pp. 556-558, 1968.
- [2] M. Milanese, J. P. Norton, H. Piet-Lahanier, and E. Walter, Eds., Bounding approaches to system identification . New York: Plenum Press, 1996.
- [3] M. Milanese and A. Vicino, 'Optimal estimation theory for dynamic systems with set membership uncertainty: an overview,' Automatica , vol. 27(6), pp. 997-1009, 1991.
- [4] E. Walter and H. Piet-Lahanier, 'Estimation of parameter bounds from bounded-error data: a survey,' Mathematics and Computers in simulation , vol. 32, pp. 449-468, 1990.
- [5] M. Milanese and G. Belforte, 'Estimation theory and uncertainty interval evaluation in presence of unknown but bounded errors: linear families of models and estimators,' IEEE Transactions on Automatic Control , vol. 27, pp. 408-414, 1982.
- [6] E. Fogel and F. Huang, 'On the value of information in system identification-bounded noise case,' Automatica , vol. 18, pp. 229-238, 1982.
- [7] R. Pearson, 'Block-sequential algorithms for set-theoretic estimation,' SIAM J. on Matrix Anal. Appl. , vol. 9, no. 4, pp. 513-527, 1988.
- [8] S. Veres, 'Polyhedron updating and relaxation for on-line parameter and state bounding,' in Proc. of IFAC Symposium on System Identification , 1994, pp. 371-376.
- [9] A. Vicino and G. Zappa, 'Sequential approximation of feasible parameter sets for identification with set membership uncertainty,' IEEE Transaction on Automatic Control , vol. 41, pp. 774-785, 1996.
- [10] L. Chisci, A. Garulli, A. Vicino, and G. Zappa, 'Block recursive parallelotopic bounding in set membership identification,' Automatica , vol. 34, no. 1, pp. 15-22, 1998.
- [11] V. Cerone, D. Piga, and D. Regruto, 'Improved parameters bounds for set-membership EIV problems,' International Journal of Adaptive Control and Signal Processing , vol. 57, no. 2, pp. 208-227, 2011.
- [12] --, 'Set-Membership Error-in-Variables Identification Through Convex Relaxation Techniques,' IEEE Transactions on Automatic Control , vol. 57, no. 2, pp. 517-522, 2012.
- [13] M. Milanese and R. Tempo, 'Optimal algorithms theory for robust estimation and prediction,' IEEE Trans. Automatic Control , vol. AC-30, no. 8, pp. 730-738, 1985.
- [14] B. Kacewicz, M. Milanese, A. Vicino, and R. Tempo, 'Optimality of central and projection algorithms for bounded uncertainty,' Systems and control letters , vol. 8, pp. 161-171, 1986.
- [15] B. Kacewicz, M. Milanese, and A. Vicino, 'Conditionally optimal algorithms and estimation of reduced order models,' J. Complexity , vol. 4, pp. 73-85, 1988.
- [16] M. Milanese, 'Properties of least-squares estimates in Set-membership identification,' Automatica , vol. 35, no. 5, pp. 767-776, 1999.
- [17] A. Garulli, 'Tight error bounds for projection algorithms in conditional set-memberhip estimation,' Systems and Control Letters , vol. 37, pp. 293-300, 1999.
- [18] A. Garulli, A. Vicino, and G. Zappa, 'Conditional central algorithms for worst-case set-membership identification and filtering,' IEEE Transactions on Automatic Control , vol. 45, no. 1, pp. 14-23, 2000.
- [19] A. Garulli, B. Kacewics, A. Vicino, and G. Zappa, 'Error Bounds for Conditional Algorithms in Restricted Complexity Set Membership Identification,' IEEE Transactions on Automatic Control , vol. 45, no. 1, pp. 160-164, 2000.

- [20] A. Ben-Tal, L. El Ghaoui, and A. Nemirovski, 'Robustness,' in Handbook of Semidefinite Programming: Theory, Algorithms, and Applications , H. Wolkowicz, R. Saigal, and L. Vandenberghe, Eds. Boston: Kluwer Academic Publishers, 2000.
- [21] A. Ben-Tal, S. Boyd, and A. Nemirovski, 'Extending scope of robust optimization: Comprehensive robust counterparts of uncertain problems,' Math. Program. Ser. B , vol. 107, pp. 63-89, 2006.
- [22] J. B. Lasserre, 'A 'joint+marginal' approach to parametric polynomial optimization,' SIAM J. Optimiz. , vol. 20, no. 4, pp. 1995-2022, 2010.
- [23] V. Cerone, D. Piga, and D. Regruto, 'Bounded error identification of Hammerstein systems through sparse polynomial optimization ,' Automatica , vol. 48, no. 10, pp. 2693-2698, 2012.
- [24] --, 'Bounding the parameters of block-structured nonlinear feedback systems,' International Journal of Robust and Nonlinear Control , vol. DOI: 10.1002/rnc.1813, 2012.
- [25] V. Cerone and D. Regruto, 'Parameter bounds evaluation of Wiener models with noninvertible polynomial nonlinearities,' Automatica , vol. 42, pp. 1775-1781, 2006.
- [26] R. T´ oth, Modeling and identification of Linear Parameter-Varying Systems . Heidelberg: Springer, Lecture Notes in Control and Information Sciences, Vol. 403, 2010.
- [27] V. Cerone, D. Piga, and D. Regruto, 'Convex relaxation techniques for set-membership identification of LPV systems,' in Proc. of American Control Conference ACC 2011 , 2011.
- [28] T. S¨ oderstr¨m, 'Errors-in-variables methods in sys tem identification,' o Automatica , vol. 43, no. 6, pp. 939-958, 2007.
- [29] A. Garulli, A. Vicino, and G. Zappa, 'Optimal induced-norm and set membership state smoothing and filtering for linear systems with bounded disturbances,' Automatica , vol. 35, no. 5, pp. 767-776, 1999.
- [30] V. Cerone, D. Piga, and D. Regruto, 'Enforcing stability constraints in set-membership identification of linear dynamic systems,' Automatica , vol. 47, no. 11, pp. 2488-2494, 2011.
- [31] L. Ljung, System Identification, Theory for the User . Upper Saddle River: Prentince Hall, 1999.
- [32] T. S¨ oderstr¨m and P. Stoica, o System Identification . Upper Saddle River: Prentice Hall, 1989.
- [33] H. Nyquist, 'The optimal l p , norm estimator in linear regression models,' Commun. Statist.-Theor. Meth. , pp. 2511-2524, 1983.
- [34] H. Akcay, H. Hjalmarsson, and L. Ljung, 'On the Choice of Norms in System Identification,' IEEE Transactions on Automatic Control , vol. 41, no. 9, pp. 1367-1372, 1996.
- [35] L. El Ghaoui and H. Lebret, 'Robust solutions to least-squares problems with uncertain data,' SIAM J. on Matrix Anal. Appl. , vol. 18, no. 4, pp. 1035-1064, 1997.
- [36] J. Nash, 'Non-cooperative games,' Annals of Math. , vol. 54, no. 2, pp. 286-295, 1951.
- [37] J. Lasserre, Moments, Positive Polynomials and Their Applications . Imperial College Press, Optimization series, Vol. 1, 2010.
- [38] R. B. Ash, Real Analysis and Probability . San Diego: Academic Press, Inc., 1972.
- [39] J. B. Lasserre, 'Global optimization with polynomials and the problem of moments,' SIAM J. Optimiz. , vol. 11, pp. 796-817, 2001.
- [40] M. Marshall, 'Representation of non-negative polynomials, degree bounds and applications to optimization,' Canad. J. Math. , vol. 61, pp. 205-221, 2009.
- [41] J. Nie, 'Optimality conditions and finite convergence of Lasserre's Hierarchy,' Dept. of Mathematics, University of California, Tech. Rep., 2012, arXiv:1206.0319, to appear in Math. Program. S´ er. A.

- [42] --, 'Certifying convergence of Lasserre's Hierarchy via Flat Truncation,' Dept. of Mathematics, University of California, Tech. Rep., 2012, arXiv:1106.2384, to appear in Math. Program. S´ er. A.
- [43] H. Waki, S. Kim, M. Kojima, and M. Muramatsu, 'Sums of squares and semidefinite programming relaxations for polynomial optimization problems with structured sparsity,' SIAM J. Optimiz. , vol. 17, no. 1, pp. 218-242, 2006.
- [44] K. Kobayashi, S. Kim, and M. Kojima, 'Correlative sparsity in primal-dual interior-point methods for lp, sdp, and socp,' Applied Mathematics and Optimization , vol. 58, no. 1, pp. 69-88, 2008.
- [45] J. B. Lasserre, 'Convergent semidefinite relaxations in polynomial optimization with sparsity,' SIAM J. Optimiz. , vol. 17, no. 1, pp. 822-843, 2006.
- [46] H. Waki, S. Kim, M. Kojima, M. Muramatsu, and H. Sugimoto, 'SparsePOP: a sparse semidefinite programming relaxation of polynomial optimization problems,' ACM Transaction on Mathematical Software , vol. 35, no. 2, 2008.
- [47] J. B. Lasserre and M. Putinar, 'Positivity and optimization for semi-algebraic functions,' SIAM J. Optimiz. , vol. 20, no. 6, pp. 3364-3383, 2010.
- [48] M. Casini, A. Garulli, and A. Vicino, 'Set-membership identification of arx models with quantized measurements,' in Proc. of 50th IEEE Conference on Decision , 2011, pp. 2806-2811.
Vito Cerone Jean-Bernard Lasserre Dario Piga Diego Regruto
A unified framework for solving a general class of conditional and robust set-membership estimation problems
In this paper we present a unified framework for solving a general class of
problems arising in the context of set-membership estimation/identification
theory. More precisely, the paper aims at providing an original approach for
the computation of optimal conditional and robust projection estimates in a
nonlinear estimation setting where the operator relating the data and the
parameter to be estimated is assumed to be a generic multivariate polynomial
function and the uncertainties affecting the data are assumed to belong to
semialgebraic sets. By noticing that the computation of both the conditional
and the robust projection optimal estimators requires the solution to min-max
optimization problems that share the same structure, we propose a unified
two-stage approach based on semidefinite-relaxation techniques for solving such
estimation problems. The key idea of the proposed procedure is to recognize
that the optimal functional of the inner optimization problems can be
approximated to any desired precision by a multivariate polynomial function by
suitably exploiting recently proposed results in the field of parametric
optimization. Two simulation examples are reported to show the effectiveness of
the proposed approach.