Alexis Bazin
email: alexis.bazin@orange.com

Bruno Jahan
email: bruno.jahan@orange.com

Maryline Hélard
email: maryline.helard@insa-rennes.fr

Impact of the Doppler Effect on the Capacity of Massive MIMO Uplink Systems: OFDM Versus FBMC/OQAM

This paper studies the impact of the Doppler effect on the capacity of an uplink communication between one or several user(s) with mobility and a base station. The OFDM and the FBMC/OQAM modulations are investigated. Moreover, massive MIMO techniques are used at the receiver side to increase the SINR and to allow for simultaneous multiuser communications. An analytic study is carried out to highlight the impact of the Doppler effect on the achievable rate of this system and to bring out benefits and drawbacks of both modulations. On the one hand, the OFDM modulation presents a loss of spectral efficiency because of its cyclic prefix and a limited robustness against the Doppler effect. On the other hand, the FBMC/OQAM performance is restricted by its intrinsic interference due to the frequency selectivity of the channel. These two modulations are compared under different mobility scenarios and the ability of the FBMC/OQAM modulation to significantly increase the achievable rate per user compared to the OFDM modulation is highlighted, especially for high speed scenarios.

I. INTRODUCTION

In [START_REF]5g white paper[END_REF], the NGMN Alliance anticipated the future needs regarding wireless networks, which should be taken into account for the development of the fifth generation of cellular networks (5G). In particular, users with high mobility such as cars and trains are expected to benefit from better quality of service in the coming years. In [START_REF]IMT Vision Framework and overall objectives of the future development of IMT for 2020 and beyond[END_REF], the Radiocommunication Sector of International Telecommunication Union (ITU-R) recommends to be able to communicate with a user up to 500 km/h under vehicle-to-vehicle (V2V) or vehicle-to-everything (V2X) scenarios. Therefore, new techniques and waveforms are investigated in order to overcome these constraints.

The filter-bank multicarrier (FBMC) modulation and in particular the FBMC with offset quadrature amplitude modulation (OQAM) [START_REF] Hirosaki | An Orthogonally Multiplexed QAM System Using the Discrete Fourier Transform[END_REF], [START_REF] Siohan | Analysis and design of OFDM/OQAM systems based on filterbank theory[END_REF] is an interesting candidate thanks to its good spectral and temporal localizations. Indeed, the FBMC/OQAM waveform does not require strict synchronization [START_REF] Fusco | Sensitivity of multi-user filterbank multicarrier systems to synchronization errors[END_REF], presents a good robustness against the Doppler effect as discussed in [START_REF] Haas | A Time-Frequency Well-localized Pulse for Multiple Carrier Transmission[END_REF]- [START_REF] Farhang-Boroujeny | OFDM Versus Filter Bank Multicarrier[END_REF] and thus is particularly suitable for high mobility applications. Furthermore, compared to the orthogonal frequency division multiplexing (OFDM) modulation, the FBMC/OQAM waveform does not need any cyclic prefix (CP) and thus achieves a better spectral efficiency.

In addition, massive MIMO techniques, which use a large number of antennas at the base station side, have gained great interest in the past years [START_REF] Larsson | Massive MIMO for next generation wireless systems[END_REF]. They can be used for uplink multiuser communication to manage multiuser interference and to improve the spectral and energy efficiency [START_REF] Ngo | Energy and Spectral Efficiency of Very Large Multiuser MIMO Systems[END_REF]. Indeed, as stated in [START_REF] Marzetta | Noncooperative Cellular Wireless with Unlimited Numbers of Base Station Antennas[END_REF] and further improved in [START_REF] Ngo | Analysis of the pilot contamination effect in very large multicell multiuser MIMO systems for physical channel models[END_REF], [START_REF] Hoydis | Massive MIMO in the UL/DL of Cellular Networks: How Many Antennas Do We Need?[END_REF], with an unlimited number of antennas, noise and multiuser effects vanished and the only remaining interference is inter-cellular interference (pilot contamination).

Only few studies combine the FBMC waveform and massive MIMO for uplink. In [START_REF] Farhang | Filter Bank Multicarrier for Massive MIMO[END_REF], [START_REF] Aminjavaheri | Frequency spreading equalization in multicarrier massive MIMO[END_REF] the authors associate multiuser massive MIMO methods for uplink and two FBMC techniques, cosine modulated multitone (CMT) and frequency spreading FBMC (FS-FBMC). They highlight the capacity of massive MIMO to reduce the frequency selectivity of the channel and called it "self equalization". According to this results, one may wonder if massive MIMO has also an effect on the time selectivity and therefore on the Doppler effect. Nevertheless, to the best of the authors' knowledge, studies on performance of massive MIMO systems rarely take into account the Doppler effect. Thus, as the FBMC/OQAM modulation is well suited for mobility applications, studying its performance is of great interest.

In this paper we compare the FBMC/OQAM and the OFDM modulations for massive MIMO uplink systems and for users with mobility in a 5G context. The system model is presented in part II. In part III, a theoretical analysis is conducted to compare these two modulations in terms of signal to interference-plus-noise ratio (SINR) and achievable rate per user. Then, in part IV, numerical results allow for conclusion concerning the best modulation as a function of the number of receive antennas and the user mobility.

The following notations are used. The symbol . indicates variables for the FBMC/OQAM modulation. Vectors and matrices are represented in boldface letters. The Euclidean norm of a matrix or vector is denoted by . , ( . ) T indicates their transpose and ( . ) H stands for their conjugate transpose. The identity matrix of size n × n is defined by I n . Finally, the conjugate value of a complex number is denoted by ( . ) * , its expected value by E[ . ], its real part by [ . ] and its imaginary ...

N u user(s)

Base station with N r antennas ...

II. SYSTEM MODEL

For this study, an uplink scenario is considered with simultaneous transmissions between N u users and a base station equipped with N r receive antennas. As in [START_REF] Farhang | Filter Bank Multicarrier for Massive MIMO[END_REF], [START_REF] Aminjavaheri | Frequency spreading equalization in multicarrier massive MIMO[END_REF] and for the sake of simplicity, we assume that all users have a single antenna, are moving at the same speed and are at the same distance from the base station. The scenario is described by Figure 1. The receive antennas are sufficiently spaced so that all propagation channels between a user and a receive antenna are considered to be independent from each other. Moreover, a perfect channel estimation at the receiver side is assumed. At the base station, a MMSE processing is carried out in order to equalize the channel and to separate users' data. Indeed, the MMSE processing performs the best performance according to [START_REF] Farhang | Filter Bank Multicarrier for Massive MIMO[END_REF]. For the rest of this paper, we set M as the number of subcarriers and F 0 as the intercarrier spacing. The symbol duration is therefore T 0 = 1/F 0 and the sampling time is T s = T 0 /M . The index for the time-domain signals and propagation channel paths is k. The propagation channel model has L h paths and the value along time of the l th path between the user n u and the antenna n r is denoted by h nu,nr l [k]. In part II-A, we define the OFDM link and in part II-B the FBMC/OQAM link.

A. OFDM link

The system model for an OFDM link is depicted in Figure 2. After the modulation process, the time domain baseband signal of user n u is:

s nu [k] = 1 √ M M -1 m=0 c nu [m]e j2π m M k , (1)
with c nu [m] the complex data on subcarrier m of variance σ 2 c . A CP is added to counteract the propagation channel delay spread. The received signal on the antenna n r is defined by: transmitted on subcarrier (m 0 + p). The received data on the antenna n r is:

