Automatic scenario generation for efficient solution of robust optimal control problems

Existing methods for nonlinear robust control often use scenario-based approaches to formulate the control problem as large nonlinear optimization problems. The optimization problems are challenging to solve due to their size, especially if the control problems include time-varying uncertainty. This paper draws from local reduction methods used in semi-infinite optimization to solve robust optimal control problems with parametric and time-varying uncertainty. By iteratively adding interim worst-case scenarios to the problem, methods based on local reduction provide a way to manage the total number of scenarios. We show that the local reduction method for optimal control problems consists of solving a series of simplified optimal control problems to find worst-case constraint violations. In particular, we present examples where local reduction methods find worst-case scenarios that are not on the boundary of the uncertainty set. We also provide bounds on the error if local solvers are used. The proposed approach is illustrated with two case studies with parametric and additive time-varying uncertainty. In the first case study, the number of scenarios obtained from local reduction is 101, smaller than in the case when all $2^{14+3\times192}$ extreme scenarios are considered. In the second case study, the number of scenarios obtained from the local reduction is two compared to 512 extreme scenarios. Our approach was able to satisfy the constraints both for parametric uncertainty and time-varying disturbances, whereas approaches from literature either violated the constraints or became computationally expensive.


Introduction
Robust optimal control problems are often solved using a scenario-based approach, where each scenario corresponds to a separate realization of uncertainty.Increasing the number of scenarios improves robustness, while increasing the size of the optimization problems.Mitigating the size of the problem by reducing the number of scenarios requires knowledge about how the uncertainty affects the system which is a challenge especially if uncertainty varies with time.This paper draws from approaches used in semi-infinite optimisation to solve robust optimal control problems by selecting scenarios in an efficient way.

Background
To ensure that the optimisation problems resulting from scenario-based approaches to robust control are tractable, the number of scenarios must be limited [5].Usually, the choice of scenarios is done from experience [11,12] and requires knowledge about both the controlled system and the uncertainty to ensure that the chosen scenarios guarantee robustness.A recent review of scenario-based methods indicated that scenario selection is highly affected by the knowledge about the uncertainty distribution [6].In practice, to tackle problems with limited knowledge about the uncertainty, it is often assumed that the worst-case scenarios lie on the boundary of the uncertainty set [21,24,36,38].The worst-case scenario in nonlinear systems may lie in the interior of the uncertainty range [19,26].In this paper we show that the worst-case scenarios may be inside the uncertainty set, even for a linear dynamic system and present a method for choosing potential worst-case scenarios assuming limited knowledge about uncertainty.
Systematic approaches to choosing scenarios for time-varying uncertainty are usually based on creating scenario trees [14,36].In these approaches, a large set of scenarios is chosen at the beginning of the time horizon.The number of elements in the set of scenarios increases combinatorially [31].However, there is no guarantee that the chosen set of scenarios includes the actual worst-case scenarios [26].To overcome this drawback, we propose to use a method derived from semi-infinite optimisation to iteratively add new scenarios to the current set, to provide more flexibility in adjusting the set of scenarios and finding worst-case scenarios.
Previous studies provided an in-depth review of semi-infinite optimization methods [9,[15][16][17]32].In particular, it has been indicated [17] that local reduction methods allow overcoming the dependence on the initial choice of scenarios.In these methods, the set of scenarios is created iteratively by alternating between solving an optimisation problem with the current set of scenarios and solving interim optimisation problems to find the maximal violation of constraints and extend the set of scenarios [4].Therefore local reduction methods enable adding scenarios that may not have been considered at the beginning, such as scenarios from the interior of the uncertainty set.
The iterative approach based on alternating between optimisation problems has been already used for robust control, which suggests potential usefulness of local reduction.The D-K iteration found in µ-synthesis consists in alternating between the synthesis of an H ∞ controller and minimisation of the singular value for the corresponding controller [33].However, the D-K iteration method depends on the optimality of each step, which indicates that the local reduction method will also be affected by how the interim optimisation problems are solved.In particular, due to the need to solve to global optimality the interim problems [9], the complexity of local reduction methods prevented their use in robust optimal control problems.An approach based on simulated annealing has been proposed to facilitate finding a global solution [25].We show that the local reduction method provides good results even if local solvers are used and provide bounds on the solution.
Existing applications of semi-infinite optimisation methods in control systems with uncertainty are limited.Semiinfinite optimization methods have been used for optimal control [13] to find optimal trajectories for robotic arms.However, the authors considered only exogenous uncertainty due to obstacles that did not affect the dynamics of the controlled systems.Parametric uncertainty for linear systems was considered [18,41] to apply semi-infinite optimisation methods for model identification.However, these works did not consider time-varying uncertainty.Similarly, semi-infinite optimisation methods has been used to solve an optimal control problem with only parametric uncertainty [26].Time-varying uncertainty was considered in works where local reduction was applied to find the interim worst-case scenarios [36].However, the authors assumed that at every time step the number of possible scenarios was finite.

Contributions
The main contributions of the current work are: • Formulation of a local reduction method as a way of automatically generating scenarios in robust nonlinear control problems if time-varying and constant uncertainties are present, • Formulation of error bounds on constraint violation if local optimization solvers are used in local reduction, • Numerical demonstration of local reduction methods for solving nonlinear robust optimal control problems with parametric uncertainty and linear robust optimal control problems with both parametric and time-varying uncertainty.
A preliminary version of the application of local reduction for scenario generation is available [40].The rest of the paper is structured as follows.Section 2 introduces robust optimal control problems.Section 3 presents the new method for solving robust optimal control problems.The numerical results are shown in Section 4 where we compared the results obtained with the scenarios from local reduction to three commonly used approaches: nominal case with no uncertainty, random case with randomly drawn realisations of uncertainty, boundary case with only extreme values.The paper ends with conclusions in Section 5.

