Chance-constrained sets approximation: A probabilistic scaling approach

Automatica 137 (2022) 110108 Contents lists available at ScienceDirect Automatica journal homepage: www.elsevier.com/locate/automatica Chance-constrained sets approximation: A probabilistic scaling approach✩ Martina Mammarella a,Victor Mirasierra b,Matthias Lorenzen c,Teodoro Alamo b, Fabrizio Dabbene a,∗ aCNR-IEIIT, c/o Politecnico di Torino, Corso Duca degli Abruzzi 24, Torino, Italy bUniversidad de Sevilla, Escuela Superior de Ingenieros, Camino de los Descubrimientos s/n, Sevilla, Spain cSystemwissenschaften, TTI GmbH, Nobelstr. 15, 70569 Stuttgart, Germany article info Article history: Received 20 January 2021 Received in revised form 17 October 2021 Accepted 19 October 2021 Available online 1 January 2022 abstract In this paper, a sample-based procedure for obtaining simple and computable approximations of chance-constrained sets is proposed. The procedure allows to control the complexity of the approximating set, by defining families of simple-approximating sets of given complexity. A probabilistic scaling procedure then scales these sets to obtain the desired probabilistic guarantees. The proposed approach is shown to be applicable in several problems in systems and control, such as the design of Stochastic Model Predictive Control schemes or the solution of probabilistic set membership estimation problems. ©2021 Published by Elsevier Ltd. 1. Introduction In real-world applications, the complexity of the phenomena encountered and the random nature of data makes dealing with uncertainty essential. In many cases, uncertainty arises in the modeling phase, in some others it is intrinsic to both the system and the operative environment, as for instance wind speed and turbulence in aircraft or wind turbine control (Prékopa, 1995). Deriving results in the presence of uncertainty is of major relevance in different areas, including, but not limited to, optimization (Sahinidis,2004) and robustness analysis (Ben-Tal & Nemirovski,1998). However, with respect to robust approaches, where the goal is to determine a feasible solution which is optimal in some sense for all possible uncertainty instances, the goal in the stochastic framework is to find a solution that is feasible for almost all possible uncertainty realizations, (Calafiore, Dabbene, & Tempo,2011;Tempo, Calafiore, & Dabbene,2012). In several applications, including engineering and finance, uncertainties in ✩This research was partially funded by the project PRIN 2017, Prot. 2017S559BB, and has received support by the Agencia Estatal de Investigación (AEI)-Spain under Grant PID2019-106212RB-C41/AEI/10.13039/501100011033. The material in this paper was partially presented at the 18th European Control Conference, June 25–28, 2019, Naples, Italy. This paper was recommended for publication in revised form by Associate Editor Marcello Farina under the direction of Editor Ian R. Petersen. ∗Corresponding author. E-mail addresses: [email protected] (M. Mammarella), [email protected] (V. Mirasierra), [email protected] (M. Lorenzen), [email protected] (T. Alamo), [email protected] (F. Dabbene). price, demand, supply, currency exchange rate, recycle and feed rate, and demographic condition are common. In these situations, it is acceptable, up to a certain safe level, to relax the inherent conservativeness of robust constraints enforcing probabilistic constraints. More recently, this probabilistic approach has been used also in unmanned autonomous vehicle navigation (Li, Zhan, Hu, & Tomizuka,2020;Mammarella, Capello, Dabbene, & Guglieri, 2018) as well as in optimal power flow (Chamanbaz, Dabbene, & Lagoa,2019;Chamanbaz, Dabbene, & Lagoa,2020). In the optimization framework, constraints involving stochastic parameters that are required to be satisfied with a prespecified probability threshold are called chance constraints (CC). In general, dealing with CC implies facing two serious challenges: the solution of difficult parametrized probability integrals and the nonconvexity of the ensuing constraints (Geng & Xie,2019). Consequently, while being attractive from a modeling viewpoint, problems involving CC are often computationally intractable, generally shown to be NP-hard, which seriously limits their applicability. However, being able to efficiently solve or approximate chance-constrained problems remains an important challenge, especially in systems and control. In the case of approximated solutions, there exists of course a fundamental trade-off between complexity of the approach and goodness of the approximation. The scientific community has devoted large research in devising computationally efficient approaches to deal with chanceconstraints. We review such techniques in Section 3, where we highlight three mainstream approaches: (i) exact techniques; (ii) robust approximations; and (iii) sample-based approximations. https://doi.org/10.1016/j.automatica.2021.110108 0005-1098/©2021 Published by Elsevier Ltd. M. Mammarella, V. Mirasierra, M. Lorenzen et al. Automatica 137 (2022) 110108 In this paper, we present what we consider an important step forward in the sample-based approach. More precisely, our developments stem from the observation that, while in the general situation one is interested in finding an optimal solution to a chance-constrained problem, there exists a significant class of practical applications in which, instead, what it is really needed is being able to construct good approximation of the chanceconstrained set. This is the case, for instance, of stochastic model predictive control (SMPC), where this approximation is necessary for post-processing in real time, see for instance (Lorenzen, Dabbene, Tempo, & Allgöwer,2017b;Mammarella et al.,2018). Motivated by these considerations, we propose a simple and efficient strategy to obtain a probabilistically guaranteed inner approximation of a chance-constrained set, with given confidence. In particular, we describe a two-step procedure that involves: (i) the preliminary approximation of the chance-constrained set by means of a so-called Simple Approximating Set (SAS); (ii) a sample-used scaling procedure that allows to properly scale the SAS so to guarantee the desired probabilistic properties. The proper selection of a low-complexity SAS allows the designer to easily tune the complexity of the approximating set, significantly reducing the sample complexity. We propose several candidate SAS shapes, grouped in two classes: (i) sampled-polytopes; and (ii) norm-based SAS. The approach we propose distinguishes itself from the previous literature on CC in the following main points. (1) It is specifically tailored towards the specific problem of approximating the chance-constrained set, as opposed to solving a specific instance of a chance-constrained problem. (2) It is very general: it applies to a very general class of uncertainty configurations. A large part of the methods available in the literature are limited to cases where the constraints depend in a ‘‘nice’’ way on the uncertainty. This is the case for instance of the solution proposed in Margellos, Goulart, and Lygeros (2014) and Nemirovski and Shapiro (2006). The reader is referred to Section 3for an overview. (3) It is highly tunable: by selecting the complexity of the approximating set, the designer has a very efficient way to control the trade-off between computational complexity and potential goodness of the approximation. The probabilistic scaling approach was presented in the conference papers (Alamo, Mirasierra, Dabbene, & Lorenzen,2019; Mammarella, Alamo, Dabbene, & Lorenzen,2020) and it is based on recent results on order statistics (Alamo, Manzano, & Camacho, 2018). The present work extends Alamo et al. (2019), Mammarella et al. (2020) in several directions. First, we perform here a thorough mathematical analysis of probabilistic scaling. Second, we provide probabilistic guarantees for a more general class of norm-based SAS. Third, we consider here joint chance constraints. This choice is motivated by the fact that enforcing joint chance constraints, which have to be satisfied simultaneously, adheres better to some applications, despite the inherent complexity. Finally, we present here a second application, besides SMPC, related to probabilistic set-membership identification. The paper is structured as follows. Section 2provides a general preamble of the problem formulation and of chance-constrained optimization, including two motivating examples. An extensive overview on methods for approximating chance-constrained sets is reported in Section 3whereas the probabilistic scaling approach has been detailed in Section 4. Section 5and Section 6are dedicated to the definition of selected candidate SAS, i.e. sampledpolytope and norm-based SAS, respectively. Last, in Section 7, we validate the proposed approach with a numerical example applying our method to a probabilistic set membership estimation problem. Main conclusions and future research directions are addressed in Section 8. Notation Given an integer N,[N]denotes the integers from 1 to N. Given x∈R,⌊x⌋denotes the greatest integer no larger than x. Given the ℓp-norm ∥·∥p, we denote with Bs pthe ℓp-norm ball of radius one in Rs, i.e. Bs p . = { z∈Rs: ∥z∥p≤1}. The Chebyshev center of a given set X, with respect to norm ∥·∥p, is denoted as Chebp(X), and it is defined as the center of the largest ℓp-norm ball inscribed in X. Given an ℓp-norm ∥·∥p, its dual norm ∥·∥p∗ is defined as ∥c∥p∗. =supz∈Bs pc⊤z,∀c∈Rs. In particular, the couples (p,p∗): (2,2), (1,∞), (∞,1) give rise to dual norms (Boyd & Vandenberghe,2004, A.1.6). Given two sets S1and S2, the notation S1⊕S2(S1⊖S2) stands for the Minkowski sum (difference) between the two sets. Given integers k,N, and parameter ε∈ [0,1], the Binomial cumulative distribution function is denoted as B(k;N, ε). = k ∑ i=0(N i)εi(1 −ε)N−i. The following definition is borrowed from the field of order statistics (Alamo et al.,2018;Mirasierra, Mammarella, Dabbene, & Alamo,2022). Definition 1 (Generalized Min).Given a collection of Nscalars Γ= {γi}N i=1∈RN, and an integer r∈ [N], we say that γ− r∈Γis the r-smallest value of Γif there is no more than r−1 elements of Γstrictly smaller than γ− r. Hence, γ− 1denotes the smallest value in Γ,γ− 2the second smallest one, and so on until γ− N, which is equal to the largest one. We also use the alternative notation min(r)(Γ). =γ− r. 2. Problem formulation Consider a robustness problem, in which the controller parameters and auxiliary variables are parametrized by means of a decision variable vector θ∈Rnθ, which is usually referred to as design parameter. Furthermore, the uncertainty vector w∈Rnwrepresents one of the admissible uncertainty realizations of a random vector with given probability distribution PrWand (possibly unbounded) support W. This paper deals with the special case where the design specifications can be decoded as a set of nℓuncertain linear inequalities F(w)θ≤g(w),(1) where F(w)=⎡ ⎢ ⎣ f⊤ 1(w) . . . f⊤ nℓ(w) ⎤ ⎥ ⎦∈Rnℓ×nθ,g(w)=⎡ ⎢ ⎣ g1(w) . . . gnℓ(w) ⎤ ⎥ ⎦∈Rnℓ, are measurable functions of the uncertainty vector w∈Rnw. In Section 8we discuss possible extensions of this approach to more general settings, in which the constraints may be nonlinear and even nonconvex. Also note that the proposed setup captures the special case of chance constraints with random right-hand side. These correspond to the choice F(w)=Fand g(w)=w. Similarly, the case of chance constraints with random technology matrix is captured by our general case. We also note that hard linear constraints on θmay be directly incorporated by introducing deterministic inequalities of the form f⊤ ℓθ≤gℓ, where fℓand gℓdo not depend on the uncertainty w. The inequality in (1) is to be interpreted component-wise, i.e. f⊤ ℓ(w)θ≤gℓ(w), ∀ℓ∈ [nℓ]. Due to the random nature of 2 M. Mammarella, V. Mirasierra, M. Lorenzen et al. Automatica 137 (2022) 110108 the uncertainty vector w, each realization of wcorresponds to a different set of linear inequalities. Consequently, each value of wgives raise to a corresponding set X(w). = { θ∈Rnθ:F(w)θ≤g(w)}.(2) In every application, one usually accepts a risk of violating the constraints. This often translates into a two-step strategy: (i) a set ˜ Wis obtained such that w∈˜ Wis satisfied with a pre-specified high probability; (ii) a robust design problem in which inequality (1) is forced to be satisfied for every w∈˜ Wis solved. This is for instance the approach proposed in Hewing, Carron, Wabersich, and Zeilinger (2018) and Margellos et al. (2014). This approach suffers from several drawbacks: (i) there is no guarantee that the ensuing robust problem is easily solvable. Indeed, this may be in general very hard, and to obtain computable solutions the authors of Margellos et al. (2014) need to make additional assumptions on the dependence of F,gon the uncertainty w. (ii) The approach in Margellos et al. (2014) does not provide a safe region (i.e. a probabilistic approximation of the chance-constrained set), but just a point satisfying the probabilistic constraint. If one uses this approach, then the result may be conservative due to the two-step procedure. In this paper, we observe that a less conservative solution can be found by choosing the set Wto encompass all possible values and characterizing the region of the design space in which the fraction of elements of W, that violate the constraints, is below a specified level. This concept is rigorously formalized by means of the notion of probability of violation (see Tempo et al.,2012). Definition 2 (Probability of Violation).Consider a probability measure PrWover Wand let θ∈Rnθbe given. The probability of violation of θrelative to inequality (1) is defined as Viol(θ). =PrW{F(w)θ≤ g(w)}. Given a constraint on the probability of violation, i.e. Viol(θ)≤ε, we denote as (joint) chance-constrained set of probability ε (shortly, ε-CCS) the region of the design space for which this probabilistic constraint is satisfied. This is formally stated in the next definition. Definition 3 (ε-CCS).Given ε∈(0,1), we define the chanceconstrained set of probability εas follows Xε . = { θ∈Rnθ:Viol(θ)≤ε}.(3) Note that the ε-CCS represents the region of the design space Rnθ for which this probabilistic constraint is satisfied and it is equivalently defined as Xε . ={θ∈Rnθ:PrW{F(w)θ≤g(w)}≥1−ε}.(4) Remark 1 (Joint vs. Individual CC).The chance-constraint θ∈Xε, with Xεdefined in (4), describes a joint chance constraint. That is, it requires that the joint probability of satisfying the inequality constraint F(w)θ≤g(w) is guaranteed to be no smaller than the probabilistic level 1−ε. We remark that this constraint is notably harder to impose than individual CC, i.e. constraints of the form θ∈Xℓ εℓ . ={θ∈Rnθ:PrW{fℓ(w)⊤θ≤gℓ(w)}≥1−εℓ}, ℓ∈ [nℓ], with εℓ∈(0,1). A discussion on the differences and implications of joint and individual chance constraints may be found in several papers, see for instance (Geng & Xie,2019;Miller & Wagner,1965) and references therein. Note that a well-known (conservative) approximation to the joint chance-constrained set is the use of multiple individual CC. Fig. 1. The ε-CCS set for ε=0.15 (smaller set), ε=0.30 (intermediate set), and ε=0.45 (larger set). We observe that all sets are nonconvex, but the nonconvexity is more evident for larger values of ε, corresponding to larger levels of accepted violation, while the set Xεappears ‘‘almost convex’’ for small values of ε. This kind of behavior is in accordance with a recent result that proves convexity of the ε-CCS for small enough values of ε, and it is usually referred to as eventual convexity (van Ackooij,2015;van Ackooij & Malick,2019). While there exist simple examples for which a closed-form evaluation of Xεis possible, see e.g. Alamo et al. (2019, Figure 1), we remark that this is not the case in general. Indeed, as pointed out in Geng and Xie (2019), typically the computation of the ε- CCS is extremely difficult, since the evaluation of the probability Viol(θ) amounts to the solution of a multivariate integral, which is NP-Hard (Khachiyan,1989). Moreover, the set ε-CCS is often nonconvex, except for very special cases. For example, Shapiro, Dentcheva, and Ruszczyński (2014, Lemma 4.60) show that the solution set of separable chance constraints can be written as the union of cones, which is nonconvex in general. Example 1 (Example of Nonconvex ε-CCS).To illustrate these inherent difficulties, we consider the following three-dimensional example (nθ=3) with w={w1, w2}, where the first uncertainty w1∈R3is a three-dimensional normal-distributed random vector with zero mean and covariance matrix Σ=[4.5 2.26 1.4 2.26 3.58 1.94 1.4 1.94 2.19], and the second uncertainty w2∈R3is a three-dimensional random vector whose elements are uniformly distributed in the interval [0,1]. The set of viable design parameters is given by nℓ=4 uncertain linear inequalities of the form F(w)θ≤[1111]⊤,(5) F(w)=[w1w2(2w1−w2)w2 1]⊤, where the square power w2 1is to be interpreted element-wise. In this case, to obtain a graphical representation of the set Xε, we resorted to gridding the design set, and, for each point θin the grid, to approximate the probability through a Monte Carlo method. This procedure is clearly unaffordable for higher dimensions. In Fig. 1 we report the plot of the obtained ε-CCS set for different values of ε. We observe that the set is indeed nonconvex. 2.1. Chance-constrained optimization Finding an optimal θ∈Xεfor a given cost function J:Rnθ→ R, leads to the chance-constrained optimization (CCO) problem min θ∈Xε J(θ),(6) 3 M. Mammarella, V. Mirasierra, M. Lorenzen et al. Automatica 137 (2022) 110108 where the cost-function J(θ) is usually assumed to be a convex, often even quadratic or linear, function. We remark that the solution of the CCO problem (6) is in general NP-hard, for the same reasons reported before. We also note that several stochastic optimization problems arising in different application contexts can be formulated as a CCO. Typical examples are for instance the reservoir system design problem proposed in Prékopa, Rapcsák, and Zsuffa (1978), where the problem is to minimize the total building and penalty costs while satisfying demands for all sites and all periods with a given probability, or the cash matching problem (Dentcheva, Lai, & Ruszczyński,2004), where one aims at maximizing the portfolio value at the end of the planning horizon while covering all scheduled payments with a prescribed probability. CCO problems also frequently arise in short-term planning problems in power systems. These optimal power flow (OPF) problems are routinely solved as part of the real-time operation of the power grid. The aim is determining minimum-cost production levels of controllable generators subject to reliably delivering electricity to customers across a large geographical area, see e.g. Chamanbaz et al. (2020) and references therein. Recently, approaches based on a probabilistic approximation of chance-constrained sets have emerged in the context of Stochastic MPC, see Lorenzen, Dabbene, Tempo, and Allgöwer (2017a), Lorenzen et al. (2017b) and Mammarella et al. (2018). These approaches exploit the sample-based results we summarize in Section 3.3 to construct offline a probabilistically guaranteed approximation of the set of all couples of control input/initial states that guarantee fulfillment of the desired input/state constraints. The possibility of constructing the approximation offline constitutes a winning feature with respect to similar approaches based on samples, since it moves all the cumbersome computation to the control design phase. In the online implementation, the only operation to be performed is to ‘‘evaluate’’ this set in correspondence to the current initial state. In this way, the original stochastic optimization program is reduced to an efficiently solvable quadratic program. This represents an undiscussed advantage, which has been demonstrated for instance in Mammarella et al. (2018). We stress that the key element of this procedure is exactly the construction of a ‘‘good’’ approximation of the ε-CCS. In the next subsection, we report an additional motivating example, which further highlights the importance of the problem at hand. 2.2. Motivating example: probabilistic set-membership estimation Consider the problem of finding ¯ θ∈Rnθsuch that |y−¯ θTϕ(x)| ≤ ρ, ∀(x,y)∈W⊆Rnx×R, where ϕ:Rnx→Rnθis a (possibly non-linear) regressor function, and ρ > 0 is a given hyperparameter accounting for modeling errors. The (deterministic) set membership estimation problem, see Bravo, Alamo, and Camacho (2006), Puig (2010) and Vicino and Zappa (1996), consists of computing the set of parameters θ that satisfy the constraint |y−θTϕ(x)| ≤ ρfor all possible values of (x,y)∈W. In the literature, this set is usually referred to as the feasible parameter set, that is FPS . = { θ∈Rnθ: |y−θTϕ(x)| ≤ ρ, ∀(x,y)∈W}. We notice that FPS could be empty if ρis chosen too small. If, for given w=(x,y), we define the set X(w)= { θ∈Rnθ: |y−θTϕ(x)| ≤ ρ}, then the feasible parameter set FPS can be rewritten as FPS = { θ∈Rnθ:θ∈X(w),∀w∈W}. The deterministic set membership problem suffers from the following limitations in real applications: (i) due to the possible non-linearity of ϕ(·), checking if a given θ∈Rnθsatisfies the constraint θ∈X(w), for every w∈W, is often a difficult problem; (ii) in many situations the robust constraint can not be checked because only samples of Ware available and therefore, only outer bounds of FPS can be computed; and (iii) there are problem instances where, because of outliers and possible nonfinite support of W, there is no point in Rnθguaranteeing the fulfillment of every possible constraint, and thus, the set FPS is empty (even for large values of ρ). To deal with this issue, one can resort to a probabilistic relaxation of the FPS. If a probability distribution is defined on W, the probabilistic set membership estimation problem is that of characterizing the set of parameters θthat satisfy PrW{|y−θTϕ(x)| ≤ ρ} ≥ 1−ε, for a given probability parameter ε∈(0,1). Hence, we can define FPSεas the set of parameters that satisfy the previous probabilistic constraint, that is, FPSε= { θ∈Rnθ:PrW{θ∈X(w)} ≥ 1−ε}. It is immediate to notice that this problem fits in the formulation proposed in this section: It suffices to define F(w)=[ϕT(x) −ϕT(x)],g(w)=[ρ+y ρ−y]. 2.3. Chance-constrained set approximations Motivated by the discussion above, we are ready to formulate the main problem studied in this paper. Problem 1 (ε-CCS Approximation).Given the set of linear inequalities (1), and a violation parameter ε, find an inner approximation of the set Xε. The approximation should be: (i) simple enough; (ii) accurate enough, (iii) easily computable. A solution to this problem is provided in the paper. In particular, regarding (i), we present a solution in which the approximating set is represented by few linear inequalities. Regarding (ii) and (iii), we propose a highly tunable and computationally efficient procedure for its construction (see Algorithm 1). Before presenting our approach, in the next section we provide a brief literature overview of different methods presented in the literature to construct approximations of the ε-CCS set. 3. Overview on different approaches to ε-CCS approximations The construction of computational efficient approximations to ε-CCS is a long-standing problem. In particular, the reader is referred to the recent reviews (Geng & Xie,2019;Lejeune & Prékopa,2018), which provide rather complete discussions on the topic, and cover the most recent results. The authors of Geng and Xie (2019) distinguish three different approaches, which we very briefly revisit here. 3.1. Exact techniques In some very special cases, the ε-CCS is convex and hence the CCO problem is efficiently computable. This is the case, for instance, of individual chance constraints with wbeing Gaussian (Kataoka,1963). Other important examples of convexity of the set Xεinvolve log-concave distribution (Prékopa,1971, 1995). General sufficient conditions on the convexity of chance constraints may be found in van Ackooij (2015), Calafiore and Ghaoui (2006), Henrion and Strugarek (2008) and Lagoa (1999). 4 M. Mammarella, V. Mirasierra, M. Lorenzen et al. Automatica 137 (2022) 110108 However, all these cases are very specific and hardly extend to the joint chance constraints considered in this work. All these previously cited references deal with continuous distributions. A different line of research concentrates instead on discrete distributions, which arise frequently in applications, either directly, or as empirical approximations of the underlying distribution (see, for example, Beraldi & Ruszczyński,2002;Prékopa, 1995). For this particular case, exact results based on the concept of p-efficiency points (Dentcheva, Prékopa, & Ruszczynski,2000) or dual methods (Prékopa,1990) have been proposed. 3.2. Deterministic approximations A second class of approaches consist in finding deterministic conditions that allow to construct a set X, which is a guaranteed inner convex approximation of the probabilistic set Xε. The classical solution consists in the applications of Chebyshevlike inequalities, see e.g. Hewing and Zeilinger (2018) and Yan, Goulart, and Cannon (2018). More recent techniques, which are proved particularly promising, involve robust optimization (Ben- Tal & Nemirovski,1998), as the convex Bernstein-based approximations introduced in Nemirovski (2012) and Nemirovski and Shapiro (2006). A particular interesting convex relaxation involves the so-called Conditional Value at Risk (CVaR), see Chen, Sim, Sun, and Teo (2010) and references therein. Finally, we point out some recent techniques based on Genz’ code for Gaussian probabilities of rectangles (Bremer, Henrion, & Möller,2015), or on polynomial moments relaxations (Jasour, Aybat, & Lagoa, 2015;Lasserre,2017). Specific solutions have been proposed for the case of discrete distributions, see the recent survey (Ahmed & Xie,2018). In particular, we point out the recent works proposing a Boolean reformulation of the feasible set of individual and joint chance constraints (see Lejeune,2012;Lejeune & Margot,2016). Nonetheless, it should be remarked that these techniques usually suffer from conservatism and computational complexity issues, especially in the case of joint chance constraints. 3.3. Sample-based techniques In recent years, a novel approach to approximate chance constraints, based on random sampling of the uncertain parameters, has gained popularity, see e.g. Alamo, Tempo, and Camacho (2009), Calafiore et al. (2011), Luedtke and Ahmed (2008), Tempo et al. (2012) and references therein. Sampling-based techniques are characterized by the use of a finite number Nof iid samples of the uncertainty {w(1), w(2), . . . , w(N)}drawn according to a probability distribution PrW. With each sample w(i),i∈ [N], we can associate the following sampled set X(w(i))= { θ∈Rnθ:F(w(i))θ≤g(w(i))},(7) sometimes referred to as scenario, since it represents an observed instance of the uncertain constraint. Then, the scenario approach considers the CCO problem (6) and approximates its solution through the following scenario problem θ∗ sc =arg min J(θ) (8) subject to θ∈X(w(i)),i∈ [N]. We note that, if the function J(θ) is convex, problem (8) becomes a linearly constrained convex program, for which very efficient solution approaches exist. Under some technical assumptions (feasibility of the problem and non-degeneracy), a fundamental result (Calafiore,2010;Calafiore & Campi,2006;Campi & Garatti, 2008,2011) provides a probabilistic certification of the constraint satisfaction for the solution to the scenario problem. In particular, it is shown that PrWN{Viol(θ∗ sc )> ε}≤B(nθ−1;N, ε),(9) where the probability in (9) is measured with respect to the samples {w(1), w(2), . . . , w(N)}. A few observations are at hand regarding the scenario approach and its relationship with Problem 1. First, if we define the sampled constraints set as XN . = N ⋂ i=1 X(w(i)),(10) we see that the scenario approach consists in approximating the constraint θ∈Xεin (6) with its sampled version θ∈XN. On the other hand, it should be remarked that the scenario approach cannot be used to derive any guarantee on the relationship existing between XNand Xε. Indeed, the nice probabilistic property in (9) holds only for the optimum of the scenario program θ∗ sc . This is a fundamental point, since the scenario results build on the so-called support constraints, which are defined for the optimum point θ∗ sc only. On the contrary, in our case we are interested in establishing a direct relation (in probabilistic terms) between the set XNand the ε-CCS Xε. This is indeed possible, but one needs to resort to results based on Statistical Learning Theory (Vapnik,1998) and in Alamo et al. (2009, Theorem 8), summarized in the following lemma. Lemma 1 (Statistical Learning Theory Bound).Given probabilistic levels δ∈(0,1) and ε∈(0,0.14), if the number of samples N is chosen so that N ≥NLT , with NLT . =4.1 ε(ln 21.64 δ+4.39nθlog2(8enℓ ε)),(11) then PrWN{XN⊆Xε}≥1−δ. The lemma, whose proof is reported in Appendix A.1, is a direct consequence of the statistical learning theory results on the so-called (α, k)-Boolean functions, given in Alamo et al. (2009, Corollary 4), where more general results are reported for cases in which εis not constrained in (0,0.14). Remark 2 (Sample-based SMPC).The learning theory-based approach discussed in this section has been applied in Lorenzen et al. (2017b) to derive offline a probabilistic inner approximation of the chance-constrained set defining the couples of input/state guaranteeing the desired input/state chance. In particular, the bound (11) is a direct extension to the case of joint CC of the result proved in Lorenzen et al. (2017b) for individual CC. To this regard, we note that the results in the previous section allow to develop a novel SMPC scheme which considers multiple constraints at the same time. These developments are omitted here for brevity, and are reported in an extended version available at Mammarella, Mirasierra, Lorenzen, Alamo, and Dabbene (2021). This work highlights the limits of Lorenzen et al. (2017b): even for a moderately sized MPC problem with nx=5 states, nu=2 inputs, prediction horizon of T=10, simple interval constraints on states and inputs, and for a reasonable choice of probabilistic parameters ε=0.05, δ=10−6, we get more than 1.6 million linear inequalities. Remark 2 motivates the approach presented in the next section, which builds upon the results presented in Alamo et al. (2019). We show how the probabilistic scaling approach directly leads to approximations of user-chosen complexity, which can be directly used in applications instead of creating the need for a post-processing step to reduce the complexity of the sampled set. 5 M. Mammarella, V. Mirasierra, M. Lorenzen et al. Automatica 137 (2022) 110108 4. The probabilistic scaling approach We propose a novel sample-based approach, alternative to the randomized procedures proposed so far. This scheme allows to maintain the nice probabilistic features of these techniques, while at the same time providing the designer with a way of tuning the complexity of the approximation. The main idea behind this approach consists of a two-step procedure: (i) first a simple initial approximation θc⊕Sof the shape of the probabilistic set Xεis obtained, and (ii) a scalable simple approximating set (Scalable SAS) of the form S(γ). =θc⊕γS(12) is considered. These sets are described by a center point θcand a lowcomplexity shape set S. The center θcand the shape Sconstitute the design parameters of the proposed approach. By appropriately selecting the shape S, the designer can control the complexity of the approximating set. The nonnegative scalar γcontrols instead the scale of the set S(γ): the larger γ, the larger will be the set. Note that we do not ask this initial set to have any guarantee of probabilistic nature. What we ask is that this set is being able to ‘‘capture" somehow the shape of the set Xε. Recipes on a possible procedure for constructing this initial set are provided in Section 5. The center θcand the set Sconstitute the starting point of a scaling procedure, which allows to derive a probabilistic guaranteed approximation of the ε-CCS, as detailed in the next subsection. 4.1. Probabilistic scaling In this section, we address the problem of how to scale the set S(γ) around its center θcto guarantee, with confidence level δ∈ (0,1), the inclusion of the scaled set into Xε. Within this samplebased procedure we assume that Nγiid samples {w(1), . . . , w(Nγ)} are obtained from PrWand, based on these, we show how to obtain a scalar ¯γ > 0 such that PrWNγ{S(¯γ)⊆Xε} ≥ 1−δ. To this end, we first define the scaling factor associated to a given realization of the uncertainty. Definition 4 (Scaling Factor).Given a Scalable SAS S(γ), with given center θc∈Rnθand shape S, and a realization w∈W, we define the scaling factor of S(γ) relative to was γ(w). ={0 if θc∈ X(w) max S(γ)⊆X(w)γotherwise,(13) with X(w)= { θ∈Rnθ:F(w)θ≤g(w)}. That is, γ(w) represents the maximal scaling that can be applied to S(γ)=θc⊕γSaround the center θcso that S(γ)⊆X(w). The following theorem, whose proof is reported in Appendix A.2, states how to obtain, by means of sampling, a scaling factor ¯γ that guarantees, with high probability, that S(¯γ)⊆Xε. Theorem 1 (Probabilistic Scaling).Given a candidate Scalable SAS S(γ), with θc∈Rnθand shape S, suppose that accuracy parameter ε∈(0,1), confidence level δ∈(0,1), integer parameter r ≥1and Nγ≥r are chosen such that B(r−1;Nγ, ε)≤δ. (14) Draw Nγiid samples {w(1), w(2), . . . , w(Nγ)}from distribution PrW, compute the corresponding scaling factors γi . =γ(w(i)),i∈ [Nγ] according to Definition 4, and let Γ. = {γi}Nγ i=1. Define ¯γ=γ− r=min(r)(Γ), i.e. ¯γis the rth smallest value of Γ(see Notation). Under these assumptions: (i) If ¯γ > 0then, with probability no smaller than 1−δ, S(¯γ)=θc⊕ ¯γS⊆Xε. (ii) If θc∈ Xεthen ¯γ=0with probability no smaller than 1−δ. We now state a Lemma, whose proof is reported in Appendix A.3, which is a generalization of a previous result presented in Mirasierra et al. (2022). The lemma provides a way to relate the choice of rand Nγby introducing the value βϵNγ. This latter represents the desired fraction of the expected number of violations, which can be interpreted as a trade-off parameter between the number of samples and the tightness of the solution. Lemma 2. Let r = ⌈βεNγ⌉, where β∈(0,1), and define κ. =(√β+√2−β √2(1 −β))2 . Then, inequality (14) is satisfied for Nγ≥κ εln 1 δ.(15) In particular, the choice β=0.5leads to r =⌈εNγ 2⌉and Nγ≥ 7.47 εln 1 δ. The above result leads to the following simple algorithm, in which we summarize the main steps for constructing the scaled set, and we provide an explicit way of determining parameter r. Algorithm 1 Probabilistic SAS Scaling 1: Given a candidate Scalable SAS S(γ), and probability levels ε and δ, choose Nγ≥7.47 εln 1 δand r=⌈εNγ 2⌉.(16) 2: Draw Nγsamples of the uncertainty w(1), . . . , w(Nγ). 3: for i=1 to Nγdo 4: Compute, according to Definition 4, the Nγscaling factors γi . =γ(w(i)),i∈ [Nγ].(17) 5: end for 6: Return ¯γ=γ− r=min(r)(Γ), the rth smallest value of Γ= {γi}Nγ i=1. A few comments are in order regarding the algorithm above. In step 4, for each uncertainty sample w(i)one has to solve an optimization problem, which amounts to finding the largest value of γsuch that S(γ) is contained in the set X(w(i)) defined in (7). If the SAS is chosen appropriately, we can show that this problem is convex and computationally very efficient: this is discussed in Section 5. Then, in step 6, one has to re-order the sequence Γ= {γ1, γ2, . . . , γNγ}so that the first element is the smallest one, the second element is the second smallest one, and so on and so forth, and then return the rth element of the reordered sequence. The properties of the output of Algorithm 1can be derived by a direct application of Theorem 1 and Lemma 2. In particular, if 6 M. Mammarella, V. Mirasierra, M. Lorenzen et al. Automatica 137 (2022) 110108 the output ¯γis larger than zero, then S(γ)⊆Xεwith probability no smaller than 1 −δ. In the next sections, we provide a ‘‘library’’ of possible candidate SAS shapes. We remind that these sets need to comply with two main requirements: (i) being a simple and low-complexity representation; and (ii) being able to capture the original shape of the ε-CCS. Moreover, in the light of the discussion after Algorithm 1, we also ask these sets to be convex. 5. Candidate SAS: Sampled-polytope First, we note that the most straightforward way to design a candidate SAS is again to recur to a sample-based procedure: we draw a fixed number NSof ‘‘design’’ uncertainty samples1 {˜w(1),..., ˜w(NS)}, and construct an initial sampled approximation by introducing the sampled-polytope XNS= NS ⋂ j=1 X(˜w(j)).(18) Note that the sampled polytope XNS, by construction, is given by the intersection of nℓNShalf-spaces. Hence, we observe that this approach provides very precise control on the final complexity of the approximation, through the choice of the number of samples NS. However, it is also clear that a choice for which NS≪ NLT implies that the probabilistic properties of XNSbefore scaling will probably not meet the required specifications (see Lemma 1). However, we emphasize again that this initial geometry does not have nor require any probabilistic guarantees, which are instead provided by the probabilistic scaling discussed in Section 4.1. It should be also remarked that this is only one possible heuristic. For instance, along this line one could as well draw many samples and then apply a clustering algorithm to boil it down to a desired number of samples. We remark that, in order to apply the scaling procedure, we need to define a center around which to apply the scaling procedure. To this end, we could compute the so-called Chebyshev center, which for a given norm ∥·∥p, is defined as the center of the largest ball inscribed in XNS, i.e. θc=Chebp(XNS). Once the center θchas been determined, the scaling procedure can be applied to the sampled-polytope SAS defined as SNS(γ). =θc⊕γ{XNS⊖θc}.(19) We note that computing the Chebyshev center of a given polytope is an easy convex optimization problem, for which efficient algorithms exist, see e.g. Boyd and Vandenberghe (2004). A possible alternative would be the analytic center of XNS, whose computation is even easier (see Boyd & Vandenberghe,2004, for further details). Note that the choice of θconly affects the goodness of the shape, but we can never know a priori if the analytic center is a better choice than any random center in the candidate SAS. Example 2 (Sample-based Approximations).To illustrate how the proposed scaling procedure works in practice in the case of sampled-polytope SAS, we revisit Example 1. To this end, a prefixed number NSof uncertainty samples were drawn, and the set of inequalities F(˜w(j))θ≤g(˜w(j)),j∈ [NS], with F(w),g(w) defined in (5), were constructed, leading to the set XNS. Then, its Chebyshev center with respect to norm ∥·∥2was 1These samples are denoted with a tilde to distinguish them from the samples used in the probabilistic scaling procedure. computed, and Algorithm 1was applied to the sampled-polytope SAS SNS(γ) defined in (19), with ε=0.05, δ=10−6, leading to Nγ=2065. We note that, in this case, the solution of the optimization problem in (13) may be obtained by bisection on γ. Indeed, for given γ, checking if SNS(γ)⊆X(w(i)) amounts to solving some simple linear programs. Two different situations were considered: a case where the number of inequalities is rather small NS=100, and a case where the complexity of the SAS is higher, i.e. NS=1000. The outcome procedure is illustrated in Fig. 2. We can observe that, for a small NS–Fig. 2(a) – the initial approximation is rather large (although it is contained in Xε, we remark that we do not have any guarantee that this will happen). In this case, the probabilistic scaling returns γ=0.8954 which is less than one. This means that, in order to obtain a set fulfilling the desired probabilistic guarantees, we need to shrink it around its center. In the second case, for a larger number of sampled inequalities – Fig. 2(b) – the initial set (the red one) is much smaller, and the scaling procedure inflates the set by returning a value of γgreater than one, i.e. γ= 1.2389. Note that choosing a larger number of samples for the computation of the initial set does not imply that the final set will be a better approximation of the ε-CCS. Finally, we compare this approach to the scenario-like ones discussed in Section 3.3. To this end, we also draw the approximation obtained by directly applying the Statistical Learning Theory bound (11). Note that in this case, since nθ=3 and nℓ=4, we need to take NLT =13,011 samples, corresponding to 52,044 linear inequalities. The resulting set is represented in Fig. 2(c). We point out that using this approximation (i) the set is much more complex, since the number of involved inequalities is much larger; (ii) the set is much smaller, hence providing a much more conservative approximation of the ε-CCS. Hence, the ensuing chance-constrained optimization problem will be computationally harder, and lead to a solution with a larger cost or even to an infeasible problem, in cases where the approximating set is too small. 6. Candidate SAS: Norm-based SAS In this section, we propose a procedure in which the shape of the scalable SAS is selected a-priori. This corresponds to situations where the designer wants to have full control over the final shape in terms of structure and complexity. The main idea is to define so-called norm-based SAS of the form Sℓp(γ). =θc⊕γHBs p,(20) where Bs pis an ℓp-ball in Rs,H∈Rnθ×s, with s≥nθ, is a design matrix (not necessarily square), and γis the scaling parameter. Note that when the matrix His square (i.e. s=nθ) and positive definite these sets belong to the class of ℓp-norm based sets originally introduced in Dabbene, Lagoa, and Shcherbakov (2010). In particular, in case of ℓ2norm, the sets are ellipsoids. This particular choice is the one studied in Mammarella et al. (2020). Here, we extend this approach to a much more general family of sets, which encompasses for instance zonotopes, obtained by letting p= ∞ and s≥nθ. Zonotopes have been widely studied in geometry, and have found several applications in systems and control, in particular for problems of state estimation and robust Model Predictive Control, see e.g. Le, Stoica, Alamo, Camacho, and Dumur (2013). Zonotopes proved to be very flexible and yield very efficient implementations. 7 M. Mammarella, V. Mirasierra, M. Lorenzen et al. Automatica 137 (2022) 110108 Fig. 2. (a–b) Probabilistic scaling approximations of the ε-CCS. Scaling procedure applied to a sampled-polytope with NS=100 (a) and NS=1000 (b). The initial sets are depicted in red, the scaled ones in green. (c) Approximation obtained by direct application of Lemma 1. Note that, in this latter case, to plot the set without out-of-memory errors a pruning procedure (Herceg, Kvasnica, Jones, & Morari,2013) of the 52,044 linear inequalities was necessary. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 6.1. Scaling factor computation for norm-based SAS We recall that the scaling factor γ(w) is defined as 0 if θc∈ X(w) and as the largest value γfor which Sℓp(γ)⊆X(w) otherwise. The following theorem, whose proof is reported in Appendix A.4, provides a direct and simple way to compute in closed form the scaling factor for a given candidate norm-based SAS. Theorem 2 (Scaling Factor for Norm-based SAS).Given a normbased SAS S(γ)=θc⊕γHBs pand a realization w∈W, define τℓ(w). =gℓ(w)−fT ℓ(w)θcand ρℓ(w). = ∥HTfℓ(w)∥p∗, with ∥·∥p∗ being the dual norm of ∥·∥p. The scaling factor γ(w)can be computed as γ(w)=min ℓ∈[nℓ]γℓ(w), with γℓ(w),ℓ∈ [nℓ], given by γℓ(w)=⎧ ⎪ ⎨ ⎪ ⎩ 0if τℓ(w)<0, ∞if τℓ(w)≥0and ρℓ(w)=0 τℓ(w) ρℓ(w)if τℓ(w)≥0and ρℓ(w)>0 . Note that γ(w) is equal to zero if and only if θcis not included in the interior of X(w). 6.2. Construction of a candidate norm-based set Similarly to Section 5, we first draw a fixed number NSof ‘‘design’’ uncertainty samples {˜w(1),..., ˜w(NS)}, and construct an initial sampled approximation (sampled-polytope SAS) XNSby means of (18). Again, we consider the Chebyshev center of XNS, or its analytical center as a possible center θcfor our approach. Analogously to what was proposed in Mammarella et al. (2020), given XNS,s≥nθand p∈ {1,2,∞}, the objective is to compute the largest set θc⊕HBs pincluded in XNS. To this end, we assume that we have a function Volp(H) that provides a measure of the size of HBs p. That is, larger values of Volp(H) are obtained for increasing sizes of HBs p. Remark 3 (On the Volume Function).The function Volp(H) may be seen as a generalization of the classical concept of Lebesgue volume of the set XNS. Indeed, when His a square positive definite matrix, some possibilities are Volp(H)=log det(H) – which is directly proportional to the classical volume definition, or Volp(H)=tr H– which for p=2 becomes the well known sum of ellipsoid semiaxes (see Boyd & Vandenberghe,2004, Chapter 8 and Dabbene, Henrion, Lagoa, & Shcherbakov,2015). These measures can be easily generalized to non square matrices. It suffices to compute the singular value decomposition. If H= UΣVT, we could use the measures Volp(H)=tr Σor Volp(H)= log det(Σ). For non square matrices H, specific results for particular values of pare known. For example, we remind that if p= ∞ and H∈Rnθ×s,s≥nθ, then θc⊕HBs ∞is a zonotope. Then, if we denote as generator each of the columns of H, the volume of a zonotope can be computed by means of a sum of terms (one for each different way of selecting nθgenerators out of the sgenerators of H); see Alamo, Bravo, and Camacho (2005) and Gover and Krikorian (2010). Another possible measure of the size of a zonotope θc⊕HBs ∞is the Frobenius norm of H(Alamo et al., 2005). Given an initial design set XNS, we elect as our candidate Scalable SAS the largest ‘‘volume’’ norm-based SAS contained in XNS. Formally, this rewrites as the following optimization problem max θc,H Volp(H) subject to θc⊕HBs p⊆XNS. (21) As it has been shown (see Appendix A.4), problem (21) is equivalent to min θc,H−Volp(H) s.t. fT ℓ(˜w(j))θc+∥HTfℓ(˜w(j))∥p∗−gℓ(˜w(j))≤0, ℓ∈ [nℓ],j∈ [NS],(22) where we have replaced the maximization of Volp(H) with the minimization of -Volp(H). We notice that the constraints are convex on the decision variables; also, the functional to minimize is convex under particular assumptions. For example when His assumed to be square and positive definite and Volp(H)=log det(H). For non square matrices, the constraints remain convex, but the convexity of the functional to be minimized is often lost. In this case, local optimization algorithms should be employed to obtain a possibly sub-optimal solution. Example 3 (Norm-based SAS).We revisit again Example 1 to show the use of norm-based SAS. We note that, in this case, the designer can control the approximation outcome by acting 8 M. Mammarella, V. Mirasierra, M. Lorenzen et al. Automatica 137 (2022) 110108 Fig. 3. Scaling procedure applied to (a) Sℓ1-SAS with NS=100, (b) Sℓ1-SAS with NS=1000, (c) Sℓ∞-SAS with NS=100, and (d) Sℓ∞-SAS with NS=1000. The initial set is depicted in red, the final one in green. The sampled design polytope SNSis represented in black. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) upon the number of design samples NSused for constructing the set SNS. In Fig. 3 we report two different norm-based SAS, respectively with p=1 and p= ∞, and for each of them we consider two different values of NS, respectively NS=100 and NS=1000. Similarly to what observed in Example 2, we see that for larger NS, the ensuing initial set becomes smaller. Consequently, we have a shrinkage process for small NSand an inflating one for large NS. However, we observe that in this case, the final number of inequalities is independent of NS(8 for Sℓ1 and 6 for Sℓ∞). 6.2.1. Relaxed computation It is worth remarking that the minimization problem of the previous subsection might be infeasible. In order to guarantee the feasibility of the problem, a soft-constrained optimization problem is proposed. With a relaxed formulation, θcis not guaranteed to satisfy all the sampled constraints. However θc∈SNSis not necessary to obtain an ε-CCS. Given ξ > 0, the relaxed version of optimization problem (22) is min θc,H,η1,...,ηNS−Volp(H)+ξ NS ∑ j=1 max{ηj,0}(23) s.t. fT ℓ(˜w(j))θc+∥HTfℓ(˜w(j))∥p∗−gℓ(˜w(j))≤ηj, ℓ∈ [nℓ],j∈ [NS]. The parameter ξserves to provide an appropriate trade off between satisfaction of the sampled constraints and the size of the obtained region. A possibility to choose ξwould be to choose it in such a way that the fraction of violations nviol/NS(where nviol is the number of elements ηjlarger than zero) is smaller than ε/2. 7. Numerical example: Probabilistic set membership estimation We now present a numerical example in which the results of the paper are applied to the probabilistic set membership estimation problem, introduced in Section 2.2. We consider the universal approximation functions given by Gaussian radial basis function networks (RBFN) (Buhmann,2000). Given the nodes [x1,x2,...,xM]and the variance parameter c, the corresponding Gaussian radial basis function network is defined as RBFN(x, θ)=θTϕ(x), where θ=[θ1. . . θM]Trepresents the weights and ϕ(x)=[exp (−∥x−x1∥2 c). . . exp (−∥x−xM∥2 c) ]T Fig. 4. Representation of the extreme values θ+and θ−and the central value θcof the FPSδ ε. is the regressor function. Given δ∈(0,1) and ε∈(0,1), the objective is to obtain, with probability no smaller than 1 −δ, an inner approximation of the probabilistic feasible parameter set FPSε, which is the set of parameters θ∈RMthat satisfies PrW{|y−θTϕ(x)| ≤ ρ} ≥ 1−ε, (24) where ρ=5, xis a random scalar with uniform distribution in [−5,5], and y=sin(3x)+σ, where σis a random scalar with a normal distribution with mean 5 and variance 1. We use the procedure detailed in Sections 4–6 to obtain a SAS of FPSε. We have taken a grid of M=20 points in the interval [−5,5]to serve as nodes for the RBFN, and a variance parameter of c=0.15. We have taken NS=350 random samples w=(x,y) to compute the initial geometry, which has been chosen to be an ℓ∞norm-based SAS of dimension 20 with a relaxation parameter of ξ=1 (see (23)). The chosen initial geometry is θc⊕HB20 ∞, where His constrained to be a diagonal matrix. When the initial geometry is obtained, we scale it around its center by means of probabilistic scaling with Algorithm 1. The number of samples required for the scaling phase to achieve ε= 0.05 and δ=10−6is Nγ=2065 and the resulting scaling factor is γ=0.3803. The scaled geometry θc⊕γHB20 ∞is, with a probability no smaller than 1 −δ, an inner approximation of FPSεwhich we will refer to as FPSδ ε. Since it is a transformation of an ℓ∞norm ball with a diagonal matrix H, we can write it as FPSδ ε= {θ:θ−≤θ≤θ+}, where the extreme values θ−, θ+∈R20 are represented in Fig. 4, along with the central value θc∈R20. 9