r nr [k] = Nu-1 nu=0 L h -1 l=0 s nu [k -l]h nu,nr l [k] + η nr [k], (2
s n u [k] s N u -1 [k ] s 0 [k ] ... ... h l nu ,n r [k ] r n r [k] r N r -1 [k ] r 0 [k ] ... ... ... ... η 0 [k ] η n r [k] η N r -1 [k ] y 0 y n r y N r -1 ... ... MMSE (W) Channel c 0 n u c 1 n u c -1 n u Frequency
y nr = 1 √ M M -1 k=0 r nr [k]e -j2π m 0 M k = Nu-1 nu=0 p H nu,nr p c nu p + b nr , (3)
with b nr the frequency domain noise component of variance σ 2 b and with:

H nu,nr p = 1 M M -1 k=0 L h -1 l=0 h nu,nr l [k]e j2π p M k e -j2π m 0 +p M l . ( 4
)
For the link between the user n u and the antenna n r , H nu,nr 0 represents the useful component of the channel when receiving the data c nu 0 on subcarrier m 0 and H nu,nr p =0

stand for the components of the channel that lead to intercarrier interference on c nu 0 coming from c nu p =0 . Let us define the channel matrices:

H nu p = H nu,0 p H nu,1 p • • • H nu,Nr-1 p T , H p = H 0 p H 1 p • • • H Nu-1 p . (5)
The received data vector is

y = y 0 y 1 • • • y Nr-1 T and the noise vector is b = b 0 b 1 • • • b Nr-1 T .
The MMSE decoding matrix is defined by:

W = H 0 H H 0 H 0 + ρ -1 I Nu -1 = W 0 W 1 • • • W Nu-1 , (6)
with ρ = σ 2 c /σ 2 b , the signal to noise ratio (SNR) which is the same for all users. The estimated received data after MMSE processing ĉ0 are defined by:

ĉ0 = W H y = ĉ0 0 ĉ1 0 • • • ĉNu-1 0 T . (7)
In a multiuser scenario, a given user n u will experience interference coming from the (N u -1) other users. The indexes Mod.

Dem. Dem. Dem.

â0,0 of this interfering users are defined by n u = n u . Thus, for the user n u , the received estimated data is therefore:

n u C.S.I. s n u [k ] s N u -1 [k] s 0 [k ] ... ... h l nu ,n r [k ] r n r [k ] r N r -1 [k] r 0 [k ] ... ... ... ... η 0 [k ] η n r [k] η N r -1 [k ] ỹ0 ỹn r ỹN r -1 ... ... MMSE ( ) + Real Part Channel a 1,0 n u a -1,0 n u a 1,1 n u a 0,1 n u a -1,1 n u a 1 ,-1 n u a 0 ,-1 n u a -1 ,-1 n u Time W a 0,0 n u Frequency
ĉnu 0 = (W nu ) H H nu 0 c nu 0 + p =0 (W nu ) H H nu p c nu p + n u =nu p (W nu ) H H n u p c n u p + (W nu ) H b, (8)
where the first, second, third and last lines of this equation respectively stand for the useful part, the Doppler interference part, the multiuser interference part and the noise part of the signal.

B. FBMC/OQAM link

In an FBMC/OQAM system, real data are sent every T 0 /2 seconds. Therefore, compared to the OFDM modulation, the orthogonality condition is only valid for the real field. Thus the symbols can be shaped by a large set of prototype filters with an increase of the complexity in return [START_REF] Farhang-Boroujeny | OFDM Versus Filter Bank Multicarrier[END_REF]. We define the discrete time formulation of this filter by g[k] and its length by L g . Figure 3 depicts the system model for an FBMC/OQAM modulation. The baseband transmitted discrete-time FBMC/OQAM signal for the user n u can be expressed as in [START_REF] Siohan | Analysis and design of OFDM/OQAM systems based on filterbank theory[END_REF] by snu [k]:

snu [k] = n∈Z M -1 m=0 a nu [m, n]g m,n [k], (9)
with a nu [m, n] the real data of variance σ 2 a = σ 2 c /2 sent at the frequency and time indexes m and n respectively and g m,n [k] defined by:

g m,n [k] = g [k -nM/2] e j2π m M k- Lg -1 2 e jφ[m,n] . (10)
In [START_REF] Ngo | Energy and Spectral Efficiency of Very Large Multiuser MIMO Systems[END_REF] the phase term φ[m, n] is defined by:

φ[m, n] = (m + n)π/2 + k φ [m, n]π, k φ [m, n] ∈ Z. ( 11
)
Similarly to (2), the received signal on the antenna n r is:

rnr [k] = Nu-1 nu=0 L h -1 l=0 snu [k -l]h nu,nr l [k] + η nr [k]. (12)
For the frequency and time indexes (m 0 , n 0 ), we define a nu p,q = a nu [m 0 + p, n 0 + q] and the received data on the antenna n r after demodulation is:

ỹnr = M -1 k=0 rnr [k]g * m0,n0 [k] = Nu-1 nu=0 p,q
Hnu,nr p,q

a nu p,q + b nr , (13)
with:

           Hnu,nr p,q = L h -1 l=0 hnu,nr p,q
[l]e -j2π m 0 +p M l e -jπ p M (Lg-1) e jΨp,q , hnu,nr

p,q [l] = Lg-1 k=0 g [k -l -qM/2] g[k]h nu,nr l [k]e j2π p M k , Ψ p,q = φ[m 0 + p, n 0 + q] -φ[m 0 , n 0 ] + πpn 0 . (14)
Compared to (4), Hnu,nr (p,q) =(0,0) represent the channel components that lead to intercarrier and intersymbol interference. Indeed, as it does not use any CP, an FBMC/OQAM modulation also suffers from intrinsic interference caused by the frequency selectivity of the channel. We now define the channel matrices:

Hnu p,q = Hnu,0 p,q Hnu,1 p,q • • • Hnu,Nr-1 p,q T , Hp,q = H0 p,q H1 p,q • • • HNu-1 p,q , (15)
and the received data vector ỹ = ỹ0 ỹ1

• • • ỹNr-1 T .
The MMSE decoding matrix is defined by:

W = H0,0 HH 0,0 H0,0 + ρ -1 I Nu -1 = W0 W1 • • • WNu-1 , (16)
so that:

â0,0 = WH ỹ = â0 0,0 â1 0,0 • • • âNu-1 0,0 T . ( 17
)
As only the real part is retained, the noise power is reduced to σ 2 b /2 and the value of the SNR is the same as in an OFDM link, meaning

ρ = σ 2 a /(σ 2 b /2) = σ 2 c /σ 2 b .
For a given user n u , the received estimated data is therefore:

ânu 0,0 = Wnu H Hnu 0,0 a nu 0,0 + (p,q) =(0,0) Wnu H Hnu p,q a nu p,q + n u =nu (p,q) Wnu H Hn u p,q a n u p,q + Wnu H b . (18)
As for [START_REF] Farhang-Boroujeny | OFDM Versus Filter Bank Multicarrier[END_REF], the first, second, third and fourth lines of ( 18) respectively stand for the useful part, the intrinsic and Doppler effect interference part, the multiuser interference part and the noise part of the receive signal.

The impact of the Doppler effect can now be analyzed through a theoretical analysis of the OFDM and FBMC/OQAM performance in the studied context.

III. THEORETICAL PERFORMANCE ANALYSIS

The aim of this part is to highlight the impact of the mobility as well as the benefits and limitations of the OFDM and FBMC/OQAM modulations. Therefore, based on their SINR, we compare the achievable rate per user for both modulations, for a finite and an infinite number of antennas at the base station.

A. Signal to interference-plus-noise ratio

For the OFDM modulation, the SINR γ nu is computed for the user n u by using [START_REF] Farhang-Boroujeny | OFDM Versus Filter Bank Multicarrier[END_REF] and for the FBMC/OQAM modulation, the SINR γnu is obtained by using [START_REF]Base Station (BS) Radio Transmission and Reception[END_REF]:

γ nu = P nu U / (P nu I + P nu M U + P nu N ) , γnu = P nu U / P nu I + P nu M U + P nu N , (19)
with:

     P nu U = (W nu ) H H nu 0 2 σ 2 c , P nu U = Wnu H Hnu 0,0 2 σ 2 a , (20)
defining the useful power,

     P nu M U = n u =nu p (W nu ) H H n u p 2 σ 2 c , P nu M U = n u =nu p,q Wnu H Hn u p,q 2 σ 2 a , (21)
the multiuser interference power and

P nu N = W nu 2 σ 2 b , P nu N = Wnu 2 σ 2 b /2, (22)
the noise power. The intrinsic interference and the Doppler effect interference powers are defined by:

     P nu I = p =0 (W nu ) H H nu p 2 σ 2 c , P nu I = (p,q) =(0,0) Wnu H Hnu p,q 2 σ 2 a , (23)
and are induced by the time and frequency selectivity of the channel. Since for the OFDM modulation the SINR does not take into account the loss of spectral efficiency, the achievable rate per user is a more fair metric and is computed in the next subsection.

B. Achievable rate per user

For the OFDM (resp. FBMC/OQAM) modulation, the achievable rate per user R (resp. R) is defined by:

R = α × E [log 2 (1 + γ nu )] , α < 1 R = E [log 2 (1 + γnu )] , (24)
with α = T 0 / (T 0 + T CP ), T CP being the length of the CP. For a low SNR, a moderate speed and if N r is low enough, the intrinsic interference and the Doppler effect interference can be neglected compared to the noise and multiuser interference and we have:

(P nu M U + P nu N ) P nu I , ( P nu M U + P nu N ) P nu I . (25)
Under this condition, one can see that the problem is simplified and is comparable to the one described in [START_REF] Farhang | Filter Bank Multicarrier for Massive MIMO[END_REF] where the authors showed that the OFDM and FBMC waveforms lead to the same SINR (γ nu = γnu ). Furthermore, according to [START_REF] Marzetta | Noncooperative Cellular Wireless with Unlimited Numbers of Base Station Antennas[END_REF], the SINR for the OFDM modulation increases as the value of N r increases. Therefore, we have R = α × R and the higher the value of N r , the higher the value of R and thus the higher the difference ( R -R). Thereby, the FBMC/OQAM modulation should outperform the OFDM modulation, thanks to the lack of CP. However, above a certain value for N r , (25) are not valid anymore and the performance of the OFDM and FBMC/OQAM modulations are dependent on P nu I and P nu I and thus on the time and frequency selectivity of the channel. Moreover, for a high SNR or under high speed conditions, (25) are never valid even with a moderate number of antennas at the base station. Therefore, in a massive MIMO system with a large number of antennas at the base station and in mobility conditions, the asymptotic achievable rate per user (when N r → +∞) could reveal to be an interesting metric as developed in the next section.

C. Asymptotic values

When N r tends towards infinity, multiuser interference (when n u = n u ) vanishes as stated in [START_REF] Marzetta | Noncooperative Cellular Wireless with Unlimited Numbers of Base Station Antennas[END_REF] and owing to the law of large numbers, we have ∀n u , n r :

                     lim Nr→∞ H n u 0 H H nu 0 = 0, lim Nr→∞ Hn u 0,0 H Hnu 0,0 = 0, lim Nr→∞ H nu 0 2 = N r E |H nu,nr 0 | 2 , lim Nr→∞ Hnu 0,0 2 = N r E Hnu,nr 0,0 2 . ( 26
)
Thus H H 0 H 0 + ρ -1 I Nu and HH 0,0 H0,0 + ρ -1 I Nu become diagonal and we have:

       lim Nr→∞ {W nu } = H nu 0 N r E |H nu,nr 0 | 2 + ρ -1 -1 , lim Nr→∞ Wnu = Hnu 0,0 N r E Hnu,nr 0,0 2 + ρ -1 -1
.

(27) Moreover, as the multiuser and noise components become negligible when N r tends towards infinity [START_REF] Marzetta | Noncooperative Cellular Wireless with Unlimited Numbers of Base Station Antennas[END_REF], the asymptotic achievable rates per user R ∞ and R∞ , for the OFDM and FBMC/OQAM modulations respectively, are defined by:

       R ∞ = α × log 2 1 + (µ 0 ) 2 p =0 (µ p ) 2 +(µ p ) 2 , R∞ = log 2 1 + (μ 0,0 ) 2 (p,q) =(0,0) (μ p,q ) 2 , (28)
with ∀n u , n r :

       µ p 2 = E (H nu,nr 0 ) * H nu,nr p 2 , µ p 2 = E (H nu,nr 0 ) * H nu,nr p 2 , μ p,q 2 = E Hnu,nr 0,0 * Hnu,nr p,q 2 . ( 29
)
R∞ is impacted by the intrinsic interference and by the Doppler effect, while R ∞ is only limited by the Doppler effect. However, as the FBMC/OQAM waveform is less sensitive to the Doppler effect (see [START_REF] Haas | A Time-Frequency Well-localized Pulse for Multiple Carrier Transmission[END_REF]- [START_REF] Farhang-Boroujeny | OFDM Versus Filter Bank Multicarrier[END_REF]), in high mobility conditions we might have R∞ ≥ R ∞ .

To conclude, the FBMC/OQAM modulation can benefit from its better spectral efficiency. However, the Doppler effect has an impact as the performances of both the OFDM and FBMC/OQAM modulations are limited by the time selectivity of the channel. Moreover, the FBMC/OQAM waveform is also impacted by the frequency selectivity of the channel because of the lack of CP. Therefore the choice between the OFDM and FBMC/OQAM modulations is not straightforward in an analytic point of view as there is a trade-off between the time and frequency selectivity of the channel. Thus, simulations are necessary to accurately compare those two waveforms.

IV. SIMULATION RESULTS

In this part, we numerically compare the achievable rate per user for the OFDM and the FBMC/OQAM modulations as a function of the number of users N u , the number of receive antennas at the base station side N r and the speed of the users. The FBMC/OQAM modulation uses the TFL1 (Time Frequency Localization) filter of length T 0 defined in [START_REF] Pinchon | Derivation of analytical expressions for flexible PR low complexity FBMC systems[END_REF], which is designed to combine a small length and a good time-frequency localization in order to counteract the time selectivity of the channel. Classical values for LTE communications are chosen, a FFT size of M = 1024 with 600 active subcarriers, a subcarrier spacing of F 0 = 15 KHz and a CP length of T CP = 4.7µs for the OFDM modulation. Moreover, the central frequency is set to 6 GHz and the SNR to ρ = 10 dB.

A. Connected car scenario

We first set the speed of the users to 50 km/h (urban configuration) and 130 km/h (highway configuration). The Extended Vehicular A (EVA) channel model is used for this simulation as it is suitable for cellular communications in a NLOS environment. Indeed, it describes propagation channel for studying LTE communications (see [START_REF]User Equipment (UE) Radio Transmission and Reception[END_REF], [START_REF]Base Station (BS) Radio Transmission and Reception[END_REF]). Due to user mobility, the coherence time limits the available duration to acquire channel state information (CSI) and thereby limits the number of users too. Thus, we set the number of users to 5 (resp. 2) for 50 km/h (resp. 130 km/h) in order to have less than 10% of the coherence time dedicated to the channel estimation. Under these conditions, the values of R ∞ and R∞ , defined in (28), are estimated through simulations. While selectivity is predominant compared to the time selectivity of the channel and should limit the performance of the FBMC/OQAM modulation with high values of N r . Figure 4 shows the achievable rate per user for the OFDM and the FBMC/OQAM modulations as a function of N r .

Even with a speed of 130 km/h, we have R R ∞ and the Doppler effect has impact neither on the OFDM performance nor on the FBMC/OQAM one. With the FBMC/OQAM modulation, we can consider that the equations in (25) are valid for N r ≤ 20. This leads to a maximum gain of approximately 0.4 b/s/Hz with speed values of 50 km/h and of 130 km/h compared to the OFDM modulation. For higher values of N r , the equations in (25) are not valid anymore and the FBMC/OQAM modulation is limited by its intrinsic interference. The difference ( R -R) thus decreases as N r increases and even becomes negative when N r ≥ 100. Regarding the central frequency, a base station with more than 100 uncorrelated antennas is actually very optimistic. Therefore, the FBMC/OQAM modulation performance is slightly above the OFDM performance for this scenario but other metrics such as complexity or spectral confinement should be considered.

B. High speed train scenario

Then, we study very high speed communications such as train connection with speed values from 150 km/h to the limit defined by [START_REF]IMT Vision Framework and overall objectives of the future development of IMT for 2020 and beyond[END_REF], 500 km/h. We choose a line-of-sight (LOS) channel model which better meets this type of high-speed scenario. The considered channel model called Rural Area is made up with 4 taps (RA4) and has a Rician K factor of 6.7 as defined in [START_REF]Digital land mobile radio communications[END_REF]. With such a mobility, the available time to estimate the propagation channel is very small and only one user can be considered. Moreover, the equations in (25) are not valid. In Figure 5, the achievable rate is computed as a function of the speed for the OFDM and FBMC/OQAM modulations and for N r = 64, N r = 128 and N r = ∞.

With this LOS propagation channel and a speed above 150 km/h, one can see that R∞ > R ∞ . Therefore the FBMC/OQAM modulation always outperforms the OFDM modulation whatever the value of the speed, of ρ and of N r . Moreover, the difference ( R -R) increases with the speed and, when the number of receive antennas tends towards infinity, the maximum value of 1.9 b/s/Hz can be reach. The FBMC/OQAM modulation is thus interesting in this high speed context. For example, with 64 receive antennas and a speed of 250 km/h, the FBMC/OQAM modulation can send 1.2 b/s/Hz more than the OFDM modulation.

V. CONCLUSION This paper presents the impact of the Doppler effect on the capacity of an uplink communication associated with multiuser massive MIMO techniques at the base station side. The OFDM and the FBMC/OQAM modulations are considered.

With a moderate speed, we demonstrate that the lack of CP for the FBMC/OQAM modulation can be an advantage when the number of antennas at the base station is low. However, when this number is higher, the intrinsic interference created by the frequency selectivity of the channel becomes non-negligible and limits the performance of the FBMC/OQAM modulation. Simulation results show that the impact of the Doppler effect is actually not significant even at 130 km/h and compared to the OFDM modulation performance, the FBMC/OQAM modulation performance is equivalent or slightly better in terms of achievable rate per user. Thus, the choice between these two modulations has to be done taking into account other metrics such as complexity or spectral confinement.

In very high mobility conditions, with a speed from 150 km/h to 500 km/h, we demonstrate that the Doppler effect limits both the OFDM and the FBMC/OQAM modulations. However, thanks to its prototype filter (TFL1), the FBMC/OQAM modulation is actually less impacted as confirmed by the simulation results. Indeed, in terms of achievable rate, the OFDM modulation performance is always outperformed by the FBMC/OQAM modulation performance in this context. Therefore, the FBMC/OQAM modulation can significantly improve the capacity of a massive MIMO uplink system in very high mobility conditions. This paper actually paves the way for deeper studies. The system model could be closer from a real system taking into account channel coding and synchronization processes and a real channel estimation for example. Similarly, a more general scenario can be considered with multiple-antenna users having heterogeneous speeds and variable distances to the base station.

Fig. 1 .

 1 Fig. 1. Scenario considered for this paper: an uplink communication between Nu users with mobility and a base station equipped with Nr antennas.

  ) with η nr [k] the noise component on this antenna. Then, a demodulation process is performed on each receive antenna. For a given subcarrier m 0 , we set c nu p = c nu [m 0 + p] the data

