domain: spi phys math info

text:
Sovanna Thai
email: <sovanna.thai@onera.fr

Clément Roos
email: clement.roos@onera.fr

Jean-Marc Biannic
email: jean-marc.biannic>@onera.fr

Probabilistic µ-Analysis for Stability and H ∞ Performance Verification

Keywords: probabilistic µ-analysis, (skewed) structured singular value, stability margin, H∞ performance level

In this paper, new µ-based algorithms are proposed in the field of probabilistic robustness analysis. The objective is to compute tight bounds on the probability for a parametrically uncertain and possibly high order system to meet some stability and performance criteria. In this approach, uncertain parameters are treated as random variables with given probability distributions. Internal stability is treated first and H∞ performance in the scalar case is considered next. The main contribution is to provide tight bounds on the probability measure in the latter case. The proposed algorithms are based on branch-and-bound techniques and tested on several benchmarks from the literature, demonstrating their efficiency and the potential of the probabilistic setting in reducing the conservatism of µ-analysis. They have been integrated in the SMART Robustness Analysis Library of the SMAC Toolbox developed by ONERA (http://w3.onera.fr/smac).

I. INTRODUCTION

µ-analysis is a standard approach to robustness analysis, typically used to identify worst-case configurations and to compute guaranteed robustness margins (see [START_REF] Doyle | Analysis of feedback systems with structured uncertainties[END_REF][START_REF] Safonov | Stability margins of diagonally perturbed multivariable feedback systems[END_REF]). Despite its recognised maturity by the control community, it often leads to pessimistic results and consequently to unnecessarily conservative designs, at the cost of some performance loss. Probabilistic methods are therefore typically preferred (see e.g. [START_REF] Chen | Risk analysis in robust control -making the case for probabilistic robust control[END_REF][START_REF] Calafiore | A survey of randomized algorithms for control synthesis and performance verification[END_REF]), with Monte Carlo simulations representing the industrial standard. But these approaches have major drawbacks: the non-guaranteed detection of worst-case configurations may lead to erroneous conclusions, and the required sample size to reach satisfactory levels of confidence and accuracy can be prohibitive.

These observations seem to indicate that the two approaches are in fact complementary: deterministic methods could be used as a first validation step, leaving simulations to be performed only on the potentially problematic configurations. But in view of the above, it seems necessary to improve the interface between the two families of methods. This is the purpose of probabilistic µ-analysis, which consists of applying µ-based techniques while considering the uncertainties as random variables with known probability distributions. It guarantees the detection of worst-case scenarios, but it also quantifies their probability of occurrence, allowing comparison with a specified threshold.

Building up on the approach of [START_REF] Khatri | Guaranteed bounds for probabilistic µ[END_REF], [START_REF] Balas | Analysis of an uav flight control system using probabilistic mu[END_REF] and [START_REF] Falcoz | Probabilistic µ-analysis for system performances assessment[END_REF], this paper deals with LTI systems subject to parametric uncertainties, and provides some efficient probabilistic µ-based algorithms to calculate bounds on the probability of internal stability and H ∞ performance satisfaction. A particular emphasis is placed on the latter case. More precisely, an original algorithm is proposed to characterise the parametric domain on which performance violation is guaranteed. It is a more difficult problem than determining a guaranteed domain of performance satisfaction, but it offers some interesting complementary information. This is a major contribution of this paper, which is organised as follows. Section II first reminds some key definitions about µ-analysis, as well as the principle of branch-and-bound techniques, which serve as a starting point for probabilistic µ-analysis. The proposed algorithms are then introduced in Section III. Section IV finally presents some numerical results obtained on several benchmarks from the literature, demonstrating the efficiency of the proposed approach.

II. PROBLEM STATEMENT

Let us consider the standard interconnection for robust stability analysis as depicted in Fig. 1 (left), where M (s) is a stable real-valued LTI plant representing the nominal closed-loop system, and ∆ is a block-diagonal LTI operator gathering all the real parametric uncertainties:

∆ = diag(δ 1 I n1 , . . . , δ N I n N ) (1)
where (δ 1 , . . . , δ N ) ∈ R N . The set of matrices with such a block-diagonal structure is denoted ∆. Let then kB ∆ = {∆ ∈ ∆ : σ(∆) < k}, where σ(∆) denotes the largest singular value of ∆. Singularity of the interconnection at frequency ω ∈ R + is characterised by det(I -M (jω)∆) = 0.

- µ-analysis relies on the structured singular value, introduced in [START_REF] Doyle | Analysis of feedback systems with structured uncertainties[END_REF] and defined as follows.

M (s) ∆ - - - u y M (s) ∆
Definition 2.1: Let ω ∈ R + be a given frequency. If no matrix ∆ ∈ ∆ makes the interconnection of Fig. 1 (left) unstable, then the structured singular value µ ∆ (M (jω)) is defined to be equal to zero. Otherwise:

µ ∆ (M (jω)) = min ∆∈∆ {σ(∆) : det(I -M (jω)∆) = 0} -1 (2)
The robust stability margin is then defined as the inverse of the maximum value of µ ∆ (M (jω)) over the whole frequency range R + :

k r = sup ω∈R+ µ ∆ (M (jω)) -1 (3)
If the structure ∆ is now split into two distinct block structures ∆ = diag(∆ 1 , ∆ 2 ), and admissible uncertainties are taken in ∆ s = diag(B ∆1 , ∆ 2 ), a skewed version of the structured singular value can be defined [START_REF] Fan | A measure of worst-case H∞ performance and of largest acceptable uncertainty[END_REF].

Definition 2.2: Let ω ∈ R + be a given frequency. If no matrix ∆ = diag(∆ 1 , ∆ 2 ) ∈ ∆ s makes the interconnection unstable, then the skewed structured singular value ν ∆s (M (jω)) is defined to be equal to zero. Otherwise:

. ν ∆s (M (jω)) = min ∆∈∆s {σ(∆2) : det(I -M (jω)∆) = 0} -1 (4)
Even putting aside the NP-hardness [START_REF] Braatz | Computational complexity of µ calculation[END_REF] of computing k r or its skewed version, which usually requires to use upper and lower bounds k r and k r instead of the real value, the robust stability margin remains inherently conservative. Indeed, it only yields the largest open ball k r B ∆ on which stability is guaranteed, regardless of the actual shape of the exact stability domain. To obtain a better approximation of the latter, branch-and-bound algorithms are used [START_REF] Newlin | Mixed mu problems and branch and bound techniques[END_REF][START_REF] Roos | Efficient computation of a guaranteed stability domain for a high-order parameter dependent plant[END_REF]. Then, without introducing much additional numerical complexity, a probabilistic framework can be considered: if the parametric uncertainties δ i are now seen as random variables with given probability density functions f i , it is of interest to compute the probability measure of the stability domain. In other words, the following problem can be stated.

Problem 2.1: Calculate the probability P ∆,f (M (s)) that the interconnection M (s) -∆ is stable when ∆ ∈ B ∆ .

Let us now consider the interconnection of Fig. 1 (right) corresponding to performance analysis. The transfer from u to y is given by the upper linear fractional transformation:

. Fu(M (s), ∆) = M22(s)+M21(s)∆(I -M11(s)∆) -1 M12(s) [START_REF] Khatri | Guaranteed bounds for probabilistic µ[END_REF] For some technical reasons to be clarified in Section III-B, it is assumed hereafter that the transfer function from u to y is single input single output. Such an assumption is not so restrictive in practice. If needed, several performance channels can be processed sequentially. The following proposition makes the link between H ∞ performance and the skewed structured singular value.

Proposition 2.1: The following statements are equivalent: 1) max

