A Posteriori Probabilistic Bounds of Convex Scenario Programs with Validation Tests

TL;DR

This paper introduces new probabilistic bounds for convex scenario programs that integrate support constraints and validation data, enabling more reliable risk assessment with fewer samples.

Contribution

It proposes a general framework for a posteriori bounds that incorporate validation tests and support constraints, improving risk evaluation in scenario optimization.

Findings

Bounds depend on support constraints and validation data

Framework encompasses many existing theoretical results

Case studies demonstrate improved risk assessment

Abstract

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.

Tables

Algorithm 1 Bisection Method for Computing A Posteriori Bounds

Inputs: Integers $\zeta$, $N$ and $M$, confidence level $\beta$, and coefficients $\left\{a_m\right\}$.

1:  Initialization: Set numerical accuracy ${\epsilon }_{\mathrm{tol}}$, and $t(k,M+1)\equiv 0,\forall k$.

2:  for $k=0,1,\cdots ,\zeta$

3:    for $l=M,M-1\cdots ,0$

4:      $LB=t(k,l+1)$, $UB=1$;

5:      while $UB-LB\ge {\epsilon }_{\mathrm{tol}}$ do

6:        $t^{\mathrm{new}}=\left(UB+LB\right)/2$;

7:        if ${\tilde{h}}_{N,M}(t^{\mathrm{new}};k,l)\ge 0$

8:          $LB=t^{\mathrm{new}}$;

9:        else

10:         $UB=t^{\mathrm{new}}$;

11:        end

12:      end

13:      $t(k,l)=\left(UB+LB\right)/2,ϵ(k,l)=1-t(k,l)$.

14:   end

15: end

16: Return $\left\{{ϵ}_{N,M}(k,l)\right\}$.

Algorithm 2 Refinement Procedure of A Posteriori Probability Bounds

Inputs: Integers $\zeta$, $N$ and $M$, confidence level $\beta$.

1: Initialization: Set coefficients $\left\{a_m^{\left(\mathrm{ite}\right)}\right\}$ and iteration counter $\mathrm{ite}=1$.

2: Do until convergence

3:   - Compute $\left\{t^{\left(\mathrm{ite}\right)}(k,l)\right\}$ based on $\left\{a_m^{\left(\mathrm{ite}\right)}\right\}$.

4:   - Solve (39) with $\left\{t^{\left(\mathrm{ite}\right)}(k,l)\right\}$ and obtain optimum $\left\{a_m^{\left(\mathrm{ite}+1\right)}\right\}$.

5:   - Update $\mathrm{ite}=\mathrm{ite}+1$.

6: end

7: Return $ϵ(k,l)=1-t^{\left(\mathrm{ite}\right)}(k,l),\forall k,\forall l$.

Table 3. TABLE I: Results of Different A Posteriori Bounds from 2000 Monte Carlo Runs ( β = 10 − 6 𝛽 superscript 10 6 \beta=10^{-6} )

$N$	$M$	$ϵ\left(s_N^{\ast }\right)$		$\eta \left(r_M^{\ast }\right)$		$ϵ(s_N^{\ast },r_M^{\ast })$
		${\mu }_{\mathrm{GAP}}$	${\sigma }_{\mathrm{GAP}}$	${\mu }_{\mathrm{GAP}}$	${\sigma }_{\mathrm{GAP}}$	${\mu }_{\mathrm{GAP}}$	${\sigma }_{\mathrm{GAP}}$
$100$	$100$	$18.18\%$	$2.74\%$	$16.92\%$	$3.63\%$	$10.25\%$	$1.83\%$
	$200$			$10.37\%$	$2.49\%$	$7.69\%$	$1.50\%$
	$500$			$5.67\%$	$1.49\%$	$4.89\%$	$1.07\%$
	$1000$			$3.67\%$	$1.03\%$	$3.36\%$	$0.80\%$
$200$	$100$	$9.82\%$	$1.50\%$	$15.32\%$	$2.89\%$	$6.74\%$	$1.14\%$
	$200$			$8.99\%$	$2.05\%$	$5.38\%$	$0.98\%$
	$500$			$4.61\%$	$1.24\%$	$3.58\%$	$0.77\%$
	$1000$			$2.90\%$	$0.83\%$	$2.51\%$	$0.60\%$
$500$	$100$	$4.13\%$	$0.62\%$	$14.00\%$	$2.06\%$	$3.40\%$	$0.53\%$
	$200$			$7.78\%$	$1.46\%$	$2.96\%$	$0.48\%$
	$500$			$3.72\%$	$0.86\%$	$2.22\%$	$0.40\%$
	$1000$			$2.23\%$	$0.57\%$	$1.65\%$	$0.34\%$

Table 4. TABLE II: Occurrences of s N ∗ superscript subscript 𝑠 𝑁 s_{N}^{*} and r M ∗ superscript subscript 𝑟 𝑀 r_{M}^{*} in 2000 Monte Carlo Runs ( M = 1000 𝑀 1000 M=1000 )

	$s_N^{\ast }$	2	3	4	5	6	7	8
	Frequency	$45.70\%$	$22.20\%$	$23.30\%$	$7.10\%$	$1.50\%$	$0.20\%$	$0.00\%$
$N=100$	$s_N^{\ast }/N$	$2.00\%$	$3.00\%$	$4.00\%$	$5.00\%$	$6.00\%$	$7.00\%$	$8.00\%$
	Conditional mean($r_M^{\ast }/M$)	$1.90\%$	$1.87\%$	$1.92\%$	$2.12\%$	$2.45\%$	$2.88\%$	-
	Frequency	$41.85\%$	$21.70\%$	$25.60\%$	$9.15\%$	$1.40\%$	$0.25\%$	$0.05\%$
$N=200$	$s_N^{\ast }/N$	$1.00\%$	$1.50\%$	$2.00\%$	$2.50\%$	$3.00\%$	$3.50\%$	$4.00\%$
	Conditional mean($r_M^{\ast }/M$)	$1.00\%$	$0.90\%$	$0.96\%$	$1.08\%$	$1.15\%$	$1.37\%$	$1.35\%$
	Frequency	$35.80\%$	$24.75\%$	$29.90\%$	$7.70\%$	$1.50\%$	$0.35\%$	$0.00\%$
$N=500$	$s_N^{\ast }/N$	$0.40\%$	$0.60\%$	$0.80\%$	$1.00\%$	$1.20\%$	$1.40\%$	$1.60\%$
	Conditional mean($r_M^{\ast }/M$)	$0.40\%$	$0.36\%$	$0.38\%$	$0.44\%$	$0.48\%$	$0.76\%$	-

Equations

