Computation of Input Disturbance Sets for Constrained Output Reachability

Linear models with additive unknown-but-bounded input disturbances are extensively used to model uncertainty in robust control systems design. Typically, the disturbance set is either assumed to be known a priori or estimated from data through set-membership identification. However, the problem of computing a suitable input disturbance set in case the set of possible output values is assigned a priori has received relatively little attention. This problem arises in many contexts, such as in supervisory control, actuator design, decentralized control, and others. In this paper, we propose a method to compute input disturbance sets (and the corresponding set of states) such that the resulting set of outputs matches as closely as possible a given set of outputs, while additionally satisfying strict (inner or outer) inclusion constraints. We formulate the problem as an optimization problem by relying on the concept of robust invariance. The effectiveness of the approach is demonstrated in numerical examples that illustrate how to solve safe reference set and input-constraint set computation problems.


I. INTRODUCTION
The theory of set invariance plays a key role in the analysis of uncertain dynamical systems, as it provides the tools for the synthesis of robust controllers that can satisfy constraints in the presence of disturbances [1]. Of particular interest are Robust Positive Invariant (RPI) sets [2], the characterization and computation of which has been a very active area of research [3]- [5]. RPI sets are used to provide robust stability and constraint satisfaction guarantees of various robust Model Predictive Control (MPC) and Reference Governor (RG) schemes [6]- [8]. These guarantees are usually established using the maximal robust positive invariant (MRPI) set [5], which is the largest RPI set included in the constraint set. The minimal RPI (mRPI) [1] set, which is the smallest RPI set for a given disturbance set [5], is used to design trajectory tubes in robust MPC [9], and to analyze the existence of MRPI sets. It was shown that an output-admissible RPI set exists for the system if and only if the mRPI set is included in the constraint set [10]. In order to enforce this inclusion, several methods were proposed in the literature to design feedback controllers that sufficiently attenuate the effects of disturbances [11], [12]. On the other hand, in applications such as fault-tolerant control [13], RPI sets that include a given set are computed and used for sensor fault isolation. All the aforementioned applications were developed under the assumption that the disturbance set is known a priori.
In many practical cases, however, while the set of admissible states can be estimated from sensor measurements or pre-specified from given constraints to be satisfied, the disturbance set is unknown, leaving the designer the task of suitably defining it, especially in case one must satisfy a given set of constraints on the system, e.g., encoding known physical limitations, or undesired states. Depending on the considered setting, one might be interested in designing a disturbance set which is either as large as possible or as small as possible, while guaranteeing that the prescribed constraints are satisfied in all circumstances. For example, in a decentralized MPC The authors are with the IMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100 Lucca, Italy. Email: s.mulagaleti@imtlucca.it application such as that presented in [12], [14], the dynamic coupling between subsystems is modeled as an additive disturbance. Then, since the disturbance is a combination of the states of the neighbors, it is desirable to obtain a large disturbance set, since this implies that the coverage of the available constraint space is maximized. Another example in which a large disturbance set is desired in presented in [15], where the disturbance set represents the set of feasible tracking references. On the other hand, if the disturbance set represents the inputs that can be applied to the system and one wants to design the actuators such that the system will be able to reach a prespecified set of outputs, computing the smallest disturbance set is of interest. Moreover, in disturbance identification techniques such as those presented in [16], [17], one is interested in obtaining a small disturbance set that can explain the data.
In this paper, we propose a method to compute a set of input disturbances acting on a dynamical system such that the resulting output set approximately matches an assigned one. Possible applications include, but are not limited to: supervisory control and decentralized MPC design, to enforce that the output of the system stays within a given set; actuator design, to size the range of the actuators so that the output covers a given range; the characterization of the robustness of a system with respect to external disturbances. This method is centered on the formulation of an optimization problem, with the input disturbance set being the unknown and the approximation error between the obtained and assigned output sets being the objective function to minimize.
We propose the formulation of the optimization problem for linear systems and polytopic sets: since the construction of the output set requires the computation of an RPI set, we extend the results of [18], [19] to encode the computation of a minimal parametrized polytopic RPI set within the optimization problem. Then, we propose to use the penalty-function method presented in [20] to solve the resulting bilevel linear program. Finally, we show the effectiveness of the approach through numerical examples related to safe reference set and input-constraint set computation problems.
The paper is organized as follows. We introduce some notation and recall basic definitions regarding set operations in Section II. Then, we introduce the problem we solve, along with relevant results to obtain a bilevel programing formulation in Section III. In Section IV we present the main results that permit the implementation of the RPI constraint. In Section V, we discuss the encoding of the inclusion constraints, following which in Section VI, we discuss the implementation of the penalty function method to solve the bilevel LP. Finally, in Section VII we present some numerical results along with some application demonstrations.