∆∈B ∆ F u (M (s), ∆) ∞ ≤ γ
2) ∀ω ∈ R + , ν ∆a (M (jω)) ≤ γ where ∆ a = diag(∆, C) In the same vein as with robust stability, branch-and-bound techniques can be used to compute a domain of guaranteed performance satisfaction, whose probability measure can then be evaluated. The problem of interest is as follows:

Problem 2.2: Given a performance level γ > 0, compute the probability

P γ ∆,f (M (s)) that F u (M (s), ∆) ∞ ≤ γ when ∆ ∈ B ∆ .

III. PROBABILISTIC µ-ANALYSIS A. Stability Analysis

The exact computation of P ∆,f (M (s)) is unsurprisingly not achievable in the general case. But branch-and-bound algorithms provide a subdomain D s of the normalised parametric domain

Θ = {δ ∈ R N : |δ i | ≤ 1} = [-1, 1]
N on which stability is guaranteed, which gives a lower bound on the sought probability. Moreover, D s conveniently takes the form of a finite union of disjoint boxes

Ds = k D (k) s where D (k) s = [a (k) 1 , b (k) 1 ] × ... × [a (k) N , b (k) N ] (6)
whose probability measures are straightforward to compute:

p(D (k) s ) = N i=1 b (k) i a (k) i f i (δ i )dδ i (7)
There is much interest in also calculating an upper bound on P ∆,f (M (s)), in order not only to get a more accurate picture of the stability domain, but also potentially to quickly invalidate the current design with respect to a specified tolerance level. In other words, to further reduce the conservatism of the proposed probabilistic approach, the computation of a domain of guaranteed instability D s is desirable, so that P ∆,f (M (s)) can be bounded as follows:

p(D s ) ≤ P ∆,f (M (s)) ≤ 1 -p(D s) (8)
This proves to be straightforward if one realises that methods to calculate the robustness margin rely on detecting when the poles of the nominally stable system M (s) reach the imaginary axis as the size of ∆ increases. Thus, the same methods applied to a nominally unstable interconnection leads to a "robustness margin for instability" [START_REF] Khatri | Guaranteed bounds for probabilistic µ[END_REF]. Integrated into a branch-and-bound scheme, this yields the desired domain of guaranteed instability.

While the branch-and-bound and the probability calculations could be kept separate, it is computationally more interesting to combine the two aspects into a single algorithm. Indeed, this allows to integrate probability-related stopping criteria into the algorithm, e.g. interrupting the algorithm if all the boxes have a probability measure less than a given threshold p min . Such a strategy is described in Algorithm 1, whose outputs are a partition of

Θ = [-1, 1] N into Θ = D s ∪ D s ∪ D su (9)
and the probability measure of each component of this partition. The third term D su corresponds to the domain of uncategorised stability. Its presence stems from the fact that, realistically, the branch-and-bound algorithm can only approximate the domains of stability and instability. Ideally, D su should be as small as possible, although depending on the validation requirements, a coarse result can be sufficient to decide if the investigated system is satisfactory. For inconclusive iterations, the partitioning of the box is classically done along its longest edge, but other strategies can be implemented. In particular, µ-sensitivities can be used to reduce computational time [START_REF] Lesprier | Improved µ upper bound computation using the µ-sensitivities[END_REF], especially for systems with numerous parametric uncertainties.

Algorithm 1 Probabilistic µ-analysis for stability

1: Inputs: M (s), ∆, f i (δ i ) for i = 1, . . . ,
p(D su ) ← p(D su ) + k p( Θk ) 29: end if B. H ∞ Performance Analysis
Let us now consider Problem 2.2. As for robust stability, the idea is to combine branch-and-bound with µ-analysis to compute both a domain of guaranteed performance satisfaction D γ and a domain of guaranteed performance violation D γ , whose probability measures can then be calculated. Proposition 2.1 provides the means to calculate D γ , but a difficulty arises in the computation of D γ . Indeed, performance violation on Θ is equivalent to

min ∆∈B ∆ sup ω∈R+ |F u (M (jω), ∆)| > γ (10)
which, as a minimax problem, is more complicated than verifying performance satisfaction. As a preliminary step to tackle this problem, a useful property on LFT inversion is reminded (see e.g. [START_REF] Doyle | Review of LFTs, LMIs, and µ[END_REF]).

Lemma 3.1: Let M = M 11 M 12 M 21 M 22
, and assume that

M 22 is nonsingular. Then F u (M, ∆) -1 = F u (N, ∆),
where

N = M 11 -M 12 M -1 22 M 21 -M 12 M -1 22 M -1 22 M 21 M -1 22 .
Remark 3.1: For a scalar performance channel, the assumption that M 22 is nonsingular simply means that it is nonzero, i.e. there is effectively a transfer from u to y. Since it is the case of interest for performance analysis, this is not restrictive. Proposition 3.1 then derives a sufficient condition for Inequality [START_REF] Newlin | Mixed mu problems and branch and bound techniques[END_REF] to hold, in the form of a skew-µ calculation involving the inverse transfer from u to y.

Proposition 3.1:

Let I = {ω ∈ R + : |M 22 (jω)| > γ}. If there exists ω 0 ∈ I such that ν ∆a (N (jω 0 )) < 1 γ (11)
where