Fig. 2 .

 2 Fig. 2. System model for an uplink communication with an OFDM modulation, Nu users and Nr antennas at the base station side.

Fig. 3 .

 3 Fig. 3. System model for an uplink communication with an FBMC/OQAM modulation, Nu users and Nr antennas at the base station side.

Fig. 4 .

 4 Fig.[START_REF] Siohan | Analysis and design of OFDM/OQAM systems based on filterbank theory[END_REF]. Achievable rate per user as a function of the number of receive antennas (Nr) with 5 users at 50 km/h and 2 users at 130 km/h and with the EVA channel model.

  Fig. 5. Achievable rate as a function of the speed with 64 and 128 receive antennas at the base station side and with the RA4 channel model.

	Achievable rate (b/s/Hz)	5 7 8 9 10 11 12 6	R∞	R ∞			Nr=64 -OFDM Nr=64 -FBMC/OQAM Nr=128 -OFDM Nr=128 -FBMC/OQAM
		150		200	250	300	350	400	450	500
						User speed (km/h)


Bernard Picinbono
email: bernard.picinbono@lss.supelec.fr

Fellow

Moments and Polyspectra of the Discrete-Time Random Telegraph Signal

Keywords:

The random telegraph signal (RTS) is an ordered signal which means that the explicit expression of the moment functions requires that the time instants appearing in these moments are put in an increasing order. This property greatly complicates the calculation of higher order polyspectra, Fourier transforms of moments. This correspondence presents an extension to the discrete-time (DT) case of a method already known in continuous time. This extension exhibits new problems which makes the calculation even more difficult. These difficulties are analyzed and the structure of trispectrum is discussed in detail.

I. INTRODUCTION

The continuous-time (CT) random telegraph signal (RTS) is one of the simplest examples of a signal constructed from a Poisson process. It plays an important role in many areas of signal processing. Its correlation function has an exponential shape and its moments of an order higher than two have a very simple mathematical expression [START_REF] Picinbono | Random Signals and Systems[END_REF]. On the other hand, the Fourier transforms (FTs) of these moments, or the polyspectra, are much more difficult to calculate, as seen in [START_REF] Picinbono | Polyspectra of ordered signals[END_REF]. But it is often more convenient to work in the frequency domain than in the time domain. This is especially the case when using linear filtering, because products are simpler to manipulate than convolutions. This property, well known for the inputoutput relationship in a linear filter, is still valid when using higher order statistics and polyspectra are often more interesting to use than time moments. This justifies the interest in the knowledge of polyspectra of various signals and there is a lot of literature on the calculation of such polyspectra of continuous-time signals (see [START_REF] Picinbono | Polyspectra of ordered signals[END_REF]) introducing sometimes various technical difficulties.

The purpose of this correspondence is to study the same problem in the discrete-time (DT) case.

The calculation of the moments is quite similar, and is briefly presented in Section II. On the other hand, the calculation of the polyspectra introduces difficulties which do not appear in the CT case. As a consequence, it is almost impossible to obtain explicitly the general expression of the polyspectra.

However, bispectrum and trispectrum can be calculated, and this is presented in Section III. The main purpose in this section is to stress the difference between the CT and DT cases, avoiding calculations that are simple transpositions from CT to DT. This allows us to understand the difficulties appearing in the determination of polyspectra of an order higher than the one of the trispectrum.