II. NOTATION AND SET OPERATIONS
Consider the sets X , Y ⊂ R n , and vectors a ∈ R na and b ∈ R n b . Given a matrix L ∈ R n×m , we denote by LX the image {y ∈ R m : y = Lx, x ∈ X } of X under the linear transformation induced by L. We denote the i-th row of matrix L by L i , the rank of L by rank(L), the image-space of L by Im(L), and the null-space of L by null(L). Given a square matrix L ∈ R n×n , ρ(L) denotes its spectral radius. The set B n p := {x : x p ≤ 1} is the unit p-norm ball in R n . A polyhedron is the intersection of a finite number of half-spaces, and a polytope is a compact polyhedron. Given two matrices L, M ∈ R n×m , L ≤ M denotes element-wise inequality. The symbols 1, 0, and I denote all-ones, all-zeros and identity matrix respectively, with dimensions specified if the context is ambiguous. The set of natural numbers between two integers m and n, m ≤ n, is denoted by I n m := {m, . . . , n}. The Minkowski set addition is defined as X ⊕ Y := {x + y : x ∈ X , y ∈ Y}. The Cartesian product is defined as X × Y := {[x y ] : x ∈ X , y ∈ Y}. The support function of a compact set X ⊆ R n for a given y ∈ R n is defined as h X (y) := max x∈X y x. Let X and Y be polytopes in R n . Then, support functions are positively homogeneous, i.e., h αX (y) = αh X (y) for any scalar α ≥ 0. Moreover, for any vector y ∈ R n , we have h We use the Hausdorff distance between polytopes X and Y defined as

III. PROBLEM DEFINITION AND APPROXIMATIONS
Consider the linear time-invariant discrete-time system with state x ∈ R nx , output y ∈ R ny and disturbance w ∈ R nw . Given a polytopic set Y := {y : Gy ≤ g} of outputs with g ∈ R m Y , our goal is to compute a disturbance set W such that Y is "reachable" by the output y, in a sense which we will define precisely later. We refer to w as a "disturbance", as it is customary in the literature of uncertain systems. Depending on the application, however, it could also represent a set of command inputs, as we will show through application examples. We work with the following standing assumption. Assumption 1: Matrix A is strictly stable, i.e., ρ(A) < 1. In this paper, we focus on the computation of a disturbance set W parametrized as the polytope W( w ) := {w : F w ≤ w } with w ∈ R m W . We further assume that the row vectors F i ∈ R nw of matrix F are given a priori, and restrict our attention to computing vector w . For simplicity, we also enforce that 0 ∈ W, which is equivalent to w ≥ 0. In the next section, we present a method to relax this restriction, i.e., permit the computation of a disturbance set W that does not contain the origin.
Given a disturbance set W( w ), the forward computation problem, which is typically tackled in the literature [5], [21], entails computing a suitable Robust Positive Invariant (RPI) set X := {x : Ax + Bw ∈ X , ∀ w ∈ W( w )}. Of particular interest is the computation of tight RPI approximations of the minimal RPI (mRPI) set Xm( w ), which is contained in every closed RPI set. It is given by the infinite Minkowski sum If w ≥ 0, i.e., W( w ) contains the origin, then Xm( w ) exists, is compact, convex and unique, and contains the origin [5]. Moreover, it is the limit of all state trajectories of (1a) under persistent disturbances w ∈ W( w ) [1]. Then, the corresponding limit set of output trajectories is Ym( w ) := CXm( w ) ⊕ DW( w ) as per (1b). The set Ym( w ) exists, is compact and convex with 0 ∈ Ym( w ) if w ≥ 0.
In this paper, we tackle the reverse computation problem, i.e., given an output polytope Y, compute the vector w such that Ym( w ) = Y. This problem, however, might not have a solution, i.e., there might not exist any w satisfying the output-set equality because of either of the following two reasons. (1) The mRPI set Xm( w ) is not finitely determined, except in a few special cases, e.g., nilpotent systems [5]. Then, depending on the set Xm( w ) and the structure of matrix C, the set Ym( w ) might also not be finitely determined. In this case, enforcing Ym( w ) = Y, with Y defined using a finite number of hyperplanes is not possible. (2) Even if the set Xm( w ) and matrix C are such that Ym( w ) is finitely determined, its shape is in general not arbitrary, but is a function of the dynamics of system (1a) and parametrization of set W( w ). This implies that enforcing Ym( w ) = Y, with Y being a user-specified arbitrarily shaped polytope, might not be possible. Hence, we instead tackle the problem This formulation includes the case Ym( w ) = Y, which holds if and only if d H (Ym( w ), Y) = 0. In the rest of this paper, we present a formulation to approximately solve Problem (3), where the approximation results from Xm( w ) not being finitely determined. Remark 1: We present a brief discussion regarding uniqueness of the set W( w ) if the equality Ym( w ) = Y holds. In case ny < nx + nw, then there exist infinitely many solutions for the equation [C D][x w ] = y for each y ∈ Y, such that the set W( w ) is nonunique. If instead ny ≥ nx + nw and rank([C D]) = nx + nw, then there exists a unique pair (x, w) corresponding to each y. By construction, this implies that W( w ) is unique. Moreover, in case D = 0, ny ≥ nx, and rank(C) = nx, the set Xm( w ) = C † Y is uniquely defined, where C † is the left-inverse of C. Then, if matrix B has a left-inverse, it can be shown that W( w ) is unique. A unified theory that includes all these cases to establish the conditions for the uniqueness of W( w ) is a subject of future research.