∆ a = diag(B ∆ , C) and N (jω) is such that (F u (M (jω), ∆)) -1 = F u (N (jω), ∆),
which guarantees that I is nonempty, and that N (jω) is well-defined for each ω ∈ I. Assume now that there exists ω 0 ∈ I such that (11) holds. From the main loop theorem ( [START_REF] Doyle | Review of LFTs, LMIs, and µ[END_REF]), this is equivalent to

   µ ∆ (N 11 (jω 0 )) < 1 max ∆∈B ∆ |F u (N (jω 0 ), ∆)| < 1 γ (13)
Since |F u (N (jω 0 ), ∆)| = |F u (M (jω 0 ), ∆)| -1 ∀∆ ∈ B ∆ , the second inequality of (13) can be rewritten as

min ∆∈B ∆ |F u (M (jω 0 ), ∆)| > γ (14)
which is a sufficient condition for Inequality [START_REF] Newlin | Mixed mu problems and branch and bound techniques[END_REF] to hold. Remark 3.2: The step from ( 13) to ( 14) strongly relies on the scalar nature of the performance channel and cannot be easily generalised to higher dimensions although this will be investigated in a near future.

A strategy to verify performance violation directly follows from Proposition 3.1, as outlined in Algorithm 2. end if 13: end if Note that the output false of Algorithm 2 merely means that the analysis was inconclusive, and does not exclude performance violation on Θ. More precisely, inconclusive cases can result from three causes:

1) the sufficient condition of Proposition 3.1 is not necessary, 2) the selected ω 0 in Algorithm 2 is not satisfactory, 3) the computation of the skew-µ upper bound is not sufficiently accurate. In particular, the heuristics used to select ω 0 can prove critical. It is readily apparent that exploring a frequency grid or frequency interval would be computationally inefficient, since only one frequency is sought, which is not even guaranteed to exist. On the other hand, the peak frequency of M 22 (s), while appearing as a natural choice for ω 0 , offers no guarantee of optimality.

The methods to validate performance satisfaction (Proposition 2.1) and performance violation (Proposition 3.1) on a given box can be integrated into a branch-and-bound scheme structurally similar to Algorithm 1. The essential difference lies in lines 8 through 24 of the algorithm, which must be replaced by a performance analysis:

• if M22 (s) ∞ < γ, then attempt to verify performance satisfaction on Θ using Proposition 2.1, • otherwise, attempt to verify performance violation on Θ using Algorithm 2. Θ is then either partitioned if the aforementioned test fails or added to the adequate domain (D γ , D γ , or the domain of uncategorised performance D γu ), whose probability measure is then updated. Note that all three issues raised on Algorithm 2 are alleviated by the branch-and-bound. Indeed, barring extreme cases of sensitivities to the uncertainties, it is realistic to assume that the transfer F u (M (s), ∆) varies less and less as ∆ is taken in smaller and smaller boxes. Thus, the sufficient condition of Proposition 3.1 becomes not only more relevant, but also less sensitive to the choice of ω 0 as a result. Moreover, the gap between the bounds of µ is known to become arbitrarily small as the domain Θ is partitioned into smaller boxes [START_REF] Roos | A set of µ-analysis based tools to evaluate the robustness properties of high-dimensional uncertain systems[END_REF], which resolves item 3. However, relying on the asymptotic behaviour of the branch-and-bound algorithm may be unsatisfactory from a computational point of view.

IV. APPLICATIONS

A. Illustrative Example

The following example, taken from [START_REF] Falcoz | Probabilistic µ-analysis for system performances assessment[END_REF], aims to illustrate the proposed algorithms on a two-dimensional problem that allows graphical representation. Consider the system given by the state-space representation

   ẋ = 0 1 -a 0 -a 1 x + 0 1 u y = 1 0 x ( 15
)
where a 0 and a 1 are two uncertain parameters described by

a 0 = 1 + q 0 , |q 0 | ≤ 2 a 1 = 0.8 + q 1 , |q 1 | ≤ 1 (16)
The parametric domain Θ is therefore the set of all possible values of the ordered pair (q 0 , q 1 ), i.e. Θ = [-2, 2] × [-1, 1]. For such an example, the Routh criterion can be used to determine that the domain of stability is exactly ] -0.5, 1]×] -0.8, 1]. Assuming a uniform distribution of the uncertainties on Θ, this represents 67.5% of the domain. The proposed algorithm for stability analysis yields P ∆,f (M (s)) ∈ [0.672, 0.679], which is a good estimation.

The same problem is now solved assuming q 0 ∼ N [-2,2] (0, (2/3) 2 ) and q 1 ∼ N [-1,1] (0, (1/3) 2 ), where N [a,b] (m, σ 2 ) denotes the truncated normal distribution. This type of distribution refers to a case where the nominal configuration is deemed more likely, and with the probability of occurrence of a configuration becoming smaller as it gets farther from the nominal point. Algorithm 1 now leads to P ∆,f (M (s)) ∈ [0.925, 0.931], which is a more optimistic result. This highlights the interest and potential of the probabilistic framework in reducing the conservatism of µanalysis. Regarding H ∞ performance analysis and setting the performance level γ = √ 2, a uniform distribution of the uncertainties leads to P γ ∆,f (M (s)) ∈ [0.344, 0.359], while the truncated normal distribution changes this estimation to P γ ∆,f (M (s)) ∈ [0.411, 0.441]. The improvement is less drastic, but this can be explained by the fact that the nominal configuration is close to the limit of performance. This case is less likely to occur for a controlled system whose controller is designed using robust synthesis methods.

Table I compares the results obtained with Algorithm 1 from [START_REF] Falcoz | Probabilistic µ-analysis for system performances assessment[END_REF] (whose Matlab code has been provided by the authors) and with the algorithms from the present paper. A dramatic reduction in the computational time can be observed. This is due to the computation of guaranteed instability or non-performance domains, which is done here but not in [START_REF] Falcoz | Probabilistic µ-analysis for system performances assessment[END_REF]. The difference is further illustrated by Fig. 2 and3. The comparison was done on a CentOS Linux 7.4 Workstation with a CPU Intel Xeon W3530 running at 2.80GHz and 6 GB of RAM, with Matlab R2017b installed. The adaptation of Algorithm 1 to H ∞ performance analysis is only applicable to boxes on which stability is already guaranteed. In the proposed implementation, a preliminary stability analysis is therefore performed, and the identified stable boxes serve as inputs for the performance analysis. Both analyses are accounted for in Table I.

