X Bombois

G Scorletti

G Mercère

Optimal identification experiment design with robust output amplitude constraints

Keywords: Optimal experiment design, Amplitude constraints, System identification

In this paper, we consider optimal experiment design with amplitude constraints for the case of a multisine excitation. The output amplitude constraint is robustified using an uncertainty region for the unknown true system. The main contribution of this paper is to treat the robust amplitude constraint without the classical gridding approximation. In particular, we provide an LMI optimization problem to verify whether the output constraint is respected for all systems in the uncertainty region. This LMI formulation can be combined with the accuracy constraint in the optimal experiment design problem. Like each optimal experiment design problem with amplitude constraints, the resulting optimization problem is non-convex, but we provide two approaches to efficiently initialize the algorithm used to solve this optimization problem. The efficiency of the proposed methodology is tested using a numerical simulation.

Introduction

In this paper, we consider optimal identification experiment design when the constraints on the input and the output signals are not expressed in term of power, but in term of time-domain amplitude.

Optimal identification experiment design has been extensively studied for prediction error identification [START_REF] Ljung | System Identification: Theory for the User, 2nd Edition[END_REF]. In most of the contributions (see, e.g., [START_REF] Ljung | System Identification: Theory for the User, 2nd Edition[END_REF][START_REF] Jansson | Input design via LMIs admitting frequency-wise model specifications in confidence regions[END_REF][START_REF] Bombois | Least costly identification experiment for control[END_REF][START_REF] Hjalmarsson | System identification of complex and structured systems[END_REF][START_REF] Annergren | Application-oriented input design in system identification: Optimal input design for control [applications of control[END_REF]), the constraints on the input and the output signals are expressed as power constraints. One of the main reasons for this is that both the inverse of the covariance matrix of the estimate θN of the parameter vector θ 0 of the true system (which is generally the measure of the model accuracy that is used in the accuracy constraint of the optimal experiment design problem) and the power of the input and output signals (which is a measure of the cost of the identification) are affine functions of the power spectrum of the excitation signal. This property allows one to formulate the optimal experiment design problem as a convex optimization problem (see, e.g., [START_REF] Jansson | Input design via LMIs admitting frequency-wise model specifications in confidence regions[END_REF][START_REF] Bombois | Least costly identification experiment for control[END_REF][START_REF] Hjalmarsson | System identification of complex and structured systems[END_REF]).

The main issue with this approach is that, in many applications (see, e.g., [START_REF] Manchester | Amplitude-constrained input design: Convex relaxation and application to clinical neurology[END_REF] for an example), the constraints are not formulated as constraints on the power of the input and output signals, but as constraints on the amplitude of the time-domain sequence of these signals. Designing optimal identification experiments with amplitude constraints is thus a very important problem in practice. Unlike in the power constraint case, such optimal experiment design problem boils down to a non-convex optimization problem. The paper [START_REF] Manchester | Input design for system identification via convex relaxation[END_REF] indeed shows that it can be reformulated as an optimization problem involving both LMI constraints and a rank constraint which makes the optimization problem non-convex. In the literature, different approaches have been proposed to address this non-convex optimization problem. An iterative approach involving different nonconvex optimization problems is proposed in [START_REF] Manchester | An algorithm for amplitude-constrained input design for system identification[END_REF] (see also [START_REF] Potters | Experiment time minimisation under parameter accuracy constraints and time-domain signal amplitude bounds[END_REF]) while the accuracy objective and the respect of the amplitude constraints are treated separately in [START_REF] Dirkx | A fast smoothingbased algorithm to generate l∞-norm constrained signals for multivariable experiment design[END_REF]. Other contributions use relaxations to approximate the original non-convex optimization problem by a convex one. A first order approximation is used in [START_REF] Suzuki | On input design for system identification in time domain[END_REF]. In [START_REF] Manchester | Input design for system identification via convex relaxation[END_REF], the non-convex rank constraint is neglected and a randomization procedure is subsequently devised to improve the results. A similar approach is used in [START_REF] Deshpande | Constrained optimal input signal design for data-centric estimation methods[END_REF][START_REF] Larsson | Experimental evaluation of model predictive control with excitation (mpc-x) on an industrial depropanizer[END_REF]. Other relaxations are proposed in [START_REF] Ebadat | Applications oriented input design for closedloop system identification: a graph-theory approach[END_REF][START_REF] Ebadat | Applications oriented input design in time-domain through cyclic methods[END_REF].

Since the output obviously depends on the unknown true parameter vector θ 0 (i.e., on the true system G(θ 0 )), we here consider an optimal experiment design problem with a robust output amplitude constraint i.e., the output amplitude constraint must be satisfied for the outputs of all systems G(θ) with θ in an uncertainty re-gion U init containing the unknown θ 0 [START_REF] Rojas | Robust optimal experiment design for system identification[END_REF]. Since θ 0 ∈ U init , such a formulation ensures that the output of the unknown true system G(θ 0 ) will respect the amplitude constraint. Robust amplitude constraints (and other types of robust constraints in optimal experiment design) are classically handled via a gridding approximation (see [START_REF] Jansson | Input design via LMIs admitting frequency-wise model specifications in confidence regions[END_REF][START_REF] Bombois | Least costly identification experiment for control[END_REF][START_REF] Rojas | Robust optimal experiment design for system identification[END_REF][START_REF] Manchester | Input design for system identification via convex relaxation[END_REF]): the uncertainty set U init (containing an infinite number of elements) is replaced by a number n of grid points of this uncertainty set U init and the robust output amplitude constraint is replaced by n constraints (one for each grid point) in the optimal experiment design problem. Since θ 0 will most probably not be among the grid points, this gridding approach entails an approximation i.e., the output of the unknown true system G(θ 0 ) is not guaranteed to satisfy the amplitude constraint. The main contribution of this paper is to develop a methodology (inspired from robustness analysis) to handle the robust output amplitude constraint without the approximation introduced by the gridding of U init . This approach will be based on the philosophy we introduced in previous contributions to handle robust output power constraints without approximation (see [START_REF] Bombois | Robust optimal identification experiment design for multisine excitation[END_REF] for multisine excitation and see [START_REF] Bombois | Identification for robust H 2 deconvolution filtering[END_REF] for filtered white noise excitation).

As mentioned above, we will develop a methodology inspired from robustness analysis to handle robust output amplitude constraints. Since the (robust) output amplitude constraint can easily be translated into a (robust) output power constraint for filtered white noise excitation, we will only consider the case of multisine excitation in this paper. More precisely, we will develop an LMI optimization problem allowing, for a given multisine excitation, to verify whether the output amplitude constraint is satisfied for the outputs of all systems G(θ) with θ ∈ U init . Since amplitude constraints must be respected at each time instant, we treat, in this LMI formulation, the time similarly as the parameter vector θ i.e., as an uncertain parameter varying in a set. This prevents the use of a gridding of the time axis (which would introduce a similar approximation as the gridding of U init ). Using the same philosophy, we also develop an LMI optimization problem to verify whether the input amplitude constraint is verified at each time instant. These LMI formulations for the input and output amplitude constraints can be combined with the accuracy constraint in the optimal experiment design problem. The resulting optimization problem obviously remains a non-convex optimization problem, but we provide two approaches to efficiently initialize the numerical algorithm used to solve this optimization problem.

Notations. Continuous-time signals will be denoted x(t) where t ∈ R is the time index while discrete-time signals will be denoted x [n] where n is the sample number. The variable s is the Laplace variable while z will denote both the Z-transform variable and the shift operator. We use j to represent √ -1. For a complex number a (i.e., a ∈ C), |a|, ∠a and Re(a) will denote, respectively, its modulus, its argument and its real part. For a real number a (i.e., a ∈ R), |a| is the abolute value of a. For a matrix A, rank(A) is the rank of A, λ min (A) (resp. λ max (A)) is the smallest (resp. largest) eigenvalue of A while A T (resp. A * ) is its transpose (resp. conjugate transpose). The matrix

     X1 0 0 0 . . . 0 0 0 Xn     
will be denoted diag(X 1 , ..., X n ) if the elements X i (i = 1, ...n) are scalar quantities, while it will be denoted bdiag(X 1 , ..., X n ) if the elements X i (i = 1, ...n) are matrices. In addition, I n represents the identity matrix of dimension n × n and ⊗, the Kronecker product.

2 Problem statement 2.1 Identification experiment and experiment design problem We consider the identification of a discrete-time model of a stable single-input single-output true system which can be described by a continuous-time transfer function G 0 (s) with input u and output y. The discrete-time data for the identification of the discrete-time model of G 0 (s) will be gathered in open loop with a sampling rate T s . The continuous-time true system will be excited by the following (continuous-time) multisine excitation signal:

u(t, A) = L i=1 (a i,s sin(ω i t) + a i,c cos(ω i t)) (1)
with A = (a 1,s , a 1,c , a 2,s , a 2,c , ..., a L,c ) T ∈ R 2L . The fundamental frequency of this multisine will be denoted ω 0 . Each of the L frequencies ω i (i = 1, ..., L) thus satisfies the following relation:

ω i = α i ω 0 (2)
for an integer α i ̸ = 0 (i = 1, ..., L). Since the discretetime data will be gathered with a sampling rate T s , we will suppose that ω L < π Ts . In the optimal experiment design problem that will be introduced in the sequel, the L frequencies ω i (i = 1, ..., L) are fixed a-priori by the user and the amplitudes a i,s and a i,c (i = 1, ..., L) will therefore be the decision variables of the optimal experiment design problem1 . These amplitudes (which are also the coefficients of the trigonometric Fourier series of u(t, A)) are gathered in the amplitude vector A ∈ R 2L defined above.

In the sequel, we will address (robust) output amplitude constraints based on a steady-state expression of the output. To ensure as much as possible that the output constraint is also satisfied during transient, the excitation signal (1) will be slowly ramped up 2 . By doing this, we prevent as much as possible that the output exceeds its steady-state expression during transient. For reference purpose, we will assume that this ramp up process is performed when t < 0 and that the system is at steady state at t = 0. At that moment, we start collecting discrete-time input-output data at a sampling rate T s applying an anti-aliasing filter to both the input and the output. This leads to a set of discrete-time data

Z N = {y[n], u[n] | n = 0, ..., N -1} which are related as follows y[n] = G(z, θ 0 )u[n] + v[n] (3) where G(z, θ 0 ) is a stable discrete-time transfer func- tion satisfying G(e jωTs , θ 0 ) = G 0 (jω) for ω ∈ [0, π
Ts ] (and therefore also for all ω i (i = 1, ..., L) in ( 1)) and where v[n] is a discrete-time signal representing the disturbance acting on the system. This disturbance is modeled as a time series v[n] = H(z, θ 0 )e[n] with e[n] a white noise of variance σ 2 e and H(z, θ 0 ) a stable, inversely stable and monic transfer function [START_REF] Ljung | System Identification: Theory for the User, 2nd Edition[END_REF]. The transfer functions G(z, θ 0 ) and H(z, θ 0 ) are parameterized with an unknown parameter vector θ 0 ∈ R k in a given model structure G(z, θ) and H(z, θ). The typical model structures (BJ, ARX, ARMAX, ...) are all rational in the vector θ [START_REF] Ljung | System Identification: Theory for the User, 2nd Edition[END_REF]. In particular, for these model structures, G(z, θ) can be written as (Z N (z)θ)/(1 + Z D (z)θ) where Z N (z) and Z D (z) are row vectors of transfer functions [START_REF] Bombois | Robustness analysis tools for an uncertainty set obtained by prediction error identification[END_REF][START_REF] Bombois | Robust optimal identification experiment design for multisine excitation[END_REF].

Under the instrumentation described above, the discrete-time signal u

[n] in (3) is given by u[n] = u(t = nT s , A) (n = 0, ..., N -1) while y[n]-v[n] = G(z, θ 0 )u[n]
(n = 0, ..., N -1) is given by y(t = nT s , A, θ 0 ) with

y(t, A, θ) = L i=1 g(ωi, θ)... ... (ai,s sin(ωit + ϕ(ωi, θ)) + ai,c cos(ωit + ϕ(ωi, θ))) (4)
where g(ω i , θ) = |G(e jωiTs , θ)| and ϕ(ω i , θ) = ∠G(e jωiTs , θ). Note also that the continuous-time signal y(t, A, θ 0 ) represents the noise-free output of the continuous-time system for t ≥ 0.

Let us now specify precisely the amplitude constraints that have to be respected when designing the multisine excitation u(t, A). Amplitude constraints. Consider the time interval T P = [0, 2π ω0 ] corresponding to a period of the multisines u(t, A) (see [START_REF] Annergren | Application-oriented input design in system identification: Optimal input design for control [applications of control[END_REF]) and y(t, A, θ 0 ) (see (4)). The multisines u(t, A) and y(t, A, θ 0 ) must satisfy

-ūmax ≤ u(t, A) ≤ ūmax ∀t ∈ T P (5) -ȳmax ≤ y(t, A, θ 0 ) ≤ ȳmax ∀t ∈ T P (6)
where ūmax > 0 and ȳmax > 0 are two user-chosen constants. In (5) and (6), we can restrict attention to t ∈ T P since u(t, A) and y(t, A, θ 0 ) are both periodic signals having a fundamental period equal to 2π ω0 .

Using the data set Z N , the prediction error estimate θN of θ 0 [START_REF] Ljung | System Identification: Theory for the User, 2nd Edition[END_REF] is given by θN = arg min θ

1 N N -1 n=0 ϵ 2 [n, θ] with ϵ[n, θ] = H -1 (z, θ) (y[n] -G(z, θ)u[n]) . (7)
Assuming that, for the chosen L in (1), the data Z N are informative [START_REF] Ljung | System Identification: Theory for the User, 2nd Edition[END_REF], the estimate θN is (asymptotically) normally distributed around θ 0 with a covariance matrix P θ > 0 that can be estimated from the data [START_REF] Ljung | System Identification: Theory for the User, 2nd Edition[END_REF]. We can also determine an expression of the inverse of the covariance matrix as a function of the true parameter vector θ 0 and of the square of the amplitudes a i,s and a i,c (i = 1, ..., L) i.e.,

P -1 θ (A, θ 0 ) = L i=1 M (ω i , θ 0 ) a 2 i,s + a 2 i,c + M v (θ 0 ) (8)
for matrices M (ω i , θ 0 ) ≥ 0 (i = 1, ..., L) and M v (θ 0 ) ≥ 0 of known expressions and which are both proportional to the number of data N (see, e.g., [START_REF] Potters | Experiment time minimisation under parameter accuracy constraints and time-domain signal amplitude bounds[END_REF]). The covariance matrix P θ > 0 is a measure of the modeling error θN -θ 0 and its inverse therefore a measure of the model accuracy.

In this paper, our objective will be to design the amplitude vector A defining the excitation signal (1) as the one solving the following (non-convex) optimization problem for a given matrix R adm ≥ 0: arg max ξ,A ξ s.t. P -1 θ (A, θ 0 ) ≥ ξ R adm and ( 5)-( 6).

Let us discuss this formulation more precisely. If you choose R adm = I k , ( 9) is an E-optimal experiment design problem whose objective is to determine the amplitude vector A of the multisine excitation (1) which leads to an estimate θN of θ 0 with the smallest3 covariance matrix P θ (A, θ 0 ) while guaranteeing that the amplitude constraints ( 5) and ( 6) are respected. Since P θ is inversely proportional to N , it is clear that higher accuracy for θN can be achieved under the amplitude constraints ( 5) and (6) if N increases (i.e., if the experiment duration increases). In this sense, the optimization problem (9) also allows to determine the minimal experiment duration N for which it is possible to find an amplitude vector A that respects (5) and ( 6) and ensures P -1 θ (A, θ 0 ) ≥ R adm for some matrix R adm representing the desired accuracy (see, e.g., [START_REF] Potters | Experiment time minimisation under parameter accuracy constraints and time-domain signal amplitude bounds[END_REF] for the details on how this can be done based on (9)).

Remark. Instead of (1), the excitation applied to the continuous-time system may also be the output of a Zero Order Hold to which a discrete-time multisine u[n] is applied. In (3), u[n] is thus here also given by L i=1 (a i,s sin(ω i nT s ) + a i,c cos(ω i nT s )). Due to the property of the Zero Order Hold, we could thus keep [START_REF] Bombois | Robustness analysis tools for an uncertainty set obtained by prediction error identification[END_REF] for the constraint on the continuous-time input or replace it by |u[n]| ≤ ūmax for n = 0, ..., 2π ω0Ts . As far as the output amplitude is concerned, we can still say that, in steady state, G(z, θ 0 )u[n] is given by y(t = nT s , A, θ 0 ) with the expression (4). However, for t ̸ = nT s , the noisefree output of the continuous-time system will slightly differ from y(t, A, θ 0 ) (this difference will be small if the gain of the continuous transfer function G 0 is sufficiently small for ω ≥ π Ts ). Consequently, we can still use the results of this paper (i.e., results aiming at satisfying ( 6)) for this alternative instrumentation, but we will then not have a formal guarantee that the output constraint is satisfied for t ̸ = nT s . Note that this was also the case in the approaches of, e.g., [START_REF] Manchester | An algorithm for amplitude-constrained input design for system identification[END_REF][START_REF] Manchester | Input design for system identification via convex relaxation[END_REF].

Remark. Since ȳmax bounds the noise-free output, if we wish to bound the noisy output of the system, ȳmax must be chosen considering the (maximal) amplitude of the disturbance acting on the system (that can be determined using an experiment with u(t, A) = 0). 2.2 Robustification of the output amplitude constraint (6) Since θ 0 is unknown, we cannot use the output amplitude constraint (6) as such. We will instead use a robustified version of this output amplitude constraint. As in, e.g., [START_REF] Bombois | Robust optimal identification experiment design for multisine excitation[END_REF], this robustified version of the constraint will be based on an initial uncertainty region U init for the unknown true parameter vector θ 0 . This initial uncertainty region U init can be obtained via an initial identification experiment that has delivered an initial estimate θinit of θ 0 with covariance matrix P θ,init . Consequently, the following ellipsoid U init is a η%-confidence region for the unknown true parameter vector θ 0 :

U init = θ ∈ R k | (θ -θinit ) T P -1 θ,init (θ -θinit ) ≤ χ (10) with χ such that P r(χ 2 (k) ≤ χ) = η (say 95 %).
For the sequel, it is important to note that U init can be equivalently rewritten as

U init = { θ ∈ R k | θT θ ≤ 1} (11)
if we perform the change of variable

θ ∆ = R(θ -θinit ) (12)
with

P -1 θ,init χ = R T R.
Let us also denote by G(z, θ) the parametrization G(z, θ) where θ has been replaced by θinit +R -1 θ. Similarly, let us define y(t, A, θ) as in ( 4) but with G(z, θ) replaced by G(z, θ). We can then robustify the output amplitude constraint (6) as follows:

-ȳ max ≤ y(t, A, θ) ≤ ȳmax ∀t ∈ T P and ∀ θ ∈ U init (13) For future reference, let us introduce the following quantities:

u(A) ∆ = max t∈T P |u(t, A)| (14) ywc (A) ∆ = max t∈T P max θ∈Uinit |y(t, A, θ)|. ( 15
)
Using these notions, ( 5) and ( 13) are respectively equivalent to

u(A) ≤ ūmax , ( 16
) ywc (A) ≤ ȳmax . ( 17
)
The robustified version of the optimal experiment design problem ( 9) is thus arg max ξ,A ξ s.t. P -1 θ (A, θinit ) ≥ ξ R adm and ( 16) and ( 17).

Note that the accuracy constraint P -1 θ ≥ ξ R adm could also be robustified using U init and the results in [START_REF] Bombois | Robust optimal identification experiment design for multisine excitation[END_REF]. However, since the focus of this paper is on the amplitude constraints, we have decided to replace the unknown θ 0 by its initial estimate θinit to make the accuracy constraint tractable.

Before proceeding, let us make the following assumption on the initial uncertainty U init . Assumption 1 The uncertainty U init defined in [START_REF] Dirkx | A fast smoothingbased algorithm to generate l∞-norm constrained signals for multivariable experiment design[END_REF] contains the unknown parameter vector θ 0 . Moreover, U init is small enough to guarantee that, like G(z, θ 0 ), G(z, θ) is a stable transfer function for all θ ∈ U init .

The first part of Assumption 1 is there to ensure that [START_REF] Larsson | Experimental evaluation of model predictive control with excitation (mpc-x) on an industrial depropanizer[END_REF] effectively implies [START_REF] Bombois | Identification for robust H 2 deconvolution filtering[END_REF]. Note that, in practice, this guarantee will only hold with a probability of η%. The second part of the assumption guarantees that (4) can be written for each θ ∈ U init . We can easily verify whether a given U init satisfies this part of Assumption 1 using the results in [START_REF] Bombois | Robustness analysis tools for an uncertainty set obtained by prediction error identification[END_REF].

Computation of u(A) and ywc (A) 3.1 Introduction

In order to address the optimal experiment design problem [START_REF] Ljung | System Identification: Theory for the User, 2nd Edition[END_REF], we will need to be able to evaluate the quantities u(A) and ywc (A) appearing in the constraints of [START_REF] Ljung | System Identification: Theory for the User, 2nd Edition[END_REF]. Note that, besides their importance for optimal experiment design, these quantities are also crucial if we wish to perform an identification experiment using an excitation (1) with an user-chosen amplitude vector A. These quantities will indeed allow the user to verify whether the chosen A satisfies the (robust) amplitude constraints ( 16) and [START_REF] Larsson | Experimental evaluation of model predictive control with excitation (mpc-x) on an industrial depropanizer[END_REF].

The quantity u(A) is the solution of the following optimization problem:

u(A) = arg minu max umax s.t. -umax ≤ u(t, A) ≤ umax ∀t ∈ TP . (19)
Similarly, the quantity ywc (A) is the solution of the following optimization problem:

ywc(A) = arg miny max ymax s.t. -ymax ≤ y(t, A, θ) ≤ ymax ∀t ∈ TP , ∀ θ ∈ Uinit. ( 20
)
In the sequel, we will show that, using robustness analysis tools, we can exactly compute, for any given A, the solution u(A) of ( 19). We will also show that we can derive a (tight) upper bound yub wc (A) for the solution ywc (A) of [START_REF] Manchester | An algorithm for amplitude-constrained input design for system identification[END_REF]. Note that, if yub wc (A) ≤ ȳmax , we have also that ywc (A) ≤ ȳmax . Consequently, the quantity yub wc (A) can be used together with u(A) to check whether a given excitation u(t, A) can be applied without risk to the considered system G 0 .

Phasors and LFT representations

Let us introduce phasor notations and transform the amplitude vector

A = (a 1,s , a 1,c , a 2,s , a 2,c , ..., a L,c ) T of dimension 2L into the complex vector A = (A 1 , A 2 , ..., A L ) T of dimension L. The i th entry A i of A is defined as A i = a i,c -j a i,s (i = 1, ..., L). (21)
Let us also define U(t) as U(t) = (e jω1t , e jω2t , ..., e jω L t ) T . We have then

u(t, A) = Re(A T U(t)). (22)
For the sequel, it is important to observe that u(t, A) depends on t via the phasors e jωit (i = 1, ..., L). Let us introduce the complex variable

τ ∆ = e jω0t (23)
with ω 0 the fundamental frequency of u(t, A). Using (2), it is clear that, for each ω i , e jωit = τ αi . Consequently, instead of a function of t, we can write u(t, A) as a function u(τ, A) of τ which has the property u(τ, A) = u(t, A) when ( 23) holds. This function u(τ, A) is given by

u(τ, A) = Re(A T U(τ )) (24
) with U(τ ) = (τ α1 , τ α2 , ..., τ α L ) T . When t takes all values in T P , the complex variable τ defined in [START_REF] Megretski | System analysis via integral quadratic constraints[END_REF] takes all values in the set

T = {τ ∈ C | τ * τ = 1}. ( 25
)
Consequently, the constraint in ( 19) is equivalent to

-u max ≤ Re(A T U(τ )) ≤ u max ∀τ ∈ T . (26)
It is also clear that the only quantity dependent on τ in [START_REF] Safonov | Stability and Robustness of Multivariable Feedback Systems[END_REF] i.e., U(τ ) = (τ α1 , τ α2 , ..., τ α L ) T can be rewritten as the following LFT in τ I α L :

p u = (τ I α L ) q u q u U(τ ) =   M 11,U M 12,U M 21,U M 22,U   M U p u 1 ( 27
)
with M 22,U = (0, ..., 0) T , M 12,U = (1, 0, ...., 0) T , M 21,U a matrix whose entries are equal to zero except the entries (i, α i ) for i = 1, ..., L which are equal to one, and with M 11,U a matrix of dimension α L × α L given by M 11,U = 0 0

I α L-1 0
.

Let us now define the quantity y(τ, A, θ) as follows:

y(τ, A, θ) = Re(A T Y(τ, θ)) (28)
with Y(τ, θ) = diag(G(e jω1Ts , θ), ...., G(e jω L Ts , θ)) U(τ ).

(29) It is clear that y(τ, A, θ) = y(t, A, θ) when [START_REF] Megretski | System analysis via integral quadratic constraints[END_REF] holds. Consequently, the constraint in ( 20) is equivalent to

-y max ≤ Re(A T Y(τ, θ)) ≤ y max ∀τ ∈ T , ∀ θ ∈ U init .
(30) Moreover, recall that the model structure G(z, θ) can be written as Z N (z)θ/(1+Z D (z)θ) (see Section 2.1) and recall also the linear change of variable [START_REF] Ebadat | Applications oriented input design in time-domain through cyclic methods[END_REF]. Consequently, diag(G(e jω1Ts , θ), ..., G(e jω L Ts , θ)) is an LFT in I L ⊗ θ (see, e.g., [START_REF] Bombois | Robust optimal identification experiment design for multisine excitation[END_REF]). Using this fact and the fact that U(τ ) is an LFT in τ I α L (see [START_REF] Suzuki | On input design for system identification in time domain[END_REF]), it is possible to find vectors of signals p τ , pθ, q τ , qθ and a matrix M Y such that Y(τ, θ) can be expressed as

p τ pθ = bdiag τ I α L , I L ⊗ θ q τ qθ     q τ qθ Y(τ, θ)     =     M τ τ 11,Y M τ θ 11,Y M τ 12,Y M θτ 11,Y M θ θ 11,Y M θ 12,Y M τ 21,Y M θ 21,Y M 22,Y     M Y     p τ pθ 1     (31)
In the sequel, we will use the symbol ⋆ to introduce shorthand notations of LFTs i.e., the shorthand notation of ( 27) is

U(τ ) = (τ I α L ) ⋆ M U and the one of (31) is Y(τ, θ) = bdiag τ I α L , I L ⊗ θ ⋆ M Y .

Multipliers

In the previous subsection, we have seen that the constraints appearing in [START_REF] Luo | Semidefinite relaxation of quadratic optimization problems[END_REF] and in [START_REF] Manchester | An algorithm for amplitude-constrained input design for system identification[END_REF] can be equivalently rewritten as ( 26) and (30), respectively. These constraints can be both considered as robust constraints since they entail the verification of inequalities for all τ ∈ T in the case of ( 26) and for all τ ∈ T and for all θ ∈ U init in the case of (30). Recall also that, in these inequalities, the quantities U(τ ) and Y(τ, θ) (i.e., the only quantities depending on τ and/or θ) are LFTs, respectively, in τ I α L and in bdiag τ I α L , I L ⊗ θ . In order to find tractable versions of ( 26) and (30), we will need to find appropriate descriptions of the uncertainty τ I α L with τ ∈ T and of the uncertainty

I L ⊗ θ with θ ∈ U init .
These appropriate descriptions are the so-called sets of multipliers for these two uncertainties (we remind the expressions of these sets in Propositions 1 and 2). In a nutshell, a set of multipliers for a given uncertainty p(t) = ∆q(t) with ∆ ∈ ∆ is an explicit and affine parametrization of the quadratic constraint satisfied by the graphs of the signals p and q when ∆ ∈ ∆ [START_REF] Safonov | Stability and Robustness of Multivariable Feedback Systems[END_REF][START_REF] Goh | Robust analysis, sectors and quadratic functionals[END_REF][START_REF] Megretski | System analysis via integral quadratic constraints[END_REF]. Proposition 1 Consider the set T defined in [START_REF] Rojas | Robust optimal experiment design for system identification[END_REF] and an arbitrary scalar α. Then,

τ I α I α * Π α τ I α I α = 0 ∀τ ∈ T (32)
when Π α = bdiag(S, -S) with S any Hermitian matrix of dimension α × α. The set of matrices Π α having this structure will be denoted Π α in the sequel.

Proof. For any Π α = bdiag(S, -S), the quadratic expression in (32) is equal to (τ * τ -1)S and this quantity is indeed equal to zero for all τ ∈ T .

Proposition 2 ([2]) Consider the set U init = { θ ∈ R k | θT θ ≤ 1}
and an arbitrary scalar L. Then,

I L ⊗ θ I L T Σ I L ⊗ θ I L ≥ 0 ∀ θ ∈ U init ( 33
)
when Σ has the following structure

Σ =   -Q ⊗ I k + B + j D P T -j Z T P + j Z Q   ( 34
)
where Q is any positive semi-definite Hermitian matrix of dimension L × L and B, D, P , Z real matrices having the structure given in Appendix A. The set of matrices Σ having the structure (34) will be denoted Σ in the sequel.

In the sequel, we will also require the following result which is a particular case of the generalized KYP lemma [START_REF] Iwasaki | Generalized kyp lemma: Unified frequency domain inequalities with design applications[END_REF]. Lemma 1 ( [START_REF] Iwasaki | Generalized kyp lemma: Unified frequency domain inequalities with design applications[END_REF]) Consider a matrix F(τ ) which depends on a complex scalar variable τ that lies in the set T defined in [START_REF] Rojas | Robust optimal experiment design for system identification[END_REF]. Assume that F(τ ) can be written as an LFT in τ I α and that this LFT is minimal (i.e., we cannot write F(τ ) as an LFT in τ I α with α < α). We thus have:

F(τ ) = (τ I α ) ⋆   F 11 F 12 F 21 F 22  
where F 11 is a matrix of dimension α × α. Consider finally an Hermitian matrix Ξ of appropriate dimension. Then, F * (τ )ΞF(τ ) ≤ 0 for all τ ∈ T is equivalent to the existence of Π α ∈ Π α (see Proposition 1) such that:

  I 0 F 11 F 12   * Πα   I 0 F 11 F 12   + F 21 F 22 * Ξ F 21 F 22 ≤ 0

LMI optimization problems

Using the LFTs and the sets of multipliers introduced in the previous subsections, we can now derive LMI optimization problems allowing to compute u(A) (see Proposition 3) and an upper bound yub wc (A) for ywc (A) (see Proposition 4).

Proposition 3 Consider a given amplitude vector A and its complex expression A (see [START_REF] Manchester | Input design for system identification via convex relaxation[END_REF]). Consider also the LFT expression [START_REF] Suzuki | On input design for system identification in time domain[END_REF] for U(τ ) and the set of multipliers Π α L corresponding to τ I α L (see Proposition 1 with α = α L ). Then, u(A) (see [START_REF] Hjalmarsson | System identification of complex and structured systems[END_REF]) is the solution u opt max of the following LMI optimization problem having as decision variables a real scalar u max ≥ 0 and two matrices

Π 1 α L and Π 2 α L in Π α L :
arg min u max s.t.

V * u,1 Π 1 α L V u,1 + V * u,2 0 0.5 0.5 -u max V u,2 ≤ 0 (35) 0 ≤ -V * u,1 Π 2 α L V u,1 + V * u,2 0 0.5 0.5 u max V u,2 (36)
with

Vu,1 =   I 0 M11,U M12,U   Vu,2 =   A T M21,U A T M22,U 0 1   .
Proof. Let us first prove that the existence of Π 1 α L ∈ Π α L such that (35) holds is equivalent to Re(A T U(τ )) ≤ u max ∀τ ∈ T . For this purpose, note that Re(A T U(τ )) ≤ u max can be rewritten as

Ū * (τ )   0 0.5 0.5 -umax   Ū(τ ) ≤ 0 with Ū(τ ) =   A T U(τ ) 1   .
Note also that, using [START_REF] Suzuki | On input design for system identification in time domain[END_REF], we have that

Ū(τ ) = (τ Iα L ) ⋆     M11,U M12,U A T M21,U A T M22,U 0 1     .
That the existence of Π 1 α L ∈ Π α L such that (35) holds is equivalent to Re(A T U(τ )) ≤ u max ∀τ ∈ T then follows from Lemma 1 (the LFT of Ū(τ ) is minimal since ( 27) is minimal). Using a similar reasoning and Lemma 1, we can also prove that the existence of Π 2 α L ∈ Π α L such that (36) holds is equivalent to -u max ≤ Re(A T U(τ )) ∀τ ∈ T . The proposition is then proven since it is shown, in Section 3, that u(A) is the solution u opt max of the optimization problem consisting in minimizing u max under the constraint [START_REF] Safonov | Stability and Robustness of Multivariable Feedback Systems[END_REF].

Proposition 4 Consider a given amplitude vector A and its complex expression A (see [START_REF] Manchester | Input design for system identification via convex relaxation[END_REF]). Consider the LFT expression (31) for Y(τ, θ), the set of multipliers Π α L corresponding to τ I α L (see Proposition 1 with α = α L ) and the set of multipliers Σ corresponding to I L ⊗ θ (see Proposition 2). Then, an upper bound for ywc (A) (see [START_REF] Iwasaki | Generalized kyp lemma: Unified frequency domain inequalities with design applications[END_REF]) is the solution y opt max of the following LMI optimization problem having as decision variables a real scalar y max ≥ 0, two matrices

Π 1 α L , Π 2 α L in Π α L and two matrices Σ 1 , Σ 2 in Σ:
arg min y max s.t.

V * y,1 Π 1 α L V y,1 +V * y,2 Σ 1 V y,2 + V * y,3 0 0.5 0.5 -y max V y,3 ≤ 0 (37) 0 ≤ -V * y,1 Π 2 α L V y,1 -V * y,2 Σ 2 V y,2 + V * y,3 0 0.5 0.5 y max V y,3 (38)
with

V y,1 =   I 0 0 M τ τ 11,Y M τ θ 11,Y M τ 12,Y   V y,2 =   0 I 0 M θτ 11,Y M θ θ 11,Y M θ 12,Y   V y,3 =   A T M τ 21,Y A T M θ 21,Y A T M 22,Y 0 0 1   .
Proof. Let us first prove that (37) and (38) imply (30). For this purpose, let us consider the signals p τ , pθ, q τ and qθ in the LFT (31) for a fixed τ ∈ T and for a fixed θ ∈ U init . Let us then pre-and post-multiply with (p * τ , p * θ , 1) and (p T τ , p T θ , 1) T the LMI constraints (37) and (38). Using (31), this yields

g * τ Π 1 α L gτ + g * θ Σ 1 gθ + Re(A T Y(τ, θ)) ≤ ymax (39) -ymax ≤ -g * τ Π 2 α L gτ -g * θ Σ 2 gθ + Re(A T Y(τ, θ)) (40)
with g τ = (p T τ , q T τ ) T and gθ = (p T θ , q T θ ) T . The above reasoning can be done for any value of τ ∈ T and for any value of θ ∈ U init . In other words, for the multipliers Π 1 α L , Π 2 α L , Σ 1 and Σ 2 found by the optimization problem, (39) and (40) hold true for all τ ∈ T and for all θ ∈ U init . Observe also that, because of (31),

g τ = τ I α L I α L q τ and gθ = I L ⊗ θ I L qθ
Consequently, due to Propositions 1 and 2, we have that 39) is positive for all τ ∈ T and for all θ ∈ U init . Similarly, -g * τ Π 2 α L g τ -g * θ Σ 2 gθ in (40) is negative for all τ ∈ T and for all θ ∈ U init . We have thus proven that (37) and ( 38) imply (30). The proposition is then proven since it is shown, in Section 3, that ywc (A) is the solution y opt max of the optimization problem consisting in minimizing y max under the constraint (30).

g * τ Π 1 α L g τ + g * θ Σ 1 gθ in (
Proposition 4 makes use of the LFT (31) for Y(τ, θ) and the sets of multipliers for the two uncertainties appearing in this LFT. The sets of multipliers given in Propositions 1 and 2 describe well the two uncertainties when they are considered separately. However, since these two uncertainties are here combined, we will generally have a tighter upper bound for ywc (A) (see, e.g., [START_REF] Bliman | A convex approach to robust stability for linear systems with uncertain scalar parameters[END_REF]) if the multiplier Σ linked to the second uncertainty θ is made dependent on τ i.e., the first uncertainty. The following proposition gives a τ -dependent version of Σ. This proposition is an extension of our earlier results in [START_REF] Bombois | Identification for robust H 2 deconvolution filtering[END_REF]. Proposition 5 Consider the set of multipliers Σ linked to I L ⊗ θ with θ ∈ U init (see Proposition 2). Consider also the variable τ defined in [START_REF] Megretski | System analysis via integral quadratic constraints[END_REF] and that varies in the set T (see [START_REF] Rojas | Robust optimal experiment design for system identification[END_REF]). Using τ and some user-chosen integer b ≥ 1, define the vector B(τ ) = (1, τ, τ 2 , ..., τ b ) T and the matrix Ψ(τ ):

Ψ(τ ) =        B(τ ) ⊗ I kL 0 0 B(τ ) ⊗ IL I kL 0 0 IL        .
Then, for each τ ∈ T , the τ -dependent matrices Ψ * (τ )PΨ(τ ) with P having the structure described below are all elements Σ of Σ. The matrix P has the following structure:

P =       -Λ ⊗ I k 0 0 Λ   P T 21 P21 P22     with Λ =      Q0 • • • Q b . . . 0 0 Q * b 0 0      P 21 =         B0 , ..., Bb + 0, D1 , ..., Db      P 0 ... P b      T -      0 Z 1 Z b      T 0, Z 1 , ..., Z b + P 0 , ..., P b 0         P22 =   j D0 -jZ T 0 jZ0 0   .
The real matrices Bi , Di , P i and Z i (i = 0, ..., b) of the parametrization of P can take any values provided that they have, respectively, the same structure as the real matrices B, D, P and Z in (34). Finally, provided that (B * (τ ) ⊗ I L )Λ(B(τ ) ⊗ I L ) ≥ 0 for all τ ∈ T , the matrix Q 0 can be any Hermitian matrix of dimension L × L and Q i (i = 1, ..., b) any complex matrices of dimension L × L .

When Q i , Bi , Di , P i and Z i are all equal to zero for i = 1, ..., b, Ψ * (τ )PΨ(τ ) becomes τ -independent and is equal to (34) with Q, B, D, P and Z replaced, respectively, by Q 0 , 2 B0 , D0 , 2P 0 and Z 0 .

Proof. See Appendix B.

In the sequel, we will use the τ -dependent multiplier Ψ * (τ )PΨ(τ ) instead of Σ in the constraints of the LMI optimization problem. The property given at the end of Proposition 5 illustrates that the τ -independent multiplier Σ is a special case of the τ -dependent multiplier Ψ * (τ )PΨ(τ ). We therefore expect that, in the vast majority of the cases, the upper bound obtained with this τ -dependent multiplier will be tighter than (or at least equal to) the upper bound obtained with Proposition 4. In order to formulate the LMI optimization problem involving the τ -dependent multiplier, we note that the (minimal) LFT expression of B(τ ) can be deduced similarly as the one of U(τ ) (see Section 3.2) and that the (minimal) LFT expressions of B(τ ) ⊗ I L and of Ψ(τ ) (that will be required in the sequel) can be easily deduced from the one of B(τ ):

B(τ ) ⊗ I L = (τ I bL ) ⋆   B 11 B 12 B 21 B 22   (41) Ψ(τ ) = τ I bL(k+1) ⋆   Ψ 11 Ψ 12 Ψ 21 Ψ 22   . ( 42
)
where B 11 (resp. Ψ 11 ) is a matrix of dimension bL × bL (resp. bL(k + 1) × bL(k + 1)). In order to address the dependence on the uncertainty τ of the multiplier Ψ * (τ )PΨ(τ ), we will need to combine the LFT (42) of Ψ(τ ) and the LFT (31) of Y(τ, θ) (i.e., the LFT that is used in Proposition 4). In particular, we will need to construct the LFT of (Y T Ψ (τ, θ), Y T (τ, θ)) T with

Y Ψ (τ, θ) = Ψ(τ )gθ
where gθ = (p T θ , q T θ ) T with pθ and qθ the internal signals in the LFT (31). That we will need this particular LFT is linked to the fact that, in (39), the multiplier Σ is used in combination with gθ. In order to derive this LFT, let us first observe that (31) gives

gθ = 0 I 0 M θτ 11,Y M θ θ 11,Y M θ 12,Y =G     p τ pθ 1     .
Then, denoting p ψ and q ψ the internal signals of the LFT (42) (i.e., p ψ = τ I bL(k+1) q ψ ), we obtain

    p ψ p τ pθ     = bdiag τ I bL(k+1) , τ I α L , I L ⊗ θ     q ψ q τ qθ              q ψ q τ qθ Y Ψ (τ, θ) Y(τ, θ)          =           M 1 0 M τ τ 11,Y M τ θ 11,Y M τ 12,Y 0 M θτ 11,Y M θ θ 11,Y M θ 12,Y M 2 0 M τ 21,Y M θ 21,Y M 22,Y                  p ψ p τ pθ 1        (43)
with

M 1 = Ψ 11 Ψ 12 G and M 2 = Ψ 21 Ψ 22 G .
We have now all the ingredients to improve the upper bound for ywc (A) using the τ -dependent multiplier Ψ * (τ )PΨ(τ ). Proposition 6 Consider a given amplitude vector A and its complex expression A (see [START_REF] Manchester | Input design for system identification via convex relaxation[END_REF]). Consider also the vector B(τ ) for a given b ≥ 1 and the corresponding matrix Ψ(τ ) (see Proposition 5). Consider finally the LFTs (41) and (43). Define α tot = α L + bL(k + 1).

Then, an upper bound for ywc (A) (see [START_REF] Iwasaki | Generalized kyp lemma: Unified frequency domain inequalities with design applications[END_REF]) is the solution y opt max of the following LMI optimization problem having as decision variables a real scalar y max ≥ 0, two matrices Π 1 αtot , Π 2 αtot in Π αtot (see Proposition 1 with α = α tot ), two matrices P 1 and P 2 having the structure of the matrix P in Proposition 5 and finally two matrices

Π 1 bL , Π 2 bL in Π bL (see Proposition 1 with α = bL):
arg min y max s.t.

V * 1 Π 1 αtot V 1 + M * 2 P 1 M 2 + V * 2 0 0.5 0.5 -y max V 2 ≤ 0 (44) 0 ≤ -V * 1 Π 2 αtot V 1 -M * 2 P 2 M 2 + V * 2 0 0.5 0.5 y max V 2 (45) 0 ≤ -V * 3 Π 1 bL V 3 + V * 4 Λ 1 V 4 (46) 0 ≤ -V * 3 Π 2 bL V 3 + V * 4 Λ 2 V 4
(47) where Λ 1 (resp. Λ 2 ) is the element Λ (see Proposition 5) of the matrix P 1 (resp. P 2 ) and where

V1 =        I 0 0 0 0 I 0 0 M1 0 M τ τ 11,Y M τ θ 11,Y M τ 12,Y        V2 =   0 A T M τ 21,Y A T M θ 21,Y A T M22,Y 0 0 0 1   V3 =   I 0 B11 B12   and V4 = B21 B22 .
Proof. Let us first note that, due to Lemma 1, ( 46) is equivalent to (B * (τ )⊗I L )Λ 1 (B(τ )⊗I L ) ≥ 0 for all τ ∈ T which is a requirement of the parametrization of Proposition 5. The LMI (47) ensures the same property for Λ 2 . Let us then prove that (44) and (45) imply (30). For this purpose, let us consider (43) for a given τ ∈ T and for a given θ ∈ U init and let us consider the corresponding signals p ψ , p τ , pθ, q ψ , q τ and qθ. Let us then pre-and postmultiply with (p * ψ , p * τ , p * θ , 1) and (p T ψ , p T τ , p T θ , 1) T the LMI constraints (44) and (45). Using (43) and the fact that Y Ψ (τ, θ) = Ψ(τ )gθ (gθ = (p T θ , q T θ ) T ), this yields

g * ψτ Π 1 α tot g ψτ +g * θ Ψ * (τ )P 1 Ψ(τ ) gθ +Re(A T Y(τ, θ)) ≤ ymax (48) -ymax ≤ -g * ψτ Π 2 α tot g ψτ -g * θ Ψ * (τ )P 2 Ψ(τ ) gθ+Re(A T Y(τ, θ)) ( 49
)
with g ψτ = (p T ψ , p T τ , q T ψ , q T τ ) T . The above reasoning can be done for any value of τ ∈ T and for any value of θ ∈ U init . In other words, for the multipliers Π 1 αtot , Π 2 αtot , P 1 and P 2 found by the optimization problem, (48) and (49) hold true for all τ ∈ T and for all θ ∈ U init . Observe also that, because of (43),

g ψτ =   τ Iα tot Iα tot   q ψτ and gθ =   IL ⊗ θ IL   qθ
with q ψτ = (q T ψ , q T τ ) T . Consequently, due to Propositions 1, 2 and 5, we have that g * ψτ Π 1 αtot g ψτ + g * θ Ψ * (τ )P 1 Ψ(τ )gθ in ( 48) is positive for all τ ∈ T and for all θ ∈ U init . Similarly, -g * ψτ Π 2 αtot g ψτg * θ Ψ * (τ )P 2 Ψ(τ )gθ in (49) is negative for all τ ∈ T and for all θ ∈ U init . We have thus proven that (44) and ( 45) imply (30). The proposition is then proven since it is shown, in Section 3, that ywc (A) is the solution y opt max of the optimization problem consisting in minimizing y max under the constraint (30).

As already mentioned, Proposition 6 will in the vast majority of the cases yield a tighter upper bound for ywc (A) than Proposition 4 and, if necessary, the conservatism can be further reduced by increasing the value of b. However, it is also clear that the LMI problem of Proposition 6 is more complex than the one of Proposition 4 and this complexity increases with the value of b. Using Proposition 6 will thus take more computation time than using Proposition 4.

Remark. The conservatism of a given upper bound yub

wc (A) for ywc (A) (computed either via Proposition 4 or via Proposition 6) can be evaluated by comparing this upper bound with a lower bound for ywc (A). This lower bound can be computed by choosing a number of grid points in U init and by determining, among those grid points θi , the one (i.e., θwc ) leading to the largest value of max t∈T grid |y(t, A, θi )| (T grid is a fine grid of T P ). The lower bound for ywc (A) is then given by max t∈T grid |y(t, A, θwc )|.

4 Optimal experiment design with robust amplitude constraints 4.1 Integration of the LMI formulations for the amplitude constraints Let us now come back to the robust optimal experiment design problem [START_REF] Ljung | System Identification: Theory for the User, 2nd Edition[END_REF] whose objective is the design of an optimal multisine excitation. This optimization problem contains the constraint [START_REF] Larsson | Experimental evaluation of model predictive control with excitation (mpc-x) on an industrial depropanizer[END_REF] where the quantity ywc (A) cannot be exactly computed. We therefore replace [START_REF] Larsson | Experimental evaluation of model predictive control with excitation (mpc-x) on an industrial depropanizer[END_REF] by the constraint yub wc (A) ≤ ȳmax that implies [START_REF] Larsson | Experimental evaluation of model predictive control with excitation (mpc-x) on an industrial depropanizer[END_REF]. This yields the following tractable optimal experiment design problem:

arg max ξ,A ξ s.t. P -1 θ (A, θinit ) ≥ ξ R adm and u(A) ≤ ūmax and yub wc (A) ≤ ȳmax ( 50
)
where u(A) can be computed via Proposition 3 and yub wc (A) via Proposition 4 or via Proposition 6. The optimization problem (50) is a constrained 4 non-convex optimization problem which can, e.g., be addressed with a SQP algorithm (this algorithm is, e.g., implemented in the Matlab function fmincon). Since (50) is non-convex, the SQP algorithm is not guaranteed to yield the global optimum and the performance of this algorithm will depend on its initialization (that must be close to the global optimum). In Section 5, we will propose two approaches to initialize the SQP algorithm for (50). These two approaches will yield an amplitude vector A app that approaches the solution of the optimization problem (50) (or the original optimization problem ( 18)). 4.2 Unconstrained optimization problem equivalent to (50) In this section, we will show that (50) can be transformed into an unconstrained optimization problem. For this purpose, let us first determine the solution of (50) in the direction of a given vector A u . Proposition 7 Let us consider the optimization problem (50) where we restrict attention to amplitude vectors A that can be written as A = κA u with a positive real scalar κ > 0 and with a given vector A u (e.g., a normalized vector such that A T u A u = 1). Then, the solutions A opt (A u ) and ξ opt (A u ) of this particular optimization problem are given by

A opt (A u ) = κ opt (A u ) A u (51) ξ opt (A u ) = 1
λ max P θ (κ opt (A u ) A u , θinit ) R adm (52
) with κ opt (A) a function defined as follows for any amplitude vector A:

κ opt (A) = min ūmax u(A) , ȳmax yub wc (A) ( 53
)
where min(x, y) is equal to x when x ≤ y and is equal to y otherwise. Proof. Let us consider two amplitude vectors A 1 = κ 1 A u and A 2 = κ 2 A u with κ 1 and κ 2 two scalars such that κ 1 > κ 2 > 0. For these two amplitude vectors, we have that P -1 θ (A 1 , θinit ) > P -1 θ (A 2 , θinit ). This means that the optimal value of the amplitude vector A = κA u will be characterized by the largest value of κ > 0 such that u(κA u ) ≤ ūmax and yub wc (κA u ) ≤ ȳmax . Since u(κA u ) = κu(A u ) and yub wc (κA u ) = κy ub wc (A u ) (see Appendix C), this largest value of κ is clearly equal to κ opt (A u ) (see ( 53)). The expressions (51) and (52) then follow directly.

As indicated in the proof of Proposition 7, κ opt (A u ) (see (53)) is the largest value of κ > 0 such that both u(κA u ) ≤ ūmax and yub wc (κA u ) ≤ ȳmax . Consequently, the scaled vector A opt (A u ) (see (51)) satisfies either u(A opt (A u )) = ūmax or yub wc (A opt (A u )) = ȳmax . Proposition 7 gives the solution of (50) in a given direction A u . In Proposition 8, we will see that the de-termination of the optimal direction can be formulated as an unconstrained optimization problem. For this purpose, we need to parametrize the set of all possible directions. Let us for this purpose introduce the vector φ = (φ 1 , φ 2 , ..., φ 2L-1 ) T containing 2L-1 angles and let us define, based on φ, the following unit vector A u (φ) ∈ R 2L :

Au(φ) =          cos(φ1)
sin(φ1)cos(φ2) ... sin(φ1)...sin(φ2L-2)cos(φ2L-1) sin(φ1)...sin(φ2L-2)sin(φ2L-1)

         . ( 54
)
The vector φ corresponds to the spherical coordinates of the vector A u (φ) and we have that the set [START_REF] Blumenson | A derivation of n-dimensional spherical coordinates[END_REF]. Using this parametrization and Proposition 7, we can now propose the following equivalent formulation of (50). Proposition 8 Consider the function κ opt (A) defined based on the computable quantities u(A) and yub wc (A) (see (53)) and the parametrization A u (φ) of the normalized amplitude vectors A u (see (54)). Then, the solutions A opt and ξ opt of the optimization problem (50) are given by

{A u (φ) | φ ∈ R 2L-1 } is equal 5 to {A u ∈ R 2L | A T u A u = 1}
A opt = κ opt (A u (φ opt )) A u (φ opt ) (55)
ξ opt = 1
λ max P θ (A opt , θinit ) R adm (56)
where φ opt is the solution of the following unconstrained optimization problem:

arg min φ∈R 2L-1 λ max P θ (κ opt (A u (φ)) A u (φ), θinit ) R adm .
(57) Moreover, the optimum A opt given in (55) satisfies either u(A opt ) = ūmax or yub wc (A opt ) = ȳmax . Proof. Each amplitude vector A is entirely characterized by its direction A u (φ) (parametrized by the vector φ) and its norm κ i.e., A(κ, φ) = κA u (φ). Recall that Proposition 7 gives the optimal value of κ for each direction A u (φ). It remains thus to determine the optimal direction (or equivalently the optimal value of φ) i.e., the value of φ yielding the maximal value of ξ opt (A u (φ)) (see ( 52)). This is exactly the objective of (57). The expressions (55) and (56) then follow from an application of Proposition 7 for A u = A u (φ opt ) and the last statement of the proposition is a consequence of the definition (53) of κ opt (A) used to solve the non-convex optimal experiment design problem

Introduction

In order to initialize the algorithm used to solve the non-convex optimal experiment design problem, we propose two approaches yielding amplitude vectors A app close to the solution of the optimal experiment design problem (50) (or equivalently close to the solution of [START_REF] Ljung | System Identification: Theory for the User, 2nd Edition[END_REF] from which (50) follows). 5.2 Semi-definite relaxation inspired from [START_REF] Manchester | Input design for system identification via convex relaxation[END_REF] In [START_REF] Manchester | Input design for system identification via convex relaxation[END_REF], an experiment design problem with amplitude constraints is considered for the case where the excitation signal is parametrized via the N elements u[n] (n = 0, ..., N -1) of the excitation sequence. Let us adapt the framework of [START_REF] Manchester | Input design for system identification via convex relaxation[END_REF] to the (more compact) parametrization of the excitation signal u(t, A) based on the amplitude vector A of dimension 2L (which is generally much smaller than N ). The optimization problem (18) can be reformulated as follows:

arg max ξ,A ξ s.t. L i=1 M (ωi, θinit) a 2 i,s + a 2 i,c + Mv( θinit) ≥ ξ R adm (58)
and

Z T u (t) AA T Zu(t) ≤ ū2 max ∀t ∈ TP and Z T y (t, θ) AA T Zy(t, θ) ≤ ȳ2 max ∀t ∈ TP and ∀ θ ∈ Uinit
where Z u (t) = (sin(ω 1 t), cos(ω 1 t), sin(ω 2 t), ..., cos(ω L t)) T and

Zy(t, θ) =          g(ω1, θ) sin(ω1t + ϕ(ω1, θ)) g(ω1, θ) cos(ω1t + ϕ(ω1, θ)) g(ω2, θ) sin(ω2t + ϕ(ω2, θ)) ... g(ωL, θ) cos(ωLt + ϕ(ωL, θ))         
.

In order to deal with the amplitude constraints that have to be satisfied for all t ∈ T P in (58), we will here for simplicity use a gridding approximation i.e., we will only impose the amplitude constraints at the time instants in a fine grid T grid of T P . A similar gridding approximation will be used (as proposed in [START_REF] Manchester | Input design for system identification via convex relaxation[END_REF]) to address the fact that the output amplitude constraint has to be satisfied for all θ ∈ U init . For this purpose, we will consider a set Θ grid containing a number of grid points θi ∈ U init .

Let us now notice as in [START_REF] Manchester | Input design for system identification via convex relaxation[END_REF] that, in the optimization problem (58), the matrix X = AA T is a positive semidefinite matrix with rank(X) = 1. Let us also notice that a 2 i,s and a 2 i,c (variables appearing in the expression of P θ ) are also entries of the matrix X. Using this fact and the approximations introduced in the previous paragraph, (58) can be reformulated as follows: In (60), we use the notation X i,j for the entry (i, j) of the matrix X. Similarly as was observed in [START_REF] Manchester | Input design for system identification via convex relaxation[END_REF], the optimization problem (59)-( 63) would be an LMI optimization problem in the absence of the rank constraint (63). The relaxation proposed in [START_REF] Manchester | Input design for system identification via convex relaxation[END_REF][START_REF] Luo | Semidefinite relaxation of quadratic optimization problems[END_REF] is to solve the optimization problem by neglecting the rank constraint (63) i.e., to solve the optimization problem (59)-(62). As shown in [START_REF] Luo | Semidefinite relaxation of quadratic optimization problems[END_REF], the optimal solution X opt of this optimization problem allows to derive, via a randomization procedure, an amplitude vector A app that approaches the solution of the original optimization problem (59)-(63). We will here slightly modify this randomization procedure to compensate for the approximations introduced by the use of T grid and Θ grid . This will be achieved by the use of the computable quantities u(A) and yub wc (A) (see Section 3). The randomization procedure is as follows: we generate a number n grid of possible A according to a normal distribution N (0, X opt ) as proposed in [START_REF] Luo | Semidefinite relaxation of quadratic optimization problems[END_REF]. We then apply Proposition 7 to each of these A i.e., for each A, we determine, in the direction of A, the solution of the optimization problem (50) (i.e., the non-convex optimization problem of which A app has to be close to the solution and that is formulated without gridding approximations). This procedure yields n grid values A scaled = κ opt (A)A and n grid values of ξ opt (A) (see ( 52)). The amplitude vector A app approaching the solution of ( 50) is then chosen as the value of A scaled corresponding to the highest value of ξ opt (A).

arg max ξ,X=X T ≥0 ξ (59)
Due to the use of Proposition 7 in the randomization procedure, we have the guarantee that the constraints ( 16) and ( 17) are satisfied for this amplitude vector A app and also that either u(A app ) = ūmax or yub wc (A app ) = ȳmax . 5.3 Use of power constraints instead of amplitude constraints The second approach consists in solving an optimal experiment design problem with power constraints instead of amplitude constraints. The considered power constraints are P(u(t, A)) ≤ ū2 max and P(y(t, A, θ)) < ȳ2 max ∀ θ ∈ U init where P(x(t)) denotes the power of the signal x(t). We therefore consider the following optimal experiment design problem:

arg max ξ,C ξ s.t. L i=1 M (ω i , θN ) c i + M v ( θN ) ≥ ξ R adm and 1 2 L i=1 c i ≤ ū2 max 1 2 L i=1 c i |G(e jωiTs , θ)| 2 ≤ ȳ2 max ∀ θ ∈ U init . ( 64
)
where c i = a 2 i,s + a 2 i,c (i = 1, ..., L) and C = (c 1 , c 2 , ..., c L ) T . Using the results in [START_REF] Bombois | Robust optimal identification experiment design for multisine excitation[END_REF], this optimization problem can be transformed into an LMI optimization problem yielding an optimal solution C opt = (c opt 1 , c opt 2 , ..., c opt L ) T . We can here also use a randomization procedure to determine, based on c opt i (i = 1, ..., L), an amplitude vector A app approaching the solution of the original optimal experiment design problem (50). Let us for this purpose generate a number n grid of possible A satisfying a 2 i,s + a 2 i,c = c opt i (i = 1, ..., L). Then, as in the previous subsection, we apply Proposition 7 to each of these A. This procedure yields n grid values A scaled = κ opt (A)A and n grid values of ξ opt (A) (see ( 52)). The amplitude vector A app approaching the solution of (50) is then here also chosen as the value of A scaled corresponding to the highest value of ξ opt (A).

Like in the previous subsection, the use of Proposition 7 in the randomization procedure gives the guarantee that the constraints ( 16) and ( 17) are satisfied for this amplitude vector A app and also that either u(A app ) = ūmax or yub wc (A app ) = ȳmax .

Numerical illustration

We consider an OE true system (3) with:

G(z, θ 0 ) = θ0,1z -1 +θ0,2z -2 1+θ0,3z -1 +θ0,4z -2 H(z, θ 0 ) = 1
where θ 0 = (θ 0,1 , θ 0,2 , θ 0,3 , θ 0,4 ) T = (0.8988, 0.1034, -0.9723, 0.8385) T and σ 2 e = 1. The sampling rate is T s = 1 s. The amplitude constraint thresholds are ūmax = ȳmax = 1. The initial uncertainty region U init in [START_REF] Dirkx | A fast smoothingbased algorithm to generate l∞-norm constrained signals for multivariable experiment design[END_REF] is characterized by θinit = (0.8 0.01 -0.9854 0.8187) T , χ = 9.49 and We wish to design a multisine excitation u(t, A) (see ( 1)) of fundamental frequency ω 0 = 0.1π rad/s (the period is thus 20 s) and containing L = 3 frequencies i.e., ω 1 = ω 0 , ω 2 = 3ω 0 and ω 3 = 5ω 0 . To design the optimal amplitude vector A of u(t, A), we consider the optimal experiment design problem (50) corresponding to an identification experiment of duration N = 1000 and where R adm is chosen equal to I 4 . Moreover, for the first part of this numerical example, we will use Proposition 4 to compute the upper bound yub wc (A) for ywc (A) (which appears in the constraints of (50)). Since R adm = I 4 , the objective function ξ of (50) is equal to λ min (P -1 θ (A, θinit )) which is clearly a measure of the accuracy of the identified parameter vector θN . Since, in this example, determining the optimal multisine with (57) takes twice less time than with (50), we will from now on consider the equivalent formulation (57) of (50). The gradient-based algorithm (implemented in the Matlab function fminunc) has been here used to solve (57). An amplitude vector A lm can easily be determined from the solution ϕ lm of this algorithm (see the remark at the end of Section 4). For all the initializations discussed in the next paragraph, the gradientbased algorithm yielded amplitude vectors A lm for which the objective function of the optimal experiment design problem 6 i.e., λ min (P -1 θ (A lm , θinit )) is equal to 187.87. These amplitude vectors A lm are thus equivalent from the point-of-view of the considered optimal experiment design problem. This equivalence can be explained by the fact that the corresponding multisines u(t, A lm ) are just time shifted versions of each other. In the sequel, we will therefore only consider one of these A lm i.e., (-0.0688, 0.2212, -0.0662, -0.0120, -0.5075, 0.4621) T .

P -1 θ,init = 10 3        0.
As mentioned in the previous paragraph, we have considered different initializations for the gradient-based algorithm. For this purpose, we have considered the procedures of Sections 5.2 and 5.3 where, in the randomization procedure, the upper bound yub wc (A) for ywc (A) is also computed via Proposition 4. Note furthermore that, in our implementation of the procedure of Section 5.2, T grid contains 201 linearly-spaced points in T P and Θ grid only contains θinit for simplicity. For both approaches, we have run the randomization procedure twice i.e., once with n grid = 100 and once with n grid = 2000. This has led to four amplitude vectors A app . The values of the objective function λ min (P -1 θ (A app , θinit )) for these four A app are given in Table 1. With respect to the four A app in Table 1, the amplitude vector A lm obtained by solving (57) using the gradient-based algorithm improves the objective function with at least 10% since λ min (P -1 θ (A lm , θinit )) = 187.87. In this example, we therefore observe that the initialization procedure gives amplitude vectors A app that are relatively close to A lm , but also that it is useful to improve the result of this initialization procedure with (57). In this example, the robustness of the gradient-based algorithm used to solve (57) is important. To show this, we have also initialized this algorithm with a vector in the direction

A u = (1/ √ 6)(1, 1, 1, 1, 1, 1)
T . This direction is far from optimal since κ opt (A u )A u yields an objective function λ min (P -1 θ (κ opt (A u )A u , θinit )) = 7.92 which is much smaller than the values in Table 1. However, initialized with this A u , the gradient-based algorithm also yields an amplitude vector A lm for which λ min (P -1 θ (A lm , θinit )) = 187.87.

Let us now analyze a bit more the amplitude vector A lm = (-0.0688, 0.2212, -0.0662, -0.0120, -0.5075, 0.4621) T delivered by the gradient-based algorithm used to solve (57). Using Proposition 3, we obtain u(A lm ) = 0.9385 which is indeed the maximal value of |u(t, A lm )| for t ∈ T P . Using Proposition 4, we obtain7 yub wc (A lm ) = 1. As expected (see the remark at the end of Section 4), for A lm , one of the amplitude constraints is active (here the output amplitude constraint). Let us now investigate the conservatism of the upper bound for ywc (A lm ) computed with Proposition 4. For this purpose, we also compute a lower bound for ywc (A lm ) using the procedure given in the remark at the end of Section 3. This leads here to 0.986477 ≤ ywc (A lm ). The conservatism linked to the LMI procedure of Proposition 4 is thus less than 1.4 % (which is small). Let us nevertheless see if a tighter upper bound can be obtained via Proposition 6 with, e.g., b = 1. This is indeed the case since Proposition 6 with b = 1 yields8 yub wc (A lm ) = 0.986550: the conservatism is now less than 0.007 % (and thus negligible).

We have tested Proposition 6 with b = 1 for a number of other A and the observed conservatism is always negligible. Consequently, one can now consider a version of the optimal experiment design problem (50) (or its equivalent (57)) where the upper bound for ywc (A) is obtained using Proposition 6 (b = 1). Due to the reduced conservatism, this new version is closer to the original optimal experiment design problem [START_REF] Ljung | System Identification: Theory for the User, 2nd Edition[END_REF]. The solution of this modified optimal experiment design problem in the direction of A lm (which was the best solution when Proposition 4 was used) can be determined via Proposition 7. In our case, this solution is A lm,bis = κ opt (A lm ) A lm with κ opt (A lm ) = min(1/0.9385, 1/0.986550) ≈ 1.014. With A lm,bis , we obtain λ min (P -1 θ (A lm,bis , θinit )) = 193.02, u(A lm,bis ) = 0.9507 (via Proposition 3) and yub wc (A lm,bis ) = 1 (via Proposition 6 with b = 1). One can wonder whether another direction than the one of A lm could be more optimal. For this purpose, we run the gradient-based algorithm to solve this new version of (57) with an initialization in the direction of A lm . This does not improve the result. Consequently, if we restrict attention to the tools developed in this paper, u(t, A lm,bis ) seems the optimal multisine.

Conclusions

We have developed a methodology to tackle the problem of optimal experiment design with robust amplitude constraints. The developed framework allows to handle these robust amplitude constraints without approximation when the identification is performed in open loop and when the considered true system has one input and one output. In future work, we will extend the framework to closed-loop identification and to multivariable systems.

A Structure of B, D, P , Z in Proposition 2

As shown in [START_REF] Barenthin | Identification for control of multivariable systems: controller validation and experiment design via LMIs[END_REF], the real matrices B and D in the expression of Σ have the following structures: with the constraints that all blocks K il (resp. R il ) satisfy

K il = -K T il ∈ R k×k (resp. R il = -R T il ∈ R k×k
). The matrix P (resp. Z) has a similar structure as B (resp. D), but with the skew-symmetric blocks replaced by row vectors of dimension k [START_REF] Barenthin | Identification for control of multivariable systems: controller validation and experiment design via LMIs[END_REF].

B Proof of Proposition 5

The proof is rather straightforward. When we perform the product Ψ * (τ )PΨ(τ ), we indeed obtain a matrix having the structure (34) for any value of τ ∈ T . More precisely, the matrix Q in (34) is given by (B * (τ )⊗ I L )Λ(B(τ

) ⊗ I L ) = Q 0 + Q 1 τ + Q * 1 τ * + ... + Q b τ + Q * b τ * .
It is clear that this matrix is Hermitian for all τ ∈ T .

The additional condition in the statement of the theorem ensures that this Hermitian matrix is also positive semi-definite for all τ ∈ T . When we perform the product Ψ * (τ )PΨ(τ ), the matrix j D in (34) is given by

C Last steps of the proof of Proposition 7

That u(κA u ) = κu(A u ) follows from the definition (14) of u(A). That ywc (κA u ) = κy wc (A u ) also follows from [START_REF] Iwasaki | Generalized kyp lemma: Unified frequency domain inequalities with design applications[END_REF]. That yub wc (κA u ) = κy ub wc (A u ) can be derived by considering the proof of Proposition 4 and, in particular, the inequalities (39) and (40). Note indeed that the existence of Π 1 α L , Π 2 α L , Σ 1 and Σ 2 such that (39) and (40) hold for all τ ∈ T and for all θ ∈ U init is equivalent to the existence of Π 3 α L , Π 4 α L , Σ 3 and Σ 4 such that the following inequalities hold for all τ ∈ T and for all θ ∈ U init :

g * τ Π 3 α L gτ + g * θ Σ 3 gθ + Re(κ A T Y(τ, θ)) ≤ κ ymax -κ ymax ≤ -g * τ Π 4 α L gτ -g * θ Σ 4 gθ + Re(κ A T Y(τ, θ)).
We can indeed choose Π 3 α L = κΠ 1 α L , Π 4 α L = κΠ 2 α L , Σ 3 = κΣ 1 and Σ 4 = κΣ 2 . An element of Π α (resp. Σ) multiplied by a positive scalar constant is indeed another element of Π α (resp. Σ). A similar reasoning can also be applied to the inequalities (48) and (49) in the proof of Proposition 6.

j

  D = j D0 + (0, D1, ..., Db ) (B ⊗ I Lk ) + (B * ⊗ I Lk )The right hand side of this expression is equal to j D0 + D1 (τ-τ * ) + ... + Db (τ b -(τ b ) * ) since Di = -DT i (i = 1, ...b).It is clear that this expression has indeed the desired structure since τ i -(τ i ) * (i = 1, ...b) is an imaginary number for all τ ∈ T . In Ψ * (τ )PΨ(τ ), the matrix jZ in (34) is given byjZ = jZ0 + (0, Z1, ..., Z b ) (B ⊗ I Lk ) -(B * ⊗ IL)The right hand side of this expression is equal tojZ 0 + Z 1 (τ -τ * ) + ... + Z b (τ b -(τ b ) *) which has thus also the desired structure. Following the same procedure, the matrix B in (34) is given by 2 B0 + B1 (τ + τ * ) + ... + Bb (τ b + (τ b ) * ) since Bi = BT i (i = 1, ...b) while the matrix P in (34) is given by 2P 0 + P 1 (τ + τ * ) + ... + P b (τ b + (τ b ) * ). These expressions have the desired structure since τ i + (τ i ) * (i = 1, ...b) is a real number for all τ ∈ T . From the above expressions, the last statement of the proposition is straightforward.Remark. In Proposition 1 of[START_REF] Bombois | Identification for robust H 2 deconvolution filtering[END_REF], we proposed a less general parametrization for P where the matrices Q i were restricted to be real matrices, the factorizations of j D was reduced to D1 (τ -τ * ) + ... + Db (τ b -(τ b ) * ) and the one of jZ to Z 1 (τ -τ * ) + ... + Z b (τ b -(τ b ) * ).

Table 1

 1

	Objective function λmin(P -1 θ (Aapp, θinit))
	Cases	λmin(P -1 θ (Aapp, θinit))
	Section 5.2 with n grid = 100	159.59
	Section 5.2 with n grid = 2000	171.38
	Section 5.3 with n grid = 100	164.69
	Section 5.3 with n grid = 2000	168.33

Preprint submitted to Automatica

Note that optimizing ai,s and ai,c (i = 1, ..., L) in (1) is the same as optimizing bi and ϕi (i = 1, ..., L) in L i=1 bi sin(ωit + ϕi).

This can, e.g., be achieved by multiplying (1) by a function γ(t) which slowly grows from zero to one.

The smallest covariance matrix P θ is here defined as the one with the smallest value of λmax(P θ ) (or equivalently the largest value of λmin(P -1 θ )).

The first constraint of (50) can be removed by considering as objective function 1/λmax P θ (A, θinit) R adm .

It is in fact sufficient to restrict the angles φi (i = 1, ..., 2L -2) to the interval [0, π] and the angle φ2L-1 to the interval [0, 2π].

Two approaches to initialize the algorithm

We here consider the objective function of (50). The objective function of (57) is its inverse.

Note that the maximal value of |y(t, A lm , θinit)| for t ∈ TP is 0.7425 ( θinit is the center of Uinit).

Note that solving the LMI problem of Proposition 6 with b = 1 takes 16 times more time than solving the one in Proposition 4. This is the reason why Proposition 4 was used in the earlier phase.

The unconstrained non-linear optimization problem (57) can be addressed using a classical gradientdescent algorithm (such algorithm is, e.g., implemented in the Matlab function fminunc). Since (57) is also nonconvex, the gradient-based algorithm is not guaranteed to yield the global optimum and the performance of this algorithm will depend on its initialization (that must be close to the global optimum). Like for (50), we can here also use the procedures of Section 5 for this initialization.

In the numerical simulations that we have performed with Matlab (see Section 6), we have observed that the gradient-based algorithm (Matlab function fminunc) used to solve (57) and the SQP algorithm (Matlab function fmincon) used to solve (50) yield the same optimum when initialized at the same point, but the gradient-based algorithm yields this optimum in less computation time.

Remark. Let us denote by φ lm the (local or global) minimum to which the gradient-based algorithm used to solve (57) converges and let denote by A lm the corresponding amplitude vector i.e., A lm = κ opt (A u (φ lm )) A u (φ lm ) (see ( 55)). Due to (53), the constraints ( 16) and ( 17) are satisfied at A lm and we have that either u(A lm ) = ūmax or yub wc (A lm ) = ȳmax . Note that, due to Assumption 1, the constraint ( 6) is also satisfied for this amplitude vector A lm .