A. Polytopic RPI set
In order to approximate the mRPI set, we consider a parametrized state-set X( x ) := {x : Ex ≤ x } with x ∈ R m X and matrix E given a priori. Then, we enforce that X( x ) is RPI for system (1a) with disturbance set W( w ), i.e., it satisfies the inclusion AX( x ) ⊕ BW( w ) ⊆ X( x ). In order for such an x to exist, however, the matrix E must satisfy some requirements, that we formulate next.
Firstly, we define the support functions . Without loss of generality, we assume that matrix E is chosen such that b(1) = 1. Then, we make the following assumption regarding the existence of an RPI set.
Assumption 2: Matrix E is chosen such that there exists anˆ x ≥ 0 satisfying the inequality c(ˆ x ) + 1 ≤ b(ˆ x ). Assumption 2 implies that there exists an RPI set X(ˆ x ) for the system x(t + 1) = Ax(t) +w(t) with disturbancesw ∈ X(1). In the following result, we show that there always exists an RPI set X( x ) for system (1a) with the disturbance set W( w ).
Proposition 1: Suppose Assumption 2 holds, then there always exists a vector x ≥ 0 for every w ≥ 0 such that the RPI condition Proof: If Assumption 2 holds, then by duality of linear programs and Farkas' lemma [1], there exist nonnegative multiplier matriceŝ Λc,Λ b ∈ R m X ×m X satisfying the relationshipŝ For a given w ≥ 0, there exists an x ≥ 0 satisfying the RPI condition c( Then, we see that setting x = d( w ) ∞ˆ x satisfies the relationships in (5) with Λc =Λc and Λ b =Λ b . Remark 2: In order to verify if Assumption 2 holds, one can solve the LP (8) in [19]: The LP is bounded if and only if Assumption 2 holds. An iterative procedure to obtain a matrix E that satisfies this requirement was presented in [22].
Remark 3: The choice of polytopic parametrizations with fixed hyperplanes for the disturbance set and the corresponding RPI set is motivated primarily by their computational convenience. In particular, the choice of matrix F is completely independent of system (1), while matrix E must satisfy Assumption 2, which depends on system (1). Conservatism introduced by this parametrization can potentially be reduced by also optimizing over the hyperplanes through the introduction of optimization variables E and F . We note that the results presented in the rest of this paper continue to hold in the presence of these additional variables. Further, alternative convex parametrizations such as ellipsoidal and zonotopic sets [1], [23] can also be considered. Embedding the computation of small RPI sets with such parametrizations within an optimization problem is a subject of future research. (3) Having established the existence of a polytopic RPI set X( x ) under Assumption 2, we will now proceed with approximating Problem (3) using this set. To that end, we first note that the inclusion Xm( w ) ⊆ X( x ) holds by the definition of the mRPI set. Then, the