The result of the calculation exhibits some interesting points concerning the trispectrum. The first important result, already verified in the CT, is that the trispectrum introduces a normal spectral density on the normal manifolds of the frequency domain. This property has various practical consequences that are discussed. The second result concerns the density of the trispectrum on the stationary manifold of the frequency domain. Indeed, it is important to verify whether or not this density is bounded, as for example for the trispectrum of linear models. This point requires a specific treatment. Indeed, we shall see in the calculations that it appears in the frequency-domain functions that are not bounded, which can introduce unbounded trispectral densities. However, despite this fact, a more detailed analysis shows that the trispectrum remains bounded. Finally, the correspondence gives the algebraic expression of the trispectrum which can be used in various problems of higher order statistics analysis [START_REF] Brillinger | An introduction to polyspectra[END_REF]

and, as a consequence, the covariance function of x(t) is

γ(τ ) = E[x(t)x(t + τ )] = a |τ | , (2)
with a = 2p -1. This shows that, as in the CT, the correlation function has an exponential shape.

However, as a can take negative values, this correlation function can have an oscillatory behavior, which is impossible in the CT case.

It results immediately from the definition that the third-order moment E[x(t 1 )x(t 2 )x(t 3 )] is zero.

In order to calculate the fourth-order moment, let us introduce the random variable

X = x(t)x(t + k 1 )x(t + k 2 )x(t + k 3 ), (3)
where 0 ≤ k 1 ≤ k 2 ≤ k 3 . We deduce from the previous properties of x(t) and u(t) that

X = u(t + 1)...u(t + k 1 )u(t + k 2 + 1)...u(t + k 3 ). (4)
and this yields

E(X) = a k2 a k3-k2 . (5)
As a consequence, we have

m 4 ({t i }) = E[x(t 1 )....x(t 4 )] = γ(θ 2 -θ 1 )γ(θ 4 -θ 3 ), (6)
where the θ i is are deduced from the t i s by the ordering permutation characterized by θ i ≤ θ i+1 . Thus the explicit expression of m 4 ({t i }) is given in terms of the ordered time instants t i , which characterizes a structure of an ordered signal [START_REF] Picinbono | Polyspectra of ordered signals[END_REF]. It is obvious that (6) can be extended to any moment of an even order and that the moments of odd orders are zero.

This shows that, as in the CT, the RTS is an ordered signal, and one can show that this property also appears for various DT Markov signals [START_REF]Picinbono Polyspectres des signaux markoviens[END_REF]. It is precisely this ordered property which greatly complicates the calculation of the polyspectra.

January 24, 2018 DRAFT

In many problems it is more interesting to work with cumulants than with moments [START_REF] Brillinger | Time Series Analysis, Data Analysis and Theory[END_REF]. At this point it is worth noting that if the expression (6) of the moments is rather simple, this is no longer the case for cumulants, because the relation between moments and cumulants makes use of a complicated transformation (see [7, p. 19]). That is why we shall reduce our calculation to moment polyspectra.

We deduce from ( 6) that fourth-order moments can be expressed in terms of the correlation function.

This also appears in the case of normal (or Gaussian) signals, for which there exists the famous relation

m 4 ({t i }) = γ(t 1 -t 2 )γ(t 3 -t 4 ) +γ(t 1 -t 3 )γ(t 2 -t 4 ) + γ(t 1 -t 4 )γ(t 2 -t 3 ). (7)
valid whatever the order of the instants ti. However, the difference between ( 6) and ( 7) is obvious and this suggests the question of the consequences of ( 6). More precisely, it is interesting to search whether or not there are other signals with moments given by ( 6). This is specified by the following result.

Theorem: Let x(t) be a stationary and ergodic DT signal of unit power and let γ(τ ) be its correlation function. If its fourth-order moment m 4 ({t i }) takes the form f (θ 2 -θ 1 )f (θ 4 -θ 3 ), where the θ i s denote the ordered permutation of the t i s, then x(t) takes only the values ±1 and

f (τ ) = γ(τ ) = a |τ | , |a| < 1. (8)
Proof: Let y(t) = x 2 (t). As x(t) is of unit power E[y(t)] = 1. As the correlation function γ y (.) of y(t) is a fourth-order moment of x(t), it results from the assumption that γ y (θ 2 -θ 1 ) = f 2 (0). This implies that y(t) is a random variable Y . As y(t) is ergodic, Y cannot be random and is equal with probability 1 to a nonrandom constant c. But as E(Y ) = 1, this constant is equal to 1 and Y = 1 with probability 1. This implies that x(t) = ±1, and therefore f 2 (0) = 1. As f (.) appears only with its product, we can assume that f (0) = 1. As a result

E[x(θ 1 )x(θ 2 )x 2 (θ 3 )] = E[x(θ 1 )x(θ 2 )] (9)
or

f (θ 2 -θ 1 ) = γ(θ 2 -θ 1 ). (10)
Similarly,

E[x(θ 1 )x 2 (θ 2 )x(θ 3 )] = E[x(θ 1 )x(θ 3 )] (11)
or

γ(θ 2 -θ 1 )γ(θ 3 -θ 2 ) = γ(θ 3 -θ 1 ). (12)
This implies that γ(τ ) = a |τ | and as x(t) is stationary, |a| < 1.

DRAFT January 24, 2018

This result is interesting because it shows that the structure of the moment ( 6) is strongly connected with the exponential correlation function.

III. TRISPECTRUM OF THE RTS

As the moments of odd order are zero, the bispectrum is also zero. The trispectrum is deduced from the fourth-order spectral moment, FT of m 4 ({t i }), defined by

M (f ) = t1 t2 t3 t4 m 4 (t) exp(-2πjf T t), (13)
where m 4 (t) is defined by ( 6), f is the frequency vector with components f i and t tt is the time vector with components t i .

It is in the calculation of this expression that the main differences between the CT and the DT cases appear. This is already true when calculating simply the power spectrum of an exponential correlation function. In the CT we have

Γ(f ) = +∞ -∞ exp(-a|τ |) exp(-2πjf τ )dτ, (14)
and the calculation can be reduced to an integral from 0 to1, with its complex conjugate. In the DT case we have

Γ(f ) = +∞ -∞ a |τ | exp(-2πjf τ ), τ integer, (15)
and in order to use two complex conjugate series, it is necessary to suppress the term corresponding to τ = 0.

The same phenomenon, though much more complicated, appears for (13), and we must distinguish between terms with distinct t i s and terms where some t i s are equal. This introduces five kinds of terms.

a) The term A corresponds to the sum (13) where all the t i s are distinct.

b) In this case there is a sum with respect to only three distinct t i s. For example, if t 1 = t 2 , the term B 12 corresponds in (13) to the sum with

[x 2 (t 1 )x(t 3 )x(t 4 )] = γ(t 3 -t 4 ).
There are, of course, six distinct terms B ij such as B 12 .

c) This corresponds to two uncoupled instants such that t 1 = t 2 , t 3 = t 4 . The corresponding moment in ( 13) is

E[x 2 (t 1 )x 2 (t 3 )] = 1,
and there are three distinct terms C ij .

January 24, 2018 DRAFT d) The terms D i correspond to two coupled time instants such that t 1 = t 2 = t 3 , t 4 . The corresponding moment in ( 13) is

E[x 3 (t 1 )x(t 4 )] = γ(t 1 -t 4 )
and there are four terms D i . e) Finally, there is one element of the sum (13) corresponding to t 1 = t 2 = t 3 = t 4 with [x 4 (t 1 )] = 1.

It now remains to calculate the sum of all these terms.

A. Distinct Time Instants

The problem is to calculate the sum (13) where all the tis are distinct and where the moment m 4 (t)

is given by ( 6). The principal difficulty in this calculation is that ( 6) is valid for ordered time instants, while (13) is independent of this ordering property.

To overcome this difficulty, the idea is to decompose (13) with distinct t i s in a sum of 24 terms corresponding to the 24 distinct permutations of the t i s. The sum associated with a given permutation can be expressed as

S P ({µ i }) = +∞ θ1=-∞ +∞ θ2=θ1+1 +∞ θ3=θ2+1 +∞ θ4=θ3+1 γ(θ 4 -θ 3 )γ(θ 2 -θ 1 ) exp -2πj i µ i θ i , (16)
where the frequencies µ i are deduced from the frequencies f i by the permutation P transforming the t i s into the θ i s. An abbreviated calculation presented in the Appendix yields

S P ({µ i }) = Γ + (µ 4 )W (µ 3 + µ 4 )Γ + (µ 2 + µ 3 + µ 4 ) ×δ(µ 1 + µ 2 + µ 3 + µ 4 ). ( 17
)
In this expression, Γ 4 (µ) is defined by

Γ + (f ) = ∞ k=1 γ(k) exp(-2πjf k), (18)
and since γ(k) = a |k| ,

Γ + (f ) = a exp(2πjf ) -a . (19)
Furthermore

W (f ) = ∞ k=1 exp(-2πjf k), (20)
W (f ) = (1/2)[V (f ) + δ(f ) -1], (21)
where δ(.) is, as in (17), the Dirac distribution, and

V (f ) = -j cot(πf ), (22)
is the FT of the signal v[k]. The function V (f ) is analogous to the principal value distribution defined as the FT of v(t) = Sign(t), where t corresponds to the CT. Thus V (f ) defined by ( 22) is also a distribution, which means that any integral containing V (f ) must be considered as a Cauchy principal value integral.

Note finally that V (f ) is an odd distribution, or that V (f ) = -V (-f ).

By using these definitions, (17) can be expressed in a form that is more convenient for what follows, or

S P ({µ i }) = T 1 -T 2 -T 3 , (23)
with

T 1 = (1/2)Γ + (µ 2 )Γ + (µ 4 )δ(µ 1 + µ 2 )δ(µ 3 + µ 4 ), T 2 = (1/2)V (µ 1 + µ 2 )Γ + (-µ 1 )Γ + (-µ 1 -µ 2 -µ 3 )δ(µ 1 + µ 2 + µ 3 + µ 4 ) T 3 = (1/2)Γ + (-µ 1 )Γ + (-µ 1 -µ 2 -µ 3 )δ(µ 1 + µ 2 + µ 3 + µ 4 ).
The two terms T 1 and T 2 are the same as those obtained in the CT (see [2, eq. ( 48)]).

Finally, the contribution of the sum over distinct t i 's in (13) is obtained by making the sum of 24 terms such that (23) corresponds to the 24 permutations of the frequencies f i '. This will be done later.

B. Terms B ij

The term B 12 is defined by