$${\mathrm{SP}}_m[{\omega }_m]:\{\begin{array}{rl}\underset{x\in {\mathbb{R}}^d}{\min }& J(x)\\ s.t.& x\in \mathcal{X}\subseteq {\mathbb{R}}^d\\ & f(x,{\delta }^{(i)})\le 0,\forall i\in {\mathbb{N}}_{1:m}\end{array}$$

$${\mathrm{SP}}_0[{\omega }_0]:\{\begin{array}{rl}\underset{x\in {\mathbb{R}}^d}{\min }& J(x)\\ s.t.& x\in \mathcal{X}\subseteq {\mathbb{R}}^d\end{array}$$

$$V(x)=\mathbb{P}\left\{\delta \in \Delta :f(x,\delta )>0\right\},$$

$${\mathbb{P}}^m\left\{V(x_m^{\ast })>ϵ\right\}\le B_m(ϵ;\zeta -1),m>\zeta ,$$

$$r_M=\sum _{l=1}^M{1}_{\{\delta :f(x,\delta )>0\}}({\tilde{\delta }}^{(l)}).$$

$${\mathbb{P}}^M\left\{V(x)>{\rho }_M(r_M)\right\}\le {\beta }^{\ast },$$

$${\mathbb{P}}^M\left\{V(x)>{\eta }_M(r_M)\right\}\le {\beta }^{\ast },$$

$$\eta_{M}(l)=\left\{\begin{array}[]{l l}\min_{\eta}\{\eta:\ B_{M}(\eta;l)\leq\beta^{*}\},&{\rm if}\ l\neq M\\ 1,&{\rm if}\ l=M\end{array}\right.$$

$${\mathbb{P}}^{N+M}\left\{V(x_N^{\ast })>ϵ(s_N^{\ast },r_M^{\ast })\right\}\le {\beta }^{\ast },$$

$$\begin{array}{rl}& \underset{g\in {\mathbb{C}}^{\zeta }[0,1]}{\inf }g(1)\\ & s.t.\left(\genfrac{}{}{0ex}{}{N}{k}\right)t^{N-k}\cdot \sum _{l=0}^M\left(\genfrac{}{}{0ex}{}{M}{l}\right)t^{M-l}(1-t{)}^l\cdot 1_{[0,1-ϵ(k,l))}(t)\\ & \le \frac{1}{k!}\frac{{\mathrm{d}}^k}{\mathrm{d}{t}^k}g(t),\forall t\in [0,1],k\in {\mathbb{N}}_{0:\zeta }\end{array}$$

$$F_k(v)={\mathbb{P}}^k\left\{V(x_k^{\ast })\le v\wedge s_k^{\ast }=k\right\},k\in {\mathbb{N}}_{0:\zeta }.$$

$$\left\{s_N^{\ast }=k\wedge r_M^{\ast }=l\right\}\bigcap \left\{s_N^{\ast }=k^{\prime }\wedge r_M^{\ast }=l^{\prime }\right\}=\varnothing ,$$

$$\begin{array}{rl}& {\mathbb{P}}^{N+M}\left\{V(x_N^{\ast })>ϵ(s_N^{\ast },r_M^{\ast })\right\}\\ =& \sum _{k=0}^{\zeta }\sum _{l=0}^M{\mathbb{P}}^{N+M}\left\{V(x_N^{\ast })>ϵ(k,l)\wedge s_N^{\ast }=k\wedge r_M^{\ast }=l\right\}.\end{array}$$

$${\mathcal{V}}_{k,l}=\left\{V(x_N^{\ast })>ϵ(k,l)\wedge s_N^{\ast }=k\wedge r_M^{\ast }=l\right\}\subseteq {\Delta }^{N+M}.$$

$${\mathbb{P}}^{N+M}\left\{V(x_N^{\ast })>ϵ(s_N^{\ast },r_M^{\ast })\right\}=\sum _{k=0}^{\zeta }\sum _{l=0}^M{\mathbb{P}}^{N+M}\left\{{\mathcal{V}}_{k,l}\right\}$$

$${\mathbb{P}}^{N+M}\{{\mathcal{V}}_{k,l}\}=\left(\genfrac{}{}{0ex}{}{N}{k}\right){\mathbb{P}}^{N+M}\{{\mathcal{A}}_1\cap {\mathcal{A}}_2\}.$$

$$\begin{array}{rl}{\mathcal{A}}_1=& \left\{V(x_N^{\ast })>ϵ(k,l)\wedge s_N^{\ast }=k\wedge I_N^{\ast }={\mathbb{N}}_{1:k}\right\},\end{array}$$

$${\mathcal{A}}_2=\left\{r_M^{\ast }=l\right\}.$$

$$\begin{array}{rl}& {\mathcal{A}}_1={\mathcal{B}}_1\triangleq \{V(x_k^{\ast })>ϵ(k,l)\wedge s_k^{\ast }=k\\ & \wedge f(x_k^{\ast },{\delta }^{(j)})\le 0,j\in {\mathbb{N}}_{k+1:N}\},\end{array}$$

$$\begin{array}{rl}& {\mathbb{P}}^{N+M}\{{\mathcal{A}}_1\cap {\mathcal{A}}_2\}\\ =& {\mathbb{P}}^{N+M}\{{\mathcal{B}}_1\cap {\mathcal{A}}_2\}\\ =& {\mathbb{P}}^{N+M}\left\{\begin{array}{rl}& V(x_k^{\ast })>ϵ(k,l)\wedge s_k^{\ast }=k\wedge r_M^{\ast }=l\\ & \wedge f(x_k^{\ast },{\delta }^{(j)})\le 0,j\in {\mathbb{N}}_{k+1:N}\end{array}\right\}\\ =& {\int }_{(ϵ(k,l),1]}\left(\genfrac{}{}{0ex}{}{M}{l}\right)(1-v{)}^{N-k+M-l}{v}^l\mathrm{d}{F}_k(v)\end{array}$$

$$\begin{array}{rl}& {\mathbb{P}}^{N+M}\left\{V(x_N^{\ast })>ϵ(s_N^{\ast },r_M^{\ast })\right\}\\ =& \sum _{k=0}^{\zeta }\sum _{l=0}^M\left(\genfrac{}{}{0ex}{}{N}{k}\right)\left(\genfrac{}{}{0ex}{}{M}{l}\right){\int }_{(ϵ(k,l),1]}(1-v{)}^{N-k+M-l}{v}^l\mathrm{d}{F}_k(v).\end{array}$$

$$\sum _{k=0}^{\min \{m,\zeta \}}\left(\genfrac{}{}{0ex}{}{m}{k}\right){\int }_{[0,1]}(1-v{)}^{m-k}\mathrm{d}{F}_k(v)=1,m\in {\mathbb{N}}_0.$$

$${\mathbb{P}}^{N+M}\left\{V(x_N^{\ast })>ϵ(s_N^{\ast },r_M^{\ast })\right\}\le {\beta }^{\prime },$$

$$\begin{array}{rl}\underset{\{F_k(\cdot )\}}{\sup }& \sum _{k=0}^{\zeta }\sum _{l=0}^M\left(\genfrac{}{}{0ex}{}{N}{k}\right)\left(\genfrac{}{}{0ex}{}{M}{l}\right){\int }_{(ϵ(k,l),1]}(1-v{)}^{N-k+M-l}{v}^l\mathrm{d}{F}_k(v)\\ s.t.& \sum _{k=0}^{\min \{m,\zeta \}}\left(\genfrac{}{}{0ex}{}{m}{k}\right){\int }_{[0,1]}(1-v{)}^{m-k}\mathrm{d}{F}_k(v)=1,m\in {\mathbb{N}}_0\\ & F_0,F_1,\cdots ,F_{\zeta }\in \mathcal{C}\end{array}$$

$$\begin{array}{rl}\underset{\{F_k(\cdot )\}}{\sup }& \sum _{k=0}^{\zeta }\sum _{l=0}^M\left(\genfrac{}{}{0ex}{}{N}{k}\right)\left(\genfrac{}{}{0ex}{}{M}{l}\right){\int }_{(ϵ(k,l),1]}(1-v{)}^{N-k+M-l}{v}^l\mathrm{d}{F}_k(v)\\ s.t.& \sum _{k=0}^{\min \{m,\zeta \}}\left(\genfrac{}{}{0ex}{}{m}{k}\right){\int }_{[0,1]}(1-v{)}^{m-k}\mathrm{d}{F}_k(v)=1,m\in {\mathbb{N}}_{0:K}\\ & F_0,F_1,\cdots ,F_{\zeta }\in \mathcal{C}\end{array}$$

$$\begin{array}{rl}\underset{{\lambda }_0,\cdots ,{\lambda }_K}{\inf }& \sum _{m=0}^K{\lambda }_m\\ s.t.& \left(\genfrac{}{}{0ex}{}{N}{k}\right)\sum _{l=0}^M\left(\genfrac{}{}{0ex}{}{M}{l}\right)(1-v{)}^{N-k+M-l}{v}^l{1}_{(ϵ(k,l),1]}(v)\\ & \le \sum _{m=k}^K{\lambda }_m\left(\genfrac{}{}{0ex}{}{m}{k}\right)(1-v{)}^{m-k},\forall t\in [0,1],k\in {\mathbb{N}}_{0:\zeta }\end{array}$$

$$g(1)=\sum _{m=0}^K{\lambda }_m,\frac{1}{k!}\frac{{\mathrm{d}}^k}{\mathrm{d}{t}^k}{t}^m=\{\begin{array}{rl}& 0,m<k\\ & \left(\genfrac{}{}{0ex}{}{m}{k}\right)t^{m-k},m\ge k\end{array}$$

$$\begin{array}{rl}\underset{g(\cdot )\in {\mathbb{P}}_K}{\inf }& g(1)\\ s.t.& \frac{1}{k!}\frac{{\mathrm{d}}^k}{\mathrm{d}{t}^k}g(t)\ge \left(\genfrac{}{}{0ex}{}{N}{k}\right)\sum _{l=0}^M\left(\genfrac{}{}{0ex}{}{M}{l}\right)t^{N-k+M-l}(1-t{)}^l\\ & \cdot 1_{[0,1-ϵ(k,l))}(t),\forall t\in [0,1],k\in {\mathbb{N}}_{0:\zeta }\end{array}$$

$${\mathbb{P}}^{N+M}\left\{V(x_N^{\ast })>ϵ(s_N^{\ast },r_M^{\ast })\right\}\le {\beta }^{\prime }\le \underset{K}{\inf }{\beta }_K^{\ast }={\beta }^{\ast }.$$

$${\mathbb{P}}^N\left\{V(x_N^{\ast })>ϵ(s_N^{\ast })\right\}\le {\beta }^{\ast },$$

$$\begin{array}{rl}& \underset{g\in {\mathbb{C}}^{\zeta }[0,1]}{\inf }g(1)\\ & s.t.\left(\genfrac{}{}{0ex}{}{N}{k}\right)t^{N-k}\cdot 1_{[0,1-ϵ(k))}(t)\le \frac{1}{k!}\frac{{\mathrm{d}}^k}{\mathrm{d}{t}^k}g(t),k\in {\mathbb{N}}_{0:\zeta }\end{array}$$

Full text

I 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 [1, 2]. Likewise, in filter design for fault detection, such effects have been taken into account in order to avoid high false alarm rates [3].

The increasing availability of data has motivated rapid development of data-driven methods in systems and control [4]. As a representative, the scenario approach has been widely adopted as a data-driven approach for reliable decision-making under uncertainty [5, 6], 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 [7], service systems [8], supply chain management [9, 10], and control design [11, 12]. 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 [5, 6]. To trade robustness for performance, extensions based on constraint discarding have been also developed [13, 14, 15].

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 [16, 17], which solves a sequence of reduced-size scenario programs that are more computationally affordable [18, 19, 20, 21]. 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 [22, 23, 24, 25]. 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 [25]. 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 [22, 23]. 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 [26]. 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 II, the technical background of scenario programs is introduced. Section III formally states the main results of this work, and discusses some practical implementation issues. An illustrative example is provided in Section IV, followed by concluding remarks in the last section.

Notations and Definitions. $\mathbb{N}_{0}$ is the set of non-negative integer numbers, and the set of consecutive non-negative integers $\{j,\cdots,k\}$ is denoted by $\mathbb{N}_{j:k}$. The $p$-norm of a matrix is denoted by $\|\cdot\|_{p}$, while $\|\cdot\|$ represents the Euclidean norm of a vector by convention. The trace operator on a square matrix is denoted by ${\rm Tr}\{\cdot\}$. The indicator function of a subset $\mathcal{S}$ is defined as $\mathbb{1}_{\mathcal{S}}(s)$. Given a positive integer $N$, a non-negative integer $m\leq N$, and $t\in(0,1)$, the Beta cumulative probability function is defined as $B_{N}(t;m):=\sum_{i=0}^{m}\binom{N}{i}t^{i}(1-t)^{N-i}$. $\mathbb{C}^{q}[0,1]$ denotes the class of $q$ times differentiable functions with continuous $q$-th derivative over the interval $[0,1]$, and $\mathbb{P}^{K}$ denotes the set of polynomials of degree $K$.

II Preliminaries

II-A Scenario Programs

Assume that $\Delta$ is a probability space, which is associated with a $\sigma$-algebra $\mathcal{F}$ and a probability measure $\mathbb{P}$. $\delta$ is a random outcome from the triplet $(\Delta,\mathcal{F},\mathbb{P})$. Moreover, denote by $(\Delta^{m},\mathcal{F}^{m},\mathbb{P}^{m})$ the Cartesian product of $\Delta$ equipped with the $\sigma$-algebra $\mathcal{F}^{m}$ and the $m$-fold probability measure $\mathbb{P}^{m}\doteq\mathbb{P}\times\cdots\times\mathbb{P}$ ($m$ times). In general, the analytic form of $\mathbb{P}$ is unknown in practice, but we are allowed to sample a set of independent scenarios $\{\delta^{(1)},\cdots,\delta^{(m)}\}$ from $\mathbb{P}^{m}$. Here, we refer to $\omega_{m}:=\{\delta^{(1)},\cdots,\delta^{(m)}\}\in\Delta^{m}$ as a multi-sample, based on which the following scenario programs can be defined:

[TABLE]

when $m>0$, and

[TABLE]

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

[TABLE]

which can also be interpreted as the risk of $\mathbb{x}$, and the corresponding reliability level is $1-V(\mathbb{x})$ [2].

Meanwhile, the following assumption is made throughout the paper, which is fairly standard in generic convex scenario programs [6, 11].

Assumption 1** (Feasibility and Uniqueness).**

For every $m\in\mathbb{N}_{0}$ and every multi-sample $\omega_{m}\in\Delta^{m}$, the optimal solution $\mathbb{x}^{*}_{m}$ to the scenario program ${\rm SP}_{m}[\omega_{m}]$ exists and is unique.

Note that in order to secure the uniqueness of $\mathbb{x}_{m}^{*}$, a suitable tie-break rule can be employed [5]. In this way, the optimal solution $\mathbb{x}_{m}^{*}$ to ${\rm SP}_{m}[\omega_{m}]$ is essentially a random variable defined over $\Delta^{m}$. Under Assumption 1 and the conventional assumption on the measurability of $V(\mathbb{x}_{m}^{*})$, a celebrated probabilistic guarantee in literature is that, the tail probability of $V(\mathbb{x}_{m}^{*})$ is dominated by the Beta cumulative probability function [6]:

[TABLE]

where $\zeta(\leq d)$ can be the Helly’s dimension [27] or its tighter bound built with structural information [28, 29]. The probabilistic guarantee (2) is closely related to the concept of support constraints, which will be made heavy use of in this paper.

Definition 1** (Support Constraint, [6]).**

A constraint of the scenario program ${\rm SP}_{m}[\omega_{m}]$ is a support constraint if its removal yields a different optimal solution of the initial problem.

Upon solving ${\rm SP}_{m}[\omega_{m}]$, the indices of induced support constraints can be described by a set $I_{m}^{*}\subseteq\mathbb{N}_{0}$, which is a set-valued mapping from the multi-sample $\omega_{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 $\omega_{m}$; hence, both $I_{m}^{*}$ and $s_{m}^{*}$ are essentially random variables on the triplet $(\Delta^{m},\mathcal{F}^{m},\mathbb{P}^{m})$. The definition of support constraints differs from that of generic active constraints, which are characterized by $\{i:\ f(\mathbb{x}_{m}^{*},\delta^{(i)})=0\}$. For convex scenario programs, support constraints are always active constraints, but the converse no longer remains true [22]. Meanwhile, under the availability of sufficient scenarios, the number of support constraints is always upper bounded by $\zeta$ [5]. For a fully-supported problem where the number of support constraints is exactly $\zeta$ with probability one (w.p.1), the inequality (2) then becomes an equality [6]; in other cases, it inevitably induces conservatism.

Another useful assumption is made as follows.

Assumption 2** (Non-Degeneracy, [22]).**

The solution to the scenario program ${\rm SP}_{m}[\omega_{m}]$ coincides w.p.1 with the solution to the program defined by support constraints only.

It is noteworthy that Assumption 2 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.

II-B Confidence Intervals for Violation Probability with Validation Tests

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

Definition 2** (Empirical Violation Frequency).**

Given a set of $M$ independent scenarios $\{\tilde{\delta}^{(1)},\cdots,\tilde{\delta}^{(M)}\}$, the empirical violation frequency of a solution $\mathbb{x}$ is calculated as:

[TABLE]

That is, the number of scenario-based validation constraints that are violated by $\mathbb{x}$.

With $\mathbb{x}$ fixed, $r_{M}$ is a random outcome on $\mathbb{P}^{M}$, which embodies useful information for estimating the binomial proportion, i.e., the true violation probability $V(\mathbb{x})$. In the context of robust design, a small violation probability is always desired, and hence an upper-bound for $V(\mathbb{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 [2].

Theorem 1** (One-Sided Chernoff Bound, [30]).**

Given a solution $\mathbb{x}$ and a pre-defined confidence level $\beta^{*}\in(0,1)$, it then holds that:

[TABLE]

where the a posteriori violation probability $\rho_{M}(r_{M})=r_{M}/M+\sqrt{\log\beta^{*}/(-2M)}$ is a random variable depending on the realization of $r_{M}$.

Theorem 2** (One-Sided Clopper-Pearson (C-P) Bound, [31, 32, 33]).**

Given a solution $\mathbb{x}$ and a pre-defined confidence level $\beta^{*}\in(0,1)$, it then holds that:

[TABLE]

where the analytic form of $\eta_{M}(\cdot)$ is expressed as:

[TABLE]

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

III Main Results

III-A A Posteriori Confidence Level Based on Violation Probability Function

Assume that we are seated in a situation where a scenario program ${\rm SP}_{N}[\omega_{N}]$ with fixed $N>\zeta$ ($N$ is the initial number of available scenarios for optimization), is first solved to pursue robustness, and then the risk of $\mathbb{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 $\mathbb{x}$, we can certainly use them to evaluate $\mathbb{x}^{*}_{N}$’s risk as a function of $r_{M}^{*}$. However, the optimal solution $\mathbb{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 $\epsilon(s_{N}^{*},r_{M}^{*})$ on the violation probability $V(\mathbb{x}^{*}_{N})$ is considered as a function of both $s_{N}^{*}$ and $r_{M}^{*}$. With a group of functions $\{\epsilon(k,l)\}$ pre-specified for $k\in\mathbb{N}_{0:\zeta},\ l\in\mathbb{N}_{0:M}$, the following finite-sample probabilistic guarantee can be established as one of our main results.

Theorem 3**.**

Assume that $\{\delta^{(1)},\cdots,\delta^{(N)},\tilde{\delta}^{(1)},\cdots,\tilde{\delta}^{(M)}\}$ are sampled from the probability space $(\Delta^{N+M},\mathcal{F}^{N+M},\mathbb{P}^{N+M})$, and $\epsilon(k,l)$ is an arbitrary $[0,1]$-valued function where $k\in\mathbb{N}_{0:\zeta},\ l\in\mathbb{N}_{0:M}$. Then it holds that

[TABLE]

where $\beta^{*}$ is the optimal value of the following variational problem, defined based on functions $\{\epsilon(k,l)\}$:

[TABLE]

Proof.

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

[TABLE]

Next, we show that the a posteriori violation probability $\mathbb{P}^{N+M}\left\{V(\mathbb{x}_{N}^{*})>\epsilon(s_{N}^{*},r_{M}^{*})\right\}$ can be calculated based on $\{F_{k}(v)\}_{k=0}^{\zeta}$, which act as the backbone of the whole machinery. Notice that the following events do not overlap with each other:

[TABLE]

if $k\neq k^{\prime}$ or $l\neq l^{\prime}$. Then the following decomposition can be made:

[TABLE]

We define the events on the probability space $\Delta^{N+M}$:

[TABLE]

which indicates that the violation probability of $\mathbb{x}_{N}^{*}$ is larger than $\epsilon(k,l)$, and there are $k$ support constraints and $l$ violated validation constraints, respectively. In this way,

[TABLE]

holds. Next we derive analytic expressions for $\mathbb{P}^{N+M}\left\{\mathcal{V}_{k,l}\right\}$. By fixing the indices of support constraints, we have:

[TABLE]

where

[TABLE]

It has been proved by [22] that, in the probability space $(\Delta^{N},\mathcal{F}^{N},\mathbb{P}^{N})$, the following events are equal up to a zero probability set:

[TABLE]

which implies that:

[TABLE]

where the last equality is due to the fact that by fixing violation probability $v$, $(1-v)^{N-k}$ is the probability that $\mathcal{X}_{\delta^{(k+1)}},\cdots,\mathcal{X}_{\delta^{(N)}}$ are satisfied by $\mathbb{x}_{k}^{*}$, and $\binom{M}{l}v^{l}(1-v)^{M-l}$ is the probability of observing exactly $l$ violated validation constraints. By using the definition of $F_{k}(v)$ in (9) and taking integral over the interval $(\epsilon(k,l),1]$, (12) can be attained. Substituting (11) and (12) into (10) yields:

[TABLE]

Next we show that the RHS of (13) is actually upper-bounded by $\beta^{*}$ given in (8). The idea is to take generalized distribution functions $\{F_{k}(v)\}_{k=0}^{\zeta}$ as decision variables subject to the following joint moment conditions [22]:

[TABLE]

Therefore, we have:

[TABLE]

where $\beta^{\prime}$ is the optimal value of the following variational problem:

[TABLE]

Here, $\mathcal{C}$ stands for the set of all positive generalized distribution functions. However, (14) 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:

[TABLE]

whose optimal value is $\beta^{\prime}_{K}$. Because only $K$ constraints are involved in (15), (15) is a relaxation of (14), so we have $\beta^{\prime}_{K}\geq\beta^{\prime},\ \forall K$. Meanwhile, the value of $\beta^{\prime}_{K}$ is non-increasing with $K$. By deriving the Lagrangian function and optimizing over $\{F_{k}(\cdot)\}_{k=0}^{\zeta}$, we obtain the following dual problem [34]:

[TABLE]

where $\lambda_{0},\cdots,\lambda_{K}$ are dual variables and $\beta^{*}_{K}$ is the optimal value. Therefore, it holds that $\beta^{\prime}\leq\beta^{\prime}_{K}\leq\beta_{K}^{*},\ \forall K\in\mathbb{N}_{0}$ due to weak duality. We define the polynomial of degree $K$ as $g(t)=\sum_{m=0}^{K}\lambda_{m}t^{m}$, which admits the following properties:

[TABLE]

By defining $t:=1-v$ and plugging (17) into (16), (16) then becomes:

[TABLE]

Denoting by $\mathcal{D}$ the feasible region of problem (8), it has been proved in [22] that the set $\left(\cup_{K\geq\zeta}\mathbb{P}_{K}\right)\cap\mathcal{D}$ is dense in $\mathcal{D}$. Therefore we can conclude that:

[TABLE]

This completes the proof. ∎

Remark 1*.*

Theorem 3 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 $\mathbb{P}$. Note that the randomness arises from sampling of both $N$ scenarios and $M$ validation samples, and $s_{N}^{*}$ is a random variable defined on $\Delta^{N}$ only, whereas $r_{M}^{*}$ is defined on the $(N+M)$-fold product probability space $\Delta^{N+M}$ since it depends on both $\mathbb{x}_{N}^{*}$ and validation data. Hence, the finite-sample probabilistic guarantee (7) 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 $\{\epsilon(k,l)\}$ on $k$ and/or $l$ imposed, Theorem 3 turns out to unify many existing theoretical results, including the wait-and-judge approach in [22] as well as the generic probabilistic guarantee (2) in [6]. First, the following result reveals itself as an obvious corollary of Theorem 3, which enforces the two-indexed bounds $\{\epsilon(k,l)\}$ to depend on $k$ only.

Corollary 1** (Theorem 1 in [22]).**

Take $\epsilon(k,l)=\epsilon(k),\forall k\in\mathbb{N}_{0:\zeta},\ l\in\mathbb{N}_{0:M}$. Under Assumptions 1 and 2, it then holds that:

[TABLE]

where $\beta^{*}$ is the optimal value of the following problem:

[TABLE]

Proof.

A straightforward substitution of the identity

[TABLE]

into (8) yields the proof. ∎

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

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

Corollary 2**.**

By taking $\epsilon(k,l)\equiv\epsilon,\ \forall k\in\mathbb{N}_{0:\zeta},\ l\in\mathbb{N}_{0:M}$, and under Assumptions 1 and 2, one obtains:

[TABLE]

Proof.

See Corollary 1 in [22]. ∎

III-B Derivation of A Posteriori Probabilistic Bounds

Theorem 3 is shown to be a generalization of many important probabilistic bounds in scenario optimization, with particular options for $\{\epsilon(k,l)\}$ used. Hence, the freedom in choosing $\{\epsilon(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(\mathbb{x}_{N}^{*})\leq\epsilon(s_{N}^{*},r_{M}^{*})$ occurs with a suitably high confidence, and in order to explicitly control the probability of a wrong declaration $V(\mathbb{x}_{N}^{*})>\epsilon(s_{N}^{*},r_{M}^{*})$, it is desirable to first assign $\beta^{*}$, which is typically small, and then determine $\{\epsilon(k,l)\}$ in reverse. In principle, there are infinitely many options of $\{\epsilon(k,l)\}$ that lead to the same confidence level $\beta^{*}$. For practical use, a certain choice of $\{\epsilon(k,l)\}$ is suggested in the following theorem. It subsumes Theorem 2 in [22] as a special case with $M=0$ (no validation tests) and specific coefficients $\{a_{m}\equiv 1/(N+1),\ \forall m\in\mathbb{N}_{0:N}\}$, where the induced bound $\epsilon(k)$ depending only on $s_{N}^{*}$ equals $\epsilon(k,0)$. However, the present proof machinery has a significant departure from that in [22], 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.

Lemma 1**.**

For fixed values of $m$ and $N$ $(m<N)$, $B_{N}(t;m)$ is a strictly decreasing function in $t\in(0,1)$.

Proof.

See [19]. ∎

Theorem 4**.**

Define the following polynomial in $t$:

[TABLE]

where $\beta\in(0,1)$ is the confidence level, and coefficients $\{a_{m}\}$ satisfy

[TABLE]

(i) $h_{N,M}(t;k,l)$ has only one solution in the interval $(0,1)$ for all $k\in\mathbb{N}_{0:\zeta}$ and $l\in\mathbb{N}_{0:M}$.

(ii) Letting the root be $t_{N,M}(k,l)$ and $\epsilon_{N,M}(k,l)=1-t_{N,M}(k,l)$, it then holds that:

[TABLE]

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

[TABLE]

Proof.

(i) First we define the following polynomial:

[TABLE]

Since $k\leq\zeta<N$ and $\sum_{m=\zeta}^{N}a_{m}>0$, $\sum_{m=k}^{N-1}a_{m}\cdot\binom{m}{k}\cdot t^{m-N}$ is a strictly decreasing function in $t\in(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)\equiv 1$ when $l=M$. Therefore, in both cases, $\tilde{h}_{N,M}(t;k,l)$ must be strictly decreasing in $(0,1)$. Meanwhile, notice that

[TABLE]

Hence, $\tilde{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 $\tilde{h}_{N,M}(t;k,l)>\tilde{h}_{N,M}(t;k,l+1)$. Since $\tilde{h}_{N,M}(t;k,l)$ is strictly decreasing in $t\in(0,1)$, we then have (22).

(iii) We first parameterize a candidate solution to (8) as:

[TABLE]

Then we will show that $g(t)$ is essentially a feasible solution to (8) with $\epsilon(k,l)$ set as $\epsilon_{N,M}(k,l)$. Its $k$-th order derivative can be computed as:

[TABLE]

Because $\epsilon_{N,M}(k,l+1)>\epsilon_{N,M}(k,l)$, the entire interval $[0,1-\epsilon_{N,M}(k,0))$ can be split as

[TABLE]

where $\epsilon_{N,M}(k,M+1)=1$ is taken for granted for notational clarity. Considering the case where $t$ falls within a certain interval $[1-\epsilon_{N,M}(k,l+1),1-\epsilon_{N,M}(k,l))$, we have:

[TABLE]

The last inequality is due to the fact that $t_{N,M}(k,l)=1-\epsilon_{N,M}(k,l)$ is the root of $h_{N,M}(t;k,l)$, and for $t\in[0,1-\epsilon_{N,M}(k,l))$, $h_{N,M}(t;k,l)>0$ holds. Therefore, if we set $\epsilon(k,l)$ in (8) to be $\epsilon_{N,M}(k,l)$, constraints in (8) hold for all $t\in[0,1]$. As such, $g(t)$ is feasible for problem (8), which indicates $\beta^{*}\leq\beta$. Finally we have:

[TABLE]

∎

The following remarks are made in order.

Remark 2*.*

The utilization of a posteriori probabilistic certificates $\{\epsilon_{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.

Remark 3*.*

The utilization of a posteriori probabilistic certificates $\{\epsilon_{N,M}(k,l)\}$ should not be oriented towards a tradeoff design with a prescribed tolerance on $V(\mathbb{x}_{N}^{*})$. Rather, they shall be used for a posteriori risk judgement of a scenario-based solution. Of common interest, a scenario program ${\rm SP}_{N}[\omega_{N}]$ may be constructed as a relaxation of the infinite-dimensional constraint $f(\mathbb{x},\delta)\leq 0,\ \forall\delta\in\Delta$, or as an approximation of the chance constraint $\mathbb{P}\{f(\mathbb{x},\delta)>0\}\leq\epsilon$ with the sample size $N$ determined based on the prior bound (19). However, the computational effort for solving ${\rm SP}_{N}[\omega_{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, $\{\epsilon_{N,M}(k,l)\}$ show greater promise than the wait-and-judge bounds $\{\epsilon_{N}(k)\}$ in addressing the practical need for risk evaluation.

Remark 4*.*

The monotonicity property (22) 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. 1, where $\{\epsilon_{N,M}(k,l)\}$ with $\beta=10^{-6}$, $\zeta=10$, $N=50$, $M=30$, and $\{a_{m}\equiv 1/(N+1)\}$ are portrayed. Meanwhile, it is easy to prove that $\{\epsilon_{N,M}(k,l)\}$ are increasing in $k$, which preserves the interpretation of generic wait-and-judge approach [22] that less support constraints indicate a potentially lower risk of randomized solutions.

Note that bounds $\{\epsilon_{N,M}(k,l)\}$ and $\{t_{N,M}(k,l)\}$ derived in Theorem 4 depend on $N$ and $M$. For notational simplicity, we will use $\{\epsilon(k,l)\}$ and $\{t(k,l)\}$ in the sequel when no confusions are caused. Meanwhile, we use $\{\epsilon(k)\}$ to denote the a posteriori bounds in the wait-and-judge approach [22] that depends on $s_{N}^{*}$ only. Now we establish a result that reveals the value of utilizing more validation information in $\{\epsilon_{N,M}(k,l)\}$.

Corollary 3**.**

With $\{a_{m}\}$ and $\beta$ fixed, it then holds that:

[TABLE]

Proof.

The second equality can be obtained by substituting the identity $B_{M}(1-t;M)=1$ into (20), and the inequality stems from (22). ∎

Corollary 3 shows that the proposed bound $\epsilon(s_{N}^{*},r_{M}^{*})$ is always no higher than $\epsilon(s_{N}^{*})$ that depends on partial information within $N$ scenarios, as can be observed from Fig. 1. 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$, $\epsilon(k,l)\leq\epsilon(k)$ still holds true, which admits a sensible explanation. It has been revealed in [22, 25] 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 $\mathbb{P}^{N}\left\{V(\mathbb{x}_{N}^{*})>\epsilon(s_{N}^{*})\right\}\leq\beta$. As a consequence, seeing an excessively large $r_{M}^{*}$ is issued with a negligible probability, so the change in $\epsilon(k,l)$ with large $l$ does not exert significant influence on the overall confidence $\beta$. Similarly, as implied by [25], the joint distribution of $(s_{N}^{*},r_{M}^{*})$ is concentrated in the region where $V(\mathbb{x}_{N}^{*})=\mathbb{E}\{r_{M}^{*}\}/M$ and $s_{N}^{*}/N$ have the same order of magnitude. This is also the region where the actual improvement of $\epsilon(k,l)$ over $\epsilon(k)$ occurs, as also confirmed by Fig. 1. In this sense, $\epsilon(k,l)$ shall be viewed as a refinement of $\epsilon(k)$ based on the new information carried by $M$ validation samples.

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

Lemma 2**.**

For fixed values of $t$ and $m$, $B_{N}(t;m)$ is strictly decreasing in $N\ (N\geq m)$.

Proof.

For $i\in\mathbb{N}_{1:N}$, it is obvious that

[TABLE]

Therefore, we have:

[TABLE]

This completes the proof of Lemma 2. ∎

Then the following relations hold between $\{\eta_{M}(\cdot)\}$ and $\{\eta_{M+1}(\cdot)\}$.

Theorem 5**.**

It holds that:

[TABLE]

Proof.

(28) immediately follows from (6). Because $\eta_{M}(l)$ is the root of $B_{M}(\eta,l)=\beta^{*}$ in $(0,1)$, and $B_{M}(\eta,l)$ is strictly decreasing in $\eta$, it remains to show that:

[TABLE]

A straightforward usage of Lemma 2 results in (29), along with the inequality $B_{M+1}(1-\eta;M-l-1)<B_{M}(1-\eta;M-l-1),\ l\in\mathbb{N}_{0:M-1}$. Hence, $1-B_{M+1}(1-\eta;M-l-1)>1-B_{M}(1-\eta;M-l-1)$ holds, which amounts to (30). ∎

In view of Theorem 5, rational refinements of $\epsilon(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 $\mathbb{x}_{N}^{*}$ is obtained, which is precisely stated as $\eta_{M+1}(l+1)>\eta_{M}(l)$. Conversely, if $\mathbb{x}_{N}^{*}$ agrees with the new validation sample, the risk level should be discounted, i.e. $\eta_{M+1}(l)<\eta_{M}(l)$. Interestingly, such properties are also possessed by $\{\epsilon_{N,M}(k,l)\}$.

Theorem 6**.**

With values of $\{a_{m}\}$ and $\beta$ fixed, the roots of $\{g_{N,M}(t;k,l)\}$ satisfy

[TABLE]

thereby indicating the following relations:

[TABLE]

The proof is similar to that for Theorem 5, and is hence omitted for brevity.

Before closing this subsection, we deal with the computational issue. Although $\{\epsilon_{N,M}(k,l)\}$ have no analytic expressions, one only needs to perform bisection to numerically search the root of $\tilde{h}_{N,M}(t;k,l)$ in $(0,1)$ to compute $\{\epsilon_{N,M}(k,l)\}$, as implied by Theorem 4(i). The pseudo code for calculating $\{\epsilon_{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 (22) is utilized, and the enumeration of $l$ is made from $M$ to [math]. 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$.

III-C Fundamental limits and a refinement procedure

Above derivations imply that lower values of $\{\epsilon(k,l)\}$ give more informative certificates on the solution’s reliability. Hence a natural question comes up: given the confidence $\beta^{*}$, what is the best possible value of $\epsilon(k,l)$ satisfying (7) in Theorem 3? Next we establish lower limits for $\{\epsilon(k,l)\}$, where $\epsilon(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.

Theorem 7**.**

Consider all admissible functions $\{\epsilon(k,l)\}$ satisfying (i) the probabilistic guarantee (7) in Theorem 3 with a prescribed $\beta^{*}$ and (ii) the monotonicity property $\epsilon(k,l)<\epsilon(k,l+1),~{}l\in\mathbb{N}_{0:M-1}$. Then for $k\geq 1$ a lower limit $\underline{\epsilon}(k,l)$ of $\epsilon(k,l)$ is given by the root of the following equation in $(0,1)$:

[TABLE]

where

[TABLE]

And, $\underline{\epsilon}(0,l)=0,\ \forall l\in\mathbb{N}_{0:M}$.

Proof.

We show by contradiction that $\epsilon(k,l)<\underline{\epsilon}(k,l)$ is impossible for $k\geq 1$. Assume that there exist admissible functions $\{\epsilon(k,l)\}$ satisfying (i) and (ii), and there exist $k^{\prime}$ and $l^{\prime}$ such that $\epsilon(k^{\prime},l^{\prime})<\underline{\epsilon}(k^{\prime},l^{\prime})$. Then we consider a fully-supported problem ${\rm SP}^{\prime}_{N}[\omega_{N}]$ where $s_{N}^{*}=k^{\prime}\leq\zeta$ occurs w.p.1. When $k^{\prime}\geq 1$, the probability density function of $V(\mathbb{x}_{N}^{*})$ in ${\rm SP}^{\prime}_{N}[\omega_{N}]$ is explicitly expressed as [6]:

[TABLE]

Hence, for this specific fully-supported problem ${\rm SP}^{\prime}_{N}[\omega_{N}]$, it follows that:

[TABLE]

Recall that $\epsilon(k^{\prime},l)<\underline{\epsilon}(k^{\prime},l)$ has been assumed, and $\{\epsilon(k,l)\}$ satisfies the monotonicity property $\epsilon(k,l)<\epsilon(k,l+1),~{}l\in\mathbb{N}_{0:M-1}$. It immediately follows that $\epsilon(k^{\prime},j)<\underline{\epsilon}(k^{\prime},l),\ \forall j\in\mathbb{N}_{0:l}$, and one obtains:

[TABLE]

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 (34). However, (36) indicates that ${\rm SP}^{\prime}_{N}[\omega_{N}]$ contradicts the distribution-free nature of $\{\epsilon(k,l)\}$ in (i), i.e., the probabilistic guarantee (7) with a given $\beta^{*}$. ∎

Fig. 2 depicts a posteriori bounds $\{\epsilon(k,l)\}$ under $\{a_{m}\equiv 1/(N+1)\}$ and the fundamental lower limits $\{\underline{\epsilon}(k,l)\}$ ($\beta=10^{-6}$, $\zeta=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 $\{\epsilon(k,l)\}$ since even knowing $k\geq 1$ beforehand does not provide too much improvement in risk evaluation. Indeed, this well inherits the merits of the wait-and-judge scenario optimization [22, 25], which can be interpreted as that the realization of $s_{N}^{*}$ is rather informative for a posteriori risk evaluation [25].

Obviously, $\{\epsilon(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 $\{\epsilon(k,l)\}$ can be improved to approach the fundamental limits $\{\underline{\epsilon}(k,l)\}$ by flexibly adjusting $\{a_{m}\}$. This can be formally cast as the following multi-objective optimization problem:

[TABLE]

where $\tau>0$ is a small positive number to ensure (21). Clearly, (37) 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 (37). The following theorem provides a sufficient and necessary condition for making improvements over the current $\{a_{m}\}$.

Theorem 8**.**

Given $N$, $M$ and $\beta$. Assume that $\{a_{m}\}$ and $\{a^{\prime}_{m}\}$ satisfying (21) are different coefficients of polynomials $\tilde{h}_{N,M}(t;k,l)$ and $\tilde{h}^{\prime}_{N,M}(t;k,l)$, whose roots in $(0,1)$ are $\{t(k,l)\}$ and $\{t^{\prime}(k,l)\}$, respectively. Then $t^{\prime}(k,l)\geq(>)t(k,l)$ holds if and only if

[TABLE]

Proof.

It is known that $\tilde{h}^{\prime}_{N,M}(t;k,l)$, whose root is $t^{\prime}(k,l)$, is strictly decreasing in $(0,1)$. Therefore, when $t^{\prime}(k,l)\geq(>)t(k,l)$, it holds that $\tilde{h}^{\prime}_{N,M}(t(k,l);k,l)\geq(>)0$, which leads to (38). The reverse is also true, which completes the proof. ∎

To derive refined coefficients $\{a^{\prime}_{m}\}$ based on present ones $\{a_{m}\}$, it suffices to ensure (38) for all $k$ and $l$ and lift the LHS of (38) as much as possible. This can be achieved, for instance, by simply solving the following linear program (LP):

[TABLE]

By successively solving LP (39), 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 $\{\epsilon(k,l)\}$ bear Pareto optimality.

Theorem 9**.**

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

[TABLE]

and $\{t^{*}(k,l)\}$ lie on the Pareto optimal boundary of (37). That is, there exist no coefficients $\{\tilde{a}^{*}_{m}\}$ that give $\{\tilde{t}^{*}(k,l)\}$ strictly dominating $\{t^{*}(k,l)\}$.

Proof.

Following the same notations in Theorem 4(i), we consider the polynomial $\tilde{h}_{N,M}(t;k,l)$ parameterized by $\{a_{m}^{({\rm ite}+1)}\}$. On one hand, $\tilde{h}_{N,M}(t^{({\rm ite}+1)}(k,l);k,l)=0$. On the other hand, $\tilde{h}_{N,M}(t^{({\rm ite})}(k,l);k,l)\geq 0$ due to (39). As a consequence, $t^{({\rm ite})}(k,l)\leq t^{({\rm ite}+1)}(k,l)$ holds in each iteration due to the monotonicity of $\tilde{h}_{N,M}(t;k,l)$ established in Theorem 4(i). Meanwhile, since $t^{({\rm ite})}(k,l)$ is upper-bounded from $1$, it must converge to a finite limit $t^{*}(k,l)$. However, the sequence $\{a_{m}^{({\rm ite})}\}$ does not necessarily converge. Because $0\leq a_{m}^{({\rm ite})}\leq 1$, there is a convergent subsequence $\{a_{m}^{({\rm ite}_{i})}\}$ with indices $1\leq{\rm ite}_{1}<{\rm ite}_{2}<\cdots$ according to the Bolzano-Weierstrass theorem, i.e.,

[TABLE]

Because $\tilde{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 (39) with $\{t(k,l)=t^{*}(k,l)\}$, and $\{t^{*}(k,l)\}$ are exactly the roots of $\tilde{h}_{N,M}(t;k,l)$ induced by $\{a_{m}^{*}\}$. Then we proceed by contradiction. Suppose that there exist $\{\tilde{a}^{*}_{m}\}$ that give $\{\tilde{t}^{*}(k,l)\}$ strictly dominating $\{t^{*}(k,l)\}$. That is, $\tilde{t}^{*}(k,l)\geq t^{*}(k,l),~{}\forall k,l$ and there exists $(k^{\prime},l^{\prime})$ such that $\tilde{t}^{*}(k^{\prime},l^{\prime})>t^{*}(k^{\prime},l^{\prime})$. Then it immediately follows from Theorem 8 that $\forall k,l$:

[TABLE]

and for the particular indices $(k^{\prime},l^{\prime})$,

[TABLE]

Consequently, $\{\tilde{a}^{*}_{m}\}$ is not only feasible for (39) 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. 3 profiles $\{\epsilon(k,l)\}$ derived with $\{a_{m}\equiv 1/(N+1)\}$ and those after refinement with $\beta=10^{-6}$, $\zeta=8$, $N=100$, $M=5$. It can be seen that values of $\{\epsilon(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 $\{\epsilon(k,l)\}$, one cannot go below the fundamental limits $\{\underline{\epsilon}(k,l)\}$ in Fig. 2.

IV Case Studies on LQR Design of Aircraft Lateral Motion

IV-A Problem Setup

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

[TABLE]

where system states $x(t)=[x_{1}(t)\ x_{2}(t)\ x_{3}(t)\ x_{4}(t)]^{\rm 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)]^{\rm T}$ include the rudder deflection and the aileron deflection. System matrices $A(\delta)$ and $B(\delta)$ are influenced by uncertain parameters $\delta$, as given by:

[TABLE]

[TABLE]

It is assumed that a $10\%$ perturbation of $\delta$ is added around its nominal value $\bar{\delta}$, 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 $-\gamma<0$. That is, the system is stabilized with a desired decay rate $\gamma>0$. In this case we choose $\gamma=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):

[TABLE]

where $P\in\mathbb{R}^{4\times 4}$ and $Y\in\mathbb{R}^{2\times 4}$. $\theta=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 $\Delta$ and only some samples of $\delta$ can be collected. To obtain an approximate solution to (42), we turn to the following scenario program with a total of $d=18$ free decision variables in matrices $P$ and $Y$:

[TABLE]

where $\{\delta^{(i)}\}_{i=1}^{N}$ are randomly collected scenarios of uncertain parameters. To solve the large-scale SDP (43), we use cvx package in MATLAB equipped with the MOSEK solver [37].

IV-B Results and Discussions

In the simulation phase, we randomly generate $N=500$ independent scenarios $\{\delta^{(i)}\}$ for solving (43), and $M=500$ validation scenarios for empirically estimating the violation frequency of the LMI in (42). The confidence is set as $\beta=10^{-6}$ (practical certainty). Since there is no knowledge about a tighter bound of Helly’s dimension, we set $\zeta=d=18$. The generic a priori bound (2) yields $\epsilon=0.0889$. The attractiveness of a posteriori bound lies in that, upon seeing $\mathbb{x}_{N}^{*}$, $s_{N}^{*}$ and $r_{M}^{*}$, the certificate of $V(\mathbb{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}\equiv 1/(N+1)$ to derive $\epsilon(s_{N}^{*},r^{*}_{M})=0.0268$, indicating that with confidence $99.9999\%$ the violation probability $V(\mathbb{x}_{N}^{*})$ is no more than $2.68\%$, which is much tighter than the a priori bound. In contrast, $\epsilon(s_{N}^{*})=0.0486$ with support constraints information used only [22]. If the one-sided C-P bound (5) for Bernoulli trials is adopted, we obtain $\eta_{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(\mathbb{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 $\beta=10^{-6}$. In each run, after deriving the solution $\mathbb{x}^{*}_{N}$, we obtain a high-fidelity estimate $\hat{V}(\mathbb{x}^{*}_{N})$ of its true violation probability $V(\mathbb{x}^{*}_{N})$ with $10^{5}$ additional Monte Carlo samples. Then the gap between a certain bound with $\hat{V}(\mathbb{x}^{*}_{N})$ can be calculated as a performance quantification, which, for instance, is ${\rm GAP}=\epsilon(s_{N}^{*},r_{M}^{*})-\hat{V}(\mathbb{x}^{*}_{N})$ for the proposed bound. By summarizing results in 2000 runs, the mean value $\mu_{\rm GAP}$ and the standard deviation $\sigma_{\rm GAP}$ of the gap are further calculated and summarized in Table I.

Many interesting observations can be attained from Table I. First, in all cases the proposed bound $\epsilon(s_{N}^{*},r_{M}^{*})$ yields the smallest mean and standard deviation of the gap, indicating its significantly reduced conservatism. This is because both $\epsilon(s_{N}^{*})$ and $\eta_{M}(r_{M}^{*})$ only rely on partial information. With fixed $N$ and small $M$, $\eta_{M}(r_{M}^{*})$ tends to be quite conservative, while the proposed bounds provides the tightest certificates by using comprehensive information. With $M$ increasing, both $\eta_{M}(r_{M}^{*})$ and $\epsilon(s_{N}^{*},r_{M}^{*})$ get lower; when $M$ is sufficiently large, the difference between $\epsilon(s_{N}^{*},r_{M}^{*})$ and $\eta_{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 $\epsilon(s_{N}^{*})$ without leveraging validation information. On the other hand, with $N$ increasing, the conservatism of both $\epsilon(s_{N}^{*})$ and $\epsilon(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 $\epsilon(s_{N}^{*},r_{M}^{*})$ and $\epsilon(s_{N}^{*})$ declines with $N$ increasing, while the difference between $\eta_{M}(r_{M}^{*})$ and $\epsilon(s_{N}^{*},r_{M}^{*})$ becomes more pronounced. Henceforth, the proposed bound $\epsilon(s_{N}^{*},r_{M}^{*})$ features a sophisticated integration of information from both support constraints and validation tests. When $N$ is small, $\epsilon(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 $\epsilon(s_{N}^{*},r_{M}^{*})$ under various choices of $N$ and $M$.

Table II 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 [18, 38]. In addition, as analyzed in Section III, the actual improvement in $\epsilon(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 II 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 $\epsilon(s_{N}^{*},r_{M}^{*})$.

Finally, we consider the practical case where independent validation scenarios are successively accumulated after $\mathbb{x}_{N}^{*}$ is obtained. For clarity, a posteriori bounds based on 20 incrementally emerging scenarios have been depicted in Fig. 4. Notice that tendencies of both $\epsilon_{N,M}(s_{N}^{*},r_{M}^{*})$ and $\eta_{M}(r_{M}^{*})$ are consistent with constraint violation results, thereby confirming the correctness of Theorem 6. Another observation is that the updating of $\epsilon_{N,M}(s_{N}^{*},r_{M}^{*})$ starts from $\epsilon_{N,0}(s_{N}^{*},0)=0.1176$, while the evolution of $\eta_{M}(r_{M}^{*})$ has to start from $M=1$. Therefore, at the beginning stage, using $\epsilon_{N,M}(s_{N}^{*},r_{M}^{*})$ tends to be much more viable than using $\eta_{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 $\eta_{M}(r_{M}^{*})$ abundant samples are entailed [2]. 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.

V 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 [22], the proposed probabilistic guarantee does not imply a conditional validity either, i.e., $\mathbb{P}^{N+M}\left\{V(\mathbb{x}_{N}^{*})>\epsilon(k,l)|s_{N}^{*}=k,r_{M}^{*}=l\right\}\leq\beta^{*}$. The recent work [39] 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.