Problem formulation 2.1 Semi-infinite optimization problem
A semi-infinite optimization problem is formulated as: where A ⊂ R n θ and B ⊂ R nρ are nonempty and compact sets, and Q and R are continuous functions of their respective arguments [4].The problem (1) has a finite number of variables θ but includes an infinite number of constraints if B has an infinite number of points.In particular, B may be uncountable.One approach to remove the infinite number of constraints consists in rewriting the constraint (1b) as: The challenge in solving the equivalent problem with constraint ( 2) is in non-differentiability of the function S(•).
The local reduction method proposed by [4] allows overcoming the non-differentiability of S(•) by sequentially solving (1) with finite subsets of constraints taken from B.
The main challenge in formulating robust optimal control problems as semi-infinite optimisation lies in inclusion of system dynamics in the form of equality constraints.In this paper, we show that optimal control problems can be formulated as semi-infinite optimization problems and solved using the local reduction method [4].

Dynamic system with uncertainty
The system to be controlled is described by a nonlinear difference equation with time-varying uncertainty where f k is continuously differentiable.The state x 0 at time zero is w.l.o.g.assumed to be equal to a given x.The control trajectory u := (u 0 , . . ., u N −1 ) is generated by a causal dynamic feedback policy u k := π k (x 0 , . . ., x k ; q 0 , . . ., q k , r) that is parameterised by q := (q 0 , q 1 , . . ., q N −1 ) ∈ R nq and r ∈ R nr .The state trajectory x := (x 0 , . . ., x N ).The time-varying uncertainty w k at time k and the constant uncertainty d affect the dynamics in both an additive and non-additive way, and take on values from compact and uncountable (infinite cardinality) sets.Uncertainty in the measured value of x k can be modelled by a suitably-defined choice of f k , π k and w k .In this work, we assume that the structure of the dynamic feedback policy, and hence the parameterisation, is already defined.A trajectory (x, u) satisfying the dynamics (3) and control policy for a given parameterization (q, r) and realisation of uncertainty (w, d), where the trajectory w := (w 0 , . . ., w

Objective function and constraints
The cost function for the optimal control problem over a horizon of length N is: ( Both the terminal cost function J f (•, •, •) and stage cost ℓ k (•, •, •, •) are continuously differentiable and depend on the uncertainty w and d.The objective of the optimal control problem is to find a feedback policy π for system (3) such that the worst-case cost in ( 5) is minimized and the constraints are satisfied for all time instants k = 0, . . ., N −1, all states x, control u, uncertainty w and d.The vector function of n g components, g k (•, •, •, •), is continuously differentiable and depends on uncertainty w and d.Note that a constraint on x N can be included by incorporating f N −1 in a suitable definition of g N −1 .
To ensure that the optimal control problem with the objective (5) and the constraints ( 6) is well-defined over the horizon N , we introduce the following assumption on the trajectories (4): Assumption 2.1.The trajectories (4) are bounded over a finite horizon N : We note that the boundedness of (4) need not imply stability of the dynamics (3).

Semi-infinite formulation
Given a set of uncertainties H ⊆ W N × D, the problem in this work is stated as: where J := {1, . . ., card H} and (x i , u i ) is the state and input trajectory associated with the i th disturbance realisation (w i , d i ) such that ) is linear jointly in all arguments, the problem (8) can often be solved using scenario-based methods for robust control [5,31], provided additional convexity assumptions are satisfied by the uncertainty set W. In this work, the dynamics from (3) are nonlinear and W is only non-empty and compact.Moreover, the set H is assumed uncountable, which is a common case if the disturbances belong to a polytope.
Proof.In contrast to (1a), the objective function in (8a) contains uncertainty.Introducing γ ∈ R, we rewrite (8) as: The problem (9) has uncertainty exclusively in the constraints.If card H is finite, then the problem ( 9) is convenient to solve numerically using tailored efficient finite-dimensional optimisation methods that exploit the sparsity in the relevant Jacobians and Hessians.However, infinite cardinality of H yields an infinite number of both constraints and variables, which means that the problem (9) needs to be further reformulated to become (1).Noticing that the constraint (9b) is equivalent to max we introduce In (11), e h is the h th column of an identity matrix I ng .Using (4) and (11), we can write (9) as: The problem ( 12) is equivalent to Taking θ := (q, r, γ) and ρ := (w, d) in (12) (similarly in (13)) we obtain the form of (1) (similarly (2)).
Theorem 2.1 makes no assumptions on the cardinality of the set H, which is uncountable in general.As a result, problem (12) has an infinite number of constraints, in general, in a similar way to (1).Using the fact that ( 12) is equivalent to (1), and ( 8) is equivalent to (12), we note that solving (12) using methods developed for semi-infinite optimization of the form (1) is equivalent to solving the optimal control problem (8).Thus, Theorem 2.1 allows one to solve the optimal control problem (8) as a semi-infinite optimization problem of the form (1) using local reduction [4].

Robust solution
We now introduce definitions of robust solutions that we are going to use in the remainder of the paper.First, we notice that the left-hand side of (13b) is equivalent to: Using ( 14) allows us to introduce the necessary definitions.
Definition 2.1 (Scenario).A scenario is a realisation of the uncertainty (w, d) ∈ W N × D.
Definition 2.4 (Solution robust to s scenarios).A triple (q * , r * , γ * ) is called a robust solution to the s scenarios in H if G max (q * , r * , γ * , w, d) ≤ 0 for all (w, d) ∈ H and card H = s.

Local reduction for optimal control
We will now extend the local reduction methods [4] to robust nonlinear optimal control.

Algorithm
The local reduction method [4] consists in iteratively solving finite-dimensional optimization problems.We use the local reduction methods for the problem ( 9) or ( 12) by iteratively solving optimal control problems parametrised by scenarios.The iterations alternate between solving minimization and maximization steps that will now be described.

Minimization step
The local reduction algorithm for robust optimal control is shown in Algorithm 1.The proposed algorithm in iteration j solves an optimal control problem of the form (9) or (12) assuming that the number of scenarios card H j at step j is finite.The algorithm needs an initial guess for the parameters of the controller.For instance, the initial guess can be obtained by solving (12) for one scenario, i.e. card H 1 = 1.Alternatively, the initial guess can be found by solving (12) for a small number of scenarios, obtained for example from a coarse discretization of the uncertainty set [30].
12 until card H j = card H j−1 ; to the current set of scenarios H j .If there exists at least one violated constraint, then a scenario corresponding to the maximum constraint violation is added to the scenario set H j+1 in the next iteration (line 7).The new set H j+1 is then used to find a new set of control parameters (line 9).The algorithm ends if no new scenarios are added, i.e. card H j = card H j−1 .
In this work, any scenario corresponding to the maximum constraint violation can be added to the set of scenarios.However, it has been shown that computational performance may be improved if multiple scenarios are added [37].

Maximization step
The maximization step consists in solving (14) with H = W N ×D.Solving ( 14) is equivalent to solving n g •(N −1)+1 optimization problems, where n g denotes the number of elements in the vector function g(•) from constraints in (6).The algorithm is presented in Algorithm 2. Without loss of generality, we assume that the first constraint to include in the maximization problem corresponds to the reformulated objective function (5).A scenario that corresponds to maximal value of this constraint is added to an auxiliary set K. The remaining n g • (N − 1) constraints are included as objectives in the respective maximization problems (line four to eight in Algorithm 2).Note that the problem corresponding to the objective (line 2) and all the problems corresponding to the constraints (line four to eight) can be solved in parallel.
All maximization problems are subject to the same equality constraints capturing the dynamics.This formulation allows us to treat the maximization problems as optimal control problems and preserve the sparsity of the relevant Jacobians and Hessians.We solve the maximization problems as optimal control problems where q, r, and γ are known parameters whereas w and d are treated as unknown inputs.Thus, the maximization problems can be solved using any off-the-shelf solver for optimal control problems.
Solving ( 14) with H = W N × D corresponds to lines four to eight in Algorithm 2 and can be done by solving a number of finite-dimensional optimization problems in parallel [41].

Convergence of Algorithm 1
The convergence of local reduction method in the case of the form (1) was shown in several previous works [4,23,27].The authors required that the sets A and B in (1) are non-empty and compact, and that the functions Q and R are continuous with respect to all their arguments.They showed that the sequence of solutions obtained for a sequence of finite and countable subsets of B converges to the solution of (1).A discussion on convergence rate of local reduction methods also follows [32].We show in Theorem 3.1 when the Algorithm 1 solves problem (8).Input: Current values of q j , r j , γ j Output: Worst case scenario (w j , d j ) in iteration j 1 Find any x * , u * , w * , d * that solves: The solution (q * , r * , γ * ) obtained from Algorithm 1 for a non-empty and compact set W N × D converges to the solution of (8 Proof.From Theorem 2.1, we have θ := (q, r, γ) and ρ := (w, d), A := F, B := H.In (12), we take Q(θ) := γ which is linear and thus continuous.Then we have R(θ, ρ) := G(q, r, γ, w, d) which is continuous because both max h,k e T h g k (•, •, •, •) and J N (•, •, •, •) are continuous.The proof follows Lemma 2.2 from a previous study [23].
Theorem 3.1 requires the constraints in G to be defined over a compact set F × W N × D. As a direct consequence, we obtain Remark 1.

Remark 1 (Boundedness of constraints
Proof.The proof follows directly from the extreme value theorem because Boundedness of G ensures that the maximization step presented in Algorithm 2 is well-posed.Section 3.5.1 will further demonstrate the impact of boundedness on the solution obtained from Algorithm 1. We also note that similarly to Assumption 2.1, the requirement of boundedness of G need not imply stability of the dynamics (3).In particular, the method can be used for solving finite horizon optimal control problems with unstable linear dynamics affected by uncertainty, as will be demonstrated in Section 4.

Constraint dropping
The method presented in Algorithm 1 assumes that the cardinality of the sets H j is increasing with j, i.e.H j ⊂ H j+1 for all j.The increasing cardinality corresponds to an increase in the size of the optimisation problem in line nine in Algorithm 1.The authors [4] provide additional convexity conditions allowing one to drop elements from the set H j .In particular, they require (1a) to be strictly convex with respect to θ and R(θ) to be convex with respect to θ for ρ ∈ B. Following Theorem 2.1, the conditions provided in previous works [4] correspond to strict convexity of (12a), and convexity of G(z(•, •, w, d), w, d, •) for any (w, d) ∈ W N × D. In the current work, we do not assume convexity of ( 11) and (12a) is only convex, not strictly convex, so dropping constraints from the set H j does not guarantee convergence of the local reduction algorithm.To enable constraint dropping, we consider a special case of problem (8) where the cost is independent of the uncertainty: Then we can adjust Algorithm 1 to Problem (20) to enable dropping constraints.First let us rewrite (11) as: and ( 14) as: where (q, r) ∈ F new ⊂ R nq × R nr are found in a subset of the whole search space.
and F new is convex, then (16) can be replaced by: where and the solution of the modified algorithm will converge to the solution of (20).

Inexact local reduction
To simplify the exact local reduction from Algorithm 1 we propose an inexact formulation of the algorithm focusing on numerical properties of the solvers used for the optimal control problems.The exact algorithm for local reduction presented in Algorithm 1 assumes that the maximisation step finds the global solutions to the maximisation problem and only one scenario obtained in this step is then added to the scenario set.These assumptions are often difficult to satisfy.In practice, there are two possible cases: • The maximisation step in a given iteration has multiple solutions in general, but only a limited number is used, • The maximisation step is solved approximately.
In particular, we will focus on analysing the case when local solvers are used.
If a global solver is used, but only a limited number of scenarios is added, the local reduction algorithm needs more iterations to find a solution than in the case of adding all the scenarios [24].Therefore there exists a trade-off between the speed of convergence of the local reduction method and the size of the problem solved in the minimisation step.We show the impact of approximate solutions by considering similarity of scenarios, i.e. when the interim worst-case scenarios from line 3 in Algorithm 1 are considered similar.Definition 3.1 (Similar scenarios).Let (w 1 , d 1 ) and (w 2 , d 2 ) be two scenarios and let ϵ w ≥ 0 and ϵ d ≥ 0 be fixed parameters.The two scenarios are similar if and Using Definition 3.1, we modify line 7 in Algorithm 1 so that the scenario (w j , d j ) in iteration j is added to the current set of scenarios if it is not similar to any of the scenarios in H j .Algorithm 3 summarizes the inexact local reduction method with the evaluation of when the scenarios are similar.
18 until card H j = card H j−1 ;

Impact of similarity of scenarios
The algorithm provides a solution that is robust to s scenarios in the sense of Definition 2.4, where s = card H * .From Remark 1 and (11), we already have that for any two scenarios (w 1 , d 1 ), (w 2 , d 2 ), the constraints g k are bounded and hence: A tighter bound can be obtained if we use Assumption 2.1 and note that for a chosen parameterization q * , r * , γ * the trajectories z are BIBO-stable w.r.t.disturbances, i.e.
Theorem 3.3 uses (31) to show the impact of the similarity of scenarios on the constraint satisfaction.
Theorem 3.3.Let (w 1 , d 1 ), (w 2 , d 2 ) be two identical scenarios with ϵ w = 1 N ϵ * w and ϵ d = ϵ * d , and ∥(w i , d i )∥ 2 ≤ ς w for i = 1, 2. Let q * , r * , γ * be the solution of (8) obtained for (w 1 , d 1 ) and let z 1 = z(q * , r * , w 1 , d 1 ), z 2 = z(q * , r * , w 2 , d 2 ).Then the constraint violation for (w 2 , d 2 ) is bounded and where z i k is the trajectory z i at time k, i = 1, 2, and L is a local Lipschitz constant for a given k.Proof.We have and We can write (33) as: where e k is the k th row of an identity matrix I N .Without loss of generality, we can assume that w 1 k − w 2 k ̸ = 0 for k = 1, . . ., N .Then from orthogonality of the set {e k (w 1 k −w 2 k )} k=1,...,N , and using (33), we get from the Pythagorean theorem: Taking into account that ∥e k (w (36), we get: for all k = 1, . . ., N .Summing up (34) and (37) gives: Then we have: where 0 denotes the origin of R n d and we used the Pythagorean theorem in R n d +1 .Then we have: where 0 denotes the origin of R nx+nu .At the same time, from the assumption that g k is continuously differentiable with respect to all its arguments, we get the local Lipschitz condition [39, Ch. 2]: where L is the local Lipschitz constant of g k .Taking a square in (41) and using (40), we obtain: which concludes the proof.
In a similar way the proof can be done for the constraint from (9d) because the cost J N is also continuously differentiable.
Theorem 3.3 shows that the satisfaction of constraints depends on: • The local Lipschitz constant of the constraints, • The response of the system to disturbances, • The choice of similarity of scenarios.
The local Lipschitz constant and the response to disturbances are inherent to the system.We note from Theorem 3.3 that the assumption about having a global solution to the maximisation problem in Algorithm 1 is crucial to ensure no constraint violation.Even if the maximisation step is solved exactly and ϵ w = ϵ d = 0, the constraint violation in (32) is defined by 4L 2 ς 2 z .At the same time, the choice of parameters in Definition 3.1 affects the constraints violations.
We also note that (31) makes no assumptions on ς z .For unstable systems, ς z may be large, thus making the bound in Theorem 3.3 uninformative.However, the value of ς z depends also on the chosen control policy π k which can be used to modify the control invariant set [3] and thus tighten the bound.
Furthermore, the impact of the similarity of scenarios provides information about the solution if local, instead of global, optimization solvers are used to solve the maximisation problems.To make this precise, let (w 1 , d 1 ) in Theorem 3.3 be a local solution from Algorithm 2, and let (w 2 , d 2 ) be a global solution.Taking ϵ * w and ϵ * z such that and ∥(w i , d i )∥ 2 ≤ ς w for i = 1, 2, the maximal constraint violation is bounded by (32).
As a direct consequence of Theorem 3.3 we obtain the following result: where we assumed that . Thus, we have: concludes the proof.

Impact of number of scenarios
The bounds obtained in Theorem 3.4 allow inference about constraint violation in a practical implementation of Algorithm 3. In particular, we will now analyse how fixing the number of scenarios H affects the constraint violation.
For simplicity, we focus on the constraint from (9b) but an analogous reasoning can be done for (9d).Let us assume that the function g k from (9b) is such that: where Q k is Lebesgue-measurable.The Lebesgue measure of the set will be called its n-dimensional volume [29, Ch. 21] and denoted Vol Q k = Ω k .We assume that the values of g k are from a uniform distribution over Q ℸ for any realization of uncertainty.Let us also introduce for every scenario (w i , d i ) ∈ H * the following subsets: where δ ≥ 0 is a constant such that S i k,δ ⊆ Q k .Then we can introduce such that S i k,δ ̸ = S j k,δ for i, j = 1, . . ., H q , i ̸ = j.The set from (47) collects H q sets S i k,δ corresponding to distinct values of constraints g k (z i k , w i k , d i ) obtained from scenarios (w i k , d i ) and trajectories z i k .We have that H δ ≤ H, with strict inequality if at least two different realizations of scenarios (w i k , d i ), (w j k , d j ) and the corresponding trajectories z i k , z j k , i ̸ = j, i, j = 1, . . ., H, lead to the same value of the constraint, i.e.
From countable subadditivity of Lebesgue measure [28], we get Vol , and then: Noticing that S i k,δ is a ball in R ng with radius δ, we have [34, p. 135]: Figure 1: Trade-off between the number of scenarios H δ and the desired threshold δ as a function of dimensionality n g of the uncertainty space where Γ is the gamma function, and then: We can now state: Theorem 3.5.Let us assume that Algorithm 3 is used to solve (8) where g k satisfies (45) for all k = 0, . . ., N − 1 with Vol Q k = Ω k , and the values of g k are uniformly distributed over Q k .Let us assume that Algorithm 3 finished with scenario set H * , card H * = H.Let us assume that a threshold δ ∈ R has been chosen so that the sets S i k,δ from (46) satisfy S i k,δ ⊂ Q k for all k = 0, . . ., N − 1, i = 1, . . ., H q .The probability that the constraint violation exceeds the threshold δ for scenario (w Proof.Dividing by Ω k > 0 all the terms in (51) we obtain k,δ into (53) and rearranging, we obtain (52).Theorem 3.5 combines the number of scenarios H δ obtained from local reduction with the threshold for constraint violation δ.Assuming that g has been normalised so that Ω k = 1, we rewrite (52) as: where ) is a constant depending on the number of constraints n g .We see in (54) that there is a trade-off between the threshold δ and the number of scenarios H δ .For a given threshold δ, increasing the number of scenarios will decrease the probability of violation.Conversely, if the number of scenarios is constant, a small threshold may be exceeded with high probability.An example of the trade-off is shown in Fig. 1 as a function of the dimensions n g .
Finally, Theorem 3.5 has a practical interpretation.Assuming that the allowed probability of constraint violation is p des , from (54) we obtain a bound on the number of scenarios H δ to satisfy this requirement for a given δ: The number of scenarios obtained from (55) can be used as an additional stopping criterion (line 18 in Algorithm 3).Finally, should the probability distributions of the uncertainty be available, our approach can be extended to explicitly include chance constraints and probability distributions a priori in the problem formulation by using an appropriate problem structure [1].

Examples
We will first show that the realisation of uncertainty leading to the maximal constraint violation can be anywhere, not necessarily on the boundary of the uncertainty set even for a linear system.We also show the performance of the local reduction applied to an unstable system.
We then show that the local reduction method described in Section 3 finds scenarios from inside the uncertainty sets and provides robust solutions to optimal control problems with uncertainty in two numerical examples: temperature control in a residential building and flow control in a centrifugal compressor.The solution provided by local reduction is then compared with the case obtained for boundary scenarios and for scenarios chosen randomly from a uniform distribution.

Scenarios not on the boundary
An example of scenarios not on the boundary for a nonlinear system was provided in previous works [19,26].We show that a linear system with parametric uncertainty may have interim worst-case scenario in the interior of the uncertainty range.The worst-case scenario in the sense of Definition 2.2 for a robust optimal control problem of the form (12) may be in the interior of the set W × D.
Let us assume that we have a system with dynamics affected by parametric uncertainty d: where k = 1, . . ., N and d ∈ [d, d], A, B are constant scalar matrices, x 0 is known.Let us assume further that the optimal control problem includes a constraint of the form: The constraint (57) must be satisfied for all k.The maximisation step in the local reduction method consists in solving a series of optimisation problems with the objective for every k: This example confirms that considering boundary scenarios would miss the actual worst-case scenario.In a similar way, selecting a priori a number of scenarios would result in adding unnecessary scenarios that may or may not be the worst-case scenario.

Unstable system
To show that the boundedness of constraints required in Theorem 3.1 need not imply stability, we analysed a system with dynamics: where a = 2.1, the uncertainty d ∈ [0.9, 1.1], and with x 1 = 0.5.The constraints to be satisfied were 0 ≤ x k ≤ 1, for all k = 1, . . ., N with N = 10.The controller u k was parameterised as an affine function of the state: and the constraints u k ∈ [−1, 1] were enforced by a smooth saturation function where β i are constants.Here β 0 = −2.0229,β 1 = 1, β 2 = 1.2963, β 3 = 1.01145.The objective was to minimise the square of the control over the whole horizon N : The trajectory obtained from maximization of the violation of the constraint x 10 ≤ 1 is shown in Fig. 2a.The time horizon is finite, N = 10, so the trajectory x is bounded for any value of d and thus the local reduction can be used.Using the local reduction algorithm resulted in three scenarios: d 1 = 1, d 2 = 0.9, d 3 = 1.1 that ensure robustness, as indicated in Fig. 2b.The trajectories in Fig. 2b were obtained for 500 randomly chosen scenarios uniformly distributed in [0.9, 1.1].
The example shows that stability of the dynamics is optional, provided the constraints are bounded.Nevertheless, it is recommended to take instability into account when numerically solving the miximization problems in Algorithm 2. To show the influence of instability on the numerical performance, we computed the violation for six different values of the horizon N and the results are collected in Table 1.For instance, the maximization problems can be solved approximately, terminating as soon as a violation has been found.Such approaches correspond to choosing large values of ϵ w and ϵ d in Algorithm 3, so the analysis in Section 3.5.1 holds.

Linear system with parametric uncertainty
This numerical example consists in a linear system with both parametric and additive time-varying uncertainty.The example describes a single zone building affected by time-varying internal heat gain, solar radiation, and external temperature [20].The objective is to follow a time-varying set-point for internal temperature x temp k .The dynamics are discrete and linear: The states x describe the indoor temperature x temp , wall temperature x wall , and the corridor temperature x corr .The control u represents the amount of heating and cooling delivered to the building.The initial condition was chosen as x 0 = 25.0 24.0 24.0 T • C.Moreover, we assume that the wall temperature and the corridor temperature can only be measured approximately, so there are two additional sources of uncertainty in the initial condition for these two states.We assume x i 0 = 24 + d i , i = wall, corr, where d i ∈ [−0.5, 0.5].We also assume that the matrices A = [a i,j ] and B = [b j ], i, j = 1, 2, 3 are affected by uncertainty: where δ i,j , η j are uncertain parameters.Two cases will be considered: Case A with δ i,j , η j ∈ [0.98, 1.02] and Case B with δ i,j , η j ∈ [0.96, 1.03].The minimal control effort is ensured by the objective function: It is assumed that the day starts at 6.00 am and lasts 12 hours.The temperature indoors must stay within limits: During the day, the indoor temperature must be kept above 23 • C and during the night can drop down to 17 • C: The maximal temperature is the same during the day and night, The optimal control problem is solved over a period of 48 hours starting at 6.00 am the first day, with N = 192.As a result, the trajectory constraints (67) impose 192 • 2 constraints corresponding to every sampling time.The three uncertain parameters, internal gain, solar radiation, and external temperature, vary with time within the limits provided in Table 2.
The control variables are parameterised as: where K and q k are decision variables.Furthermore, we include saturation of the control inputs: The saturation was approximated by a smooth function: where β i are constants.Here β 0 = −5030, β 1 = 2.937, β 2 = 0.003, β 3 = 1207.In total, there are 14 uncertain parameters affecting the matrices A, B, and the initial condition for the wall and corridor temperatures.We assume no knowledge about the scenarios, except the ranges of uncertainty.

Results -Case A
In case A, we ran Algorithm 1 with the time-varying uncertainties from Table 2 and parametric uncertainties δ i,j , η j ∈ [0.98, 1.02].

Overall performance
The local reduction method in Case A reduced the number of scenarios to five.The resulting controller obtained for the interim worst-case scenarios was then validated for 500 random realisations from a uniform distribution of uncertainty.The validation of the controller is shown in top left plot in Fig. 3.The black curves stay within the green bounds corresponding to constraints (67).The results suggest that local reduction was able to find a robust solution despite using a local solver for maximisations.
The results also indicate that the local reduction method handles time-varying uncertainty without specifying the scenarios over the whole time horizon.This is because there is no need to specify time-varying scenarios as they Figure 4: An example of a time-varying scenario found by the local reduction in Case A will be found in the maximisation step in Algorithm 2.Moreover, Algorithm 2 treats time-varying uncertainty as one realisation over the whole horizon, thus overcoming the limitations of separate robust horizon [36].
An example of an interim worst-case scenarios obtained in the maximisation step is shown in Fig. 4.

Comparison with other approaches
The results obtained from local reduction are then compared with three scenario-based approaches from the literature [6]: • Nominal approach, with a controller obtained assuming there is no uncertainty, i.e. (w, d) = 0 (further denoted as "Nominal") • Randomised approach, with a controller obtained for a number of randomly chosen scenarios (further denoted as "Random") • Extreme approach, with a controller obtained for three scenarios: nominal, lower bound, and upper bound for all uncertainties [36] (further denoted as "Nominal+two extreme") Validating the nominal controller with 500 random scenarios shows that the approach based on nominal values leads to violation of constraints as shown in the plot 'Nominal' in the left column of Fig. 3.
The second set of controllers we used was derived using three sets of random scenarios: five scenarios because five scenarios were found in local reduction, 100 scenarios, and 250 scenarios.The results are shown in the plot 'Random' in the left column of Fig. 3, with black corresponding to the controller obtained from five scenarios, yellow to the controller with 100 scenarios, and blue to the controller with 250 scenarios.In all the cases the controller violated at least one of the bounds (100 scenarios gave 0.2 • C, 250 scenarios gave 0.1 • C), with the controller based on five scenarios violating both the lower and upper bound (1.1 • C).Even though the violation decreased with increasing the number of scenarios, further increasing the number of random scenarios to 600 proved unsuccessful in avoiding the violation.Larger problems could not be solved on the computer.
A possible reason for the random controller being unable to satisfy the constraints is due to not including extreme scenarios in the scenario set.If we were to take only extreme values for every uncertainty and consider all extreme scenarios, we would need to solve a problem with 2 14+3×192 scenarios, which is intractable.To reduce the number of scenarios, we chose to use the nominal scenario, combined with two extreme scenarios.The extreme scenarios were taken as all uncertainties on their lower or upper bound simultaneously.The results of validating the controller for 500 scenarios are shown in the plot 'Nominal+two extreme' in the left column of Fig. 3.The controller based on the extreme scenarios was also able to avoid constraint violations with three scenarios.
The results of the comparison show that the local reduction method provides better results than the nominal control or a controller based on a random choice of scenarios.At the same time, the number of scenarios found in the maximisation step from Algorithm 2 is comparable to the controller based only on nominal and two extreme scenarios.Further analysis of the performance comparison will be done for Case B in Section 4.3.2 to show that other controllers fail if parametric uncertainty is more significant.

Time performance
The results from Fig. 3 indicate that the local reduction method enables reducing the number of scenarios compared to approaches based on random choice or on time-varying extreme scenarios.Figure 5 shows in the left column the time necessary to solve each step of the local reduction method.The plot in the top left shows the time to solve the minimisation problem as a function of iterations of the local reduction.The iterations correspond to the number of scenarios included in the minimisation problem.The number of scenarios is relatively small (five scenarios), so the time for each iteration is below 3.5 s.
The second plot in the left column shows the time for solving the maximisation problem if the objective is considered (line 2 in Algorithm 2).The time to solve the maximisation problem in every iteration is comparable to the time to solve the minimisation problem.The subsequent two plots show the time for the maximisation of the constraint violation corresponding to the lower and upper bounds over the overall time horizon (one maximisation per time step, lines 4-8 in Algorithm 2).In both cases, the average time to solve a single maximisation problem was 0.35 s.The total time for finding the five scenarios was 10 min 48 s.The algorithm can be parallelised so that the maximisation problems are solved simultaneously [41].Therefore, it can be expected that the time to find a solution in a single iteration of the local reduction method will be equivalent to the solution of the minimisation problem and the maximal time needed to solve the maximisation problems.

Results -Case B
In case B, we ran Algorithm 1 with the time-varying uncertainties from Table 2 and parametric uncertainties δ i,j , η j ∈ [0.96, 1.03].The range for parameters δ i,j , η j was chosen to increase the parametric uncertainty in the dynamics while still ensuring that a controller of the form (69) exists.
Overall performance In contrast to Case A, which has given only five scenarios, the local reduction method in Case B found 101 scenarios.A validation for 500 scenarios is shown in the top right plot in Fig. 3.The plot shows that the controller obtained for 101 scenarios from the local reduction avoided violating constraints.This result indicates that local reduction can handle parametric uncertainty.

Comparison with other approaches
The controller in Case B has also been compared with the same set of controllers as in Case A, obtained for the new range of uncertainty.As expected, the nominal controller and the random controller were unable to satisfy the constraints (middle plots in the right column in Fig. 3).In contrast to Case A, the controller based on nominal and two extreme scenarios was also unable to satisfy the constraints, as shown in the bottom right plot in Fig. 3.The black lines after 24 hours cross the green lines so that the lower bound on the temperature is violated (0.5 • C).Therefore, taking extreme scenarios may be insufficient, as shown in Section 4.1.
The performance of local reduction in handling parametric uncertainty will be also confirmed in the nonlinear case study in Section 4.4.
Computational time performance Figure 5 shows in the right column the time performance of the elements of the local reduction method in terms of time needed to solve them.As expected, the time to obtain the solution to the minimisation problems increases with iterations.The increase is due to the fact that the number of scenarios  In this work, we also assume that the structure of the dynamic feedback policy in (3) is known.The improved computational time of the local reduction can be used to validate whether the chosen control parametrisation is suitable for robustness, because it enables obtaining a solution more quickly.Thus, if the results for a given parametrisation are unsatisfactory, a different parametrisation can be evaluated.

Results -properties of interim scenarios
Finally, we show the impact of the choice of when two scenarios are considered similar in Algorithm 3. The results of varying ϵ are collected in Table 3.The time was obtained using BenchmarkTools.jl[7].
As expected, a high threshold for similarity of scenarios leads to fewer scenarios added to the problem.This is visible in particular in Case B with more significant parametric uncertainty, where the high threshold ϵ = 0.1 led to two scenarios, whereas a lower threshold ϵ = 0.001 led to 101 scenarios.The middle column in Table 3 shows that robustness to the three scenarios in Case B for ϵ = 0.1 is insufficient to robustify the system against random realisations of uncertainties.Conversely, both ϵ = 0.01 and ϵ = 0.001 seem to robustify the system against the random realisations.Unless explicitly stated, the paper considered ϵ = 0.001.
The number of scenarios also affects the time necessary to solve the resulting optimisation problem corresponding to all the scenarios (right column of Table 3).

Dynamics
A further case study is presented to show how our proposed method can be used in nonlinear systems.We want to design a flow controller for a centrifugal compressor.The dynamics for a compressor are nonlinear [8]: where p s and p d are the suction and discharge suction pressures, a 01 , V s , A 1 , L c , J are constant parameters defining the geometry of the compressor, the piping, and the shaft, m SP is the controller for the recycle valve, m r is the mass flow through the recycle valve, m is the mass flow through the compressor, ω is the speed of the shaft of the compressor in rad s −1 , τ is torque provided by a flow controller, τ c is the reaction torque of the compressor.The function Π(•, •) gives the pressure ratio across a compressor as a function of compressor mass flow and speed: The coefficients α i , i = 0, . . ., 5 are usually estimated from operating data.Here we assume α 0 = 2.691, α 1 = −0.014,α 2 = −0.041,α 3 = 0.0009, α 4 = 0.0002, α 5 = 0.00002.The uncertainty in α i is described in Section 4.4.3.
The value of m in and m out captures the external mass flows on the suction and discharge side, respectively.The mass flows depend on the pressures p s and p d , and external pressures p in and p out : where A in , A out , A rec represent the inlet, outlet and recycle valve orifice areas and k in , k out , k rec the respective valve gains.The values of constant parameters were taken from [22].The value of u rec is obtained from an auxiliary PI controller and can take values between 0 and 1, ensured by a smoothed saturation function of the form (71) with β 0 = 0.072, β 1 = 0.071, β 2 = 5.279, β 3 = −0.001.

Optimal control
The objective is to reach the desired flow level m d = 100 kg s −1 without violating speed and flow constraints, imposed due to safety.The objective function was formulated as: where t f = 100 s.The constraints on the mass flow and the speed are: The parametrisation from (77) is typical for centrifugal compressors [22].The torque that can be applied to the compressor must be between zero and 1000 Nm.The bounds on the torque were ensured by a smoothed saturation function of the form (71) with β 0 = 73.324,β 1 = 0.072, β 2 = 0.005, β 3 = 0.

Uncertainties
The uncertainties we considered in this case study are in the valve gains k in , k out , k rec , and correspond to ±5%, and in the parameters α i in the polynomial compressor map (73), and correspond to ±2%.Thus, there are nine uncertain parameters.

Results
To find the flow controller from (77), the dynamics were discretised using the trapezoidal collocation method with time step 0.5 s.

Overall performance
The local reduction method applied to the compressor case study resulted in two scenarios.Figure 6 shows the results of the validation for the two scenarios obtained from the local reduction method with ϵ = 10 −6 .The controller obtained from the local reduction did not violate the constraints on either the mass flow (top left) or the speed of the compressor (top right).Comparison with other approaches If we were to consider all extreme realisations of the uncertainties, we would obtain 2 9 = 512 scenarios.The local reduction method we propose in this paper reduced the number of scenarios to two.Conversely, the nominal controller was not able to satisfy the constraints and both the mass flow and the speed violated their upper limits (second row of plots in Fig. 6).The controllers based on randomly chosen scenarios (from a uniform distribution) are shown in the third row, with black indicating a controller based on two random scenarios.The controller based on two random scenarios was insufficient to ensure constraint satisfaction and both the mass flow and the speed of the compressor violated their upper limits.The minimal number of random scenarios needed to ensure constraint satisfaction was nine (yellow).Finally, the bottom row in Fig. 6 shows the performance of the nominal+extreme controller obtained for three scenarios (one scenario on the lower bound, one scenario on the upper bound, and one scenario with no uncertainty).The controller based on the nominal and two extreme scenarios was insufficient to satisfy constraints.Thus, the comparison with other scenario-based approaches confirms the potential of local reduction for solving robust optimal control problems with parametric uncertainty.

Conclusions
Solving robust nonlinear optimal control problems is challenging, especially if the knowledge about the uncertainty is limited.Scenario-based approaches provide a way of reformulating the optimal control problems as nonlinear optimization problems.The choice of scenarios and their number affects the robustness of the solution as well as computational complexity of the resulting optimisation problems.In this work, we formulated robust optimal control problems with time-varying and parametric uncertainty as semi-infinite optimisation problems to facilitate the choice of scenarios.The new formulation enabled usage of semi-infinite optimisation algorithms, such as local reduction methods.By adding interim worst-case scenarios, the local reduction method enables finding a trade-off between the size of the resulting optimization problem and robustness of the solution to the original optimal control problem.We overcome the dependence on global solvers in the original local reduction formulation by proposing inexact local reduction and providing theoretical bounds on possible constraint violation.The new method consists in solving multiple optimal control problems of reduced size compared to the full scenario-based optimisation.In particular, the small control problems can be solved in parallel, further improving the computational speed.The performance of our approach was evaluated in two case studies with both additive and parametric uncertainty: thermal comfort control in a residential building and mass flow control in a centrifugal compressor.A comparison with common approaches based on a random choice of scenarios and on extreme scenarios indicates that local reduction allows solving robust optimal control problems in an efficient way while ensuring robustness.In particular, the case studies confirm that the proposed inexact local reduction method allows finding worst-case scenarios in the interior of the uncertainty sets.As a result, the new method was able to handle larger parametric uncertainty than other scenario-based approaches.
In this work we required that the constraints must be satisfied for all realisations of the uncertainty.In the future, it would be advisable to look at the conservatism of the obtained solutions and possible relaxations of this requirement.In particular, tighter bounds on constraint violation can be derived if a distribution of the uncertainty is available.Future work could include numerical improvements of approximate local reduction, including warmstarting and use of custom nonlinear optimization solvers, as well as explicit parallelisation of the optimal control problems.

4 j=0u 4
k corresponds to a different optimisation problem of the form (59). Let us now take k = 4, x 0 = 0, A = −0.5,B = 1, and we are looking for the maximal constraint violation for a constant d ∈ [−0.5, 0.5].Let us assume that the current optimal control input u k , k = 0, . . ., 4 is u 0 = u 2 = u 3 = −1 and u 1 = u 4 = 1.The maximisation problem from (59) becomes a maximisation of a fourth order polynomial of d: max d −j (−0.5 + d) j (60) Looking for the maximum of the polynomial (60) yields d ≈ 0.2 which is not on the boundary of the interval [−0.5, 0.5].
(a) Trajectory x corresponding to maximization of upper bound (black) with the respective bounds (green) (b) Validation of the scenarios from local reduction for 500 randomly chosen scenarios uniformly distributed in [0.9, 1.1]

Figure 2 :
Figure 2: Application of local reduction for the unstable system (61)

Figure 3 :
Figure 3: Comparison of local reduction with scenario based approaches in Case A (left) and Case B (right)

Figure 5 :
Figure 5: Impact of significant parametric uncertainty on local reduction for Case A (left) and Case B (right).The first row shows the time necessary for solving the minimization problem (line 3 in Algorithm 3) as a function of iteration of local reduction, corresponding to the number of scenarios in the current set H. The second row shows the time necessary to solve the first maximization problem related to the objective as a function of iteration of local reduction (line 1 in Algorithm 2).The two bottom rows show the time necessary to solve the maximization problems corresponding to n g = 2 trajectory constraints (67) as a function of samples over the horizon of 48 h, N = 192, (lines 4-8 in Algorithm 2)

Table 1 :
Constraint violation obtained for the unstable system (61) with different time horizon N

Table 2 :
Ranges of uncertain parameters throughout the day

Table 3 :
Influence of the tolerance for checking the similarity of scenarios on the resulting number of scenarios, maximal constraint violation over 500 random scenarios, and the time to obtain a solution for the scenarios obtained iteration is greater than in the previous one.At the same time, the time necessary to solve a single maximisation problem remained similar across the iterations.This result indicates the potential for parallelisation to improve performance.