B 12 (f ) = t1 t2 t3 γ(t 4 -t 3 ) × exp{-2πj[(f 1 + f 2 )t 1 + f 3 t 3 + f 4 t 4 ]}, (24)
in which the t i s are distinct. The short calculations presented in the Appendix show that

B 12 (f ) = [Γ(f 4 ) -1]δ(f 1 + f 2 )δ(f 3 + f 4 ) -[Γ(f 3 ) + Γ(f 4 ) -2]δ( f i ), (25)
January 24, 2018 DRAFT wherre Γ(f ) is the power spectrum of the signal defined by

Γ(f ) = γ(k) exp(-2πjkf ). (26)
There are six distinct such terms: B 12 , B 13 , B 14 , B 23 , B 24 , B 34 .

C. Terms C ij

The term C 12 is defined by

C 12 (f ) = t1 t3 =t1 exp{-2πj[(f 1 + f 2 )t 1 + (f 3 + f 4 )t 2 ]} (27)
and its value is, of course,

D 12 (f ) = δ(f 1 + f 2 )δ(f 3 + f 4 ) -δ( f i ). (28)
There are three different such terms ; C 12 , C 13 , C 14 .

D. Terms D i

The term D 4 is defined by

D 4 (f ) = t1 t4 =t1 γ(t 1 -t 4 ) exp{-2πj[(f 1 + f 2 + f 3 )t 1 + f 4 t 4 ]} (29)
and its value is

D 4 (f ) = [Γ(f 4 ) -1]δ( f i ). (30)
There are four distinct such terms :

D 1 , D 2 , D 3 ; D 4 .

E. Terms E

This term corresponding to equal t i s is, of course,

E(f ) = δ( f i ). ( 31
)
and there is only one such term.

The trispectrum is obtained by the sum of all the previous terms. In order to better understand its structure, it is necessary to make this sum by conveniently grouping the terms with similar forms.

DRAFT January 24, 2018

IV. STRUCTURE OF THE TRISPECTRUM All the previous equations show that the fourth-order spectral moment M (f ) given by ( 13) is zero outside the stationary manifold defined by f 1 + f 2 + f 3 + f 4 = 0. This is obviously a consequence of the stationarity of the RTS x(t), and can be expressed as

M 4 (f ) = Γ 3 (f 1 , f 2 , f 3 )δ(f 1 + f 2 + f 3 + f 4 ), (32)
which defines the trispectrum Γ 3 (f 1 , f 2 , f 3 ) as the density of the spectral moment on the stationary manifold. However, this density is not necessarily bounded because some terms of M 4 (f ) contain two Dirac terms δ(.), which introduces distributions on submanifolds of the stationary manifold. In our case, these submanifolds are only the normal manifolds (see [START_REF] Picinbono | Random Signals and Systems[END_REF] and [START_REF] Picinbono | Polyspectra of ordered signals[END_REF]) and we shall first calculate the corresponding density.

A. Distribution on the Normal Manifolds

The first normal manifold (FNM) is defined by

f 1 + f 2 = 0 ; f 3 + f 4 = 0, (33)
and there are two other distinct similar manifolds defined by replacing the first equation of (33) either by f 1 + f 3 = 0 or by f 1 + f 4 = 0. By using the results obained above a simple calcuation shows that the contribution to M 4 (f ) of the FNM is

T 1N M = Γ(f 1 )Γ(f 3 )δ(f 1 + f 2 )δ(f 3 + f 4 ), (34)
where Γ(f ) is

Γ(f ) = 1 -a 2 1 + a 2 -2a cos(2πf ) . ( 35
)
This shows that the density on the FNM is simply the product Γ(f 1 )Γ(f 3 ) of power spectra, and this is called a normal density. The same result is, of course, valid for the other two normal manifolds. This is a very important point and it is related to various properties of the RTS such as ergodicity [START_REF] Picinbono | Ergodicity and fourth-order spectral moments[END_REF] and asymptotic normality (the Central Limit Theorem) [START_REF] Picinbono | Polyspectra of ordered signals[END_REF]. Furthermore, it has a direct consequence on the calculation of fourth-order spectral cumulants. Indeed, this cumulant is obtained by subtracting from M 4 (f ) the contribution of normal manifolds (NMs) with normal density given by (34). As a consequence, there is no contribution of the NMs to the spectral cumulant, which is one of the most important advantages of cumulants with respect to moments.

It is worth pointing out that the same result is valid for the continuous-time RTS, but the calculation is much simpler, because the terms B, C, and D do not appear in the CT case.

January 24, 2018 DRAFT

  is the FT of the signal w[k] = 1 for k > 0 and zero otherwise. Let v[k] be the signal Sign of k, equal to +1 for k > 0 and even. It is obvious that w[k] = (1/2){v[k] + 1 -δ[k]}. By Fourier transformation this yields

  [START_REF] Dalle Molle | Trispectral analysis of stationary random times series[END_REF].

	II. MOMENTS OF THE RTS
	The RTS is a real stationary signal x(t), t integer, taking only the values ±1 and satisfying the
	recursion x(t) = x(t 1 )u(t), where u(t) is a strictly white noise, or a sequence of independent and
	identically distributed (i.i.d.) random variables taking only the values +1 and -1 with probabilities p
	and 1 -p, respectively. As a consequence, x 2k (t) = 1, and the same relations are valid for u(t). The
	stationarity assumption implies that E[x(t)] = 0. It results from the recursion that
	x(t)x(t + τ ) = x 2 (t)u(t + 1)u(t + 2)....u(t + τ ),
	DRAFT	January 24, 2018

B. Bounded Density on the Stationary Manifold

Following the same procedure as for the normal manifolds we find that the bounded density on the stationary manifold can be expressed as

where

This sum is completely different from the one obtained in the CT case, and this is quite natural. This fact already appears for the power spectrum which is given by (35) in the DT case while the FT of the

Finally it results from our previous calculations that the trispectrum remains bounded on the stationary manifold of the frequency domain. This situation also appears in the CT domain, but with entirely different values of the trispectrum.

V. CONCLUSION

It results from the previous discussion that the fourth-order spectral moment (13) of the DT random telegraph signal is distributed on the normal manifold with a normal density, product of power spectra.

The remaining part of the spectral moment is bounded on the stationary manifold. It is clear that the main difficulty of this calculation is due to the necessity to distinguish in (13) the terms with distinct tis and terms where some tis are equal. The calculation remains possible for the trispectrum, but it appears extremely tedious for polyspectra of orders higher than 3, and this difficulty does not appear for the continuous time, where polyspectra of any order of the random telegraph signal have been calculated [START_REF] Picinbono | Polyspectra of ordered signals[END_REF].


Ngoc Anh Nguyen

Sorin Olaru

Robust output feedback control: convex lifting based approach

This paper presents a design method for robust output feedback control of linear, discrete-time, invariant systems affected by both state and output additive disturbances. The method relies on the so-called convex lifting, which is defined on the N -step controllable set. It is proven that the proposed method guarantees the recursive feasibility and robust stability in the sense that the closed loop converges to a given robust positively invariant set as time tends to infinity. Moreover, the method only requires the resolution of a linear programming problem at each sampling instant. Finally, a numerical example is considered to illustrate the methodology.

I. INTRODUCTION

Robust control is still an active topic in control theory, although extensive studies have been dedicated to over several decades, see among the others [START_REF] Morari | Robust process control[END_REF], [START_REF] Zhou | Essentials of robust control[END_REF]. In these robust control problems, full knowledge of the state is required, however, in practice measurement is usually affected by external noises. Therefore, suitable observer is required to estimate the state and to integrate in control design. As a consequence, control strategy becomes output feedback, since it is a function of outputs. Among the available methods, a basic idea is to transform a robust output feedback control problem into a robust state feedback counterpart, see [START_REF] Findeisen | State and output feedback nonlinear model predictive control: An overview[END_REF]. Accordingly, robust control methodologies can be applied. Note that if the underlying system is linear and subject to linear constraints, a linear Luenberger observer can be deployed. On this direction, model predictive control (MPC) strategies were proposed in [START_REF] Mayne | Robust output feedback model predictive control of constrained linear systems[END_REF], [START_REF] Sui | Robust output feedback model predictive control for linear systems via moving horizon estimation[END_REF]. A common point of MPC strategies is that suitable terminal constraints, describing a robust positively invariant set, are imposed at the end of the prediction horizon, leading to robust stability [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF]. However, this also leads to an exponential increase of the number of constraints along the prediction horizon and subsequently to a demanding online evaluation.

To reduce the computational complexity, set-theoretic methods can also be of use by determining a priori a suitable control Lyapunov function, as employed in [START_REF] Blanchini | Ultimate boundedness control for uncertain discretetime systems via set-induced Lyapunov functions[END_REF]. This method only requires to solve a linear program at each sampling instant. However, it also requires the contractivity of the feasible region over which this control Lyapunov function is defined. The computation of such a contractive set is usually based on an iterative procedure, known to be demanding. Recall that the vertex control presented in [START_REF] Gutman | An algorithm to find maximal state constraint sets for discrete-time linear dynamical systems with bounded controls and states[END_REF] also requires the contractivity of the feasible region. As an extension, the Ngoc Anh Nguyen, Sorin Olaru are with the Laboratory of Signals and Systems (L2S, UMR CNRS 8506), CentraleSupélec-CNRS-UPS, Université Paris Saclay, 3 rue Joliot Curie, Plateau de Moulon, 91192, Gif-sur-Yvette, France. Ngocanh.Nguyen.rs@gmail.com, Sorin.Olaru@centralesupelec.fr interpolation based approach was put forward in [START_REF] Nguyen | Constrained control of uncertain, time-varying linear discrete-time systems subject to bounded disturbances[END_REF], in which the interpolation gain serves as a control Lyapunov function. Accordingly, robust stability is ensured by the nonincrease of this gain along the dynamics. Moreover, the socalled convex lifting approach has recently been presented in [START_REF]Robust control design based on convex liftings[END_REF], [START_REF] Nguyen | Convex liftings-based robust control design[END_REF]. Unlike the above set-theoretic methods, this concept is not a control Lyapunov function, however the contractivity of the feasible region is still required.