Remark 4.2: The limiting threshold on the probability measure of the analysed boxes is used slightly differently between the algorithm from [START_REF] Falcoz | Probabilistic µ-analysis for system performances assessment[END_REF] and Algorithm 1. The former bisects all uncategorised boxes, then checks the probability measure of the resulting boxes to determine if they are to be analysed or not. The latter uses the probability measure of uncategorised boxes to determine if partitioning should be done. Thus, for the same limiting threshold p min , Algorithm 1 analyses more boxes. Fig. 2.

Stability (left) and performance (right) analyses without instability/non-performance verification [START_REF] Falcoz | Probabilistic µ-analysis for system performances assessment[END_REF]: the green domain satisfies the tested criterion, the red domain is uncategorised Fig. 3. Stability (left) and performance (right) analyses with instability/nonperformance verification: the green domain satisfies the tested criterion, the blue domain invalidates it, and the red domain is uncategorised

B. Numerical Results

Algorithm 1 is now applied to several benchmarks described in the literature. For these tests, all normalised uncertainties follow the truncated normal distribution N [-1,1] (0, (1/3) 2 ), and two partitioning strategies are used: bisection along the longest edge, and bisection along the highest µ-sensitivity. The limiting threshold for bisecting boxes is set to p min = 10 -4 . The results are summarised in Table II. The notation m × p in the description of ∆ means that ∆ contains m blocks of size p × p, all blocks being real diagonal. R denotes the ratio of the probability measure of the categorised boxes (stable or unstable) over the total probability measure of the operating domain Θ, i.e. R = p(Ds)+p(Ds) p(Ds)+p(Ds)+p(Ds u ) = p(D s ) + p(D s). The obtained ratios are close to 1 for all benchmarks, meaning that the domains of stability and instability are well approximated by D s and D s respectively, leaving only a small domain D su as uncategorised (in the sense of probability measure).

Regarding the partitioning strategies, while the use of µsensitivities tends to improve both the domains computation and the running time, this is not always the case, benchmark 10 being a counterexample. However, this system aside, the use of µ-sensitivities leads to significant time reduction as the dimension of the problem increases, as observed on benchmarks 19, 21 and 26, whose parametric domains have respective dimensions 9, 9 and 23. This shows that, in spite of the effectiveness of naive branching schemes [START_REF] Zhu | Soft vs. hard bounds in probabilistic robustness analysis[END_REF], there is still an interest in investigating more sophisticated ones.

Table III shows test results for probabilistic performance analysis, where R now denotes the ratio R = p(Dγ )+p(Dγ ) p(Dγ )+p(Dγ )+p(Dγ u ) (which differs from p(D γ ) + p(D γ ) if stability is not achieved on the complete domain Θ). As mentioned in Section III-B, the peak frequency of M 22 (s) is selected as the frequency ω 0 . The relevance of µ-sensitivities for partitioning is less consistent, leading to a slight decrease of the probability measure of D γu at the cost of increased execution time for all but benchmark B5. The latter benefits the most from the use of µ-sensitivities, with a drastic reduction of p(D γu ) (more than 25%) and an execution time reduced by a factor of approximately 3.

To assess the relevance of the chosen ω 0 , sample configurations have been drawn in each uncategorised box to see whether a configuration satisfying the performance level could be found. This would mean that a suitable ω 0 does not exist, i.e. there is nothing to gain from improving the choice of ω 0 . For each nominally non-performing N -dimensional box, the samples drawn are the corners of the box, for a total of 2 N samples per box. The reason behind this sampling method is that the farthest points from the centre of the box are empirically the most likely to exhibit the most substantial changes on their transfer function. Also, the relatively small number of dimensions in the considered benchmarks keeps the sample size from being prohibitive. Results show that the loss on p(D γ ) remains reasonable on the tested benchmarks, not exceeding 3%. This means that no change in the heuristics can lead to dramatic improvements of the algorithm. However, if improvements are sorely needed, then a possible solution could be to consider as additional candidates for ω 0 the peak frequencies of the 2N configurations located at the centre of the (N -1)dimensional faces of the currently studied N -dimensional box. This choice constitutes a compromise between taking configurations that are far from the centre of the box, and limiting the number of candidate frequencies, since taking the corners as was done for the sampling would lead to 2 N additional frequencies.

V. CONCLUSION In this paper, new contributions are proposed to improve the efficiency of robustness analysis by combining deterministic µ-based methods with a probabilistic framework. Both robust stability and H ∞ performance are addressed. The proposed algorithms are implemented in Matlab and tested on several benchmarks taken from the literature. The resulting routines are available in the SMART Library of the SMAC Toolbox at http://w3.onera.fr/smac. Beside the efficiency of the algorithms, test results showcase that probabilistic µ-analysis can be an adequate response to reduce the conservatism of a purely deterministic approach. In the overall validation process of a controlled system, probabilistic µ-analysis complements Monte Carlobased methods nicely. It provides relevant inputs for the sampling strategy of the simulation phase, and could ultimately lighten the tedium of control laws development. In light of the above, probabilistic µ-analysis appears as an attractive method for the validation of controlled systems with parametric uncertainties, provided that reliable probability density functions are available. Future work will focus on extending the presented H ∞ performance analysis method to the MIMO case.

Fig. 1 .

 1 Fig. 1. Standard interconnections for stability analysis (left) and performance analysis (right)

Algorithm 2 1 : 8 :

 218 Verification of performance violation on Θ Inputs: M (s), ∆, γ 2: Output: true if performance violation is guaranteed on Θ, false otherwise 3: if M 22 (s) ∞ ≤ γ then ν ∆a (N (jω 0 )), where N (jω 0 ) is given by Lemma 3.1 if ν ∆a (N (jω 0 )) < 1/γ then 9:

  N , p min , p max 2: Outputs: D s , D s, D su , p(D s ), p(D s), p(D su ) 3: L ← {Θ} 4: D s , D s, D su ← ∅ 5: p(D s ), p(D s), p(D su ) ← 0 6: while L = ∅ and p(D s ) + p(D s) < p max do

	7:	select and remove the box Θ from L with the highest
		probability measure and compute the interconnection
		M (s) -∆ normalised w.r.t.	Θ
	8:	compute a guaranteed robust stability margin k r for
		M (s) -	∆
	9:	if k r ≥ 1 then	analysis is conclusive
	10: 11: 12:	if F u ( M (s), ∆) is nominally stable then add Θ to D s p(D s ) ← p(D s ) + p( Θ)
	13:	else
	14:		add Θ to D s
	15:
	16:	end if
	17:	else		analysis is inconclusive
	23:	end if
	24:	end if
	25: end while
	26: if L = ∅ then
	27:	add all remaining boxes Θk of L in D su
	28:

p(D s) ← p(D s) + p( Θ) 18: if p( Θ) > p min then 19: partition Θ and add the boxes obtained in L 20:

else 21: add Θ to D su 22: p(D su ) ← p(D su ) + p( Θ)

  then Inequality[START_REF] Newlin | Mixed mu problems and branch and bound techniques[END_REF] holds, and performance violation is guaranteed on Θ, i.e. ∀∆ ∈ B ∆ .

	Proof: First, note that nominal performance violation is
	equivalent to
	sup

ω∈R+ |M 22 (jω)| > γ


Mohamed Elobaid
email: mohamed.elobaid@uniroma1.it

Mattia Mattioni
email: mattia.mattioni@uniroma1.it

Salvatore Monaco
email: salvatore.monaco@uniroma1.it

Dorothée Normand-Cyrot

On unconstrained MPC through multirate sampling

Keywords: Control of Sampled-Data Systems, Model-Predictive, Optimal Control

This paper formally highlights how the multirate sampled data equivalent model can be exploited for prediction in an MPC formulation in order to mitigate the possible instability arising from an MPC design while ensuring prefixed boundedness of the control amplitude. This last aspect is in particular addressed and solved with reference to the class of systems which admit, under feedback, a computable sampled model.

INTRODUCTION

Model Predictive Control (MPC) has become a widely investigated research area in linear and nonlinear control and their applications (e.g., [START_REF] Boucher | La commande prédictive[END_REF]; [START_REF] Camacho | Model predictive control[END_REF]; [START_REF] Borrelli | Predictive control for linear and hybrid systems[END_REF]; [START_REF] Kwon | Stabilizing state-feedback design via the moving horizon method[END_REF]). Roughly speaking the control action is designed by solving a constrained optimization problem subject to the system dynamics and possibly additional requirements and bounds. The simplest way to design and implement the resulting control law makes use of a sampled data model of the plant for prediction that is exploited for solving the optimization problem over a finite prediction horizon of length n p . In this context, one implicitly assumes the dependence of the predicted future values of the output over n c future controls actions, generally referred to as the control horizon. Driving the output trajectory to a desired reference at the sampling instants via discrete time model predictive control while preserving stability in closed loop may induce difficulties especially when the plant (and the model used for prediction) is not minimum phase. In that case, the choice of the prediction and control horizons plays a key role. In fact, typically, one sets n c much smaller than n p to address this fact as well as minimizing the computations required [START_REF] Clarke | Generalized predictive control-part i. the basic algorithm[END_REF]). As proven in [START_REF] Monaco | Zero dynamics of sampled nonlinear systems[END_REF], singlerate (or standard) sampling generally induces unstable and extra zero-dynamics so that minimum-phaseness is lost independently of the continuous-time plant properties; to overcome such a pathology, a multirate sampling procedure has been properly introduced to preserve the continuous-time internal properties. With this in mind, it is proposed and shown in the sequel that the use of a mul-tirate (MR) sampled equivalent model at the prediction and implementation level overcomes the aforementioned problems in an MPC control scheme. Recalling moreover the peculiar properties of MR sampling when applied to nonholonomic systems [START_REF] Monaco | An introduction to motion planning under multirate digital control[END_REF]), it is shown that a MR based MPC control scheme can be fruitfully employed to limit the amplitude and increase the robustness of MR solution to the steering problem. This result relies on the fact that MPC restitutes the MR steering solution. The work is organised as follows. In Section 2 recalls on single and multi-rate sampling are given and the problem under investigation is stated. Section 3 is devoted to the proposed MR-MPC control scheme and prove its effectiveness. Section 4 investigates the relation between MPC and MR controllers in the sampled data context with reference to the steering control of nonholonomic dynamics; the example of the unicycle dynamics is used to verify the effectiveness of the proposed control scheme. Concluding remarks end the paper.

RECALLS AND STATEMENT OF THE PROBLEM

Notation and definitions

