A Posteriori Probabilistic Bounds of Convex Scenario Programs with Validation Tests
Chao Shang Fengqi You







































































A Posteriori Probabilistic Bounds of Convex Scenario Programs with Validation Tests
    Chao Shang, Member, IEEE
        and Fengqi You, Senior Member, IEEE





C. Shang is with Department of Automation, Beijing National Research Center for Information Science and Technology, Tsinghua University, Beijing 100084, China (e-mail: c-shang@tsinghua.edu.cn).

F. You is with College of Engineering, Cornell University, Ithaca, New York 14853, USA (e-mail: fengqi.you@cornell.edu).

    December 30, 2023
==========================================================================================================================================================================================================================================================================================================================================================================================================






Scenario programs have established themselves as efficient tools towards decision-making under uncertainty. To assess the quality of scenario-based solutions a posteriori, validation tests based on Bernoulli trials have been widely adopted in practice. However, to reach a theoretically reliable judgement of risk, one typically needs to collect massive validation samples. In this work, we propose new a posteriori bounds for convex scenario programs with validation tests, which are dependent on both realizations of support constraints and performance on out-of-sample validation data. The proposed bounds enjoy wide generality in that many existing theoretical results can be incorporated as particular cases. To facilitate practical use, a systematic approach for parameterizing a posteriori probability bounds is also developed, which is shown to possess a variety of desirable properties allowing for easy implementations and clear interpretations. By synthesizing comprehensive information about support constraints and validation tests, improved risk evaluation can be achieved for randomized solutions in comparison with existing a posteriori bounds. Case studies on controller design of aircraft lateral motion are presented to validate the effectiveness of the proposed a posteriori bounds.




Scenario approach, stochastic programming, Bernoulli trials, data-driven decision-making.







§ INTRODUCTION



THE widespread presence of uncertainty has been invariably a crucial issue in design, analysis and optimization of complex systems, and the ignorance of uncertainty can lead to decisions that are fragile in the real-world uncertain environment. In the past decades, decision-making under uncertainty has raised immense research efforts across various communities. Typical examples include robust optimal control, where control inputs are designed to yield invariably satisfactory performance under model mismatch, unmeasured disturbance, and measurement noise <cit.>. Likewise, in filter design for fault detection, such effects have been taken into account in order to avoid high false alarm rates <cit.>.



The increasing availability of data has motivated rapid development of data-driven methods in systems and control <cit.>. As a representative, the scenario approach has been widely adopted as a data-driven approach for reliable decision-making under uncertainty <cit.>, where a finite number of constraints defined on past samples of uncertainty are enforced to seek robustness. It has found wide applications in operations of smart grid <cit.>, service systems <cit.>, supply chain management <cit.>, and control design <cit.>. The scenario approach is known to generalize well, which has been revealed theoretically that the probability of constraint violation of the solution can be safely controlled with high confidence by choosing an adequately large sample size <cit.>. To trade robustness for performance, extensions based on constraint discarding have been also developed <cit.>.

Due to the inherent randomness of scenario sampling, the optimal solution of a scenario program is itself a random variable; hence, the risk of the randomized solution, especially the violation probability, is a fortiori uncertain. Therefore, for robust and secure decision-making, a fair evaluation of the risk of the randomized solution must be performed before its implementation in real-world situations. If the risk is considered to be too high, the solution will not be accepted, and further refinement is necessitated. Towards this goal, the simplest approach is to carry out validation test on a collection of new instances of the uncertainty. If few violations are seen on the validation set, the solution is believed to offer a high protection level, and vice versa. Technically, the probabilistic guarantee is established based on finite-sample bounds for Bernoulli trials. The validation strategy is also adopted in the so-called sequential scenario approach originally developed by <cit.>, which solves a sequence of reduced-size scenario programs that are more computationally affordable <cit.>. After obtaining the candidate solution in each iteration, validation test is performed to verify whether its violation probability is smaller than a pre-specified threshold. If so, the sequential algorithm terminates and a final solution is returned.



Recently, an emerging line of research has concentrated on the so-called “wait-and-judge" scenario approach <cit.>. To attain a robust design, sufficient data samples are first incorporated for optimization, and a posteriori assessment of the reliability of scenario-based solutions can be made by counting observed decisive constraints, which can be interpreted as the complexity of the scenario program <cit.>. In comparison with a priori violation probability bounds, significant improvement can be achieved by the wait-and-judge approach, and it is even unnecessary to upper-bound the number of support constraints in advance <cit.>. However, this scheme still has some limitations in practical use. First, there may be computational limits that prohibit the utilization of all available data in the scenario program to seek robustness <cit.>. To alleviate the computational burden, one has to use only a fraction of scenarios for optimization and give up valuable information within unused scenarios. Second, after implementing a scenario-based decision new instances of uncertainty could keep arriving, and one may want to use these out-of-sample data to refine the risk judgement of the decision in hand. In both cases, additional information in scenarios that are not involved in the design problem cannot be effectively utilized for risk evaluation by the present wait-and-judge theory.

In this article, we seek to fill the knowledge gap by establishing a general class of a posteriori probabilistic bounds for scenario-based solutions, where a decision-maker has to solve a scenario program with a fixed sample size, and then further evaluates the solution's risk on out-of-sample realizations of uncertainty. The proposed probabilistic bounds depend simultaneously on the number of decisive support constraints and the outcome of Bernoulli trials. In this way, refined evaluation of solution's risk can be attained by leveraging an increased amount of information from both decisive support constraints and out-of-sample performance on validation data. The established probabilistic bounds turn out to subsume many existing a priori and a posteriori probabilistic bounds as special cases, thereby enjoying widespread generality. To make practical use of the proposed bounds, an efficient parameterization approach is developed. By incorporating information from validation tests, the proposed bounds turn out to be no higher than existing wait-and-judge bounds that only use partial information. We also show that they bear clear interpretations because of their consistency with the validation performance, thereby greatly enhancing the practicability of wait-and-judge scenario optimization. In addition, we establish fundamental lower limits for the proposed bounds, and develop a refinement procedure to reduce the conservatism. Finally, we illustrate our main results with case studies on linear quadratic regulator (LQR) design of aircraft lateral motion, where the proposed bounds provide empirically better judgements than both the wait-and-judge approach and the outcome of Bernoulli trials. Comprehensive comparisons indicate that the efficacy of proposed bounds owes to a sophisticated synthesis of support constraint information and validation test information. For small validation sample size, it mainly relies on support constraint information and largely improves upon bounds based on Bernoulli trials, while under large validation sample size, validation test information tends to have a dominating effect and contribute to reliable judgements.

The layout of this paper is organized as follows. In Section <ref>, the technical background of scenario programs is introduced. Section <ref> formally states the main results of this work, and discusses some practical implementation issues. An illustrative example is provided in Section <ref>, followed by concluding remarks in the last section.

Notations and Definitions. ℕ_0 is the set of non-negative integer numbers, and the set of consecutive non-negative integers {j,⋯,k} is denoted by ℕ_j:k. The p-norm of a matrix is denoted by ·_p, while · represents the Euclidean norm of a vector by convention. The trace operator on a square matrix is denoted by Tr{·}. The indicator function of a subset 𝒮 is defined as 1_𝒮(s). Given a positive integer N, a non-negative integer m ≤ N, and t ∈ (0,1), the Beta cumulative probability function is defined as B_N(t;m) := ∑_i=0^m Nit^i(1-t)^N-i. C^q[0,1] denotes the class of q times differentiable functions with continuous q-th derivative over the interval [0,1], and P^K denotes the set of polynomials of degree K.




§ PRELIMINARIES





 §.§ Scenario Programs