This paper presents an alternative method, relying on a suitable convex lifting as an extension of the method in [START_REF]Robust control design based on convex liftings[END_REF], [START_REF] Nguyen | Convex liftings-based robust control design[END_REF]. More precisely, this convex lifting is defined over an N -step controllable set, known not to be contractive. Similar to these set-theoretic methods, terminal constraints are not imposed. However, unlike a control Lyapunov function, this function is equal to 0 over a given robust positively invariant set and strictly positive outside this region. In addition, this function is shown to be strictly decreasing along the closed-loop dynamics outside this invariant set, leading to its convergence to 0 as time tends to infinity. Accordingly, robust stability is guaranteed by means of the convergence of the closed-loop dynamics to a given robust positively invariant set. Meanwhile, this method also requires to solve a linear program at each sampling instant, it thus could be useful for systems with fast dynamics, e.g. cantilever beam system [START_REF] Gulan | Embedded linear model predictive control for 8-bit microcontrollers via convex lifting[END_REF].

Nomenclature

Throughout the paper, R, N, N >0 denote the field of real numbers, the set of nonnegative integers and the positive integer set, respectively. The index set I N is also defined as I N = {1, 2, . . . , N } with respect to a given N ∈ N >0 . A polyhedron is defined as the intersection of finitely many closed halfspaces. A polytope is defined as a bounded polyhedron. Also, V(P ) denotes the set of vertices of polytope P. We use int(S) to denote the interior of a full-dimensional set S and conv(S) denotes the convex hull of S. Also, given a set S and a matrix A of suitable dimension, we use AS to denote the following set:

AS = {Ax : x ∈ S} . Given two sets S 1 , S 2 ⊂ R d , by S 1 \S 2 , we denote the follow- ing set: S 1 \S 2 := x ∈ R d : x ∈ S 1 , x / ∈ S 2 .
Also, the Minkowski sum of these two sets, denoted by S 1 ⊕S 2 , is defined as follows:

S 1 ⊕ S 2 := {x 1 + x 2 : x 1 ∈ S 1 , x 2 ∈ S 2 } .
Finally, the Potryagin difference S 1 S 2 is defined by: S 1 S 2 := {x : x ⊕ S 2 ⊆ S 1 } .

II. PROBLEM SETTINGS

In this paper, we consider a linear time-invariant system, affected by both additive state and output disturbances:

x k+1 = Ax k + Bu k + w k y k = Cx k + v k (1)
where x k , u k denote the state and control variable at time k, while w k , v k represent the additive state and output disturbances at time k, respectively. These variables and disturbances satisfy:

x k ∈ X, u k ∈ U, w k ∈ W, v k ∈ V, (2)
where the constraint sets X ⊂ R dx , U ⊂ R du , W ⊂ R dx , V ⊂ R dy are assumed to be polytopes, containing the origin in their interior, with given d x , d u , d y ∈ N >0 . This system is assumed to satisfy the following assumption. Assumption 1: The pair (A, B) is controllable and the pair (A, C) is observable. This paper aims to design robust output feedback controller which is able to cope with the state and output disturbances, and to stabilize system (1), while guaranteeing the constraints in [START_REF]Positively invariant polyhedral sets of discrete-time linear systems[END_REF]. To tackle this problem, a common approach is to make use of a Luenberger observer. More precisely, this observer is written in the following form:

xk+1 = Ax k + Bu k + L(y k -ŷk ), ŷk = C xk . (3)
If one defines e k = x k -xk , then we obtain the following autonomous system, according to (1) and ( 3),

e k+1 = (A -LC)e k + w k -Lv k . (4)
Note that since w k ∈ W and v k ∈ V, then w k -Lv k , considered as additive disturbances of system (4), is also bounded in W ⊕ (-LV), i.e., w k -Lv k ∈ W ⊕ (-LV). Therefore, if one chooses a gain L such that A -LC is strictly stable, then it is proven in [START_REF] Gilbert | Linear systems with state and control constraints: The theory and application of maximal output admissible sets[END_REF] that a robust positively invariant set exists (this concept will be recalled in Sect. III). Such a gain L can be determined by means of the Riccati equation for the system (A T , C T ).

For ease of presentation, let Ω e denote a robust positively invariant set for the autonomous system (4). Accordingly, for any initial error state e 0 ∈ Ω e , then e k ∈ Ω e for all k ∈ N >0 . The observer system (3) can be written in the following form:

xk+1 = Ax k + Bu k + Lv k + LCe k , ŷk = C xk . (5)
We remark that if the initial error state e 0 ∈ Ω e , then Lv k + LCe k is also considered as additive disturbance for the observer system [START_REF] Blanchini | Set-theoretic methods in control[END_REF]. Furthermore, such an additive disturbance is also bounded in LV⊕LCΩ e . Note that X, V, W are assumed to be polytopes, therefore Ω e is bounded, leading to the boundedness of LV ⊕ LCΩ e . If we compute a polytopic robust positively invariant set Ω e for the autonomous system (4), then we can ensure that the additive disturbance Lv k + LCe k is bounded in polytope LV⊕LCΩ e . Also, it can easily be observed that 0 ∈ int(LV⊕LCΩ e ), since 0 ∈ int(V) and 0 ∈ int(Ω e ). As a consequence, we now return to a robust control design problem for the linear system (5) affected by bounded additive state disturbances, whose variables and parameters are bounded as below:

xk ∈ X Ω e , u k ∈ U, Lv k + LCe k ∈ LV ⊕ LCΩ e , (6)
under the assumption that X Ω e is non-empty. In the sequel, we will present a robust control design method for system [START_REF] Blanchini | Set-theoretic methods in control[END_REF], which relies on a suitable convex lifting. Unlike the counterpart presented in [START_REF]Robust control design based on convex liftings[END_REF], [START_REF] Nguyen | Convex liftings-based robust control design[END_REF], a convex lifting employed in this paper is defined over the N -step controllable set, usually known not to be contractive. Note that the contractivity property is always required in most of existing methods related to piecewise linear control Lyapunov functions [START_REF] Blanchini | Ultimate boundedness control for uncertain discretetime systems via set-induced Lyapunov functions[END_REF], [START_REF] Gutman | An algorithm to find maximal state constraint sets for discrete-time linear dynamical systems with bounded controls and states[END_REF], [START_REF] Nguyen | Implicit improved vertex control for uncertain, time-varying linear discretetime systems with state and control constraints[END_REF]. It is worth emphasizing that robust stability in this paper is not in the sense of Lyapunov. Instead, we will prove that such a constructive convex lifting is strictly decreasing along the closed-loop dynamics outside a given robust positively invariant set and consequently the closed loop of the observer system ( 5) is convergent to this set, leading to robust stability of system (1).

III. CONSTRUCTION OF SUITABLE CONVEX LIFTING

Before going to the main result, we recall some important concepts which are of use later in the proposed control design method. Positive invariance concept has been investigated in many studies [START_REF] Bitsoris | On the positive invariance of polyhedral sets for discretetime systems[END_REF]- [START_REF] Bitsoris | Constrained regulation of linear systems[END_REF], [START_REF] Blanchini | Set-theoretic methods in control[END_REF] and deployed in different control design strategies. In case the underlying system is affected by disturbances, the robust positive invariance concept is of use instead.

Definition 3.1: Given an admissible control law u = K x ∈ U, a set Ω x ⊆ X Ω e is called robust positively invariant with respect to system (5) subject to constraints [START_REF] Findeisen | State and output feedback nonlinear model predictive control: An overview[END_REF]

iff (A + BK)Ω x ⊕ LV ⊕ LCΩ e ⊆ Ω x.
Similar to determining an observer gain, such a local controller u = K x ∈ U in Definition 3.1 can also be obtained via the Riccati equation for the pair (A, B). According to this local controller, one can use existing algorithms to compute a robust positively invariant set Ω x, see for instance [START_REF] Gilbert | Linear systems with state and control constraints: The theory and application of maximal output admissible sets[END_REF], [START_REF] Kolmanovsky | Theory and computation of disturbance invariant sets for discrete-time linear systems[END_REF], [START_REF] Nguyen | Implicit improved vertex control for uncertain, time-varying linear discretetime systems with state and control constraints[END_REF], [START_REF] Rakovic | Invariant approximations of the minimal robust positively invariant set[END_REF]. Hereafter, Ω x is assumed to be a full-dimensional polytope.

Another important concept is the feasible region. Unlike the one deployed in [START_REF] Nguyen | Convex liftings-based robust control design[END_REF], in this paper we make use of the N -step controllable set as the feasible region with respect to a given N ∈ N >0 . Its definition is recalled below. Definition 3.2: Consider system (5) subject to constraints [START_REF] Findeisen | State and output feedback nonlinear model predictive control: An overview[END_REF]. Let a robust positively invariant set Ω x and N ∈ N >0 be given. A set denoted by K N (Ω x) ⊆ X Ω e is called the N -step controllable set if any point, belonging to this set, can reach Ω x in N steps, while staying inside X Ω e despite any disturbances in LV ⊕ LCΩ e , i.e.,

K 0 (Ω x) = Ω x, K N (Ω x) = x0 ∈ X Ω e : ∃u 0 , . . . , u N -1 ∈ U s.t. xN ∈ Ω x, ∀v i ∈ V, e i ∈ Ω e , i ∈ {0} ∪ I N -1 .
For brevity, the computation of K N (Ω x) is referred to [START_REF] Kerrigan | Robust constraint satisfaction: Invariant sets and predictive control[END_REF] for further details. We now present the construction of a suitable function according to the above ingredients, which will be deployed in the proposed control design. For ease of presentation, let N (x) denote such a function defined over K N (Ω x). According to its construction presented in the sequel, N (x) was proven in [START_REF] Nguyen | Convex lifting: Theory and control applications[END_REF], [START_REF] Nguyen | Convex liftings-based robust control design[END_REF] to be a convex lifting, therefore, we will call this function a convex lifting throughout the rest of the paper. For simplicity, the formal definition of this concept and related ones are referred to [START_REF] Nguyen | Any discontinuous PWA function is optimal solution to a parametric linear programming problem[END_REF], [START_REF]Inverse parametric linear/quadratic programming problem for continuous PWA functions defined on polyhedral partitions of polyhedra[END_REF], [START_REF] Nguyen | Constructive solution of inverse parametric linear/quadratic programming problems[END_REF]. This function should be convex, nonnegative, equal to 0 over Ω x, strictly positive over