B. Approximating Problem
by the same definition. Based on this set, we approximate Problem (3) as the bilevel programing problem in which the disturbance set is computed by the upper-level problem, for which a corresponding RPI set X( x ) is computed by the lowerlevel problem in (6b). The lower-level problem is formulated in such a way that X( x ) is the tightest RPI approximation of the mRPI set, as seen in the objective function d H (X( x ), Xm( w )). The constraintset of this problem is nonempty according to Proposition 1, and all feasible¯ The rationale for formulating this problem follows from the triangle inequality: For a given w ≥ 0 and x ≥ 0, the inequality holds. The second part of the inequality is minimized by the upperlevel problem in (6). With respect to the first part, let us define From basic properties of Minkowski algebra, it follows that Since d x H is minimized by the lower-level problem (6b), Problem (6) minimizes an upper bound to Problem (3) as Finally, in order to eliminate the mRPI set from Problem (6b) we use the following results from [18], which state that the solution of Problem (6b) can be obtained using fixed-point iterations. In recalling these results, we denote d( w ) by d for ease of notation.
1) The sequence generated by the iterative procedure x . This fixed-point satisfies the equalities 2) The fixed-point reached from the initial-condition for all¯ x ∈ E(d( w )). Since this solution also has the smallest 1norm value over all¯ x ∈ E(d( w )), Problem (6) is equivalent to , there exist problem settings which have stronger requirements with respect to inclusions of the output-set Y.
Example 1: Consider the case in which (1a) represents a linear system equipped with a stabilizing feedback controller, such that w represents the feedforward reference signal. The system is subject to constraints y ∈ Y. Then, an inner-approximation version of Problem (8) enforces Ym( w ) ⊆ Y to compute the set of references W( w ) that satisfy the system constraints.
Example 2: Consider the case in which system (1) is used to model a dynamic disturbance y(t), e.g., pedestrian behavior [24] or an asset price, possibly estimated from linear time-series analysis. Let Y represent experimental data from the disturbance-generating system. Then, an outer-approximation version of Problem (8) enforces Ym( w ) to compute a disturbance set W( w ), using which a realistic simulator can be designed.
In order to formulate the inner-approximation version of Problem (8), we note that enforcing In order to formulate the outer-approximation version of Problem (8), we note that ⊕ N t=0 A t BW( w ) ⊆ Xm( w ) for all N ≥ 0. Then, the inclusion Y ⊆ Ym( w ) can be enforced by choosing an index N ≥ 0, and appending the constraint Y ⊆ ⊕ N t=0 CA t BW( w ) ⊕ DW( w ). Hence, we formulate We now formulate the assumptions that the output-set Y must satisfy in order to guarantee feasibility of Problems (9) Under Assumption 3-Inner, vector g ≥ 0, and ( x , w ) = 0 are feasible solutions of Problem (9): this condition is necessary and sufficient for the existence of a set W( w ) solving the innerapproximation problem.
Under Assumption 3-Outer, all y ∈ Y can be reached from the origin with feasible inputs w. Then, Problem (10) is feasible for all N ≥ nx: this condition is necessary and sufficient for the existence of a set W( w ) solving the outer-approximation problem. Moreover, if rank([CB CAB · · · CA nx−1 B D]) = ny, then Problem (10) is feasible for every nonempty Y and N ≥ nx since system (1) is then output-controllable.
In the rest of this paper, we transform Problems (9)-(10) into implementable forms. To that end, in the next section, we discuss the characterization of the RPI constraints.