All functions and vector fields defining the dynamics are assumed smooth and complete over the respective definition spaces. M U (resp. M I U ) denotes the space of measurable and locally bounded functions u : R → U (u :

I → U , I ⊂ R) with U ⊆ R. U δ ⊆ M U denotes the set of piecewise constant functions over time intervals of fixed length δ ∈]0, T * [; i.e. U δ = {u ∈ M U s.t. u(t) = u k , ∀t ∈ [kδ, (k + 1)δ[; k ≥ 0}. When u(t) ∈ R m then u j k , (u k ) j , u (j)
k are the jth component of u for t ∈ [kδ, (k + 1)δ[; k ≥ 0}, u k raised to the power j and the jth derivative respectively. Given a vector field f , L f denotes the Lie

derivative operator, L f = n i=1 f i (•)∇ xi with ∇ xi := ∂ ∂xi while ∇ = (∇ x1 , . . . , ∇ xn ).
The Lie exponential operator is denoted as e L f and defined as e L f := I + i≥1

L i f i! . A function R(x, δ) = O(δ p
) is said to be of order δ p (p ≥ 1) if whenever it is defined can be written as R(x, δ) = δ p-1 R(x, δ) and there exist function θ ∈ K ∞ and δ * > 0 such that ∀δ ≤ δ * , | R(x, δ)| ≤ θ(δ * ).

Sampled data systems and multirate sampling

The following recalls on sampled-data systems are given (see [START_REF] Monaco | Issues on nonlinear digital control[END_REF] and the references therein). Given a SISO system ẋ = f (x) + g(x)u, y = h(x)

(1) and considering u(t) ∈ U δ and y(t) = y(kδ) for t ∈ [kδ, (k+ 1)δ[ (δ the sampling period), the dynamics of (1) at the sampling instants is described by the single-rate sampleddata equivalent model

x k+1 = F δ (x k , u k ), y k =h(x k )
(2) with x k := x(kδ), y k := y(kδ), u k := u(kδ). The mapping F δ (•, u) : R n × R → R n gets the form of formal series expansion in powers of δ that is (dropping the time subscript for clarity)

F δ (x, u) = e δ(L f +uLg) x = x + i>0 δ i i! (L f + uL g ) i x. (3)
It is a matter of computations to verify that if (1) has well defined relative degree r ≤ n, the relative degree of the sampled-data equivalent (2) always falls to r d = 1; namely, one has

y k+1 = r i=0 δ i i! L i f h(x) x k + δ r r! u k L g L r f h(x) x k + O(δ r+1 ) so that ∇ u k y k+1 = δ r r! L g L r f h(x) x k + O(δ r+1 ) = 0.
As a consequence, whenever r > 1, the sampling process induces a further zero-dynamics of dimension r -1 (the socalled sampling zero-dynamics) that is in general unstable for small values of δ when r > 1. As a consequence, dynamics-inverting controllers via single-rate sampling do not guarantee internal stability. Multirate sampling has been developed in a nonlinear context to overcome those issues. Namely, by setting

u(t) = u i k for t ∈ [(k + i -1) δ, (k + i) δ[ for i = 1, . . . , r and y(t) = y k for t ∈ [kδ, (k + 1)δ[,
the multirate equivalent model of order r of (1) gets the form

x k+1 =F δ m (x k , u k ), Y k = H(x k ) (4) with δ = δ r , u ∈ R r = u 1 . . . u r , a dummy output vector H(x) = h(x) L f h(x) . . . L r-1 f h(x) and F δ m (x k , u k ) =e δ(L f +u 1 k Lg) . . . e δ(L f +u r k Lg) x x k = F δ m (•, u r k ) • • • • • F δ (x k , u 1 k ).
One gets so far a MIMO system possessing vector relative degree r δ = (1, . . . , 1) and zero dynamics which inherits the zero-dynamics stability properties of ( 1). It must be recalled that exact computation of the sampled data equivalent model cannot in general be achieved, so that approximated models are usually computed by truncation in different ways of the expansion [START_REF] Mpc-Mr | = 0.2I, Q = 10I Fig. 4. MR-MPC Vs MR tracking with R = 0.5I, Q = 2I[END_REF]. Related to the properties of ( 3), the notion of exact and finite computability of the sampled data model, possibly under preliminary feedback, has also been introduced in Monaco and Normand-Cyrot (1992).

Multirate digital steering of chained forms

Multirate sampling has been shown to be of interest in the digital design of nonholonomic and under actuated mechanical systems (see [START_REF] Monaco | An introduction to motion planning under multirate digital control[END_REF]). Under preliminary continuous-time feedback, a nonholonomic mechanical system admits the so-called single chained dynamics

ξ1 = u 1 , ξ2 = u 2 , ξi = ξ i-1 u 1 i = 3, . . . , n (5)
admitting a finitely computable sampled dynamics; moreover, the one-step ahead dynamics can be easily inverted with respect to the control input to deduce a multirate feedback ensuring exact deadbeat steering. For, consider ( 5) and (ξ 0 , ξ f ) then there exists a multirate control of order (n -1) on u 2 and 1 on u 1 (i.e. δ = δ n-1 ) such that system (5) is exactly steered in one step δ from ξ 0 to ξ f . The feedback ensuring exact steering for (5) can be easily deduced by inverting the multirate sampled data equivalent model of (5) provided by

ξ 1,k+1 = ξ 1,k + δu 1,k , ξ 2,k+1 = ξ 2,k + δ i u i 2,k ξ 3,k+1 = ξ 3,k + δu 1,k ξ 2,k + η 1 ( δ2 , u 1,k )u 2,k (6)
. . .

ξ n,k+1 = ξ n,k + G(δ, ξ k , u 1,k ) + η n-2 ( δn-1 , u 1,k )u 2,k with η 1 (•), G(•) and η n (•) being η 1 (u 1,k ) = 1 2! c 1 1 u 1,k c 2 1 u 1,k . . . c n-1 1 u 1,k G(δ, ξ i,k , u 1,k ) = ξ n,k + δu 1,k ξ n-1,k + δ 2 2! (u 1,k ) 2 ξ n-2,k + . . . + δ n-2 (n -2)! (u 1,k ) n-2 ξ 2,k η n-2 (u 1,k ) = 1 (n -1)! c 1 n-2 u 1,k c 2 n-2 u 1,k . . . c n-1 n-2 u 1,k
with some suitable constants c j i . For all fixed ξ 0 = ξ[kδ], ξ f = ξ[kδ + δ], directly solving the above system in the unknowns (u 1,k , u 2,k ), one gets

u 1,k = 1 δ (ξ 1,k+1 ) -ξ 1,k ) u 2,k = η -1 (ξ k+1 -G(•)) (7)
with η(•), G(•) and ξ(•) being the compact forms of the corresponding elements in (6). Taking into account the preliminary feedback, steering is achieved under piecewise continuous control designed on the basis of the multirate sampled model. Roughly speaking, the multirate sampling is used as a trajectory planner.

A key thing to note on the control solution above is that, while it does steer the system to the desired final state, it is indeed an inverting controller and the control effort might grow unboundedly so making the feedback not implementable in practice. To overcome this issue, we shall improve such a feedback via MPC.

Problem statement

Consider the continuous-time system (1) under sampling, with relative degree r ≤ n being minimum phase . Hereinafter we shall address the problem of driving the output trajectory, to a desired reference ν(t) at the sampling instants t = kδ, k ≥ 0 via discrete time MPC [START_REF] Camacho | Model predictive control[END_REF]) while preserving stability in closed loop; that is

y k = ν k , k ≥ k * with ν k = ν(kδ) by minimizing the cost functional J = np i=1 e k+i 2 Q + u k+i-1 2 R = np i=1 L(x k+i , ν k+i , u k+i-1 ) (8)
with Q > 0, R ≥ 0 being appropriate penalizing weights on the tracking error and input magnitude and n p being the prediction horizon; moreover, e is a suitably defined error map, such that e k = 0 iff y k = ν k . MPC induces a constrained optimization problem subject to the dynamics (4) and possibly additional requirements and bounds. To solve this problem several methods are available, the simplest to implement of which is the so called direct single shooting [START_REF] Hicks | Approximation methods for optimal control synthesis[END_REF]) by plugging ( 4) into (8) so getting

J = np i=1 L(•, u k+i-1 ) • (F δ m (•, u k+i-1 ) • . . . • (F δ m (x k , u k )).
Hence, an optimal solution u e = u k . . . u k+nc-1 is computed by solving ∇ ue J = 0 with n c being the so-called control horizon. For our purposes it is interesting to note that in its usual implementation MPC makes use of a single-rate sampled data model of the plant of the form (2) for prediction. This induces the loss of the minimum-phaseness, so forcing the designer to set n c < n p to recover internal stability (or using terminal penalties and/or constraints sets) while also defining a dynamical controller (in the sense of using feedback on the states and also on the previous controls) to ensure off-set free tracking. In the approach we are proposing, the use of multirate sampled data model will provide a static feedback overcoming both issues.

PREDICTIVE MULTIRATE DIGITAL CONTROL OF NONLINEAR SYSTEMS

With reference to the problem statement in the previous section, and the augmented output vector, we set out to state our main result, however to do so one needs the following assumption; Assumption 1. Measures of ν and its derivatives ν (i) for i = 1 . . . r -1 are available at all t = kδ, k ≥ 0.

For the sake of compactness, we shall define the extended output vector dynamics for n p = n c future values as

Y e k+1 = A e Y k + B e (x k )u e k + Θ(x k , u e k )
with

Y e k+1 = Y k+1 Y k+2 . . . Y k+np and u e k = u k u k+1 . . . u k+nc-1 B e (•) = L g L r-1 f h(•)     B 0 . . . 0 AB B . . . 0 . . . A np-1 B A np-2 B . . . B     , B = D∆(r) D = diag( δr /r!, . . . , δ) ∆(r + j) =    (r + j) r+j -(r + j -1) r+j . . . (j + 1) r+j -(j) r+j . . . (r + j) j -(r + j -1) j . . . (j + 1) j -(j) j    A =    1 δ . . . δ r-1 /(r -1)! . . . 0 0 . . . 1    A e =     A A 2 . . . A np    
with j ≥ 0 and Θ(•) containing all higher order terms in O( δr+1 ) (?). The next result shows that the problem in Section 2.4 is always solvable with n p = n c > 1 under multirate feedback, provided the relative degree is well defined.

Theorem 2. Let ( 1) possess relative degree r ≤ n, and (4) be its multirate equivalent model of order r. Consider the MPC problem with cost functional ( 8) and

e k = Y k -ν k with ν k = ν k . . . ν (r-1) k
. Then, there exists δ * > 0 such that for all δ ∈ [0, δ * [, the MPC problem is solvable with internal stability for all n p = n c ≥ 1. The feedback is defined as the unique solution to the equality

(K(x, u e )Q e B e (x) + R e )u e = B e (•)(ν e -A e Y -Θ(x, u e )) (9)
with

Q e = I ⊗ Q, R e = I ⊗ R K(x, u e ) = (B e (x) + ∇ ue Θ(x, u e ))
ν e k = ν k+1 . . . ν k+np .

Sketch of proof. To prove that ( 9) is optimal, we first rewrite (8) as J = Y e k+1 -ν e k Qe + u e k Re whose jacobian is clearly annihilated by the solution to (9). Existence of a feedback solution can be deduced by rewriting (9) a formal series expansion in powers of δ and applying the implicit function theorem. Indeed, the term (D -1 K(x, u e )Q e B e (x) + R e ) is invertible as δ → 0 [START_REF] Monaco | Issues on nonlinear digital control[END_REF]; [START_REF] Mattioni | Immersion and invariance stabilization of strictfeedback dynamics under sampling[END_REF]). Internal stability is ensured by the minimum phaseness of the continuous-time plant which is consequently preserved under multirate sampling [START_REF] Monaco | Zero dynamics of sampled nonlinear systems[END_REF]). Remark 3. The solution obtained in Theorem 2 is implicitly defined by the above equality and is a formal series in powers of δ. Such a solution cannot be generally exactly computed in practice although several procedures are available for deducing approximation up to any desired order so to guarantee the required performances (see [START_REF] Monaco | Issues on nonlinear digital control[END_REF] for further details). Remark 4. As R → 0, the feedback defined by (9) coincides with the deadbeat inverting control that steers the output to the desired ν in one step of length δ. Such a feedback comes with an effort that is in general, inversely proportional to δ, thus by suitably setting R one can reduce the effort while still guaranteeing off-set free tracking in finite time. Roughly speaking, with reference to a minimum-phase plant, Theorem 2 suggests the use of MR-MPC control law of order equal to the relative degree. In what follows, we explicitly define this solution for nonholonomic systems that are feedback-equivalent to chained forms [START_REF] Brockett | Differential geometric control theory: proceedings of the conference held at Michigan Technological University[END_REF]). In doing so, we formally show that as R → 0 one recovers the standard deadbeat control. As a byproduct, we also provide an extension of Theorem 2 to the case of MIMO systems for which the relative degree might not be defined.

Y e k+1 =      C d A d C d A 2 d . . . C d A np d      x k +      C d B d 0 . . . 0 C d A d B d C d B d . . . 0 . . . C d A np-1 d B d C d A np-2 d B d . . . C d A np-nc d B d      u e k

PREDICTIVE MULTIRATE STEERING FOR CHAINED FORMS

For illustrative purposes, the following discussion will consider the chained form (5) with n = 3 albeit the arguments extend to the general case as highlighted. For, suppose one wants the state of the system to converge to a desired trajectory ν ∈ R n . Considering u 2 (t) = u i 2,k = u 2 (kδ + (i -1) δ) for i = 1, 2, one can write the output prediction over n p = 1 as

x 1,k+1 x k+1 = 1 0 0 A 2 m (•) x 1,k x k + 1 0 0 R m (•) u 1,k u 2,k (10) with u 2 = (u 1 2 u 2 2 ) and A m (•) = 1 0 δu 1,k 1 , b =   δ δ2 2! u 1,k   , R m (•) = [A m (•)b b]
which compactly rewrites as Y = F (δ, u 1 )X+G(δ, u 1 , u 2 )u 2 . One can then proceed in a similar fashion to the previous section, using the cost index ( 8) with e = col(e 1 , e 2 , e 3 ) = x 1 x -ν and setting ∇ u J = 0, so getting (when Q = I) that the optimal control u = u is solution to (dropping time subscript for clarity)

u 1 = -2δe 1 -e 3 (3 δ2 u 1 2 + δ2 u 2 2 + 4 δx 2 ) 2δ 2 + (3 δ2 u 1 2 + δ2 u 2 2 + 4 δx 2 ) u 1 2 = - 2 δ(e 2 + δu 2 2 ) -3 δ2 u 1 (e 3 + 2 δu 1 x 2 -0.5 δ2 u 1 u 2 2 ) 2 δ2 + 4.5 δ4 u 2 1 u 2 2 = - 2 δ(e 2 + δu 1 2 ) -δ2 u 1 (e 3 -2 δu 1 x 2 -1.5 δ2 u 1 u 1 2 ) 2 δ2 -0.5 δ4 u 2 1 (11) with u i 2 (i = 1
, 2) being the two controls resulting from multirate of order 2 over u 2 and ν i , i = 1, 2, 3 being the reference values over the single step prediction horizon. To show that the solution to this system of equations coincides with that of the multirate inverting controller, it is sufficient to show that the multirate inverting solution, is indeed a solution of this system of equations. Proposition 6. The multirate inverting controller ( 7) is a solution of the optimal control problem with the cost function (8) and n p = n c = 1, Q = I, R = 0.

Proof: Starting from ( 7), one has that

u 1 = ν 1 -x 1 δ , u 1 2 = - ν 2 -x 2 2 δ - x 3 -ν 3 + 2 δ u 1 x 2 δ2 u 1 u 2 2 = 3 (ν 2 -x 2 ) 2 δ + x 3 -ν 3 + 2 δ u 1 x 2 δ2 u 1
solves (11). As a matter of fact, (11) admit two solutions for (u 1 , u1 2 , u 2 2 ), one of which corresponds to the solution u 1 = 0 which is discarded 1 whereas the other one is u 1 = e1 δ and

u 1 2 = -2 δ e 3 -δ ν 1 ν 2 -3 δ ν 1 x 2 + δ ν 2 x 1 + 3 δ x 1 x 2 2 δ2 (ν 1 -x 1 )
u 2 2 = --2 δ e 3 -3 δ ν 1 ν 2 -δ ν 1 x 2 + 3 δ ν 2 x 1 + δ x 1 x 2 2 δ2 (ν 1 -x 1 ) clearly coinciding with the expression above from (7). A( δ, u 1,k+np-i-1 ) (13) and A(•) = A 2 m (•). We can then substitute this expression of prediction in our cost function and take the partial derivatives with respect to each u i and, then prove that whenever n p = n c the multirate inversion control is an optimum control with respect to our cost function. The following statement highlights this fact.

  Remark 5. It is rather straightforward to show that when (1) is linear (i.e. f (x) = F x, g(x) = G and h(x) = Cx) one recovers the known output trajectory prediction of the discrete time model with A d = e F δ , B d = [ Āδ (r-1) B . . . B], Ā = e F δ , B = δ 0 e F s dsG and C d = C; i.e,

  which in compact form can be written as Y e k+1 = A e x k + B e u e k . Along the lines of[START_REF] Borrelli | Predictive control for linear and hybrid systems[END_REF], the optimal control is u * e = (B e Q e B e + R e ) -1 B e Q e (ν e -A e x). When implementing the above optimal control trajectory a receding horizon algorithm (i.e. selecting the first m components of u e [k] and discarding the rest and repeating at each sampling instant) one hasu k = (I m 0 . . . 0) (B e Q e B e + R e ) -1 B e Q e (ν e -A e x).Thus, when n p = n c , Q e = I and R e = 0 the MR-MPC feedback reduces tou k = (C d B d ) -1 (ν -C d A d x), which is the classical dynamics inverting feedback.

4. 1

 1 The case of n p = n c > 1We write the prediction model for the two components of the states vector as followsx 1,k+np = x 1,k + δ np-1 i=0 u 1,k+i x k+np = φ(k + n p , k)x k + + n p -1, k + i + 1)R(•, u 1,k+i )u 2,k+i where φ(k + n p , k) = np-1 i=0

Since, this solution doesn't bring the error on the state x 1 to zero.

Partially funded by Université Franco-Italienne/Università Italo-Francese (UFI/UIF) through the Vinci Program.

Proposition 7. As R → 0, the control minimizing (8) computed over the prediction model ( 12) reduces to the multirate plant inversion solution (7) if n p = n c and Q = I.

Proof: The proof follows from induction starting with Proposition 6 which proves n p = n c = 1. By assuming that the statement holds for some n p = N , we show it holds for n p = N + 1. Let us split the cost functional (8) along the prediction model ( 12) as follows

Denoting by u * mr the multirate inverse controller which satisfies by assumption ∇ u J 1 (u * mr ) = 0, it remains to prove that that ∇ u J 2 (u * mr ) = 0. For, notice that

and recalling that u * mr is of the form (7), the proof proceeds as follows; one gets that as ∇ u2 J 2 = 0 (for compactness we omit the time variable k and we write i for k + i)

Inspecting the second equation above, namely ∇ u2 J 2 = 0 gives two possibilities, either the term -2nφ(N, i + 1)R(•, u 1,i ) = 0 or the term between the large brackets is 0, which upon inspection is 0 exactly when we set for i

with Ξ( δ, u 1 ) collecting the product terms of φ(N, i + 1)R(•, u 1,i ) which coincides with the solution obtained from (7). We then write only ∇ u 1,N J 2 , since by assumption and from the expression above for u

By substituting u 1 2,N , u 2 2,N as in ( 14) one recovers u 1,N = ν 1,N +1 -x 1,N δ which possesses the same form as in Proposition 6.

It is rather intuitive to see that the discussion above holds for general chained forms, and the statements can be extended along the same lines, albeit the notations and algebraic manipulations will get rather long and cumbersome, the following statement summarizes this. Theorem 8. Consider the dynamics of the form (5) admitting multirate equivalent model ( 6). Then, the multirate plant-inverting feedback (7) solves the MPC problem with (8) under prediction model ( 6) with perfect steering, whenever

SIMULATION RESULTS AND COMMENTS

Simulations are performed to compare the proposed control scheme (MPC-MR) with respect to the standard MPC implementation (Figs 1,2) and the classical MR control (Figs 3,4). In all the cases MPC-MR is implemented with n p = n c , δ = 1. Figure 1 clearly emphasises the pathology motivating this work; the MPC may fail when n p = n c and no stability constraints are incorporated, even with no penalty on the control. To prevent this, as suggested in the literature, in Fig 2 MPC works with n p > n c ; the comparison with the proposed MPC-MR with R > 0 in this case shows the better performance of our solution. A deeper comparison is proposed in Figs 3 and4 where the proposed MPC-MR and the standard MR solutions are shown for steering and tracking maneuvers under the penalty R > 0. The proposed MPC-MR scheme appears to be the natural context to be adopted to account for the control amplitude in standard MR.

CONCLUSION

We establish an intuitive interpretation of MR inverting controllers, by highlighting the roles of the prediction and control horizons, and their relations to the relative degree. We then motivate the use of MPC with a MR prediction model through penalizing the controls, and obtaining comparable performance to the plant inversion controller, while maintaining low control effort. Future works concern the application of this improved MPC scheme to several case studies as in power systems or automotive control [START_REF] Giuseppi | Risk adverse virtual power plant control in unsecure power systems[END_REF]; [START_REF] Gionfra | Combined feedback linearization and mpc for wind turbine power tracking[END_REF]. Ongoing work is addressing the extension to other classes of systems, possibly non-minimum phase [START_REF] Mattioni | On partially minimum-phase systems and disturbance decoupling with stability[END_REF].