K N (Ω x)\Ω x. Furthermore, it should satisfy N (v) > max x∈Ki-1(Ω x) N (x) for all v ∈ V(K i (Ω x))\K i-1 (Ω x).
All these requirements are of use later to prove the convergence of N (x) to 0 along the trajectories, as it will be clear later. To construct N (x), we start with the following initial point, for h 0 = 0 and a given constant h 1 > 0:

V 0 := v T h 0 T : v ∈ V(Ω x) ⊂ R dx+1 V 1 := v T h 1 T : v ∈ V(K 1 (Ω x)) ∪ V 0 Π 1 := conv(V 1 ) 1 (x) := min z z s.t. xT z T ∈ Π 1 . (7)
The above function 1 (x) represents the following properties. Lemma 3.1: Function 1 (x) is continuous, convex and piecewise affine. The proof follows as a direct consequence of Theorems IV-3 and IV-4 in [START_REF]Postoptimal analyses, parametric programming and related topics[END_REF]. The following important properties of 1 (x) are of use later.

Lemma 3.2: Function 1 (x) satisfies:

1) 0 ≤ 1 (x) ≤ h 1 for all x ∈ K 1 (Ω x); 2) 1 (x) = 0 for x ∈ Ω x; 3) 1 (x) > 0 for x ∈ K 1 (Ω x)\Ω x; 4) 1 (x) = h 1 for x ∈ V(K 1 (Ω x))\Ω x.
Proof: Any point xT z T ∈ Π 1 can be expressed as:

α 0 (v), α 1 (v) ≥ 0, v∈V(Ω x ) α 0 (v) + v∈V(K1(Ω x )) α 1 (v) = 1, x z = v∈V(Ω x ) α 0 (v) v 0 + v∈V(K1(Ω x)) α 1 (v) v h 1 . (8)
Accordingly, one can easily observe that

0 ≤ z = v∈V(Ω x ) α 0 (v)0 + v∈V(K1(Ω x )) α 1 (v)h 1 ≤ h 1 ,
leading to claim 1). Note that the left-hand inequality becomes equality only if α 1 (v) = 0 for all v ∈ V(K 1 (Ω x)), leading to x ∈ Ω x, as stated in claim 2). Moreover, for any point x ∈ K 1 (Ω x)\Ω x expressed above, there exists at least one vertex v ∈ V(K 1 (Ω x)) such that α 1 (v) > 0, leading to z > 0 and subsequently to claim 3). Finally, since 1 (x) is convex over K 1 (Ω x), this function attains its maximal value h 1 at vertices of K 1 (Ω x). Note however that any x ∈ V(K 1 (Ω x)) ∩ Ω x can be expressed as in [START_REF] Gilbert | Linear systems with state and control constraints: The theory and application of maximal output admissible sets[END_REF] with

α 1 (v) = 0 for all v ∈ V(K 1 (Ω x)), yielding 1 (x) = 0 < h 1 .
In other words,

1 (x) = h 1 for x ∈ V(K 1 (Ω x))\Ω x.
In order to construct a convex lifting N (x) as required above, we construct intermediate functions i (x) for i ∈ N >0 such that N (x) = i (x) for all x ∈ K i (Ω x) and for all i ∈ I N . Starting from [START_REF]Postoptimal analyses, parametric programming and related topics[END_REF], the construction of N (x) is presented step-by-step in the sequel. Suppose function i (x), associated with the polytopic partition

X (i) j j∈I M (i) of K i (Ω x)
, is written in the form:

i (x) = (a (i) j ) T x + b (i) j for x ∈ X (i) j . (9)
We also define the following function:

σ i (x) = max j∈I M (i) (a (i) j ) T x + b (i) j for x ∈ R dx . (10)
Accordingly, i+1 (x) is constructed as follows, with respect to a given constant > 0:

h i+1 := min h h s.t. σ i (v) + ≤ h, ∀v ∈ V(K i+1 (Ω x))\K i (Ω x) V i+1 := v T h i+1 T : v ∈ V(K i+1 (Ω x)) ∪ V i Π i+1 := conv(V i+1 ) i+1 (x) := min z z s.t. xT z T ∈ Π i+1 . (11)
The constructions ( 9)-( 11) are repeated until i = N -1.

Similar to 1 (x) shown in Lemma 3.1, i+1 (x) for all i ∈ I N -1 possess the following properties: Lemma 3.3: Functions i+1 (x) are convex, continuous, piecewise affine, for all i ∈ I N -1 .

Lemma 3.4:

σ i (x) = i (x) for all x ∈ K i (Ω x), ∀i ∈ I N .
Proof: Since i (x) is a convex, piecewise affine function, proven in Lemma 3.3, it can be written as follows:

i (x) = max j∈I M (i) (a (i) j ) T x + b (i) j for x ∈ K i (Ω x).
Accordingly, the representation of σ i (x) in (10) leads to σ i (x) = i (x) for all x ∈ K i (Ω x). We now prove that function N (x) satisfies the aforementioned requirements.

Lemma 3.5: For i ∈ I N , we obtain 1)

0 < h i < h i+1 2) i (x) = h i for x ∈ V(K i (Ω x))\K i-1 (Ω x) 3) 0 ≤ i (x) ≤ h i 4) i+1 (x) = i (x) for x ∈ K i (Ω x).
Proof: We will prove claims 1), 2) and 3) in the same time. As proven in Lemma 3.2, for i = 1, claims 2) and 3) hold true and 0 < h 1 , we now prove that h 1 < h 2 . In fact, as σ 1 (x) is convex over K 2 (Ω x), it thus attains its maximal value at vertices of K 2 (Ω x). Accordingly, there exists a

vertex v ∈ V(K 2 (Ω x))\K 1 (Ω x) such that σ 1 (v) ≥ σ 1 (x) for all x ∈ K 2 (Ω x). Therefore, if one chooses x ∈ V(K 1 (Ω x)), then (11) yields h 2 ≥ + σ 1 (x) = + 1 (x) = + h 1 > h 1 .
To prove claim 3) for i = 2, we consider xT z T ∈ Π 2 , this point can be described as below

α j (v) ≥ 0, 2 j=0 v∈V(Kj (Ω x )) α j (v) = 1 xT z T = 2 j=0 v∈V(Kj (Ω x)) α j (v) v T h j T . ( 12
)
As a consequence of claim 1) for i = 1, we obtain

0 ≤ z = 2 j=0 v∈V(Kj (Ω x )) α j (v)h j ≤ h 2 ,
leading to claim 3) for i = 2.

To prove claim 2) holds true for i = 2, consider any point xT z T ∈ Π 2 , this can be written as in [START_REF] Kerrigan | Robust constraint satisfaction: Invariant sets and predictive control[END_REF]. Note that 2 (x) is convex over K 2 (Ω x), therefore 2 (x) attains its maximal value h 2 at vertices of K 2 (Ω x). However, any point x ∈ V(K 2 (Ω x)) ∩ K 1 (Ω x) can be expressed as in form [START_REF] Kerrigan | Robust constraint satisfaction: Invariant sets and predictive control[END_REF] with

α 2 (v) = 0, ∀v ∈ V(K 2 (Ω x)), leading to 2 (x) ≤ h 1 < h 2 . Roughly speaking, 2 (x) = h 2 for x ∈ V(K 2 (Ω x))\K 1 (Ω x).
The proof of claims 1), 2) and 3) for N ≥ i ≥ 3 follows the same arguments as deployed above.

For claim 4), we first prove that it holds for i = 1. In fact, as proven in claim 1) that h 2 ≥ h 1 + , then the construction in [START_REF] Herceg | Multi-parametric toolbox 3.0[END_REF]

leads to h 2 ≥ σ 1 (v) + for all v ∈ V(K 2 (Ω x)).
Accordingly, any point xT z T ∈ Π 2 described as in [START_REF] Kerrigan | Robust constraint satisfaction: Invariant sets and predictive control[END_REF] yields:

z = 2 j=0 v∈V(Kj (Ω x )) α j (v)h j (13a) ≥ 1 j=0 v∈V(Kj (Ω x )) α j (v)σ 1 (v) + v∈V(K2(Ω x )) α 2 (v)(σ 1 (v) + ) (13b) ≥ σ 1 (x) + v∈V(K2(Ω x)) α 2 (v) ≥ σ 1 (x). (13c)
Note that inclusion (13c) is obtained due to the convexity of σ 1 (x). Also, the right-hand inequality of (13c) becomes equality only if α 2 (v) = 0 for all v ∈ V(K 2 (Ω x)), leading to x ∈ K 1 (Ω x). Therefore, 2 (x) = σ 1 (x) = 1 (x) for all x ∈ K 1 (Ω x). Suppose claim 4) holds until i, leading to h j ≥ σ i (v) for all v ∈ V(K j (Ω x)) and all j ∈ I i . We prove that it holds for i + 1. In fact, following the same argument of use to prove h 2 ≥ h 1 + , we can obtain h i+1 ≥ + h i . Accordingly, any point xT z T ∈ Π i+1 can be written in the following form:

α j (v) ≥ 0, i+1 j=0 v∈V(Kj (Ω x )) α j (v) = 1 xT z T = i+1 j=0 v∈V(Kj (Ω x )) α j (v) v T h j T . ( 14
)
As a consequence, the definition of h i+1 in (11) yields:

z = i+1 j=0 v∈V(Kj (Ω x )) α j (v)h j (15a) ≥ i+1 j=0 v∈V(Kj (Ω x )) α j (v)σ i (v) + v∈V(Ki+1(Ω x)) α i+1 (v) ≥ σ i (x) + v∈V(Ki+1(Ω x )) α i+1 (v) ≥ σ i (x). ( 15b
)
Note that the left-hand inequality of (15b) is due to the convexity of σ i (x). The equality holds true only if

α i+1 (v) = 0 for all v ∈ V(K i+1 (Ω x)), leading to x ∈ K i (Ω x). In other words, i+1 (x) = σ i (x) = i (x) for x ∈ K i (Ω x).

IV. ROBUST OUTPUT FEEDBACK CONTROL DESIGN

This section presents the main result of this paper, i.e., a control design procedure based on the convex lifting N (x), constructed in the preceding section. This procedure is summarized in Algorithm 1.

Algorithm 1 Control design procedure