IV. CHARACTERIZATION OF RPI CONSTRAINTS
In this section, we show that the lower-level problems in (9),(10) that characterize the minimal parametrized RPI set can be replaced by the equality c( x ) + d( w ) = x , i.e., the equivalence holds. For ease of notation, we denote d( w ) by d in the rest of this section, since the results are presented for a fixed w ≥ 0. Firstly, we recall from Lemma 1 that the fixed-point solution exists, and satisfies the equality c( x the equivalence in (11) holds with x = x # (d). In the following result from [19], the uniqueness of the fixed-point was shown under a slightly restrictive assumption.
. We now present a brief discussion regarding the restrictions imposed by the assumption d > 0: recalling the definition , and w > 0. While the positivity condition w > 0 can be enforced easily through a linear constraint in Problems (9)-(10), the former condition holds only if the additional assumption E i / ∈ null(B ) (or the stronger assumption rank(B) = nx) is satisfied: these assumptions restrict the class of systems and RPI set parametrizations that are often encountered. Moreover, they lead to excessively conservative RPI set parametrizations. For example, an uncontrollable system would require an RPI set parametrization that always includes the origin within its interior.
In the following result, we use continuity properties of support functions to show that the uniqueness of fixed-point holds without these additional assumptions. To that end, we introduce the perturbed disturbance vector d δ := d + δ1 defined for a scalar δ > 0.
Theorem 1: Suppose Assumption 2 holds and d ≥ 0, then the fixed-point starting from the same initial point x . We first show by induction that at all iterations k, the limit At the first fixed-point iteration we have x Hence, by definition of d δ , (12) holds for k = 1, i.e., lim We will now proceed by induction by exploiting the fact that if (12) holds at iteration index k, then At iteration index k + 1 we have = lim Since c is a continuous function of d, for all d ≥ 0, we have Hence, (13) implies Thus, by induction, (12) holds for all iteration indices k. Since Assumptions 2 holds, then according to Lemma 1, the iterates , d) and x * ( x [0] , d δ ) respectively. Therefore, (12) implies where we used the fact that x # (d δ ) is unique and independent of x [0] , as proven in Lemma 2. Finally, we note that the function x * ( x [0] , d) is well defined for all d ≥ 0 from Lemma 1, and is continuous in d ≥ 0 since it is the composition of continuous functions. This implies that the limit lim Having shown that there exists a unique fixed-point x # (d), we proceed by replacing the lower-level optimization problems in (9) and (10) by the equivalent condition c( x ) + d = x .
Remark 4: While all the results presented assume that 0 ∈ W( w ), there exist cases where it is not known a priori if the origin belongs to the disturbance set. Such cases can be accommodated in the formulation of Problems (9)-(10) by considering the disturbance set parametrization {w} ⊕ W( w ), where 0 ∈ W( w ) if w ≥ 0, andw represents the origin offset. Then, an RPI set parametrized as {x} ⊕ From basic properties of support functions, this inclusion is equivalent to . The first part of this inequality can then be eliminated by using the steady-state offset valuex = (I−A) −1 Bw. Hence, by appending the optimization variablesw andx, along with the equality constraintx = (I−A) −1 Bw, Problems (9),(10) can be modified to accommodate disturbance sets not including the origin. Note that the corresponding output-set is then {Cx + Dw} ⊕ Y( x , w ): the inclusions with respect to Y should be modeled by considering this offset. Since this modification is relatively straightforward, we skip further details because of space constraints.

V. CHARACTERIZING HAUSDORFF DISTANCE AND ENCODING INCLUSION CONSTRAINTS
In this section, we use the inclusion encoding formulation presented in [25] under the following assumption on the output-set Y = {y : Gy ≤ g} that is slightly stronger than Assumption 3.
Assumption 4: The set Y is full-dimensional, i.e., there exists somê y ∈ R ny and a scalarˆ > 0 such that {ŷ} ⊕ˆ B ny 2 ⊂ Y. For the inner-approximation problem, this assumption implies that vector g > 0. For the outer-approximation problem, this assumption implies system (1) must be output-controllable for feasibility. (9) We use the Hausdorff distance given by