Assume that Δ is a probability space, which is associated with a σ-algebra ℱ and a probability measure ℙ. δ is a random outcome from the triplet (Δ,ℱ,ℙ). Moreover, denote by (Δ^m,ℱ^m,ℙ^m) the Cartesian product of Δ equipped with the σ-algebra ℱ^m and the m-fold probability measure ℙ^m ≐ℙ×⋯×ℙ (m times). In general, the analytic form of ℙ is unknown in practice, but we are allowed to sample a set of independent scenarios {δ^(1),⋯,δ^(m)} from ℙ^m. Here, we refer to ω_m := {δ^(1),⋯,δ^(m)}∈Δ^m as a multi-sample, based on which the following scenario programs can be defined:

    SP_m[ω_m]: {min_x∈ℝ^d    J(x) 
     s.t.    x∈𝒳⊆ℝ^d 
        f(x, δ^(i)) ≤ 0, ∀ i ∈ℕ_1:m .

when m > 0, and

    SP_0[ω_0]: {min_x∈ℝ^d    J(x) 
     s.t.    x∈𝒳⊆ℝ^d
     .

with m = 0 as a degenerate case. In this work, we focus on convex scenario programs, where as regularity conditions, the set 𝒳 is assumed to be convex and closed, the objective function J(x) is convex, and f(x, δ) is lower semi-continuous and convex in x∈ℝ^d with the value of δ fixed. Given a candidate decision x∈𝒳, its violation probability is defined as:

    V(x) = ℙ{δ∈Δ: f(x, δ) > 0
    },

which can also be interpreted as the risk of x, and the corresponding reliability level is 1 - V(x) <cit.>.

Meanwhile, the following assumption is made throughout the paper, which is fairly standard in generic convex scenario programs <cit.>.
[Feasibility and Uniqueness]

For every m ∈ℕ_0 and every multi-sample ω_m ∈Δ^m, the optimal solution x^*_m to the scenario program SP_m[ω_m] exists and is unique.


Note that in order to secure the uniqueness of x_m^*, a suitable tie-break rule can be employed <cit.>. In this way, the optimal solution x_m^* to SP_m[ω_m] is essentially a random variable defined over Δ^m. Under Assumption <ref> and the conventional assumption on the measurability of V(x_m^*), a celebrated probabilistic guarantee in literature is that, the tail probability of V(x_m^*) is dominated by the Beta cumulative probability function <cit.>:

    ℙ^m { V(x_m^*) > ϵ}≤ B_m(ϵ;ζ-1), m > ζ,

where ζ (≤ d) can be the Helly's dimension <cit.> or its tighter bound built with structural information <cit.>. The probabilistic guarantee (<ref>) is closely related to the concept of support constraints, which will be made heavy use of in this paper. 

[Support Constraint, <cit.>]
A constraint of the scenario program SP_m[ω_m] is a support constraint if its removal yields a different optimal solution of the initial problem.


Upon solving SP_m[ω_m], the indices of induced support constraints can be described by a set I_m^* ⊆ℕ_0, which is a set-valued mapping from the multi-sample ω_m. The cardinality of I_m^*, i.e. the associated number of support constraints, is given by s_m^* = |I_m^*|. Likewise, s_m^* is a function of ω_m; hence, both I_m^* and s_m^* are essentially random variables on the triplet (Δ^m,ℱ^m,ℙ^m). The definition of support constraints differs from that of generic active constraints, which are characterized by {i: f(x_m^*, δ^(i)) = 0 }. For convex scenario programs, support constraints are always active constraints, but the converse no longer remains true <cit.>. Meanwhile, under the availability of sufficient scenarios, the number of support constraints is always upper bounded by ζ <cit.>. For a fully-supported problem where the number of support constraints is exactly ζ with probability one (w.p.1), the inequality (<ref>) then becomes an equality <cit.>; in other cases, it inevitably induces conservatism.

Another useful assumption is made as follows.

[Non-Degeneracy, <cit.>]

The solution to the scenario program SP_m[ω_m] coincides w.p.1 with the solution to the program defined by support constraints only.


It is noteworthy that Assumption <ref> is quite mild for general convex scenario programs, because it excludes anomalous situations where the solution defined by support constraints lies exactly on boundaries of other constraints with a nonzero probability.



 §.§ Confidence Intervals for Violation Probability with Validation Tests


Given a candidate solution x, a usual approach to evaluate its violation probability V(x) is to perform Bernoulli trials with some independent validation scenarios {δ̃^(1),⋯,δ̃^(M)}. There are two possible outcomes in each trial, i.e. “success" and “failure", which are characterized by f(x, δ̃^(l)) ≤ 0 and f(x, δ̃^(l)) > 0, respectively. By summarizing results of Bernoulli trials, the empirical violation frequency on the validation dataset can be obtained.

[Empirical Violation Frequency]
Given a set of M independent scenarios {δ̃^(1),⋯,δ̃^(M)}, the empirical violation frequency of a solution x is calculated as:

    r_M = ∑_l=1^M 1_{δ:f(x, δ) > 0 } (δ̃^(l) ).

That is, the number of scenario-based validation constraints that are violated by x.


With x fixed, r_M is a random outcome on ℙ^M, which embodies useful information for estimating the binomial proportion, i.e., the true violation probability V(x). In the context of robust design, a small violation probability is always desired, and hence an upper-bound for V(x) based on r_M shall be constructed to quantify the risk of an unreliable solution. This can be achieved by adopting the following one-sided bounds for evaluating the binomial proportion from a series of Bernoulli trials <cit.>.



Given a solution x and a pre-defined confidence level β^* ∈ (0,1), it then holds that:

    ℙ^M{ V(x) > ρ_M(r_M)  }≤β^*,

where the a posteriori violation probability ρ_M(r_M) = r_M/M + √(logβ^*/(-2M)) is a random variable depending on the realization of r_M.




Given a solution x and a pre-defined confidence level β^* ∈ (0,1), it then holds that:

    ℙ^M{ V(x) > η_M(r_M) }≤β^*,

where the analytic form of η_M(·) is expressed as:

    η_M(l) = {[ min_η{η: B_M(η;l) ≤β^* },                   if l  M;                        1,                    if l=M ] .



The above probabilistic bounds are referred to as a posteriori bounds, in the sense that they provide certificates on V(x) after the value of r_M is revealed <cit.>. Their use in validating scenario-based solutions can be illustrated as follows. Considering the case where there are M=100 samples for validating x, and r_M turns out to be 10, setting β^* = 10^-6, one then obtains the one-sided Chernoff bound ρ_100(10) = 0.3628 and the one-sided C-P bound η_100(10) = 0.3045, respectively, indicating that the violation probability V(x_N^*) does not exceed 36.28% ( or 30.45%) with a very high confidence of 1 - 10^-6 (i.e. 99.9999%). Note that the calculation of ρ_M(r_M) is easier than that of η_M(r_M). However, it is known that ρ_M(r_M) is generally more conservative than η_M(r_M), and it is possible that ρ_M(r_M) have values greater than one, especially when the value of r_M is large; by contrast, η_M(r_M) is always between zero and one, thereby enjoying better practicability.




§ MAIN RESULTS





 §.§ A Posteriori Confidence Level Based on Violation Probability Function


Assume that we are seated in a situation where a scenario program SP_N[ω_N] with fixed N>ζ (N is the initial number of available scenarios for optimization), is first solved to pursue robustness, and then the risk of x^*_N is further evaluated by computing the empirical violation frequency r_M^* on M extra validation samples. Since both the one-sided Chernoff bound and the one-sided C-P bound are applicable to any candidate solution x, we can certainly use them to evaluate x^*_N's risk as a function of r_M^*. However, the optimal solution x_N^* as well as s_N^* have some inherent randomness, whose realizations embody useful information about reliability. In this section, we present a more general theory to concurrently tackle two classes of randomness. To be more specific, the a posteriori bound ϵ(s_N^*,r_M^*) on the violation probability V(x^*_N) is considered as a function of both s_N^* and r_M^*. With a group of functions {ϵ(k,l) } pre-specified for k ∈ℕ_0:ζ, l ∈ℕ_0:M, the following finite-sample probabilistic guarantee can be established as one of our main results.



Assume that {δ^(1), ⋯, δ^(N), δ̃^(1), ⋯, δ̃^(M)} are sampled from the probability space (Δ^N+M,ℱ^N+M,ℙ^N+M), and ϵ(k,l) is an arbitrary [0,1]-valued function where k ∈ℕ_0:ζ, l ∈ℕ_0:M. Then it holds that

    ℙ^N+M{ V(x_N^*) > ϵ(s_N^*,r_M^*) }≤β^*,

where β^* is the optimal value of the following variational problem, defined based on functions {ϵ(k,l) }:

    inf_g ∈C^ζ[0,1] g(1) 
             s.t.   Nk t^N-k·∑_l=0^M  Ml t^M-l (1-t)^l·1_[0,1-ϵ(k,l)) (t) 
                          ≤1/k! d^k/ dt^kg(t), ∀ t ∈ [0,1], k ∈ℕ_0:ζ



 The proof bears resemblance to Theorem 1 in <cit.> so we only provide details for key steps herein. First, the following group of generalized distribution functions is defined:

    F_k(v) = ℙ^k { V(x_k^*) ≤ v  s_k^* = k }, k∈ℕ_0:ζ.

Next, we show that the a posteriori violation probability ℙ^N+M{ V(x_N^*) > ϵ(s_N^*,r_M^*) } can be calculated based on { F_k(v) }_k=0^ζ, which act as the backbone of the whole machinery. Notice that the following events do not overlap with each other:

    { s_N^* = k  r_M^* = l }⋂{ s_N^* = k'  r_M^* = l' } = ∅,

if k  k' or l  l'. Then the following decomposition can be made:

    ℙ^N+M{ V(x_N^*) > ϵ(s_N^*,r_M^*) }
    
    =     ∑_k=0^ζ∑_l=0^M ℙ^N+M{ V(x_N^*) > ϵ(k,l)  s_N^* = k  r_M^* = l }.

We define the events on the probability space Δ^N+M:

    𝒱_k,l = { V(x_N^*) > ϵ(k,l)  s_N^* = k  r_M^* = l }⊆Δ^N+M.

which indicates that the violation probability of x_N^* is larger than ϵ(k,l), and there are k support constraints and l violated validation constraints, respectively. In this way,

    ℙ^N+M{ V(x_N^*) > ϵ(s_N^*,r_M^*) } = ∑_k=0^ζ∑_l=0^M ℙ^N+M{𝒱_k,l}

holds. Next we derive analytic expressions for ℙ^N+M{𝒱_k,l}. By fixing the indices of support constraints, we have:

    ℙ^N+M{𝒱_k,l} = Nkℙ^N+M{𝒜_1 ∩𝒜_2}.

where

    𝒜_1 =    { V(x_N^*) > ϵ(k,l)  s_N^* = k  I_N^* = ℕ_1:k},
    
    𝒜_2 = { r_M^* = l }.

It has been proved by <cit.> that, in the probability space (Δ^N,ℱ^N,ℙ^N), the following events are equal up to a zero probability set:

    𝒜_1 = ℬ_1 ≜{ V(x_k^*) > ϵ(k,l)  s_k^* = k 
                         f(x_k^*,δ^(j)) ≤ 0, j ∈ℕ_k+1:N},

which implies that:

    ℙ^N+M{𝒜_1 ∩𝒜_2}
    
    =     ℙ^N+M{ℬ_1 ∩𝒜_2}
    
    =     ℙ^N+M{    V(x_k^*) > ϵ(k,l)  s_k^* = k   r_M^* = l 
         f(x_k^*,δ^(j)) ≤ 0, j ∈ℕ_k+1:N}
    
    =     ∫_(ϵ(k,l),1]Ml (1-v)^N-k+M-l v^l dF_k(v)

where the last equality is due to the fact that by fixing violation probability v, (1-v)^N-k is the probability that 𝒳_δ^(k+1), ⋯, 𝒳_δ^(N) are satisfied by x_k^*, and Mlv^l(1-v)^M-l is the probability of observing exactly l violated validation constraints. By using the definition of F_k(v) in (<ref>) and taking integral over the interval (ϵ(k,l),1], (<ref>) can be attained. Substituting (<ref>) and (<ref>) into (<ref>) yields:

    ℙ^N+M{ V(x_N^*) > ϵ(s_N^*,r_M^*) }
    
    =     ∑_k=0^ζ∑_l=0^M NkMl∫_(ϵ(k,l),1] (1-v)^N-k+M-l v^l dF_k(v).


Next we show that the RHS of (<ref>) is actually upper-bounded by β^* given in (<ref>). The idea is to take generalized distribution functions { F_k(v) }_k=0^ζ as decision variables subject to the following joint moment conditions <cit.>:

    ∑_k=0^min{m,ζ}mk∫_[0,1] (1-v)^m-k dF_k(v) = 1, m ∈ℕ_0.

Therefore, we have:

    ℙ^N+M{ V(x_N^*) > ϵ(s_N^*,r_M^*) }≤β',

where β' is the optimal value of the following variational problem:

    sup_{ F_k(·) }   ∑_k=0^ζ∑_l=0^M NkMl∫_(ϵ(k,l),1] (1-v)^N-k+M-l v^l dF_k(v) 
     s.t.    ∑_k=0^min{m,ζ}mk∫_[0,1] (1-v)^m-k dF_k(v) = 1, m ∈ℕ_0
         F_0, F_1, ⋯, F_ζ∈𝒞

Here, 𝒞 stands for the set of all positive generalized distribution functions. However, (<ref>) is a generalized moment problem with infinite moment constraints, making the optimization problem difficult to solve analytically. Hence we turn to the following truncated problem:

    sup_{ F_k(·) }   ∑_k=0^ζ∑_l=0^M NkMl∫_(ϵ(k,l),1] (1-v)^N-k+M-l v^l dF_k(v) 
     s.t.    ∑_k=0^min{m,ζ}mk∫_[0,1] (1-v)^m-k dF_k(v) = 1, m ∈ℕ_0:K
        F_0, F_1, ⋯, F_ζ∈𝒞

whose optimal value is β'_K. Because only K constraints are involved in (<ref>), (<ref>) is a relaxation of (<ref>), so we have β'_K ≥β', ∀ K. Meanwhile, the value of β'_K is non-increasing with K. By deriving the Lagrangian function and optimizing over { F_k(·) }_k=0^ζ, we obtain the following dual problem <cit.>:

    inf_λ_0,⋯,λ_K    ∑_m=0^K λ_m 
     s.t.     Nk∑_l=0^M Ml (1-v)^N-k+M-l v^l1_(ϵ(k,l),1] (v) 
        ≤∑_m=k^K λ_m mk (1-v)^m-k, ∀ t ∈ [0,1], k ∈ℕ_0:ζ

where λ_0,⋯,λ_K are dual variables and β^*_K is the optimal value. Therefore, it holds that β' ≤β'_K ≤β_K^*, ∀ K ∈ℕ_0 due to weak duality. We define the polynomial of degree K as g(t) = ∑_m=0^K λ_m t^m, which admits the following properties:

    g(1) = ∑_m=0^K λ_m, 1/k! d^k/ dt^k t^m = {    0,            m < k 
       mkt^m-k,  m ≥ k
     .

By defining t := 1 - v and plugging (<ref>) into (<ref>), (<ref>) then becomes:

    inf_g(·) ∈P_K    g(1) 
     s.t.     1/k! d^k/ dt^kg(t) ≥Nk∑_l=0^M Ml t^N-k+M-l (1-t)^l
                    ·1_[0,1-ϵ(k,l)) (t), ∀ t ∈ [0,1], k ∈ℕ_0:ζ

Denoting by 𝒟 the feasible region of problem (<ref>), it has been proved in <cit.> that the set ( ∪_K ≥ζP_K  ) ∩𝒟 is dense in 𝒟. Therefore we can conclude that:

    ℙ^N+M{ V(x_N^*) > ϵ(s_N^*,r_M^*) }≤β' ≤inf_Kβ_K^* = β^*.

This completes the proof.



Theorem <ref> enjoys wide generality in assessing a posteriori violation probabilities of scenario-based solutions, since it holds for any convex scenario program irrespective of the probability distribution ℙ. Note that the randomness arises from sampling of both N scenarios and M validation samples, and s_N^* is a random variable defined on Δ^N only, whereas r_M^* is defined on the (N+M)-fold product probability space Δ^N+M since it depends on both x_N^* and validation data. Hence, the finite-sample probabilistic guarantee (<ref>) only applies to a design-and-validate procedure in one-shot, and one is not allowed to repeat the validation phase over and over (with the same or a different solution) until a satisfactory performance is obtained.


Next we show that with independence of {ϵ(k,l)} on k and/or l imposed, Theorem <ref> turns out to unify many existing theoretical results, including the wait-and-judge approach in <cit.> as well as the generic probabilistic guarantee (<ref>) in <cit.>. First, the following result reveals itself as an obvious corollary of Theorem <ref>, which enforces the two-indexed bounds {ϵ(k,l) } to depend on k only.

[Theorem 1 in <cit.>]

Take ϵ(k,l) = ϵ(k), ∀ k ∈ℕ_0:ζ, l ∈ℕ_0:M. Under Assumptions 1 and 2, it then holds that:

    ℙ^N{ V(x_N^*) > ϵ(s_N^*) }≤β^*,

where β^* is the optimal value of the following problem:

    inf_g ∈C^ζ[0,1] g(1) 
             s.t.   Nk t^N-k·1_[0,1-ϵ(k)) (t) ≤1/k! d^k/ dt^kg(t), k ∈ℕ_0:ζ



A straightforward substitution of the identity

    ∑_l=0^M Ml (1-v)^M-l v^l = 1

into (<ref>) yields the proof.


As an established bound in the wait-and-judge approach, Corollary <ref> implies that only information embedded in training data can be utilized to evaluate V(x_N^*). As aforesaid, one may want to carry out Bernoulli trials on extra scenarios in conjunction with (<ref>) to update the risk judgement. In this sense, Theorem <ref> has wider applicability than Corollary <ref> because the information within support constraints and empirical validation tests can be synthesized.


Based on Corollary <ref>, one can easily cast the generic probabilistic guarantee (<ref>) as a corollary of Theorem <ref> by trivially assuming constant violation probabilities.


By taking ϵ(k,l) ≡ϵ, ∀ k ∈ℕ_0:ζ, l ∈ℕ_0:M, and under Assumptions 1 and 2, one obtains:

    ℙ^N { V(x_N^*) > ϵ}≤ B_N(ϵ;ζ-1).



See Corollary 1 in <cit.>.




 §.§ Derivation of A Posteriori Probabilistic Bounds

Theorem <ref> is shown to be a generalization of many important probabilistic bounds in scenario optimization, with particular options for {ϵ(k,l) } used. Hence, the freedom in choosing {ϵ(k,l) } endows the proposed scheme with capability in evaluating the solution's reliability based on all possible realizations of s_N^* and r_M^*. One typically wishes to ensure that the event V(x_N^*) ≤ϵ(s_N^*,r_M^*) occurs with a suitably high confidence, and in order to explicitly control the probability of a wrong declaration V(x_N^*) > ϵ(s_N^*,r_M^*), it is desirable to first assign β^*, which is typically small, and then determine {ϵ(k,l) } in reverse. In principle, there are infinitely many options of {ϵ(k,l) } that lead to the same confidence level β^*. For practical use, a certain choice of {ϵ(k,l) } is suggested in the following theorem. It subsumes Theorem 2 in <cit.> as a special case with M = 0 (no validation tests) and specific coefficients {a_m ≡ 1/(N+1), ∀ m ∈ℕ_0:N}, where the induced bound ϵ(k) depending only on s_N^* equals ϵ(k,0). However, the present proof machinery has a significant departure from that in <cit.>, which is no longer applicable to the present case where validation test results have been included and general non-negative coefficients { a_m } are considered. Before proceeding, we present the following preparatory lemma. 

For fixed values of m and N (m<N), B_N(t;m) is a strictly decreasing function in t ∈ (0,1).


See <cit.>.




Define the following polynomial in t:

    h_N,M(t;k,l) 
    
    =     β∑_m=k^N a_m mk t^m-k - Nk t^N-k B_M(1 - t; l),

where β∈ (0,1) is the confidence level, and coefficients {a_m} satisfy

    a_m≥ 0, m ∈ℕ_0:N, ∑_m=ζ^N-1 a_m > 0, ∑_m=0^N a_m = 1.


(i) h_N,M(t;k,l) has only one solution in the interval (0,1) for all k ∈ℕ_0:ζ and l ∈ℕ_0:M.

(ii) Letting the root be t_N,M(k,l) and ϵ_N,M(k,l) = 1 - t_N,M(k,l), it then holds that:

    t_N,M(k,l) > t_N,M(k,l+1),  ϵ_N,M(k,l) < ϵ_N,M(k,l+1),    
    
    l ∈ℕ_0:M-1.


(iii) The choice of above functions {ϵ_N,M(k,l) } gives rise to the following finite-sample probabilistic guarantee:

    ℙ^N+M{ V(x_N^*) > ϵ_N,M(s_N^*,r_M^*) }≤β.



 (i) First we define the following polynomial:

    h̃_N,M(t;k,l) 
    ≜    β∑_m=k^N a_m ·mk· t^m-N - Nk B_M(1 - t; l) 
    
                       =    h_N,M(t;k,l) · t^k-N

Since k≤ζ<N and ∑_m=ζ^N a_m > 0, ∑_m=k^N-1 a_m ·mk· t^m-N is a strictly decreasing function in t ∈ (0,1). In addition, according to Lemma 1, -B_M(1 - t; l) is also strictly decreasing in t when l < M, and B_M(1 - t; l) ≡ 1 when l = M. Therefore, in both cases, h̃_N,M(t;k,l) must be strictly decreasing in (0,1). Meanwhile, notice that

    lim_t → 0^+h̃_N,M(t;k,l) = + ∞, 
    h̃_N,M(1;k,l)     = β∑_m=k^N a_m ·mk - NkB_M(0;l) 
       ≤β∑_m=0^N a_m ·Nk - Nk
        = (β - 1)Nk < 0

Hence, h̃_N,M(t;k,l) has exactly one root in (0,1), and so does h_N,M(t;k,l).

(ii) From the definition it is obvious that B_M(1 - t; l) < B_M(1 - t; l+1), which implies h̃_N,M(t;k,l) > h̃_N,M(t;k,l+1). Since h̃_N,M(t;k,l) is strictly decreasing in t ∈ (0,1), we then have (<ref>).

(iii) We first parameterize a candidate solution to (<ref>) as:

    g(t) = β∑_m=0^N a_m · t^m.

Then we will show that g(t) is essentially a feasible solution to (<ref>) with ϵ(k,l) set as ϵ_N,M(k,l). Its k-th order derivative can be computed as:

    1/k! d^k/ dt^kg(t) = β∑_m=k^N a_m ·mkt^m-k, k ∈ℕ_1:ζ.

Because ϵ_N,M(k,l+1) > ϵ_N,M(k,l), the entire interval [0,1 - ϵ_N,M(k,0)) can be split as

    [0,1 - ϵ_N,M(k,0)) 
    
    =     ⋃_l=0^M [1 - ϵ_N,M(k,l+1),1 - ϵ_N,M(k,l)),

where ϵ_N,M(k,M+1) = 1 is taken for granted for notational clarity. Considering the case where t falls within a certain interval [1 - ϵ_N,M(k,l+1),1 - ϵ_N,M(k,l)), we have:

    1/k! d^k/ dt^kg(t) 
        - Nk∑_j=0^M Mj t^N-k+M-j (1-t)^j1_[0,1-ϵ_N,M(k,j)) (t) 
    
    =     1/k! d^k/ dt^kg(t) - Nk∑_j=0^l Mj t^N-k+M-j (1-t)^j
    
    =     β∑_m=k^N a_m ·mkt^m-k - Nk t^N-k B_M(1-t;l) 
    
    =     h_N,M(t;k,l) 
    
    >     0

The last inequality is due to the fact that t_N,M(k,l) = 1 - ϵ_N,M(k,l) is the root of h_N,M(t;k,l), and for t ∈ [0, 1 - ϵ_N,M(k,l)), h_N,M(t;k,l)>0 holds. Therefore, if we set ϵ(k,l) in (<ref>) to be ϵ_N,M(k,l), constraints in (<ref>) hold for all t ∈ [0,1]. As such, g(t) is feasible for problem (<ref>), which indicates β^* ≤β. Finally we have:

    ℙ^N+M{ V(x_N^*) > ϵ_N,M(s_N^*,r_M^*) }≤β^* ≤β.



The following remarks are made in order.

The utilization of a posteriori probabilistic certificates {ϵ_N,M(k,l) } does not require any distributional information of uncertainty. Such a distribution-free characteristic makes the proposed bounds particularly suitable for data-driven decision-making.



The utilization of a posteriori probabilistic certificates {ϵ_N,M(k,l) } should not be oriented towards a tradeoff design with a prescribed tolerance on V(x_N^*). Rather, they shall be used for a posteriori risk judgement of a scenario-based solution. Of common interest, a scenario program SP_N[ω_N] may be constructed as a relaxation of the infinite-dimensional constraint f(x, δ) ≤ 0, ∀δ∈Δ, or as an approximation of the chance constraint ℙ{f(x, δ) > 0 }≤ϵ with the sample size N determined based on the prior bound (<ref>). However, the computational effort for solving SP_N[ω_N] may be unaffordable, and one has to settle for a reduced sample size but leverages unused samples to evaluate the risk. Another situation is that after implementing the randomized solution, one may be interested in further evaluating its risk with new samples successively accumulated. In such situations, {ϵ_N,M(k,l) } show greater promise than the wait-and-judge bounds {ϵ_N(k) } in addressing the practical need for risk evaluation.



The monotonicity property (<ref>) is consistent with our intuition that the lower the empirical violation frequency, the lower the a posteriori violation probability. This can be also evidenced from the example in Fig. <ref>, where {ϵ_N,M(k,l)} with β = 10^-6, ζ = 10, N = 50, M = 30, and { a_m ≡ 1/(N+1) } are portrayed. Meanwhile, it is easy to prove that {ϵ_N,M(k,l)} are increasing in k, which preserves the interpretation of generic wait-and-judge approach <cit.> that less support constraints indicate a potentially lower risk of randomized solutions.




Note that bounds {ϵ_N,M(k,l) } and { t_N,M(k,l) } derived in Theorem <ref> depend on N and M. For notational simplicity, we will use {ϵ(k,l) } and { t(k,l) } in the sequel when no confusions are caused. Meanwhile, we use {ϵ(k) } to denote the a posteriori bounds in the wait-and-judge approach <cit.> that depends on s_N^* only. Now we establish a result that reveals the value of utilizing more validation information in {ϵ_N,M(k,l) }.



With { a_m } and β fixed, it then holds that:

    ϵ(k) =     ϵ_N,0(k,0) 
    
    =     ϵ_N,M(k,M) 
    
    >     ϵ_N,M(k,l), k ∈ℕ_0:ζ, l ∈ℕ_0:M-1




The second equality can be obtained by substituting the identity B_M(1-t;M) = 1 into (<ref>), and the inequality stems from (<ref>).


Corollary <ref> shows that the proposed bound ϵ(s_N^*, r_M^*) is always no higher than ϵ(s_N^*) that depends on partial information within N scenarios, as can be observed from Fig. <ref>. This is an expected outcome because having more information implies having a more accurate evaluation of a solution's risk. Moreover, even for a large l, ϵ(k,l) ≤ϵ(k) still holds true, which admits a sensible explanation. It has been revealed in <cit.> that with high confidence r_M^* / M cannot be much higher than s_N^* / N due to the correctness of the wait-and-judge bound ℙ^N{ V(x_N^*) > ϵ(s_N^*) }≤β. As a consequence, seeing an excessively large r_M^* is issued with a negligible probability, so the change in ϵ(k,l) with large l does not exert significant influence on the overall confidence β. Similarly, as implied by <cit.>, the joint distribution of (s_N^*,r_M^*) is concentrated in the region where V(x_N^*) = 𝔼{ r_M^* } / M and s_N^* / N have the same order of magnitude. This is also the region where the actual improvement of ϵ(k,l) over ϵ(k) occurs, as also confirmed by Fig. <ref>. In this sense, ϵ(k,l) shall be viewed as a refinement of ϵ(k) based on the new information carried by M validation samples.



Next we investigate the incremental monotonicity of generic one-sided C-P bounds {η_M(l)}, which turns out to be well inherited by {ϵ_N,M(k,l) }. This can be precisely stated as follows. Suppose that upon obtaining x_N^*, k support constraints have been found, and l violations have occurred on M validation samples {δ̃^(j)}, giving rise to the one-sided C-P bound η_M(l). Afterwards, a validation sample δ̃^(M+1) carrying new information arrives, based on which the new probabilistic bound η_M+1(·) is derived. More exactly, the new bound shall be η_M+1(l+1) if f(x^*_N,δ̃^(M+1)) > 0, and η_M+1(l) if f(x^*_N,δ̃^(M+1)) ≤ 0. Before proceeding, the following lemma is presented.


For fixed values of t and m, B_N(t;m) is strictly decreasing in N (N ≥ m).


For i ∈ℕ_1:N, it is obvious that 
    N+1i = Ni + Ni-1.
 Therefore, we have:

    B_N+1(t;m) 
    
    =     ∑_i=0^m N+1it^i(1-t)^N+1-i
    
    =     (1 - t)^N+1 + ∑_i=1^m Nit^i(1-t)^N+1-i_=(1-t)B_N(t;m)
        + ∑_i=1^m Ni - 1t^i(1-t)^N+1-i_=t B_N(t; m - 1)
    
    =     (1-t)B_N(t;m) + t B_N(t; m - 1) 
    
    <     (1-t)B_N(t;m) + t B_N(t;m) 
    
    =     B_N(t;m)

This completes the proof of Lemma <ref>.


Then the following relations hold between {η_M(·) } and {η_M+1(·) }.


It holds that:

    η_M(l) > η_M+1(l), l ∈ℕ_0:M, 
    η_M+1(l+1) > η_M(l), l ∈ℕ_0:M-1, 
    η_M+1(M+1) = η_M(M) = 1.



(<ref>) immediately follows from (<ref>). Because η_M(l) is the root of B_M(η,l) = β^* in (0,1), and B_M(η,l) is strictly decreasing in η, it remains to show that:

    B_M(η; l) > B_M + 1(η; l), l ∈ℕ_0:M, 
    
    B_M+1(η; l+1) > B_M(η; l), l ∈ℕ_0:M-1.

A straightforward usage of Lemma <ref> results in (<ref>), along with the inequality B_M+1(1 - η; M-l-1) < B_M(1 - η; M-l-1), l ∈ℕ_0:M-1. Hence, 1 - B_M + 1(1 - η; M-l-1) > 1 - B_M(1 - η;M-l-1) holds, which amounts to (<ref>).


In view of Theorem <ref>, rational refinements of ϵ(s_N^*, r_M^*) can be made once more validation samples are successively accumulated. If constraint violation is witnessed on a new validation sample, a higher risk level of x_N^* is obtained, which is precisely stated as η_M+1(l+1) > η_M(l). Conversely, if x_N^* agrees with the new validation sample, the risk level should be discounted, i.e. η_M+1(l) < η_M(l). Interestingly, such properties are also possessed by {ϵ_N,M(k,l) }.



With values of { a_m } and β fixed, the roots of { g_N,M(t;k,l) } satisfy

    t_N,M(k,l) < t_N,M+1(k,l), l ∈ℕ_0:M, 
    
    t_N,M+1(k,l+1) < t_N,M(k,l), l ∈ℕ_0:M-1, 
    
    t_N,M+1(k,M+1) = t_N,M(k,M),

thereby indicating the following relations:

    ϵ_N,M(k,l) > ϵ_N,M+1(k,l), l ∈ℕ_0:M, 
    ϵ_N,M+1(k,l+1) > ϵ_N,M(k,l), l ∈ℕ_0:M-1, 
    ϵ_N,M+1(k,M+1) = ϵ_N,M(k,M).



The proof is similar to that for Theorem <ref>, and is hence omitted for brevity.

Before closing this subsection, we deal with the computational issue. Although {ϵ_N,M(k,l)} have no analytic expressions, one only needs to perform bisection to numerically search the root of h̃_N,M(t;k,l) in (0,1) to compute {ϵ_N,M(k,l)}, as implied by Theorem <ref>(i). The pseudo code for calculating {ϵ_N,M(k,l)} for all k and l is outlined in Algorithm 1, which can be readily implemented with numerical softwares. In order to reduce the search space, the monotonicity property (<ref>) is utilized, and the enumeration of l is made from M to 0. For a particular choice of k and l, one only needs to execute the inner-loop in Algorithm 1 with LB = 0 and UB = 1.





 §.§ Fundamental limits and a refinement procedure

Above derivations imply that lower values of {ϵ(k,l) } give more informative certificates on the solution's reliability. Hence a natural question comes up: given the confidence β^*, what is the best possible value of ϵ(k,l) satisfying (<ref>) in Theorem <ref>? Next we establish lower limits for {ϵ(k,l) }, where ϵ(k,l) is restricted to be strictly increasing in l, which is in agreement with the intuition that seeing more violations leads to a discounted reliability level.



Consider all admissible functions {ϵ(k,l) } satisfying (i) the probabilistic guarantee (<ref>) in Theorem <ref> with a prescribed β^* and (ii) the monotonicity property ϵ(k,l) < ϵ(k,l+1),  l ∈ℕ_0:M-1. Then for k ≥ 1 a lower limit ϵ(k,l) of ϵ(k,l) is given by the root of the following equation in (0,1):

    ∑_j=0^l z_j B_N + M(ϵ; k + j - 1) = β^*,

where

    z_j = NkMj/N+Mk+j·k/k + j.

And, ϵ(0,l) = 0, ∀ l ∈ℕ_0:M.



We show by contradiction that ϵ(k,l) < ϵ(k,l) is impossible for k ≥ 1. Assume that there exist admissible functions {ϵ(k,l) } satisfying (i) and (ii), and there exist k' and l' such that ϵ(k',l') < ϵ(k',l'). Then we consider a fully-supported problem SP'_N[ω_N] where s_N^* = k' ≤ζ occurs w.p.1. When k' ≥ 1, the probability density function of V(x_N^*) in SP'_N[ω_N] is explicitly expressed as <cit.>:

    p(α) = Nk'α^k'-1· k' · (1-α)^N-k'.

Hence, for this specific fully-supported problem SP'_N[ω_N], it follows that:

    ℙ^N+M{ V(x_N^*) > ϵ(s_N^*,r_M^*) }
    
    =     ∑_j=0^M ℙ^N+M{ V(x_N^*) > ϵ(k',j)  r_M^* = j }. 
    
    =     ∑_j=0^M ∫_ϵ(k',j)^1 Mj p(α) ·α^j · (1-α)^M-j dα
    
    =     ∑_j=0^M ∫_ϵ(k',j)^1 Nk'Mjk' α^k'+j-1·  (1-α)^N+M-k'-j dα
    
    =     [Integration by parts]
    
    =     ∑_j=0^M z_j ∑_i=0^k'+j-1N+Miϵ(k',j)^i [1 - ϵ(k',j)]^N+M-i
    
    =     ∑_j=0^M z_j B_N + M(ϵ(k',j); k' + j - 1)

Recall that ϵ(k',l) < ϵ(k',l) has been assumed, and {ϵ(k,l) } satisfies the monotonicity property ϵ(k,l) < ϵ(k,l+1),  l ∈ℕ_0:M-1. It immediately follows that ϵ(k',j) < ϵ(k',l), ∀ j ∈ℕ_0:l, and one obtains:

    ℙ^N+M{ V(x_N^*) > ϵ(s_N^*,r_M^*) }
    
    =     ∑_j=0^M z_j B_N + M(ϵ(k',j); k' + j - 1) 
    ≥    ∑_j=0^l z_j B_N + M(ϵ(k',j); k' + j - 1) 
    
    >     ∑_j=0^l z_j B_N + M(ϵ(k',l); k' + j - 1) 
    
    =     β^*

where the last inequality stems from Lemma 1 that B_M(t;m) is strictly decreasing in t, and the last equality is due to (<ref>). However, (<ref>) indicates that SP'_N[ω_N] contradicts the distribution-free nature of {ϵ(k,l) } in (i), i.e., the probabilistic guarantee (<ref>) with a given β^*.



Fig. <ref> depicts a posteriori bounds {ϵ(k,l) } under {a_m ≡ 1 / (N+1)} and the fundamental lower limits {ϵ(k,l) } (β = 10^-6, ζ = 8, N = 100, M = 5), where two surfaces are not far away from each other except for k=0. Hence, we pay a small price for distribution-free certificates {ϵ(k,l) } since even knowing k ≥ 1 beforehand does not provide too much improvement in risk evaluation. Indeed, this well inherits the merits of the wait-and-judge scenario optimization <cit.>, which can be interpreted as that the realization of s_N^* is rather informative for a posteriori risk evaluation <cit.>.



Obviously, {ϵ(k,l) } are decided by coefficients {a_m}, and different choice of {a_m} lead to different probabilistic certificates. A natural question is that whether {ϵ(k,l) } can be improved to approach the fundamental limits {ϵ(k,l) } by flexibly adjusting {a_m}. This can be formally cast as the following multi-objective optimization problem:

    max     t(k,l), k∈ℕ_0:ζ, l∈ℕ_0:M
     s.t.     0< t(k,l) <1, k∈ℕ_0:ζ, l∈ℕ_0:M
          β∑_m=k^N a_m mk t(k,l)^m-N
        = Nk B_M(1 - t(k,l); l), k∈ℕ_0:ζ, l∈ℕ_0:M
        ∑_m=ζ^N-1 a_m ≥τ, ∑_m=0^N a_m = 1, a_m ≥ 0, m ∈ℕ_0:N

where τ > 0 is a small positive number to ensure (<ref>). Clearly, (<ref>) is non-convex and is hence difficult to solve. Next we seek to devise a tailored algorithm. Suppose that a certain choice of coefficients {a_m} is feasible for (<ref>). The following theorem provides a sufficient and necessary condition for making improvements over the current {a_m}.



Given N, M and β. Assume that {a_m} and {a'_m} satisfying (<ref>) are different coefficients of polynomials h̃_N,M(t;k,l) and h̃'_N,M(t;k,l), whose roots in (0,1) are {t(k,l)} and {t'(k,l)}, respectively. Then t'(k,l) ≥ (>) t(k,l) holds if and only if

    β∑_m=k^N a'_m mk t(k,l)^m-N≥ (>) Nk B_M(1 - t(k,l);l).




It is known that h̃'_N,M(t;k,l), whose root is t'(k,l), is strictly decreasing in (0,1). Therefore, when t'(k,l) ≥ (>) t(k,l), it holds that h̃'_N,M(t(k,l);k,l) ≥ (>) 0, which leads to (<ref>). The reverse is also true, which completes the proof.


To derive refined coefficients { a'_m } based on present ones { a_m }, it suffices to ensure (<ref>) for all k and l and lift the LHS of (<ref>) as much as possible. This can be achieved, for instance, by simply solving the following linear program (LP):

    max_{ a_m' }    ∑_k=0^ζ∑_l=0^M ∑_m=k^N a_m' mk t(k,l)^m-N
     s.t.       β∑_m=k^N a'_m mk t(k,l)^m-N
        ≥Nk B_M(1 - t(k,l);l), ∀ k ∈ℕ_0:ζ,  l ∈ℕ_0:M
        ∑_m=ζ^N-1 a_m' ≥τ, ∑_m=0^N a_m' = 1, a_m' ≥ 0, m ∈ℕ_0:N


By successively solving LP (<ref>), coefficients { a_m } and the associated {t(k,l)} can be refined. The implementation details are summarized in Algorithm 2. It turns out that the resultant roots {t(k,l)} and a posteriori bounds {ϵ(k,l) } bear Pareto optimality.






In Algorithm 2, the sequence {t^( ite)(k,l)} converges, i.e.,

    lim_ ite→ +∞ t^( ite)(k,l) = t^*(k,l) < + ∞, ∀ k, l,

and {t^*(k,l)} lie on the Pareto optimal boundary of (<ref>). That is, there exist no coefficients {ã^*_m} that give {t̃^*(k,l) } strictly dominating { t^*(k,l) }.



Following the same notations in Theorem <ref>(i), we consider the polynomial h̃_N,M(t;k,l) parameterized by { a_m^( ite+1)}. On one hand, h̃_N,M(t^( ite + 1)(k,l);k,l) = 0. On the other hand, h̃_N,M(t^( ite)(k,l);k,l) ≥ 0 due to (<ref>). As a consequence, t^( ite)(k,l) ≤ t^( ite + 1)(k,l) holds in each iteration due to the monotonicity of h̃_N,M(t;k,l) established in Theorem <ref>(i). Meanwhile, since t^( ite)(k,l) is upper-bounded from 1, it must converge to a finite limit t^*(k,l). However, the sequence { a_m^( ite)} does not necessarily converge. Because 0 ≤ a_m^( ite)≤ 1, there is a convergent subsequence { a_m^( ite_i)} with indices 1 ≤ ite_1 <  ite_2 < ⋯ according to the Bolzano-Weierstrass theorem, i.e.,

    lim_i→ +∞ a_m^( ite_i) = a_m^* < +∞, ∀ m.


Because h̃_N,M(t;k,l) is continuous, by a continuity argument it is not difficult to show that { a_m^* } is an optimal solution to (<ref>) with {t(k,l) = t^*(k,l)}, and {t^*(k,l)} are exactly the roots of h̃_N,M(t;k,l) induced by { a_m^* }. Then we proceed by contradiction. Suppose that there exist {ã^*_m} that give {t̃^*(k,l) } strictly dominating { t^*(k,l) }. That is, t̃^*(k,l) ≥ t^*(k,l), ∀ k,l and there exists (k',l') such that t̃^*(k',l') > t^*(k',l'). Then it immediately follows from Theorem <ref> that ∀ k,l:

    β∑_m=k^N ã^*_m mk t^*(k,l)^m-N≥   Nk B_M(1 - t^*(k,l);l) 
    
    =    β∑_m=k^N a^*_m mk t^*(k,l)^m-N,

and for the particular indices (k',l'),

    β∑_m=k'^N ã^*_m mk' t^*(k',l')^m-N >    Nk' B_M(1 - t^*(k',l');l') 
    
    =    β∑_m=k'^N a^*_m mk' t^*(k',l')^m-N.

Consequently, {ã^*_m} is not only feasible for (<ref>) with {t(k,l) = t^*(k,l)}, but also leads to an objective value strictly higher than { a_m^* } does. This contradicts the optimality of { a_m^* }.


To obtain different Pareto optimal certificates, one could use different initializations of { a_m }. Fig. <ref> profiles {ϵ(k,l) } derived with {a_m ≡ 1 / (N+1)} and those after refinement with β = 10^-6, ζ = 8, N = 100, M = 5. It can be seen that values of {ϵ(k,l) } for all k and l become lower after refinement. While doing so one shall keep in mind that no matter how much improvement can be made over {ϵ(k,l) }, one cannot go below the fundamental limits {ϵ(k,l) } in Fig. <ref>.











§ CASE STUDIES ON LQR DESIGN OF AIRCRAFT LATERAL MOTION




 §.§ Problem Setup

In this section, we adopt a simulation example of controller design of aircraft lateral motion from <cit.>, which has been widely used as a testbed for probabilistic control design <cit.>, to demonstrate the developed theory. The system is described by the following state-space model:

    ẋ(t) = A(δ)x(t) + B(δ)u(t),

where system states x(t) = [x_1(t) x_2(t) x_3(t) x_4(t)]^ T are the blank angle and its derivative, the sideslip angle, and the yaw rate, respectively. Control inputs u(t) = [u_1(t) u_2(t)]^ T include the rudder deflection and the aileron deflection. System matrices A(δ) and B(δ) are influenced by uncertain parameters δ, as given by:

    A(δ) =  [
    [             0             1             0             0;             0           L_p           L_β           L_r;           g/V             0           Y_β            -1;          N_β̇           N_p N_β + N_β̇Y_β    N_r - N_β̇;               ] ], 
    
    B(δ) =  [
    [                 0                 0;                 0             L_δ_a;             Y_δ_r                 0; N_δ_r + N_β̇Y_δ_r             N_δ_a;                   ] ],


    δ = [L_p L_β L_r g/V Y_β N_β̇ N_p N_r L_δ_a Y_δ_r N_δ_r N_δ_a]^ T.

It is assumed that a 10% perturbation of δ is added around its nominal value δ̅, which is uniformly distributed.

The control goal is to design a state feedback controller u=Kx such that real parts of eigenvalues of the closed-loop system are smaller than -γ < 0. That is, the system is stabilized with a desired decay rate γ > 0. In this case we choose γ = 0.5. A possible formulation of the controller design problem is given by the following semi-definite program (SDP), which involves an infinite number of linear matrix inequalities (LMIs):

    min_P ≻ 0,Y     Tr{ P }
     s.t.     P ≽θ I 
        A(δ)P + PA(δ)^ T + B(δ)Y + Y^ TB(δ) + 2γ P ≼ 0, 
                                             ∀δ∈Δ

where P ∈ℝ^4 × 4 and Y ∈ℝ^2 × 4. θ = 0.01 is a small positive number to ensure the positive definiteness of P. The control gain K can be finally computed as K = YP^-1.



Assume that there is no complete knowledge about Δ and only some samples of δ can be collected. To obtain an approximate solution to (<ref>), we turn to the following scenario program with a total of d = 18 free decision variables in matrices P and Y:

    min_P ≻ 0,Y     Tr{ P }
     s.t.     P ≽θ I 
        A(δ^(i))P + PA(δ^(i))^ T + B(δ^(i))Y + Y^ TB(δ^(i)) 
                            + 2γ P ≼ 0, i ∈ℕ_0:N

where {δ^(i)}_i=1^N are randomly collected scenarios of uncertain parameters. To solve the large-scale SDP (<ref>), we use cvx package in MATLAB equipped with the MOSEK solver <cit.>.





 §.§ Results and Discussions

In the simulation phase, we randomly generate N = 500 independent scenarios {δ^(i)} for solving (<ref>), and M = 500 validation scenarios for empirically estimating the violation frequency of the LMI in (<ref>). The confidence is set as β = 10^-6 (practical certainty). Since there is no knowledge about a tighter bound of Helly's dimension, we set ζ = d = 18. The generic a priori bound (<ref>) yields ϵ = 0.0889. The attractiveness of a posteriori bound lies in that, upon seeing x_N^*, s_N^* and r_M^*, the certificate of V(x_N^*) can be refined. For example, in a particular simulation run, s_N^* = 3 support constraints and r^*_M = 2 times of violations have been revealed. Accordingly, we use a_m ≡ 1/(N+1) to derive ϵ(s_N^*,r^*_M) = 0.0268, indicating that with confidence 99.9999% the violation probability V(x_N^*) is no more than 2.68%, which is much tighter than the a priori bound. In contrast, ϵ(s_N^*) = 0.0486 with support constraints information used only <cit.>. If the one-sided C-P bound (<ref>) for Bernoulli trials is adopted, we obtain η_M(r_M^*) = 0.0376. Hence in this instance, by synthesizing comprehensive information from both decisive support constraints and validation tests, a tighter certificate on the violation probability V(x_N^*) can be obtained.

Next, we carry out 2000 Monte Carlo simulation runs for N=100,200,500 and M=100,200,500,1000 with β = 10^-6. In each run, after deriving the solution x^*_N, we obtain a high-fidelity estimate V̂(x^*_N) of its true violation probability V(x^*_N) with 10^5 additional Monte Carlo samples. Then the gap between a certain bound with V̂(x^*_N) can be calculated as a performance quantification, which, for instance, is GAP = ϵ(s_N^*,r_M^*) - V̂(x^*_N) for the proposed bound. By summarizing results in 2000 runs, the mean value μ_ GAP and the standard deviation σ_ GAP of the gap are further calculated and summarized in Table <ref>.

Many interesting observations can be attained from Table <ref>. First, in all cases the proposed bound ϵ(s_N^*,r_M^*) yields the smallest mean and standard deviation of the gap, indicating its significantly reduced conservatism. This is because both ϵ(s_N^*) and η_M(r_M^*) only rely on partial information. With fixed N and small M, η_M(r_M^*) tends to be quite conservative, while the proposed bounds provides the tightest certificates by using comprehensive information. With M increasing, both η_M(r_M^*) and ϵ(s_N^*,r_M^*) get lower; when M is sufficiently large, the difference between ϵ(s_N^*,r_M^*) and η_M(r_M^*) vanishes, because the effect of validation tests becomes increasingly dominant. Note that in this case the wait-and-judge approach gives a constant bound ϵ(s_N^*) without leveraging validation information. On the other hand, with N increasing, the conservatism of both ϵ(s_N^*) and ϵ(s_N^*,r_M^*) is reduced, since a larger N leads to improved robustness of randomized solutions. This can also be evidenced from the fact that, under the same M, the difference between ϵ(s_N^*,r_M^*) and ϵ(s_N^*) declines with N increasing, while the difference between η_M(r_M^*) and ϵ(s_N^*,r_M^*) becomes more pronounced. Henceforth, the proposed bound ϵ(s_N^*,r_M^*) features a sophisticated integration of information from both support constraints and validation tests. When N is small, ϵ(s_N^*,r_M^*) tends to depend more on validation test, while when M is small, information carried by s_N^* becomes dominant. This yields an explanation for empirically better performances of ϵ(s_N^*,r_M^*) under various choices of N and M.


  





Table <ref> reports frequencies of s_N^* and statistics of r_M^* with M = 1000 during 2000 Monte Carlo runs. The realization of s_N^* is always between 2 and 8, which is much smaller than d = 18. Hence, the information carried by s_N^* helps give a non-conservative certificate due to the non-fully-supportedness, which is not rare in practical situations <cit.>. In addition, as analyzed in Section III, the actual improvement in ϵ(s_N^*, r_M^*) occurs primarily in the region where the joint distribution of (s_N^*, r_M^*) is concentrated. This can also be seen from Table <ref> where the mean value of r_M^*/M, which is computed conditionally to the realization of s_N^*, is always similar to or lower than that of s_N^*/N, and excessively large r_M^* seldom occurs. Hence, these “positive" validation outcomes from Bernoulli trials also take effect in reducing the conservatism of ϵ(s_N^*, r_M^*).





Finally, we consider the practical case where independent validation scenarios are successively accumulated after x_N^* is obtained. For clarity, a posteriori bounds based on 20 incrementally emerging scenarios have been depicted in Fig. <ref>. Notice that tendencies of both ϵ_N,M(s_N^*,r_M^*) and η_M(r_M^*) are consistent with constraint violation results, thereby confirming the correctness of Theorem <ref>. Another observation is that the updating of ϵ_N,M(s_N^*,r_M^*) starts from ϵ_N,0(s_N^*,0) = 0.1176, while the evolution of η_M(r_M^*) has to start from M=1. Therefore, at the beginning stage, using ϵ_N,M(s_N^*,r_M^*) tends to be much more viable than using η_M(r_M^*). This fundamentally owes to the usage of an increased level of information in s_N^*. It also justifies that, to reach a reliable certificate with η_M(r_M^*) abundant samples are entailed <cit.>. In this sense, the proposed method enjoys desirable applicability in practical situations, especially when validation scenarios arrive in an incremental manner, or M is not too large.




§ CONCLUDING REMARKS

In this work, we proposed a general class of a posteriori probabilistic bounds for convex scenario programs with violation probability of scenario-based solutions assessed on additional validation scenarios. The a posteriori bounds turned out to be a function of both the number of support constraints and violation frequency on validation datasets. It can assess the feasibility of scenario-based solutions with randomness overall considered, which arises from both sampling N scenarios for formulating the problem and sampling M scenarios for validation. Thanks to the utilization of more information, refined evaluation of solution's risk can be made by the proposed bounds. It has been shown the established result contains the existing bounds in scenario optimization as special cases, thereby enjoying wide generality. For practical use, we developed a class of a posteriori violation probabilities under prespecified confidence levels, which are shown to possess a variety of desirable properties allowing for easy implementations and clear interpretations. Comprehensive case studies demonstrated the efficacy of the proposed a posteriori bounds in judging robustness of solutions to convex scenario programs.

Finally, it is noteworthy that similar to the wait-and-judge optimization <cit.>, the proposed probabilistic guarantee does not imply a conditional validity either, i.e., ℙ^N+M{ V(x_N^*) > ϵ(k,l) | s_N^* = k, r_M^* = l }≤β^*. The recent work <cit.> provides an interesting idea for attaining such a conditional validity, which is worth further investigating. Another promising direction is to generalize these results in a repetitive design scheme to achieve a “prescribed" risk level by deliberately choosing N and M.






















§ ACKNOWLEDGMENTS

The first author, C. Shang, was supported by National Science and Technology Innovation 2030 Major Project (No. 2018AAA0101604) of the Ministry of Science and Technology of China, and National Natural Science Foundation of China (Nos. 61873142, 61673236). F. You acknowledges financial support from the National Science Foundation (NSF) CAREER Award (CBET-1643244). The authors would like to thank the editor and anonymous reviewers for their valuable and constructive comments, which helped significantly improve this paper.





10

url@samestyle








zhou1996robust
K. Zhou, J. C. Doyle, K. Glover et al., Robust and Optimal
  Control.1em plus 0.5em minus 0.4emPrentice hall New Jersey,
  1996, vol. 40.

tempo2012randomized
R. Tempo, G. Calafiore, and F. Dabbene, Randomized Algorithms for
  Analysis and Control of Uncertain Systems: With Applications.1em plus
  0.5em minus 0.4emSpringer Science & Business Media, 2012.

ding2008model
S. X. Ding, Model-Based Fault Diagnosis Techniques: Design Schemes,
  Algorithms, and Tools.1em plus 0.5em minus 0.4emSpringer
  Science & Business Media, 2008.

van2020data
H. J. Van Waarde, J. Eising, H. L. Trentelman, and M. K. Camlibel, “Data
  informativity: A new perspective on data-driven analysis and control,”
  IEEE Transactions on Automatic Control, 2020.

calafiore2005uncertain
G. Calafiore and M. C. Campi, “Uncertain convex programs: Randomized solutions
  and confidence levels,” Mathematical Programming, vol. 102, no. 1,
  pp. 25–46, 2005.

campi2008exact
M. C. Campi and S. Garatti, “The exact feasibility of randomized solutions of
  uncertain convex programs,” SIAM Journal on Optimization, vol. 19,
  no. 3, pp. 1211–1230, 2008.

wang2012chance
Q. Wang, Y. Guan, and J. Wang, “A chance-constrained two-stage stochastic
  program for unit commitment with uncertain wind power output,” IEEE
  Transactions on Power Systems, vol. 27, no. 1, pp. 206–215, 2012.

deng2016decomposition
Y. Deng and S. Shen, “Decomposition algorithms for optimizing multi-server
  appointment scheduling with chance constraints,” Mathematical
  Programming, vol. 157, no. 1, pp. 245–276, 2016.

santoso2005stochastic
T. Santoso, S. Ahmed, M. Goetschalckx, and A. Shapiro, “A stochastic
  programming approach for supply chain network design under uncertainty,”
  European Journal of Operational Research, vol. 167, no. 1, pp.
  96–115, 2005.

you2013multicut
F. You and I. E. Grossmann, “Multicut Benders decomposition algorithm for
  process supply chain planning under uncertainty,” Annals of Operations
  Research, vol. 210, no. 1, pp. 191–211, 2013.

calafiore2006scenario
G. C. Calafiore and M. C. Campi, “The scenario approach to robust control
  design,” IEEE Transactions on Automatic Control, vol. 51, no. 5, pp.
  742–753, 2006.

calafiore2011research
G. C. Calafiore, F. Dabbene, and R. Tempo, “Research on probabilistic methods
  for control system design,” Automatica, vol. 47, no. 7, pp.
  1279–1293, 2011.

luedtke2008sample
J. Luedtke and S. Ahmed, “A sample approximation approach for optimization
  with probabilistic constraints,” SIAM Journal on Optimization,
  vol. 19, no. 2, pp. 674–699, 2008.

campi2011sampling
M. C. Campi and S. Garatti, “A sampling-and-discarding approach to
  chance-constrained optimization: Feasibility and optimality,” Journal
  of Optimization Theory and Applications, vol. 148, no. 2, pp. 257–280,
  2011.

care2015scenario
A. Carè, S. Garatti, and M. C. Campi, “Scenario min-max optimization and the
  risk of empirical costs,” SIAM Journal on Optimization, vol. 25,
  no. 4, pp. 2061–2080, 2015.

oishi2007polynomial
Y. Oishi, “Polynomial-time algorithms for probabilistic solutions of
  parameter-dependent linear matrix inequalities,” Automatica, vol. 43,
  no. 3, pp. 538–545, 2007.

dabbene2010randomized
F. Dabbene, P. S. Shcherbakov, and B. T. Polyak, “A randomized cutting plane
  method with probabilistic geometric convergence,” SIAM Journal on
  Optimization, vol. 20, no. 6, pp. 3185–3207, 2010.

care2014fast
A. Carè, S. Garatti, and M. C. Campi, “FAST–fast algorithm for the
  scenario technique,” Operations Research, vol. 62, no. 3, pp.
  662–671, 2014.

alamo2015randomized
T. Alamo, R. Tempo, A. Luque, and D. R. Ramirez, “Randomized methods for
  design of uncertain systems: Sample complexity and sequential algorithms,”
  Automatica, vol. 52, pp. 160–172, 2015.

chamanbaz2016sequential
M. Chamanbaz, F. Dabbene, R. Tempo, V. Venkataramanan, and Q.-G. Wang,
  “Sequential randomized algorithms for convex optimization in the presence of
  uncertainty,” IEEE Transactions on Automatic Control, vol. 61, no. 9,
  pp. 2565–2571, 2016.

calafiore2017repetitive
G. C. Calafiore, “Repetitive scenario design,” IEEE Transactions on
  Automatic Control, vol. 62, no. 3, pp. 1125–1137, 2017.

campi2018wait
M. C. Campi and S. Garatti, “Wait-and-judge scenario optimization,”
  Mathematical Programming, vol. 167, no. 1, pp. 155–189, 2018.

care2019wait
A. Carè, S. Garatti, and M. C. Campi, “The wait-and-judge scenario
  approach applied to antenna array design,” Computational Management
  Science, pp. 1–19, 2019.

campi2018general
M. C. Campi, S. Garatti, and F. A. Ramponi, “A general scenario theory for
  non-convex optimization and decision making,” IEEE Transactions on
  Automatic Control, vol. 63, no. 12, pp. 4067–4078, 2018.

garatti2019risk
S. Garatti and M. C. Campi, “Risk and complexity in scenario optimization,”
  Mathematical Programming, pp. 1–37, 2019.

you2018distributed
K. You, R. Tempo, and P. Xie, “Distributed algorithms for robust convex
  optimization via the scenario approach,” IEEE Transactions on
  Automatic Control, vol. 64, no. 3, pp. 880–895, 2018.

calafiore2010random
G. C. Calafiore, “Random convex programs,” SIAM Journal on
  Optimization, vol. 20, no. 6, pp. 3427–3464, 2010.

schildbach2013randomized
G. Schildbach, L. Fagiano, and M. Morari, “Randomized solutions to convex
  programs with multiple chance constraints,” SIAM Journal on
  Optimization, vol. 23, no. 4, pp. 2479–2501, 2013.

zhang2015sample
X. Zhang, S. Grammatico, G. Schildbach, P. Goulart, and J. Lygeros, “On the
  sample size of random convex programs with structured dependence on the
  uncertainty,” Automatica, vol. 60, pp. 182–188, 2015.

chernoff1952measure
H. Chernoff, “A measure of asymptotic efficiency for tests of a hypothesis
  based on the sum of observations,” The Annals of Mathematical
  Statistics, vol. 23, no. 4, pp. 493–507, 1952.

clopper1934use
C. J. Clopper and E. S. Pearson, “The use of confidence or fiducial limits
  illustrated in the case of the binomial,” Biometrika, vol. 26, no. 4,
  pp. 404–413, 1934.

conover1980practical
W. J. Conover, Practical Nonparametric Statistics.1em plus 0.5em
  minus 0.4emWiley New York, 1980.

thulin2014cost
M. Thulin, “The cost of using exact confidence intervals for a binomial
  proportion,” Electronic Journal of Statistics, vol. 8, no. 1, pp.
  817–840, 2014.

anderson1987linear
E. J. Anderson and P. Nash, Linear Programming in Infinite-Dimensional
  Spaces: Theory and Applications.1em plus 0.5em minus 0.4emJohn
  Wiley & Sons, 1987.

calafiore2006probabilistic
G. Calafiore and F. Dabbene, Probabilistic and Randomized Methods for
  Design under Uncertainty.1em plus 0.5em minus 0.4emSpringer,
  2006.

levine2018probabilistic
W. S. Levine, “Probabilistic and randomized tools for control design,” in
  The Control Systems Handbook.1em plus 0.5em minus 0.4emCRC Press, 2018, pp. 1497–1520.

grant2008cvx
M. Grant, S. Boyd, and Y. Ye, “CVX: Matlab software for disciplined convex
  programming,” 2008.

welsh2011robust
J. S. Welsh and H. Kong, “Robust experiment design through randomisation with
  chance constraints,” IFAC Proceedings Volumes, vol. 44, no. 1, pp.
  13 197–13 202, 2011.

garatti2019learning
S. Garatti and M. C. Campi, “Learning for control: A Bayesian scenario
  approach,” in 2019 IEEE 58th Conference on Decision and Control
  (CDC).1em plus 0.5em minus 0.4emIEEE, 2019, pp. 1772–1777.