Input: A local controller u = K x associated with Ω x, N (x) and K N (Ω x). Output: optimal controller u * (x k ) at each instant.

1: Compute N (x k ). 2: If xk ∈ Ω x, then u * (x k ) = K xk . 3: Else Solve the following problem: γ * (u * k ) T T = arg min γ, u k γ s.t. N (Ax k + Bu k + w) ≤ γ N (x k ), (Ax k + Bu k ) ⊕ (LV ⊕ LCΩ e ) ⊆ K N (Ω x), u k ∈ U, γ ≥ 0, ∀w ∈ V(LV ⊕ LCΩ e ). (16) 4
: u * (x k ) = u * k . 5: End 6: k ← k + 1. Return to step 1.
We now need to prove that the controller designed in Algorithm 1 can both guarantee the recursive feasibility and closed-loop stability of system (5), while still satisfying constraints [START_REF] Findeisen | State and output feedback nonlinear model predictive control: An overview[END_REF]. For ease of presentation, we first prove that the above controller makes function N (x k ) strictly decreasing along the closed-loop dynamics outside Ω x. This observation is formally stated in the following result.

Proposition 4.1: Consider controller designed in Algorithm 1 and function N (x), then for any x ∈ K N (Ω x)\Ω x, it follows that

N (Ax + Bu * (x) + w) < N (x) for all w ∈ LV ⊕ LCΩ e .
Proof: First, we will prove that any v ∈

i∈I N V(K i (Ω x))\Ω x satisfies N (Av + Bu * (v) + w) < N (v) for all w ∈ LV ⊕ LCΩ e .
In fact, for ease of presentation, define the following index, with respect to

x ∈ K N (Ω x), i * (x) = min j∈I N ∪{0} j s.t. x ∈ K j (Ω x).
Obviously, for an x ∈ K N (Ω x), the definition of i * (x) is unique. According to claim 2) of Lemma 3.5, for any

v ∈ i∈I N V(K i (Ω x))\Ω x, it yields N (v) = h i * (v) > 0, i * (v) ≥ 1. (17)
Also, there exists a control law u

(v) ∈ U such that Av + Bu(v) + w ∈ K i * (v)-1 (Ω x), leading to N (Av +Bu(v)+w) ≤ h i * (v)-1 , ∀w ∈ LV⊕LCΩ e , (18)
as shown in claim 3) of Lemma 3.5. As a consequence of claim 1) in Lemma 3.5, inclusions [START_REF] Nguyen | Implicit improved vertex control for uncertain, time-varying linear discretetime systems with state and control constraints[END_REF] and ( 18) lead to

N (Av + Bu(v) + w) ≤ h i * (v)-1 < h i * (v) = N (v). (19)
Otherwise, if v ∈ V(Ω x), then there also exists a controller u(v) = Kv such that Av + Bu(v) + w ∈ Ω x for all w ∈ LV ⊕ LCΩ e , because Ω x is a robust positively invariant set. Therefore, claim 4) of Lemma 3.5 and claim 2) of Lemma 3.2 lead to

N (Av + Bu(v) + w) = N (v) = 0. (20)
Consider any point x ∈ K N (Ω x)\Ω x, there exists a region

X (N ) j in the polytopic partition X (N ) j j∈I M (N ) of K N (Ω x)
such that x ∈ X (N ) j

. Accordingly, x can be expressed as:

x = v∈V(X (N ) j ) α(v)v, α(v) ≥ 0, v∈V(X (N ) j ) α(v) = 1. (21)
This yields

N (x) = v∈V(X (N ) j ) α(v) N (v) (22a) > v∈V(X (N ) j ) α(v) N (Av + Bu(v) + w) (22b)
≥ N (Ax + B v∈V(X (N ) j ) α(v)u(v) + w) (22c) ≥ N (Ax + Bu * (x) + w). (22d)
Note that equality (22a) holds because N (x) is affine over X (N ) j

. Since x ∈ K N (Ω x)\Ω x, there exists at least one vertex v ∈ V(X (N ) j

)\Ω x such that α(v) > 0. Accordingly, inclusions [START_REF] Nguyen | Convex lifting: Theory and control applications[END_REF] and ( 20) lead to the strict inequality (22b). Moreover, inclusion (22c) holds due to the convexity of N (x). Finally, inequality (22d) is obtained because of u(v) ∈ U for all v ∈ V(X (N ) j

) and the convexity of U. The proof is complete.

Proposition 4.2: Consider controller designed in Algorithm 1 and function N (x), then for any x ∈ Ω x, it yields

N (Ax + Bu * (x) + w) = N (x) = 0, ∀w ∈ LV ⊕ LCΩ e .
Proof: Due to the robust positive invariance of Ω x, for any x ∈ Ω x, it yields for all w ∈ LV ⊕ LCΩ e , Ax + Bu * (x) + w = (A + BK)x + w ∈ Ω x. [START_REF] Nguyen | Constructive solution of inverse parametric linear/quadratic programming problems[END_REF] Accordingly, claim 2) of Lemma 3.2 leads to

N (Ax + Bu * (x) + w) = 0, ∀w ∈ LV ⊕ LCΩ e . ( 24
)
The proof is complete.

The main results of the paper are presented in the sequel. Theorem 4.3: Consider system (5) subject to constraints [START_REF] Findeisen | State and output feedback nonlinear model predictive control: An overview[END_REF]. Then, for any initial point x0 ∈ K N (Ω x), controller designed in Algorithm 1 guarantees the recursive feasibility and robust stability.

Proof: The recursive feasibility is guaranteed by suitable γ < 1 according to Propositions 4.1 and 4.2. Moreover, Proposition 4.1 shows that N (x k ) is strictly decreasing along the closed-loop dynamics outside Ω x, it is also bounded in [0, h N ] , therefore { N (x k )} ∞ k=0 is convergent to 0. As a consequence, xk converges to Ω x as time tends to infinity, leading to robust stability of system (5).

Theorem 4.4: Consider system (1) subject to constraints [START_REF]Positively invariant polyhedral sets of discrete-time linear systems[END_REF]. If the initial points e 0 ∈ Ω e and x0 ∈ K N (Ω x), then controller designed in Algorithm 1 also guarantees the recursive feasibility and robust stability.

Proof:

As Ω e is a robust positively invariant set of system (4), then e k ∈ Ω e , when e 0 ∈ Ω e . Also, as shown in Theorem 4.3, controller designed in Algorithm 1 also guarantees the recursive feasibility for system (5), therefore xk ∈ X Ω e . As a consequence, x k = e k + xk ∈ X. Moreover, this controller satisfies the constraints in (2), leading to the recursive feasibility for system [START_REF] Bitsoris | On the positive invariance of polyhedral sets for discretetime systems[END_REF].

As for robust stability of system (1) and controller designed in Algorithm 1, it can be observed that xk converges to Ω x, while e k is bounded in Ω e , therefore x k = e k + xk is convergent to Ω x ⊕ Ω e as time tends to infinity, leading to its robust stability. Recall that the vertex control in [START_REF] Gutman | An algorithm to find maximal state constraint sets for discrete-time linear dynamical systems with bounded controls and states[END_REF] and the one in [START_REF] Blanchini | Ultimate boundedness control for uncertain discretetime systems via set-induced Lyapunov functions[END_REF] cannot guarantee that the closed-loop dynamics converge to a given robust positively invariant set, i.e., the result in Proposition 4.1. Note also that the construction of N (x) becomes more demanding when the dimension of the state space increases, since it requires vertex enumeration. However, as this construction is performed offline, it is reasonable to assume that powerful computational resources are available.

V. NUMERICAL EXAMPLE

In this section, a numerical example is considered to illustrate the proposed robust output feedback control design procedure. To this end, a DC-DC converter model is employed and given below:

x k+1 = 1 0.0075 -0.143 0.996 x k + 4.798 0.115 u k + w k y k = -0.5 1 x k + v k . (25)
Also, the state, control variable and disturbances satisfy:

-[10 10] T ≤ x k ≤ [3 3] T , u k ∞ ≤ 0.5, w k ∞ ≤ 0.01, v k ∞ ≤ 0.005.
We choose the observer gain L = [-0.5640 0.4838] and a local controller u = [-0.2313 0.1781] x. Accordingly, the minimal robust positively invariant set Ω e is computed, using the algorithm presented in [START_REF] Rakovic | Invariant approximations of the minimal robust positively invariant set[END_REF], to enlarge the feasible region for the observer system (3). Also, Ω x is chosen as the maximal output admissible set associated with the local controller u = [-0.2313 0.1781] x. The set Ω x ⊕Ω e and the 20-step controllable set K 20 (Ω x) ⊕ Ω e are both represented in Fig. 1. Subsequently, a convex lifting 20 (x) is constructed following the construction in [START_REF]Postoptimal analyses, parametric programming and related topics[END_REF] and [START_REF] Herceg | Multi-parametric toolbox 3.0[END_REF] with h 1 = 10 -2 and = 10 -3 . This convex lifting is presented in Fig. 2. Also, it is shown in this figure the evolution of 20 (x k ) along the closed-loop dynamics of system [START_REF] Blanchini | Set-theoretic methods in control[END_REF]. It verifies that the convex lifting 20 (x k ) is strictly decreasing along the dynamics outside Ω x, leading to the convergence of 20 (x k ) to 0. Roughly speaking, x k converges to Ω x ⊕ Ω e as shown more clearly in Fig. 1. Finally, the numerical example of this paper is simulated in the environment of MPT 3.0 [START_REF] Herceg | Multi-parametric toolbox 3.0[END_REF].

VI. CONCLUSIONS

This paper presented a method for robust output feedback control of linear systems affected by both state and output additive disturbances. The method made use of a suitable convex lifting, which was shown to be strictly decreasing along the closed-loop dynamics outside a given robust positively invariant set. Accordingly, the closed-loop dynamics were shown to be convergent to this invariant set, leading to robust stability. The method was shown to only require solving a linear program at each sampling instant. Finally, a numerical example was considered to illustrate the method.

K 20 ( 1 x 2 Fig. 1 .Fig. 2 .

 201212 Fig. 1. Maximal output admissible set Ω x ⊕ Ωe, the 20-step controllable set K 20 (Ω x) ⊕ Ωe and the closed-loop dynamics of system (1).