A. Inner-approximation Problem
for Problem (9) In order to encode the inclusion Y ⊆ Y( x , w ) ⊕ B( ), we use the sufficient conditions presented in [25,Theorem 1], which states that the inclusion holds if there exist variables z I := {Σ I , Θ I , Π I } with Σ I ∈ R (nx+nw+ny )×ny , Θ I ∈ R nx+nw+ny and Π I ∈ R (m X +m W +m B )×m Y satisfying ( x , w , , z I ) ∈ Ξ I , where is a set of linear equality and inequality constraints, and Π I is a matrix of nonnegative elements. Since we enforce w ≥ 0, and x ≥ 0 by construction, the set Y( x , w ) is always nonempty. Then, there always exists some ≥ 0 such that the inclusion Y ⊆ Y( x , w ) ⊕ B( ) holds. Under Assumption 4, it then follows from [25, Theorem 1] that there always exist variables z I such that Ξ I is nonempty. Since this inclusion encoding is only sufficient, the computed value of Hausdorff distance is an upper-bound to the actual value.
In order to encode the inclusion Y( x , w ) ⊆ Y, we use the support functions defined as to enforce the inequality l( x ) + m( w ) ≤ g. Hence, using the definition in (17), RPI set equivalence in (11), and the proposed inclusion encodings, we write Problem (9) as B. Outer-approximation Problem (10) Similar to (17), we use the Hausdorff distance given by for Problem (10). Then, we approximately encode the inclusion Y( x , w ) ⊆ Y ⊕ B( ) using the support functions defined as The approximation results from the fact that this condition is only necessary for the inclusion to hold. Hence, the computed value of Hausdorff distance is a lower-bound to the actual value.
In order to encode the inclusion Y ⊆ ⊕ N t=0 CA t BW( w ) ⊕ DW( w ), we again use the sufficient conditions presented in [25,Theorem 1], which states that the inclusion holds if there exist variables is a set of linear equality and inequality constraints. Here, ⊗ denotes the Kronecker product. The set Ξ O is nonempty under Assumption 4 under the same reasoning as that for set Ξ I .
Hence, using the definition in (19), RPI set equivalence in (11), and the proposed inclusion encodings, we write Problem (10) as

VI. NUMERICAL OPTIMIZATION
In this section, we adopt a penalty function approach to solve Problems (18)- (20). To that end, we first note that w might be unbounded above in both these problems in case of a nonminimal representation of W( w ). We tackle this issue by introducing the support function , such that q( w ) = w if and only if W( w ) is in minimal representation. Then, we modify the objective function of Problems (18)- (20) as where σ > 0 is some scalar tuning parameter. This modification ensures that (a) the solution w is such that W( w ) is in a minimal representation; (b) the solution is not perturbed, since W(q( w )) = W( w ). This modification is not required for x , since the RPI constraint enforces uniqueness of x for a given value of w . Then, we propose to use the penalty function approach presented in [20] to solve the problems with the modified objective function. In the rest of this section, we present the approach for the innerapproximation problem. Since a very similar method follows for the outer-approximation problems, we skip further details because of space constraints. Considering Problem (18) along with objective function (21), we note that all the LPs formulating the support functions are feasible and bounded for every bounded x ≥ 0 and w ≥ 0. Hence, they satisfy strong duality [26]. This property is exploited in the penalty function algorithm to compute local optima. Introducing the optimal primal and dual variables strong duality of the LPs implies that these values satisfy the primal and dual feasibility conditions, and have a zero duality gap. For the support function c i ( x ), these conditions for a given x ≥ 0 are Introducing these variables along with the optimality conditions, Problem (18) is reformulated to a single-level problem. Within this reformulation, a penalty function approach is followed to penalize the duality gap using a constant K > 0 to obtain where z * ∈ R m X (nx+nw)+m Y (nx+nw)+m W nw denotes the optimal primal variables, λ * ∈ R m X (m X +m W )+m Y (m X +m W )+m 2 W denotes the optimal dual variables, and the penalty function penalizes the duality gap of the support function LPs. We denote Problem (22) as F( x , w , , z I , z * , λ * , K). The main idea behind the approach presented in [20] is that there exists a penalty parameter K * such that, if Problem (22) is solved with K > K * , then the duality gap P( x , w , z * , λ * ) = 0 at the optimal solution, and this solution also solves the original Problem (18) with objective function (21). Then in order to solve Problem (22), an iterative algorithm is proposed, with each iteration composed of solving two LPs.
Denoting an iteration by the subscript {l}, the first LP solved is F( x , w , , z I , z * , λ  (22) and duality gap is zero, the algorithm is terminated. Else, the procedure is repeated with K {l} ≥ K {l−1} . This algorithm was shown to converge to a local optimal solution of Problem (18) in [20]. The approach was further extended in [27] to obtain the global minimizer, and future work focuses on an efficient implementation of this method.
Remark 5: We propose to initialize the optimization algorithm using the scaling ζ ≥ 0 as w . This value can be computed using the one-step procedure in [19], and ζ can be selected by solving an LP that enforces desired inclusions with respect to the output-set Y. The dual variables corresponding to these LPs can be set as λ * {0} . We skip further details because of space constraints.
Remark 6: Alternative procedures to compute the disturbance sets can be derived from the methods presented in [28]- [30], by formulating their optimization problems with w as an optimization variable and enforcing a fixed feedback gain. While the formulation in [28], [29] involves solving LPs, the reduced-complexity polytopes can be excessively conservative. The formulation in [30] accommodates full-complexity polytopes, with the nonlinear terms in the resulting optimization problem dealt with through a Newton-type procedure. Comparison with this method is a subject of future investigation.

VII. NUMERICAL EXAMPLES
We now present two examples, with the first related to the innerapproximation problem, and the second to the outer-approximation problem. These results are obtained using the penalty function algorithm discussed in the previous section. CPLEX LP solver [31] was used to solve the LPs. The sets are plotted using plotting tools from the Multi-Parametric Toolbox [32]. For both the examples, we use the set B( ) = {y : I −I y ≤ } to define the Hausdorff distance in (17)- (19), and choose σ = 1 in (21).

A. Computation of safe reference-sets for supervisory control
We consider the system with input-constraints u ∈Û := {u : |u| ≤ 2 3 }. We assume that the system is equipped with an LQI-tracking controller such that z tracks a reference signal w: an integral-action state q with dynamics q(t + 1) = q(t) + z(t) − w(t) is appended, and the state x = [z q ] is introduced. Then, an LQI feedback  . Bottom-right plot shows resulting closed-loop input response. Observe that the input bounds are respected.
For this system, we aim to design a supervisory controller that saturates the reference signal such that input-constraints are respected: we compute the largest reference saturation limits w = [w 1w2 ] such that for all w ∈ W( w ) = {w : |w 1 | ≤w 1 , |w 2 | ≤w 2 }, we have u ∈Û. Moreover, the supervisory controller does not have access to the state x(t) of the system, which implies the set W( w ) should guarantee input-constraint satisfaction for all reachable x.
In order to compute these bounds, we note that if w ∈ W( w ), the state of the closed-loop system always belongs to the mRPI set as x ∈ Xm( w ) (provided x(0) ∈ Xm( w )). Then, the condition u ∈Û in equivalent to the inclusion KXm( w ) ⊆Û. Finally, we assume that the references are always bounded as w ∈Ŵ := {w : |w| ≤ 5 5 }. Hence, we compute the vector w such that the inclusions KXm( w ) ⊆Û and W( w ) ⊆Ŵ hold. We write , and the output-set Y =Û ×Ŵ, based on which we solve the inner-approximation problem (18): We approximate the mRPI set using the RPI set where the matrix E is composed of hyperplanes defining the set ⊕ 5 t=0 A t BŴ (A,B denote the matrices of the closed-loop system). This choice results in m X = 240. The result of solving this problem using the methods presented in this paper is shown in Figure 1. The computed saturation bounds arew 1 = 1.6172,w 2 = 4.0125.
We also plot the set KX γ,µ , where X γ,µ is a tight RPI approximation of the mRPI set Xm( w ) od presented in [33]: This method considers the system x(t + 1) = Ax(t) +w(t) with w ∈ BW( w ) ⊕ γB 4 ∞ . Labeling the mRPI set for this system as X γ m ( w ), the tightly approximating RPI set satisfies X γ,µ ⊆ X γ m ( w ) ⊕ µB 4 ∞ . We observe that KX( x ) characterizes a fairly tight approximation of the set KXm( w ). Closed-loop trajectories are plotted with references w sampled from the vertices of W( w ), for which the input response satisfies the input-constraints. Hence, if x(0) ∈ X( x ), the supervisory controller can command any reference w ∈ W( w ) with guaranteed input-constraint satisfaction.
Remark 7: The mRPI set is suitable to formulate the problem in Example A since we do not have access to the state x(t). If Parameters m X = 30, m W = 6, γ = 10 −5 , µ = 10 −6 Fig. 2: Results of solving the outer-approximation problem. Input, state and output trajectories are plotted with w(0), x(0), y(0) denoted by black * , w(100), x(100), y(100) denoted by red * . Observe that the vertices of Y are reachable from x(0) = 0 with w ∈ W( w ).
this limitation is overcome, then a reference governor scheme [8] is more suitable to design the supervisory controller, which uses control invariant sets to guarantee constraint satisfaction.

B. Computation of input-constraint sets for output reachability
We consider system (1) with initial-state x(0) = 0, for which we compute the smallest input-constraint set W( w ) = {w : F w ≤ w } with rows F i = [sin(2π(i − 1)/m W ) cos(2π(i − 1)/m W )] for each i ∈ I m W 1 , such that all y ∈ Y can be reached with control inputs w ∈ W( w ). To that end, we use f ( w ) as a measure of the set W( w ), and formulate the optimization problem P N defined as w,N := arg min s.t. y = Σ N −1 t=0 CA t Bwy(t) + Dwy(N ), wy(t) ∈ W( w ), ∀ y ∈ Y, ∀ t ∈ I N 0 , such that W( w,N ) is the smallest input-constraint set in which there exist inputs driving the output of system (1) to all y ∈ Y from the origin in N -steps. If Assumption 3-Outer holds, then P N is feasible for all N ≥ nx. It can then be shown that the sequence of optimal values {f ( w,N )} N is non-increasing, and converges to the optimal value of the problem w * := arg min where Xm( w ) is the mRPI set corresponding to W( w ). This follows from the idea that the mRPI set is the closure of the largest 0reachable set [1]. Hence, computing the smallest input-constraint set entails solving Problem (23). We choose f ( w ) = d H (Ym( w ), Y), such that Problem (23) is equivalent to Problem (3) along with the output-set inclusion constraint. This choice ensures that we compute an input-constraint set W( w ) whose 0-reachable set in the output space tightly includes the target output-set Y. We approximately solve Problem (23) based on the outer-approximation formulation in Problem (10): we approximate the mRPI set using the polytopic RPI set X( x ) = {x : Ex ≤ x } with rows E i = [sin(2π(i − 1)/m X ) cos(2π(i − 1)/m X )] for each i ∈ I m X 1 . Using this set, we formulate Problem (20). The results of solving this problem using the methods presented in this paper are shown in Figure 2 (20). We see that the computed set W( w ) is such all y ∈ Y are reachable from the origin. We also plot tight approximation RPI set X γ,µ of the mRPI set Xm( w ) using the methods presented in [33], in a manner similar to the previous example. We observe through the set Y γ,µ := CX γ,µ ⊕ DW( w ) that Y ⊆ Y γ,µ ⊆ Y( x , w ) holds, thus ensuring the desired reachability.
In conclusion, one can design feedback controllers to select inputs w from the input-constraint set W( w ), with the guarantee that for any x(0) ∈ Xm( w ), there always exist feasible inputs to reach every target output y ∈ Ym( w ) ⊃ Y. In Figure 2, we also plot state, input and output trajectories with x(0) = 0 and y(100) ∈ Y to demonstrate the reachability.

VIII. CONCLUSIONS
We have presented a method for computing an input disturbance set for discrete-time linear time-invariant systems such that the reachable set of outputs approximates an assigned set. To that end, we formulated an optimization problem in order to minimize the approximation error. Finally, we presented some numerical results to demonstrate the feasibility of the approach and two possible practical applications. Future research will further develop the solution algorithm by considering: (a) alternative solution methods such as, e.g., value function approaches [26]; (b) optimizing also over matrices E and F . Finally, the potential of this technique when applied to feedback controller synthesis and to system identification problems will be investigated.