2022

[1]Bowen Li

[1]Sven Leyffer

1]Mathematics and Computer Science Division, Argonne National Laboratory,
Lemont, IL, USA

# Polytopic Superset Algorithm for Nonlinear Robust Optimization

<EMAIL_ADDRESS>Kibaek Kim<EMAIL_ADDRESS><EMAIL_ADDRESS>[

###### Abstract

Nonlinear robust optimization (NRO) is widely used in different applications,
including energy, control, and economics, to make robust decisions under
uncertainty. One of the classical solution methods in NRO is an outer
approximation method that iteratively solves a sample-based nonlinear problem
and updates the sample set by solving an auxiliary problem subject to the
uncertainty set. Although it guarantees convergence under certain assumptions,
its solution iterates are generally infeasible in the original NRO problem,
and it provides only a lower bound on the optimal objective value. We propose
a new algorithm for a class of NRO problems that generates feasible solution
iterates and provides both lower and upper bounds to the optimal objective
value. In each iteration, the algorithm solves the reformulation of an NRO
subproblem with respect to the polytopic supersets of the original uncertainty
set and uses a cutting plane method to improve the supersets over iteration.
If the NRO subproblem is infeasible, we provide a feasibility restoration step
to detect whether the original NRO problem is infeasible or construct a new
superset to restore the feasibility of the NRO subproblem. Further, we prove
that our superset algorithm converges to the optimal solution of the original
NRO problem. In numerical studies, we use application instances from portfolio
optimization and production cost minimization and compare the performance
between the superset algorithm and an outer approximation method called
Polak’s algorithm. Our result shows that the superset algorithm is more
advantageous than Polak’s algorithm when the number of robust constraints is
large.

###### keywords:

Robust optimization, Supersets, Cutting plane method

##  1 Introduction

We consider the following structured nonlinear robust optimization (NRO)
problem,

$\displaystyle\min_{x}\ f(x)\quad\mbox{s.t. }u_{i}^{\top}h_{i}(x)\leq
b_{i}(x),\ \forall u_{i}\in\mathcal{U}_{i},\ i=1,2,...,I,$ (1)

where $x\in\mathbb{R}^{n}$ are the decision variables,
$f(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$ is the objective function, and $I$
denotes the number of robust constraints. For each constraint $i$,
$u_{i}\in\mathbb{R}^{p_{i}}$ are the uncertain parameters of dimension
$p_{i}$; $\mathcal{U}_{i}$ are the uncertainty sets; and
$h_{i}(x):\mathbb{R}^{n}\rightarrow\mathbb{R}^{p_{i}}$ and
$b_{i}(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$ are the nonlinear functions on
$x$.

For clarity and simplicity, we consider only the following single-constrained
NRO for the remainder of the paper:

$\displaystyle\min_{x}\ f(x)\quad\mbox{s.t. }u^{\top}h(x)\leq b(x),\ \forall
u\in\mathcal{U},$ (2)

where $\mathcal{U}:=\\{u\in\mathbb{R}^{p}:\ g(u)\leq 0\\}$ is the uncertainty
set governed by function $g:\mathbb{R}^{p}\rightarrow\mathbb{R}^{m}$, where
$m$ represents the number of constraints in uncertainty set $\mathcal{U}$. We
denote the components of $g$ by $g_{j}:\mathbb{R}^{p}\rightarrow\mathbb{R}$
with $j\in J=\\{1,...,m\\}$. We also define the feasible set
$\mathcal{F}:=\\{x\in\mathbb{R}^{n}:\ u^{\top}h(x)\leq b(x),\ \forall
u\in\mathcal{U}\\}$. It is straightforward to extend the developed results of
(2) to the multiple-constrained NRO (1) or to include general nonlinear
constraints.

######  Assumption 1.

We make the following assumptions for the NRO problem (2):

  1. 1.

For all $j\in J$, $g_{j}(u)$ are convex and continuously differentiable in
$u$.

  2. 2.

The functions $f(x)$, $h(x)$, and $b(x)$ are continuously differentiable in
$x$.

  3. 3.

The uncertainty set $\mathcal{U}$ is compact, nonempty, and the feasible set
$\mathcal{F}$ is bounded.

######  Remark:.

From Assumption 1, because $g_{j}(u)$ are convex, we can efficiently check if
set $\mathcal{U}$ is empty or not. If $\mathcal{U}=\emptyset$, then we have no
robust constraint in (2). Otherwise, we conclude that $\mathcal{U}$ is a
convex set; the gradient functions $\nabla g_{j}(u)$ are Lipschitz continuous
in $u$ for all $j=1,...,m$. The set $\mathcal{F}$ is closed because it is an
intersection of an infinite number of closed sets. From closedness and
boundedness of $\mathcal{F}$, we further conclude that $\mathcal{F}$ is a
compact set, which implies that (2) has a solution if $\mathcal{F}$ is
nonempty. We will discuss more about the feasibility issue in Section 3.5.

NRO has many important applications, of which we select two as our motivation
for this paper. One instance of NRO is in portfolio optimization port1 ; port2
. This problem maximizes the risk-adjusted expected return under the
uncertainty of the mean and covariance of the asset returns by finding the
optimal asset allocation. Different models of the uncertainty set are
discussed in port2 . The dimension of the uncertain parameter is quadratic in
the number of assets. Another example of NRO arises in production cost
minimization prod_cost . Under the uncertainty of price, this problem
minimizes the production cost as well as the cost associated with production
ramping while satisfying the daily demand requirements and ramping
limitations. The number of robust constraints is proportional to the
scheduling horizon. All the conditions in Assumption 1 are satisfied by these
numerical instances.

Different iterative methods have been developed to solve NRO and can be
broadly categorized into discretization methods dis1 ; tut4 , exchange methods
outer1 ; polak1997 ; polak_boyd ; polak3 , and local reduction methods tut1 ;
tut3 ; tut5 . In this work, our benchmark is an outer approximation method
called Polak’s algorithm polak1997 ; polak_boyd  in the category of the
exchange methods. Polak’s algorithm iteratively solves a sample-based
subproblem that approximates the uncertainty set $\mathcal{U}$ of the original
NRO problem (2) by a finite sample set. Because the finite sample set is a
subset of the uncertainty set $\mathcal{U}$, the sample-based subproblem can
be seen as an outer approximation of (2). As the iteration progresses, new
sample points are obtained by solving a worst-case constraint violation
problem and are accumulated to be new constraints of the next-iteration
subproblem. Before termination, the resulting solution iterates of Polak’s
algorithm are not feasible in the original NRO problem (2) and provide only a
lower bound on the optimal objective value of (2).

The major benefit of the affine relationship between the constraint $h(x)$ and
uncertainty parameter $u$ in (2) is the existence of a reformulation approach.
If we approximate the uncertainty set $\mathcal{U}$ with a polytopic superset
and apply linear programming duality to the robust constraints, then the
resulting NRO subproblem can be exactly reformulated into a finite-dimensional
nonlinear problem similar to nonlin1 ; rob_book , which can be solved by off-
the-shelf nonlinear optimization solvers. Hence, we propose a superset
algorithm that iteratively solves the reformulation of the NRO subproblem with
polytopic supersets and aims to improve the solution quality by gradually
constructing better supersets. We show that our proposed superset algorithm
follows an iterative structure similar to Polak’s algorithm polak_boyd  but
generates feasible iterates.

Our contribution is the development of a superset algorithm that iteratively
solves the reformulation of an NRO subproblem with polytopic supersets of the
original uncertainty set. We propose different cutting plane methods to
improve the supersets iteratively. Compared with Polak’s algorithm polak1997 ;
polak_boyd , the solution iterates of our algorithm are feasible in the
original NRO problem and provide both upper and lower bounds to the optimal
objective value. We show that with the proposed cutting plane methods, the
solution of the NRO subproblem with polytopic supersets converges to the
optimal solution of the original NRO problem. In addition, we provide a
feasibility restoration algorithm to restore the feasibility if the initial
supersets are overly conservative or detect whether the original NRO problem
is infeasible. To demonstrate the computation performance, we compare the
superset algorithm with Polak’s algorithm with test instances from portfolio
optimization and production cost minimization. We show that the superset
algorithm outperforms Polak’s algorithm in most of the test cases, especially
when the number of robust constraints is large.

The outline of this paper is as follows. Section 2 discusses theoretical
background about general NRO problems. Section 3 discusses the benefits of
using polytopic supersets and gives a detailed description and theoretical
analysis of our proposed superset algorithm. Section 4 gives the convergence
analysis, and Section 5 shows the numerical results of comparing the superset
algorithm with Polak’s algorithm in portfolio optimization and production cost
minimization. Section 6 summarizes our work and briefly discusses future
plans.

##  2 Background

In this section we give some theoretical results from the literature as
background. For simplicity, we define the following constants and notations
for the rest of the paper. First, for the convergence analysis in Section 4,
we define a unified Lipschitz constant $L$ for all the gradient functions
$\nabla g_{j}(u)$ and radius $R$ for the bounded set $\mathcal{U}$. Second, we
define the following notation rules: For a finite set $U$, $\lvert U\rvert$
denotes its number of elements; and for a closed set $\mathcal{U}$,
$\partial\mathcal{U}$ denotes its boundary. We define $[x]^{+}=\max\\{x,0\\}$.
For a vector $x$, we define $x\nleq 0$ if there exists a component $x_{i}$
such that $x_{i}>0$.

###  2.1 General Properties of Robust Optimization

Here we recall some of the main theoretical results regarding the optimal
solution of (2). Given $x\in\mathcal{F}$, we define the active set
$\mathcal{A}(x):=\\{u\in\mathcal{U}:\ u^{\top}h(x)=b(x)\\}$ and Jacobian of
$h(x)$ as $\nabla h(x)=[\nabla h_{1}(x),...,\nabla h_{p}(x)]$. Next, we give
the definition of the Mangasarian–Fromovitz constraint qualification (MFCQ) of
(2):

######  Definition 1 (Mangasarian–Fromovitz Constraint Qualification emfcq1 ;
emfcq2 ; emfcq3 ).

We say that $x^{*}$ satisfies the MFCQ if there exists $s\in\mathbb{R}^{n}$
such that

$\displaystyle u^{\top}h(x^{*})-b(x^{*})+s^{\top}\left(\nabla h(x^{*})u-\nabla
b(x^{*})\right)<0,\ \forall u\in\mathcal{U}.$

Now with MFCQ we give the first-order condition tut1 ; tut2 ; sven_report  as
follows.

######  Theorem 1.

Suppose $x^{*}\in\mathcal{F}$ satisfies MFCQ. If $x^{*}$ is a local minimizer
of (2), then there exist a finite subset $A^{*}\subseteq\mathcal{A}(x^{*})$
and multipliers $\lambda^{*}_{u}\geq 0$ for each $u\in A^{*}$ such that

$\displaystyle\nabla f(x^{*})+\sum_{u\in A^{*}}\lambda^{*}_{u}\left(\nabla
h(x^{*})u-\nabla b(x^{*})\right)=0.$ (3)

Theorem 1 motivates the error measure $errm(x,\lambda_{u},A)$ to evaluate the
first-order condition error given solution estimate $x$, a finite set $A$, and
the corresponding multipliers $\lambda_{u}$ for all $u\in A$:

$\displaystyle errm(x,\lambda_{u},A)$ $\displaystyle=\left\|\nabla
f(x)+\sum_{u\in{A}}\lambda_{u}\left(\nabla h(x)u-\nabla
b(x)\right)\right\|_{2}$ (4a)
$\displaystyle+\left[\max_{u\in\mathcal{U}}u^{\top}h(x)-b(x)\right]^{+}$ (4b)
$\displaystyle+\sum_{u\in{A}}\ [-\lambda_{u}]^{+}$ (4c)
$\displaystyle+\sum_{u\in{A}}\ \lvert u^{\top}h(x)-b(x)\rvert\
+\sum_{u\in{A}}\ \|u-\text{proj}_{\mathcal{U}}(u)\|_{2}.$ (4d)

The function $\text{proj}_{\mathcal{U}}()$ can be any projection operator to
the set $\mathcal{U}$. Without loss of generality, we use the Euclidean
projection and define

$\displaystyle\text{proj}_{\mathcal{U}}(x)=\text{argmin}_{u\in\mathcal{U}}\
\|u-x\|_{2}^{2}.$ (5)

Comparing (4) with Theorem 1, we observe that (4a) corresponds to the first-
order condition (3), (4b) corresponds to primal feasibility, (4c) corresponds
to the non-negativity of the Lagrangian multipliers, and (4d) checks the
activity condition $A\subseteq\mathcal{A}(x)$. We note that
$A\subseteq\mathcal{A}(x)$ if and only if (4d) equals zero.

###  2.2 Polak’s Algorithm

As a benchmark, we also consider an outer approximation method, namely,
Polak’s algorithm sven_report ; polak1997 ; polak_boyd . The algorithm is
based on a sample-based outer approximation of (2). At iteration $k$, given a
finite sample set $U_{k}\subseteq\mathcal{U}$, the algorithm solves the
following finite-dimensional nonlinear problem:

$\displaystyle\min_{x}\ f(x)\quad\mbox{s.t. }u^{\top}h(x)\leq b(x),\ \forall
u\in U_{k}.$ (6)

Because we approximate $\mathcal{U}$ with a finite sample set $U_{k}$, (6) is
a relaxation of (2). The solution of (6) is not guaranteed to be feasible in
(2) and provides only a lower bound on the optimal objective value of (2). In
each iteration $k$, to improve the current solution $x_{k}$, we solve the
following problem:

$\displaystyle\max_{t,u\in\mathcal{U}}\ t\quad\mbox{s.t.
}u^{\top}h(x_{k})-b(x_{k})\geq t.$ (7)

Let $\hat{u}_{k}$ and $t_{k}$ be the solution and optimal objective value of
(7), respectively. Clearly, $t_{k}\leq 0$ if and only if
$x_{k}\in\mathcal{F}$. When $t_{k}>0$, it corresponds to the worst-case
constraint violation at uncertainty realization $\hat{u}_{k}$. In this case,
$\hat{u}_{k}$ will be added to the sample set $U_{k}$ and generates a new
constraint in the subproblem of the next iteration. As the sample size
increases, the algorithm terminates when $t_{k}\leq\epsilon$ for some
tolerance $\epsilon>0$. The algorithmic steps are summarized in Algorithm 1.

initialization: Set $k\leftarrow 0$; initialize $U_{0}\leftarrow\emptyset$ and
some $\epsilon>0$.

repeat

sample-based subproblem: Let $x_{k}$ be the solution of (6).

update sample set: Let $\hat{u}_{k}$ and $t_{k}$ be the solution of (7).

Set $U_{k+1}\leftarrow U_{k}\cup\\{\hat{u}_{k}\\}$.

Set $k\leftarrow k+1$.

until $t_{k-1}\leq\epsilon$.

return $x_{k-1}$.

Algorithm 1 Polak’s Algorithm

Note that when $k=0$, (6) can be unbounded. In practical applications,
however, deterministic constraints or the inclusion of a nominal sample can
resolve this issue without loss of generality. By applying the proof from
polak_boyd  to (2), Algorithm 1 is guaranteed to converge.

##  3 Polytopic Superset Algorithm

In this section we propose a new polytopic superset algorithm and discuss its
details. Specifically, in Section 3.1 we define a polytopic superset NRO; and,
based on this, we give an algorithm sketch to the superset algorithm. In
Section 3.2 we develop a reformulation approach given the affine relationship
between $u$ and $h(x)$ in (2) to efficiently solve the polytopic superset NRO.
In Section 3.3, we discuss the cutting plane methods that we use to improve
the polytopic superset NRO. In Section 3.4 we extend the error measure (4) to
the proposed polytopic superset algorithm and use it as the termination
criteria. In Section 3.5 we develop a feasibility restoration step to either
detect whether (2) is infeasible or construct an appropriate superset that
guarantees the feasibility of the polytopic superset NRO. In Section 3.6 we
give the full algorithmic steps of the polytopic superset algorithm and
provide some theoretical analysis.

###  3.1 Sketch of a Polytopic Superset Algorithm

In each iteration $k$, we define the following NRO subproblem using a
polytopic superset $\mathcal{S}_{k}\supseteq\mathcal{U}$ of the current
iteration:

$\displaystyle\min_{x}\ f(x)\quad\mbox{s.t. }u^{\top}h(x)\leq b(x),\ \forall
u\in\mathcal{S}_{k}.$ (8)

Based on (8), we propose the following superset algorithm.

initialization: Set $k\leftarrow 0$; initialize a box superset
$\hat{\mathcal{S}}\supseteq\mathcal{U}$ and some $\epsilon>0$.

feasibility restoration: Detect if (2) is infeasible, or construct a new
superset $\mathcal{S}_{0}$ from $\hat{\mathcal{S}}$ that guarantees the
feasibility of (8).

repeat

superset subproblem: Let $x_{k}$ be the solution of (8).

update superset: Update $\mathcal{S}_{k}$ by adding cutting planes.

Set $k\leftarrow k+1$.

until $x_{k-1}$ is optimal to (2) for some tolerance $\epsilon$.

return $x_{k-1}$.

Algorithm 2 Sketch of the Superset Algorithm

We initialize the algorithm using a box superset denoted $\hat{\mathcal{S}}$
that contains the compact uncertainty set $\mathcal{U}$. By construction, the
set $\hat{\mathcal{S}}$ is compact. As the first step, we use feasibility
restoration to either provide a certificate of infeasibility to (2) or
construct a superset $\mathcal{S}_{0}$ based on $\hat{\mathcal{S}}$ such that
(8) is guaranteed to be feasible when $k=0$ (see Section 3.5). Over the
iteration, we keep solving the NRO subproblem (8) and update the superset with
cutting planes generated using $x_{k}$ (see Section 3.3). These cutting planes
help reduce the gap between the superset $\mathcal{S}_{k}$ and the uncertainty
set $\mathcal{U}$. The algorithm terminates when $x_{k}$ is optimal to (2)
with some given tolerance $\epsilon$.

Because $\mathcal{S}_{k}\supseteq\mathcal{U}$ for all $k$, all iterates of (8)
are also feasible in (2) and provide an upper bound on the optimal objective
value of (2). In other words, the NRO subproblem (8) can be seen as an inner
approximation of (2), which is distinct from the classical outer approximation
methods polak1997 ; polak_boyd ; polak3 . In Section 3.6 we will discuss how
the superset algorithm also provides a lower bound on the optimal objective
value of (2) through the Euclidean projection used in the termination
criteria.

###  3.2 Reformulation of the Polytopic Superset NRO

In this section we develop a reformulation approach to efficiently solve (8).
First, we model the polytopic superset for each iteration $k$ using linear
inequalities. We define $\mathcal{S}_{k}:=\\{u:\ B_{k}u\leq d_{k}\\}$ with
linear coefficients $B_{k}\in\mathbb{R}^{b_{k}\times p}$ and
$d_{k}\in\mathbb{R}^{b_{k}}$, where $b_{k}$ denotes the number of linear
inequality constraints in $\mathcal{S}_{k}$. In our algorithm, $B_{k}u\leq
d_{k}$ corresponds to a set of supporting hyperplanes for the convex
uncertainty set $\mathcal{U}$, and we show in Section 3.3 how to update
$\mathcal{S}_{k}$. Because of the affine relationship between $h(x)$ and $u$
in (8), we have the following proposition by applying the linear programming
duality to the robust constraint.

######  Proposition 1.

Given polytopic supersets $\mathcal{S}_{k}:=\\{u\in\mathbb{R}^{p}:\ B_{k}u\leq
d_{k}\\}$, problem (8) can be reformulated into the following finite-
dimensional nonlinear optimization problem:

$\displaystyle\min_{x,\gamma_{k}}\ $ $\displaystyle f(x)$ (9a) s.t.
$\displaystyle\gamma_{k}^{\top}d_{k}\leq b(x),$ (9b)
$\displaystyle(B_{k})^{\top}\gamma_{k}=h(x),$ (9c)
$\displaystyle\gamma_{k}\geq 0.$ (9d)

Proof: The proof can be obtained by generalizing the reformulation approach
used for a robust linear program with a polytopic uncertainty set nonlin1 ;
rob_book  to our structured nonlinear problem (8). First, because
$\mathcal{S}_{k}$ is compact, we have

$\displaystyle\max_{u\in\mathcal{S}_{k}}u^{\top}h(x_{k})\leq
b(x_{k})\Leftrightarrow\ u^{\top}h(x_{k})\leq b(x_{k}),\ \forall
u\in\mathcal{S}_{k}.$ (10)

Next, by linear programming duality, we have

$\displaystyle\min_{\gamma\geq
0,(B_{k})^{\top}\gamma=h(x_{k})}\gamma^{\top}d_{k}=\max_{u\in\mathcal{S}_{k}}u^{\top}h(x_{k}).$
(11)

Substituting the robust constraint in (8) in the left-hand side of (11), we
obtain reformulation (9). Our superset algorithm iteratively solves
reformulation (9) of the NRO subproblem (8) with polytopic supersets. The
reformulation is a finite-dimensional nonconvex problem that can be solved by
off-the-shelf nonlinear optimization solvers. From Proposition 1, (8) and (9)
are equivalent in the sense that they have the same global optimal solution.
Unfortunately, even if (8) is convex, $h(x)$ is convex and $b(x)$ is concave,
the reformulation (8) of (9) is nonconvex. However, we show next, that we can
at least recover the equivalence of stationary points in this reformulation.
In the next theorem we show that the stationary points of (9) satisfy the
first-order condition (3) of the infinite-dimensional NRO subproblem (8).

######  Theorem 2.

If $(x_{k},\gamma_{k})$ satisfies the KKT conditions of (9), then $x_{k}$ also
satisfies the first-order condition (3) (applied to (8)).

Proof: First, we analyze the KKT conditions for (9) where $\mu_{k}$, $v_{k}$,
and $\nu_{k}$ are the corresponding multipliers for (9b), (9c), and (9d),
respectively:

$\displaystyle\nabla f(x_{k})+\nabla h(x_{k})v_{k}-\mu_{k}\nabla b(x_{k})=0,$
(12a) $\displaystyle-B_{k}v_{k}+\mu_{k}d_{k}-\nu_{k}=0,$ (12b)
$\displaystyle\mu_{k}\geq 0,\ \nu_{k}\geq 0,\
\mu_{k}\left(\gamma_{k}^{\top}d_{k}-b(x_{k})\right)=0,\
\nu_{k}^{\top}\gamma_{k}=0,$ (12c) $\displaystyle\gamma_{k}^{\top}d_{k}\leq
b(x_{k}),\ \gamma_{k}\geq 0,\ B_{k}^{\top}\gamma_{k}=h(x_{k}).$ (12d)

When $\mu_{k}>0$, we have that

$\displaystyle\eqref{stt3}$ $\displaystyle\Rightarrow\
\gamma_{k}^{\top}d_{k}=b(x_{k}),$ (13a) $\displaystyle\eqref{stt2}\mbox{ and
}\eqref{stt4}$ $\displaystyle\Rightarrow\
B_{k}(\frac{v_{k}}{\mu_{k}})=d_{k}-\frac{\nu_{k}}{\mu_{k}}\Rightarrow
B_{k}(\frac{v_{k}}{\mu_{k}})\leq
d_{k}\Leftrightarrow\frac{v_{k}}{\mu_{k}}\in\mathcal{S}_{k}.$ (13b)

Further, we have

$\displaystyle\eqref{eq:13b}\Rightarrow\gamma_{k}^{\top}B_{k}(\frac{v_{k}}{\mu_{k}})=\gamma_{k}^{\top}d_{k}-\gamma_{k}^{\top}(\frac{\nu_{k}}{\mu_{k}})\Rightarrow(\frac{\nu_{k}}{\mu_{k}})^{\top}h(x_{k})=b(x_{k}),$
(14)

where the second implication comes from (12c) and (12d). Further, from (10),
(13), and (14), we conclude that $u^{\top}h(x_{k})=b(x_{k})$ is a supporting
hyperplane to $\mathcal{S}_{k}$ and
$\frac{v_{k}}{\mu_{k}}\in\partial\mathcal{S}_{k}$.

Next, we show that $\mu_{k}=0$ implies $v_{k}=0$ by contradiction. Suppose
$v_{k}\neq 0$. Because $\mathcal{S}_{k}$ is compact by construction, we have
that $B_{k}$ has full-column rank. From (12b) and (12c), we then have
$B_{k}v_{k}=-\nu_{k}\leq 0$ (since $\mu_{k}=0$). Because $v_{k}\neq 0$ and
$B_{k}$ is full column rank, we conclude that $\nu_{k}\neq 0$. This implies
that for any $\hat{u}\in\mathcal{S}_{k}$, we can pick $\alpha>0$ arbitrarily
large such that $\hat{u}+\alpha v_{k}\in\mathcal{S}_{k}$ because
$B_{k}(\hat{u}+\alpha v_{k})=B_{k}\hat{u}-\alpha\nu_{k}\leq d_{k}$, which
contradicts the compactness of $\mathcal{S}_{k}$. Thus, we have proved that
$\mu_{k}=0$ implies $\nu_{k}=0$.

Now, we can rewrite (12a) as follows:

$\displaystyle\nabla f(x_{k})+\mu_{k}\left(\nabla
h(x_{k})\left(\frac{v_{k}}{\mu_{k}}\right)-\nabla
b(x_{k})\right)=0,\quad\mbox{if $\mu_{k}>0$},$ (15a) $\displaystyle\nabla
f(x_{k})=0,\quad\mbox{if $\mu_{k}=0$}.$ (15b)

From (15) we conclude that $x_{k}$ satisfies the first-order condition of (8).

We observe that by substituting the optimal multipliers of (9) into (15), we
get the first-order condition of (8). Depending on the strict positivity of
$\mu_{k}$, however, we denote the solution pair of (9) as
$(x_{k},\lambda_{k},u_{k})$, where

$\displaystyle\lambda_{k}=\begin{cases}\mu_{k}>0,&\text{if }\mu_{k}>0\\\
0,&\text{otherwise}\end{cases},\quad
u_{k}=\begin{cases}\frac{v_{k}}{\mu_{k}},&\text{if }\mu_{k}>0\\\ \text{any
}u\in\mathcal{U},&\text{otherwise.}\end{cases}$ (16)

###  3.3 Projection and Cuts

To improve the solution of (9) over iteration, we use cutting planes to reduce
the size of $\mathcal{S}_{k}$ and remove the violated uncertainty realization
in the current solution pair. At iteration $k$, we have the solution pair
$(x_{k},\lambda_{k},u_{k})$ from the reformulation (9). If
$u_{k}\notin\mathcal{U}$, we add cutting planes that are valid for the
uncertainty set $\mathcal{U}$, to remove the violated point $u_{k}$ from
$\mathcal{S}_{k}$ and reduce the gap between the superset $\mathcal{S}_{k}$
and the uncertainty set $\mathcal{U}$. Next, we will provide different ways to
generate cutting planes that exclude such points $u_{k}$. We define the
Jacobian of $g(u)$ as $\nabla g(u)=[\nabla g_{1}(u),...,\nabla g_{m}(u)]$.

####  3.3.1 Kelley’s Cutting Plane

Given a point $u_{k}$ with $g(u_{k})\nleq 0$, we first define Kelley’s cutting
plane on $u$ as from Kelley

$\displaystyle\nabla g(u_{k})^{\top}u\leq\nabla
g(u_{k})^{\top}u_{k}-g(u_{k}).$ (17)

Because $g$ is convex, $\nabla g(u_{k})^{\top}(u-u_{k})+g(u_{k})\leq g(u)\leq
0$ for all $u\in\mathcal{U}$. We note that $u_{k}$ violates the cut because
$g(u_{k})\nleq 0$.

We also note that the cutting plane (17) may not support the uncertainty set
$\mathcal{U}$. For example, with uncertainty set
$\mathcal{U}:=\\{(u_{1},u_{2})\ \lvert\ u_{1}^{2}+u_{2}^{2}\leq 1\\}$ with
violated $u=(2,0)\notin\mathcal{U}$, the resulting cutting plane is $u_{1}\leq
1.25$, which does not support $\mathcal{U}$. Next, we provide two alternatives
to address this issue.

####  3.3.2 Euclidean Projection Cut

We can generate an alternative cutting plane at Euclidean projection
$z_{k}=\text{proj}_{\mathcal{U}}(u_{k})$ that removes $u_{k}$ and supports
$\mathcal{U}$. The cutting plane is given by

$\displaystyle\nabla g(z_{k})^{\top}u\leq\nabla
g(z_{k})^{\top}z_{k}-g(z_{k}).$ (18)

Cutting plane (18) is also valid for $\mathcal{U}$ and
$z_{k}\in\partial\mathcal{U}$ because of the Euclidean projection. Hence, (18)
supports $\mathcal{U}$. The exclusion of $u_{k}$ when $u_{k}\notin\mathcal{U}$
was proved in Theorem 1 of eudcut .

####  3.3.3 Gradient-Free Cut

We replace the Euclidean projection cut with its gradient-free version.
Instead of using $\nabla g$, we use the Euclidean projection $z_{k}$ and
define the gradient-free cut as follows:

$\displaystyle(u_{k}-z_{k})^{\top}u\leq(u_{k}-z_{k})^{\top}z_{k}.$ (19)

Inequality (19) is valid for $\mathcal{U}$ because of the variational
characteristic of projection proj_prop . Because
$z_{k}\in\partial\mathcal{U}$, (19) supports $\mathcal{U}$. The exclusion of
$u_{k}$ can also be seen by plugging it in because $u_{k}\neq z_{k}$.

######  Example:.

Here we consider the following NRO problem to demonstrate the difference
between the cutting planes:

$\displaystyle\min_{x_{1},x_{2}}\quad$ $\displaystyle-x_{1}-x_{2}$ (20) s.t.
$\displaystyle x_{1}^{2}u_{1}+x_{2}^{2}u_{2}\leq 6,\
\forall(u_{1},u_{2})^{\top}\in\mathcal{U},$

where $\mathcal{U}=\\{u_{1}^{2}+u_{2}^{2}\leq 1,u_{1}\geq 0,u_{2}\geq 0\\}$.
The solution of (20) is $x^{*}=(\sqrt{3\sqrt{2}},\sqrt{3\sqrt{2}})$ with
optimal multiplier $\lambda^{*}=\sqrt{3\sqrt{2}}/6$ and
$u^{*}=(\sqrt{2}/2,\sqrt{2}/2)$.

(a) $\mathcal{S}_{0}$ and first cut.

(b) $\mathcal{S}_{1}$ and second cut.

(c) $\mathcal{S}_{2}$ and third cut.

Figure 1: Illustration of the superset update in Algorithm 4 using Kelley’s
cutting plane (top) and Euclidean projection cut (bottom).

In Figure 1 we illustrate the superset update in Algorithm 2 with all proposed
types of cutting planes for the first three iterations (i.e., $k=0,1,2$). The
algorithm starts with initial superset $\mathcal{S}_{0}=\\{0\leq u_{1}\leq 1,\
0\leq u_{2}\leq 2\\}$ (rectangular). In this example, gradient-free cuts
coincide with the Euclidean projection cuts and result in the same algorithm
progress because of the simplicity of the uncertainty set. We observe that the
Euclidean projection cut and gradient-free cut are more effective than
Kelley’s cutting plane in terms of improving the superset to better
approximate the uncertainty set. The second cut and third cut from Kelley’s
cutting method are close to each other and can be less effective. The detailed
iteration steps are summarized in Table 1(b).

Table 1: Solution pair $(x_{k},u_{k},\lambda_{k})$ and cutting planes obtained
in the first three iterations of Algorithm 2 (i.e., $k=0,1,2$) where
$a_{1}=\frac{4\sqrt{6}+\sqrt{12}-\sqrt{5}-2\sqrt{10}}{7}$ and
$a_{2}=\frac{(\sqrt{10}-\sqrt{3})(1+2\sqrt{2})}{7}$.

(a) Kelley’s Cutting Plane. | Iteration
---
$k$
$x_{k}$ | $u_{k}$ | $\lambda_{k}$ |  $k$th Cutting Plane
$0$ | $(2,1)$ | $(1,2)$ | $1/4$ | $u_{1}+2u_{2}\leq 3$
$1$ | $(\sqrt{3},\sqrt{3})$ | $(1,1)$ | $\sqrt{3}/6$ | $u_{1}+u_{2}\leq 1.5$
$2$ | $(2,2)$ | $(0.75,0.75)$ | $1/3$ | $u_{1}+u_{2}\leq 443/300$

(b) Euclidean Projection Cut.

| Iteration
---
$k$
$x_{k}$ | $u_{k}$
0 | $(2,1)$ | $(1,2)$
1 | $\left(\frac{\sqrt{30\sqrt{5}}}{5},\frac{2\sqrt{15\sqrt{5}}}{5}\right)$ | $\left(\sqrt{5}(\sqrt{2}-1),\frac{\sqrt{10}(\sqrt{2}-1)}{2}\right)$
2 | $\left(\sqrt{\frac{6a_{2}}{a_{1}a_{2}+a_{1}^{2}}},\sqrt{\frac{6a_{1}}{a_{1}a_{2}+a_{2}^{2}}}\right)$ | $(a_{1},a_{2})$

| Iteration
---
$k$
$\lambda_{k}$ |  $k$th Cutting Plane
0 | $1/4$ | $\sqrt{5}u_{1}+2\sqrt{5}u_{2}\leq 5$
1 | $\frac{\sqrt{6}(\sqrt{2}+1)}{12}$ | $\sqrt{6}u_{1}+\sqrt{3}u_{2}\leq 3$
2 | $\frac{\sqrt{a_{1}+a_{2}}}{2\sqrt{6a_{1}a_{2}}}$ | $\sqrt{\frac{a_{1}^{2}}{a_{1}^{2}+a_{2}^{2}}}u_{1}+\sqrt{\frac{a_{2}^{2}}{a_{1}^{2}+a_{2}^{2}}}u_{2}\leq 1$

###  3.4 Termination Criteria

We consider the NRO subproblem (8) using $\mathcal{S}$ and denote its solution
pair by $(\bar{x},\bar{\lambda},\bar{u})$ (from its reformulation (9)). To
evaluate the solution quality to (2), we use (4) as follows:

$\displaystyle
errm(\bar{x},\bar{\lambda},\\{\bar{u}\\})=\|\bar{u}-\text{proj}_{\mathcal{U}}(\bar{u})\|_{2}.$
(21)

Because $\bar{u}\in\mathcal{S}$ and $\mathcal{S}\supseteq\mathcal{U}$, (21)
can be nonzero. All the remaining components in (4) are zero because of the
proof of Theorem 2. For brevity, we will use a short notation $errm(x)$ in the
future discussion of the superset algorithm.

##### Comparison with Polak’s Subproblem

As comparison, we consider the sample-based finite nonlinear problem (6) using
$U_{k}$ and denote its solution as $\hat{x}$. Define a finite set
$\hat{A}:=\\{u\in U_{k}\ \lvert\ u^{\top}h(\hat{x})=b(\hat{x})\\}$, and let
$\hat{\lambda}_{u}$ be the corresponding optimal multipliers of (6) for all
$u\in\hat{A}$. We use (4) to evaluate the solution quality of $\hat{x}$ to (2)
as follows:

$\displaystyle
errm(\hat{x},\hat{\lambda}_{u},\hat{A})=\left[\max_{u\in\mathcal{U}}u^{\top}h(\hat{x})-b(\hat{x})\right]^{+}.$
(22)

Because $\hat{x}$ can guarantee feasibility only for the given sample set
$U_{k}$ but not for the original uncertainty set $\mathcal{U}$, (22) can be
nonzero. All the remaining components (4a), (4c), and (4d) are zero because of
the KKT condition of (6).

From the comparison above, we observe that the solutions of the polytopic
superset NRO and sample-based subproblem of Polak’s algorithm achieve
different KKT guarantees. The iterates from Polak’s algorithm violate primal
feasibility. On the other hand, the iterates from the superset algorithm
violate the activity check but are guaranteed to satisfy primal feasibility.

###  3.5 Feasibility Restoration

If (8) is infeasible for $k=0$, there are two possibilities: either the
original problem (2) is infeasible or the initial $\mathcal{S}_{0}$ is too
conservative. To detect whether (2) is infeasible or restore the feasibility
of (8), we define the following phase-I problem to (2):

$\displaystyle\min_{x,p\geq 0}\ p\quad\mbox{s.t. }u^{\top}h(x)\leq b(x)+p,\
\forall u\in\mathcal{U}.$ (23)

One can easily see that (23) is always feasible.

######  Proposition 2.

(23) has optimal objective value $p^{*}=0$ if and only if (2) is feasible.

The proof can be obtained by direct checking. To solve (23), we will use the
superset algorithm starting with the initial box superset $\hat{\mathcal{S}}$.
We define the NRO subproblem to (23) with polytopic superset $\mathcal{S}_{k}$
at iteration $k$ as

$\displaystyle\min_{x,p\geq 0}\ p\quad\mbox{s.t. }u^{\top}h(x)\leq b(x)+p,\
\forall u\in\mathcal{S}_{k}.$ (24)

Similar to (23) being a phase-I problem to (2), (24) is a phase-I problem to
(8) and can be solved with its reformulation from Proposition 1. Based on
(24), we design the following feasibility restoration algorithm.

1:initialization: Set $k\leftarrow 0$, and initialize superset
$\mathcal{S}_{0}\leftarrow\hat{\mathcal{S}}$ and some $\epsilon>0$.

2:repeat

3: superset subproblem: Let $({x}_{k},p_{k})$ solve (24).

4: if $p_{k}=0$ then

5: return ${x}_{k}$, $p_{k}$, $\mathcal{S}_{k}$.
$\ignorespaces\triangleright$ $x_{k}$ is feasible in (2).

6: end if

7: update superset: Update $\mathcal{S}_{k}$ with the following cutting
planes:

8: $\bullet$ Kelley’s cutting plane: $\nabla g(u_{k})^{\top}u\leq\nabla
g(u_{k})^{\top}u_{k}-g(u_{k})$, or

9: $\bullet$ Euclidean Projection cut: $\nabla g(z_{k})^{\top}u\leq\nabla
g(z_{k})^{\top}z_{k}-g(z_{k})$, or

10: $\bullet$ Gradient-free cut:
$(u_{k}-z_{k})^{\top}u\leq(u_{k}-z_{k})^{\top}z_{k}.$

11: Set $k\leftarrow k+1$.

12:until $errm({x}_{k-1})\leq\epsilon$.

13:return ${x}_{k-1}$, $p_{k-1}$, $\mathcal{S}_{k-1}$.
$\ignorespaces\triangleright$ (2) is infeasible with $p_{k-1}>0$.

Algorithm 3 Feasibility Restoration:
$(x,p,\mathcal{S})\leftarrow\text{FR}(\hat{\mathcal{S}},\epsilon).$

Here the error measure $errm$ in line 12 follows the same definition to (21).
Algorithm 3 has two termination criteria. When $p_{k}=0$ in line 4, we obtain
${x}_{k}$, which is feasible in (2). In this case, superset $\mathcal{S}_{k}$
is returned and will be used in the main loop of Algorithm 2. Otherwise, in
line 13 we terminate with an $\epsilon$-stationary point $x_{k-1}$ to the
phase-I problem (23) with a certificate of infeasibility $p_{k-1}>0$. In
Section 4 we show that the Algorithm 3 results in the convergence to
stationarity.

###  3.6 Full Algorithm

The complete steps of the superset algorithm are described in Algorithm 4.
First, we recall the Euclidean projection
$z_{k}=\text{proj}_{\mathcal{U}}(u_{k})$.

initialization: Set $k\leftarrow 0$; initialize a box superset
$\hat{\mathcal{S}}\supseteq\mathcal{U}$ and some $\epsilon>0$.

feasibility restoration:
$(x_{0},p_{0},\mathcal{S}_{0})\leftarrow\text{FR}(\hat{\mathcal{S}},\epsilon)$.

if $p_{0}>0$ then

return $x_{0}$ and the infeasibility of (2), $p_{0}$.

end if

repeat

superset subproblem: Let $x_{k}$ be the solution of (9).

update supersets: Update $\mathcal{S}_{k}$ with cutting planes:

$\bullet$ Kelley’s cutting plane: $\nabla g(u_{k})^{\top}u\leq\nabla
g(u_{k})^{\top}u_{k}-g(u_{k})$, or

$\bullet$ Euclidean Projection cut: $\nabla g(z_{k})^{\top}u\leq\nabla
g(z_{k})^{\top}z_{k}-g(z_{k}),$ or

$\bullet$ Gradient-free cut:
$(u_{k}-z_{k})^{\top}u\leq(u_{k}-z_{k})^{\top}z_{k}.$

Set $k\leftarrow k+1$.

until $errm(x_{k-1})\leq\epsilon.$

return $x_{k-1}$.

Algorithm 4 Superset Algorithm

If (2) is infeasible or if the initial box superset $\hat{\mathcal{S}}$ is
overly conservative, then (8) can be infeasible. Here we use the feasibility
restoration step (see Section 3.5) to either restore the feasibility by
constructing a new superset $\mathcal{S}_{0}$ for the NRO subproblem (8) or
return a stationary point $x_{0}$ to a phase-I problem (23) with an
infeasibility certificate $p_{0}>0$.

After exiting feasibility restoration with $p_{0}=0$, we iteratively solve the
NRO subproblem (8) and update the superset $\mathcal{S}_{k}$ with cutting
planes. We define the feasible set of (8) at iteration $k$ as
$\mathcal{F}_{k}$. As we keep adding cutting planes, we conclude that
$\mathcal{F}_{k}\subseteq\mathcal{F}_{k+1}\subseteq\mathcal{F}$ for all $k$
(i.e., the feasible set of (8) gets larger) and the sequence of the resulting
supersets $\\{\mathcal{S}_{k}\\}$ is all compact sets. We also achieve a
monotonic convergent upper bound on the optimal objective value of (2) (i.e.,
$f(x_{k})\geq f(x_{k+1})\geq f^{*}$, where $f^{*}$ is the optimal objective
value of (2)). In Section 3.3 we will show that Algorithm 4 generates solution
iterates that converge to the solution of (2) with the proposed cutting plane
types.

In addition to the upper bound, a lower bound on the optimal objective value
can be computed by solving a problem like (6), whose sample sets are collected
from the Euclidean projection $z_{k}$ in the error measure $errm(x_{k})$ (21).
Alternatively, we can combine the superset algorithm with any effective outer
approximation methods to obtain faster lower bound convergence (e.g., dis1 ;
tut4 ; outer1 ). This combined algorithm will have a decreasing objective gap
to terminate with.

Compared with Polak’s algorithm, the superset algorithm is capable of
providing both lower and upper bounds to the optimal objective value of (2),
whereas Polak’s algorithm provides only a lower bound. In terms of
feasibility, the superset algorithm generates iterates that are feasible in
(2) whereas the iterates of Polak’s algorithm are generally infeasible in (2)
before convergence. This allows us to define a superset algorithm with early
terminations in practical applications if needed. In terms of assumption
requirements, the superset algorithm additionally requires convexity of the
uncertainty set $\mathcal{U}$, so that the cutting plane method can be used to
update supersets over iteration.

##  4 Convergence Analysis

We prove the convergence for the sequence of the NRO subproblem solutions
generated by Algorithm 4. We develop a unified convergence proof for Algorithm
4 with each of the proposed cutting plane types in Section 3.3. This implies
that a hybrid cutting plane strategy also converges.

######  Proposition 3.

In Algorithm 4 with Kelley’s cutting plane, (17), Euclidean projection cut
(18), or gradient-free cut (19), there exists a convergent subsequence indexed
by $\\{k_{p}\\}$ and $u^{*}\in\mathcal{U}$ such that
$\\{u_{k_{p}}\\}\rightarrow u^{*}$.

Proof: From Algorithm 4, we obtain a sequence of solutions $\\{x_{k}\\}$ and
uncertainty realizations $\\{u_{k}\\}$. We define the corresponding Euclidean
projection sequence $\\{z_{k}\\}$, where
$z_{k}=\text{proj}_{\mathcal{U}}(u_{k})$ for all $k$. We claim that there
exists a convergent subsequence, indexed by $\\{k_{p}\\}$, and
$u^{*}\in\mathcal{U}$ such that $\\{u_{{k_{p}}}\\}\rightarrow u^{*}$.

Next, we consider each cutting plane in turn and prove the claim by
contradiction. If the convergence does not occur, then there exists $r>0$
(independent of $k$) such that for any given $k$ and $t$ and $0\leq t\leq
k-1$, we have the following:

  1. 1.

Kelley’s cutting plane (17):

$\displaystyle r$ $\displaystyle\leq\max_{j\in J}g_{j}(u_{t})\leq\max_{j\in
J}\nabla g_{j}(u_{t})^{\top}(u_{t}-u_{k})\leq L\|u_{t}-u_{k}\|_{2},$

where the second inequality comes from the fact that $u_{k}$ satisfies
Kelley’s cutting plane at iteration $t$ and the following inequality comes
from the Lipschitz continuity of $\nabla g_{j}$. Then, it follows that for
every subsequence (i.e., $k_{l}<k_{m}$), we have

$\displaystyle\|u_{{k_{l}}}-u_{k_{m}}\|_{2}\geq r/L.$ (25)

  2. 2.

Euclidean projection cut (18):

$\displaystyle r\leq\|u_{t}-z_{t}\|_{2}^{2}$
$\displaystyle=(u_{t}-u_{k})^{\top}(u_{t}-z_{t})+(u_{k}-z_{t})^{\top}(u_{t}-z_{t})$
(26a) $\displaystyle\leq(u_{t}-u_{k})^{\top}(u_{t}-z_{t})$ (26b)
$\displaystyle\leq\|u_{t}-u_{k}\|_{2}\|u_{t}-z_{t}\|_{2}$ (26c)
$\displaystyle\leq 2R\|u_{t}-u_{k}\|_{2},$ (26d)

where (26d) comes from the fact that $u_{t},z_{t}\in\mathcal{S}_{0}$. By
construction, $\mathcal{S}_{0}$ is a compact set and hence bounded with radius
$R$. Next, we show (26b) by proving $(u_{t}-z_{t})^{\top}(u_{k}-z_{t})\leq 0$.
For all $j\in\hat{A}_{t}=\\{j\in J:g_{j}(z_{t})=0\\}$ , we have

$\displaystyle g_{j}(z_{t})=0\ \Rightarrow\ \nabla
g_{j}(z_{t})^{\top}(u_{k}-z_{t})\leq 0,$ (27)

because $u_{k}$ satisfies the Euclidean projection cut (18) at iteration $t$.
Next, based on the KKT condition of (5), we have

$\displaystyle 2(z_{t}-u_{t})+\sum_{j\in\hat{A}_{t}}\lambda_{jt}\nabla
g_{j}(z_{t})=0,$ (28)

where $\lambda_{jt}\geq 0$ are the optimal multipliers for all
$j\in\hat{A}_{t}$. Then, from (27) and (28), we get

$\displaystyle
2(u_{t}-z_{t})^{\top}(u_{k}-z_{t})=\sum_{j\in\hat{A}_{t}}\lambda_{jt}\nabla
g_{j}(z_{t})^{\top}(u_{k}-z_{t})\leq 0.$

Next, from (26d), it follows that for every subsequence (i.e., $k_{l}<k_{m}$),
we have

$\displaystyle\|u_{{k_{l}}}-u_{{k_{m}}}\|_{2}\geq r/2R.$ (29)

  3. 3.

Gradient-free cut (19):

$\displaystyle r\leq\|u_{t}-z_{t}\|_{2}^{2}$
$\displaystyle=(u_{t}-u_{k})^{\top}(u_{t}-z_{t})+(u_{k}-z_{t})^{\top}(u_{t}-z_{t})$
(30a) $\displaystyle\leq(u_{t}-u_{k})^{\top}(u_{t}-z_{t})$ (30b)
$\displaystyle\leq\|u_{t}-u_{k}\|_{2}\|u_{t}-z_{t}\|_{2}$ (30c)
$\displaystyle\leq 2R\|u_{t}-u_{k}\|_{2},$ (30d)

where (30b) comes from the fact that $u_{k}$ satisfies the cutting plane
generated at iteration $t$, that is,

$\displaystyle(u_{k}-z_{t})^{\top}(u_{t}-z_{t})\leq 0,$

and (30d) comes from the fact that $u_{t},z_{t}\in\mathcal{S}_{0}$ and
$\mathcal{S}_{0}$ is compact and hence bounded by construction. Then, it
follows that for every subsequence (i.e., $k_{l}<k_{m}$), we have

$\displaystyle\|u_{{k_{l}}}-u_{{k_{m}}}\|_{2}\geq r/2R.$ (31)

Inequalities (25), (29), and (31) imply that $\\{u_{k}\\}$ does not contain a
Cauchy subsequence, which contradicts the fact that
$\\{u_{k}\\}\in\mathcal{S}_{0}$ because the compactness of $\mathcal{S}_{0}$
implies that any sequence in $\mathcal{S}_{0}$ contains a Cauchy subsequence.
Therefore, $\\{u_{k}\\}$ contains a convergent subsequence with indices
$\\{k_{p}\\}$ that converges to a point $u^{*}\in\mathcal{U}$.

######  Theorem 3.

With the convergent subsequence (with indices $\\{k_{p}\\}$) generated by
Algorithm 4, we obtain the following three conclusions:

  1. 1.

The sequence of $\\{\|z_{{k_{p}}}-u_{{k_{p}}}\|_{2}\\}$ converges to $0$.

  2. 2.

The sequence of $\\{z_{{k_{p}}}\\}$ converges to $u^{*}$.

  3. 3.

There exist $x^{*}\in\mathcal{F}$ and a subsequence of $\\{x_{k_{p}}\\}$ with
indices $\\{k_{l}\\}$ such that $\\{x_{k_{l}}\\}$ converges to $x^{*}$ and
$\\{errm(x_{k_{l}})\\}$ converges to 0.

Proof: We have $\\{\|z_{{k_{p}}}-u_{{k_{p}}}\|_{2}\\}\rightarrow 0$ because

$\displaystyle
0\leq\|z_{{k_{p}}}-u_{{k_{p}}}\|_{2}\leq\|u_{{k_{p}}}-u^{*}\|_{2}.$

The inequality holds because $z_{k}$ is the Euclidean projection of $u_{k}$
onto the uncertainty set $\mathcal{U}$, and we have
$\\{z_{{k_{p}}}\\}\rightarrow u^{*}$ because

$\displaystyle 0$
$\displaystyle\leq\|z_{{k_{p}}}-u^{*}\|_{2}\leq\|z_{{k_{p}}}-u_{{k_{p}}}\|_{2}+\|u_{{k_{p}}}-u^{*}\|_{2}.$

Because $\mathcal{F}$ is compact and $x_{k_{p}}\in\mathcal{F}$ for all
$k_{p}$, there exists a subsequence with indices $\\{k_{l}\\}$ and
$x^{*}\in\mathcal{F}$ such that $\\{x_{k_{l}}\\}\rightarrow x^{*}$. Further,
we have $\\{errm(x_{k_{l}})\\}\rightarrow 0$ because of conclusion 1.

Following the third conclusion in Theorem 3, we conclude that with any cutting
planes (17), (18), or (19), Algorithm 4 converges to a stationary solution of
(2).

######  Remark:.

Comparing the convergence analysis of Polak’s algorithm polak_boyd  and our
superset algorithm, we see that the main difference is that Polak’s algorithm
uses an assumption that problem (7) is solved to global optimality. In our
superset algorithm, because we aim to improve the supersets by removing the
violated points, we assume that $\mathcal{U}$ is a convex set to take
advantage of the cutting plane methods. Note that convexity is a sufficient
condition for global optimality, which means that the two sets of assumptions
are not significantly different.

##  5 Numerical Results

We compare the computational performance of the superset algorithm and Polak’s
algorithm using the following applications: portfolio optimization and
production cost minimization.

###  5.1 Simulation Setup

We implement Polak’s and the superset algorithms using Julia 1.6.5 and run the
numerical studies on a Linux system with a single-thread, 2.20 GHz Intel(R)
Xeon(R) CPU E7-8890 and 16 GB RAM. We use the filterSQP solver filtersqp  for
all the smooth nonlinear optimization problems and Mosek 9.3 mosek  for the
semi-definite programming problem in portfolio optimization. In the superset
algorithm, we use the Euclidean projection cut and the gradient-free cut. For
termination tolerance, we use $\epsilon=10^{-5}$.

###  5.2 Robust Asset Allocation

We consider the following multiperiod asset allocation problem with
transaction costs port1 ; port2 ; port_horizon1 ; port_horizon2 . We define
$n$ to be the number of assets and the decision horizon as
$t\in\\{1,...,T\\}$. Coefficient $U$ denotes the time-invariant linear
transaction cost that is independent of the asset type. For each period $t$,
$x_{t}\in\mathbb{R}^{n}$ represents the allocation decision for each asset,
and $\lambda_{t}>0$ is a predefined risk aversion parameter. Parameters
$\mu_{t}\in\mathbb{R}^{n}$ and $Q_{t}\in\mathbb{R}^{n\times n}$ represent the
uncertain mean and covariance of the asset return with uncertainty sets
$\mathcal{U}_{\mu t}$ and $\mathcal{U}_{Qt}$. Problem (32) maximizes the total
profit from total risk-adjusted expected return $s_{t}\in\mathbb{R}$ excluding
the transaction cost $c_{t}\in\mathbb{R}$ under the uncertainty of mean and
covariance. The left-hand side of (32b) gives the risk-adjusted expected
return $s_{t}$ in period $t$. The left-hand side of (32d) gives the
transaction cost $c_{t}$ associated with the changes on allocation decision
$x_{t}$. The problem contains $T$ robust constraints with $n(n+1)$ uncertainty
variables per robust constraint.

$\displaystyle\max_{x_{t}\geq 0,\ s_{t},\ c_{t}\geq 0}\quad$
$\displaystyle\sum_{t=1}^{T}(s_{t}-c_{t})$ (32a) s.t
$\displaystyle\mu_{t}^{\top}x_{t}-\lambda x_{t}^{\top}Q_{t}x_{t}\geq s_{t},\
\forall\ \mu_{t}\in\mathcal{U}_{\mu_{t}},\ \forall
Q_{t}\in\mathcal{U}_{Q_{t}},\ \forall t\in\\{1,...,T\\},$ (32b)
$\displaystyle\sum_{i=1}^{n}x_{it}=1,\ \forall t\in\\{1,...,T\\},$ (32c)
$\displaystyle U\|x_{t+1}-x_{t}\|_{1}\leq c_{t},\ \forall t\in\\{1,...,T\\},.$
(32d)

Further, we define uncertainty sets $\mathcal{U}_{\mu
t}:=\\{\underline{\mu}_{it}\leq\mu_{it}\leq\overline{\mu}_{it},\ i=1,...,n\\}$
and $\mathcal{U}_{Qt}:=\\{{\color[rgb]{0,0,0}Q\succeq 0},\
\sum_{i,j}\frac{(Q_{ijt}-C_{ijt})^{2}}{r_{ijt}^{2}}\leq 1,\ i=1,...,n,\
j=1,...,n,\\}$. $C_{ijt}$ and $r_{ijt}$ are the corresponding center and
radius for entry $Q_{ijt}$. The uncertainty set $\mathcal{U}_{\mu t}$ and
$\mathcal{U}_{Qt}$ are constructed by using the Australian stock price dataset
port_data . In our numerical study we choose three cases among
$T\in\\{7,14,21\\}$ and $n=2$ (i.e., a total uncertainty dimension of $6T$)
with transaction cost selection $U\in\\{0.05,0.15,0.25,0.35\\}$. Because
$\lambda_{t}$ does not interfere with the problem dimension, we use
$\lambda_{t}=1$ for simplicity in these computational tests. For each given
horizon and asset setting, we randomly select $3$ groups of assets and
demonstrate the average computational performance across the groups.

Table 2: Average iteration and runtime of Polak’s and the superset algorithms
for $T=\\{7,14,21\\}$.

(a) Runtime (Sec). | Horizon
---|---
| $T=7$ | $T=14$ | $T=21$
| Transaction
---
Cost $U$
| Polak’s
---
Algorithm
| Superset
---
Algorithm
| Polak’s
---
Algorithm
| Superset
---
Algorithm
| Polak’s
---
Algorithm
| Superset
---
Algorithm
0.05 | 1.23 | 1.55 | 2.22 | 1.82 | 4.48 | 1.83
0.15 | 0.20 | 0.13 | 1.08 | 0.31 | 2.63 | 0.35
0.25 | 0.20 | 0.11 | 1.08 | 0.22 | 2.54 | 0.58
0.35 | 0.08 | 0.06 | 0.27 | 0.13 | 0.68 | 0.28

(b) Iteration Count. | Horizon
---|---
| $T=7$ | $T=14$ | $T=21$
| Transaction
---
Cost $U$
| Polak’s
---
Algorithm
| Superset
---
Algorithm
| Polak’s
---
Algorithm
| Superset
---
Algorithm
| Polak’s
---
Algorithm
| Superset
---
Algorithm
0.05 | 3.0 | 7.0 | 4.3 | 5.7 | 7.0 | 4.0
0.15 | 3.0 | 3.3 | 3.0 | 3.7 | 3.7 | 2.3
0.25 | 3.0 | 3.3 | 4.7 | 2.3 | 4.7 | 4.0
0.35 | 3.0 | 1.3 | 3.0 | 1.3 | 3.0 | 2.7

The results are summarized in Table 2(b). We observe that Polak’s algorithm
requires more iterations and runtime compared to the superset algorithm when
the horizon increases (i.e., more robust constraints in (32)). As the number
of robust constraints increases, more nonlinear constraints are added in
Polak’s subproblem (6) as iteration goes, while the constraint dimension of
the superset algorithm remains the same. As a result, when the horizon
increases, the superset algorithm demonstrates a steady computational
performance and an increasing advantage over the Polak’s algorithm. In terms
of the transaction cost change, when $U$ increases, we observe reductions on
the iteration count and runtime for both algorithms.

To summarize the computational difference between Polak’s algorithm and the
superset algorithm, we provide the performance profile performprofile  of the
CPU runtime in Fig. 2). We conclude that the superset algorithm outperforms
Polak’s algorithm in $92\%$ of the instances. Specifically, for $52\%$ of the
cases, the superset algorithm is more than 2.0 times faster than the Polak’s
algorithm; for $27\%$ of the cases, the superset algorithm is more than 3
times faster; for $17\%$ of the cases, the superset algorithm is more than 4
times faster.

Figure 2: Performance profile of CPU runtime comparison: the superset
algorithm wins $92\%$ of the multiperiod portfolio optimization instances.

###  5.3 Power System Application: Production Cost

We consider a production cost minimization problem for a given horizon
$t\in\\{1,...,T\\}$ between two generators with well-forecast load $d_{t}$ and
uncertain cost prod_cost . We assume the first generator is expensive with
ramping costs and limits, whereas the second generator is cheap and has no
ramping limitation. For each period $t$, we denote the productions of the
first generator as $x_{t}$ with cost $c_{1,t}$ and the ramp rate $u_{t}$ with
control cost $c_{3,t}$. Because of power balance at each period $t$, we denote
the production of the second generator as $d_{t}-x_{t}$ with cost $c_{2,t}$.
For each $t$, we further define the uncertainty set as $\mathcal{C}_{t}$. The
full problem is formulated as follows:

$\displaystyle\min_{x_{t},\ u_{t},\ s_{t}}\quad$
$\displaystyle\sum_{t=1}^{T}s_{t}$ (33a) s.t $\displaystyle
c_{1,t}(x_{t}-d_{t})^{2}+c_{2,t}x_{t}^{2}$
$\displaystyle\quad\quad+c_{3,t}u_{t}^{2}\leq s_{t},\
\forall(c_{1,t},c_{2,t},c_{3,t})\in\mathcal{C}_{t},\forall t\in\\{1,...,T\\}$
(33b) $\displaystyle x_{t+1}=x_{t}+u_{t},\ \forall t\in\\{1,...,T\\},$ (33c)
$\displaystyle-U\leq u_{t}\leq U,\ \forall t\in\\{1,...,T\\},$ (33d)

where the left-hand side of (33b) represents the cost at time $t$ as $s_{t}$
and (33c) and (33d) represent the ramp dynamic and limit $U$, respectively. In
the simulation we fix $T=24$, and the problem contains $73$ variables with
$24$ robust constraints and $48$ deterministic constraints.

For $d_{t}$, we use the seven-day load forecast dataset from PJM pjm_data  and
select three daily load-forecast patterns from Dec. 2021 (see Fig. 3). For
uncertainty sets $\mathcal{C}_{t}$, we first define the nominal costs
$o_{1}=10$, $o_{2}=5$, $o_{3}=1$ and then construct ellipsoidal uncertainty
sets as
$\frac{(c_{1,t}-o_{1})^{2}}{r_{1}^{2}}+\frac{(c_{2,t}-o_{2})^{2}}{r_{2}^{2}}+\frac{(c_{3,t}-o_{3})^{2}}{r_{3}^{2}}\leq
1$ with $r_{i}=po_{i},\ i=1,2,3$, where $p\in(0,1]$ is a percentage parameter
controlling the size of the ellipsoid. As the comparison setup, we use 3 load
patterns, 3 different values of $U\in\\{80,120,160\\}$, and 3 different values
of $p\in\\{0.7,0.8,0.9\\}$ and in total run 27 instances with every
combination.

Figure 3: Three daily load patterns used in the testing.

First, we show the iteration count comparison. The iteration count of all
instances is summarized in Table 3. In all cases, the superset algorithm takes
fewer iterations than Polak’s algorithm does. Both methods have stable
iteration count regardless of the changes on load pattern, ramp limits, and
the size of the uncertainty set.

Table 3: Iteration count comparison for all testing instances.

|  | Ramp Rate $U$
---|---|---
|  | $80$ | $120$ | $160$
| Load
---
Pattern
| Uncertainty
---
Set Size $p$
| Polak’s
---
Algorithm
| Superset
---
Algorithm
| Polak’s
---
Algorithm
| Superset
---
Algorithm
| Polak’s
---
Algorithm
| Superset
---
Algorithm
1 | $0.7$ | 10 | 6 | 10 | 6 | 10 | 6
1 | $0.8$ | 10 | 6 | 10 | 6 | 10 | 6
1 | $0.9$ | 10 | 6 | 10 | 6 | 10 | 6
2 | $0.7$ | 10 | 6 | 10 | 6 | 10 | 6
2 | $0.8$ | 10 | 6 | 10 | 6 | 10 | 6
2 | $0.9$ | 10 | 6 | 10 | 6 | 10 | 6
3 | $0.7$ | 10 | 6 | 10 | 6 | 10 | 6
3 | $0.8$ | 10 | 6 | 10 | 6 | 10 | 6
3 | $0.9$ | 10 | 6 | 10 | 6 | 10 | 6

Next, we show the runtime comparison. The runtimes for all instances are
summarized in Table 4. In all cases, superset algorithm takes less CPU time
than Polak’s algorithm does. Given the stable iteration count in Table 3,
superset algorithm also has stable CPU time regardless of the changes on test
settings. On the other hand, the CPU time of Polak’s algorithm is greatly
dependent on the load patterns, ramp limits and the size of uncertainty set
but no general monotonic relation can be observed.

Table 4: Runtime (sec) comparison for all testing instances.

|  | Ramp Rate $U$
---|---|---
|  | $80$ | $120$ | $160$
| Load
---
Pattern
| Uncertainty
---
Set Size $p$
| Polak’s
---
Algorithm
| Superset
---
Algorithm
| Polak’s
---
Algorithm
| Superset
---
Algorithm
| Polak’s
---
Algorithm
| Superset
---
Algorithm
1 | $0.7$ | 0.767 | 0.748 | 1.653 | 0.744 | 1.430 | 0.742
1 | $0.8$ | 0.873 | 0.715 | 1.913 | 0.725 | 0.972 | 0.716
1 | $0.9$ | 0.806 | 0.731 | 1.271 | 0.760 | 1.000 | 0.750
2 | $0.7$ | 1.230 | 0.715 | 1.095 | 0.743 | 1.704 | 0.737
2 | $0.8$ | 1.192 | 0.703 | 1.244 | 0.711 | 1.480 | 0.706
2 | $0.9$ | 0.766 | 0.755 | 1.313 | 0.749 | 1.323 | 0.750
3 | $0.7$ | 1.703 | 0.721 | 1.139 | 0.721 | 1.126 | 0.727
3 | $0.8$ | 0.908 | 0.708 | 1.617 | 0.696 | 1.301 | 0.684
3 | $0.9$ | 1.112 | 0.748 | 1.019 | 0.745 | 1.729 | 0.740

To summarize the computational difference between Polak’s algorithm and the
superset algorithm, we provide the performance profile performprofile  of the
CPU runtime in Fig. 4. We observe that, in all cases, the superset algorithm
outperforms Polak’s algorithm. Specifically, for $59\%$ of the cases, superset
algorithm is more than 1.5 times faster than the Polak’s algorithm; for $22\%$
of the cases, superset algorithm is more than 2 times faster than the Polak’s
algorithm.

Figure 4: Performance profile of CPU runtime comparison: the superset
algorithm is superior in all the production cost minimization instances.

##  6 Conclusion

We have developed a superset algorithm for a class of structured NRO problems.
The algorithm iteratively solves the reformulation of an NRO subproblem with
the polytopic supersets of the uncertainty set. Different cutting plane
methods are proposed to improve the supersets over iteration. We showed that
the solution iterates from the superset algorithm are feasible in the original
NRO problem and provide both lower and upper bounds to the optimal objective
value. We proved the convergence of the superset algorithm under the
assumption that uncertainty sets are convex. We also provided a feasibility
restoration algorithm to detect whether the NRO is infeasible or restore the
feasibility of the NRO subproblem of the superset algorithm by constructing a
new superset.

To evaluate the computational performance, we compared the superset algorithm
with Polak’s algorithm in applications including portfolio optimization and
production cost minimization. We demonstrated that the superset algorithm is
more advantageous than Polak’s algorithm when the number of robust constraints
is large.

For future work, we plan to extend the superset algorithm to more general NRO
formulations by relaxing the current structural assumptions of the affine
relationship between the constraint and uncertainty parameter.

Acknowledgments

This work was supported by the Applied Mathematics activity within the U.S.
Department of Energy, Office of Science, Advanced Scientific Computing
Research, under Contract DE AC02-06CH11357.

## References

  *   * (1) Tütüncü, R.H., Koenig, M.: Robust asset allocation. Annals of Operations Research 132, 157–187 (2004)
  * (2) Lobo, M.S., Boyd, S.P.: The worst-case risk of a portfolio. Technical report, Stanford University (2000)
  * (3) Xu, W., Anitescu, M.: Exponentially accurate temporal decomposition for long-horizon linear-quadratic dynamic optimization. SIAM Journal on Optimization 28(3), 2541–2573 (2018)
  * (4) Polak, E., He, L.: Rate-preserving discretization strategies for semi-infinite programming and optimal control. SIAM Journal on Control and Optimization 30(3), 548–572 (1992)
  * (5) Goberna, M.Á., López, M.A.: Semi-infinite programming: recent advances. Nonconvex Optimization and its Applications (2001)
  * (6) Reemtsen, R.: Some outer approximation methods for semi-infinite optimization problems. Journal of Computational and Applied Mathematics 53(1), 87–108 (1994)
  * (7) Polak, E.: Optimization. Springer, New York (1997)
  * (8) Mutapcic, A., Boyd, S.: Cutting-set methods for robust convex optimization with pessimizing oracles. Optimization Methods and Software 24(3), 381–406 (2009)
  * (9) Pätzold, J., Schöbel, A.: Approximate cutting plane approaches for exact solutions to robust optimization problems. European Journal of Operational Research 284(1), 20–30 (2020)
  * (10) Hettich, R., Kortanek, K.O.: Semi-infinite programming: Theory, methods, and applications. SIAM Review 35(3), 380–429 (1993)
  * (11) Fiacco, A.V., Kortanek, K.O.: Semi-Infinite Programming and Applications: An International Symposium Austin, Texas, September 8–10, 1981 vol. 215. Springer, Berlin (2012)
  * (12) Reemtsen, R., Rückmann, J.-J.: Semi-infinite Programming vol. 25. Springer, New York (1998)
  * (13) Ben-Tal, A., Den Hertog, D., Vial, J.-P.: Deriving robust counterparts of nonlinear uncertain inequalities. Mathematical Programming 149(1), 265–299 (2015)
  * (14) Ben-Tal, A., El Ghaoui, L., Nemirovski, A.: Robust Optimization. Princeton university press, New Jersey (2009)
  * (15) Stein, O.: How to solve a semi-infinite optimization problem. European Journal of Operational Research 223(2), 312–320 (2012)
  * (16) Jongen, H.T., Twilt, F., Weber, G.W.: Semi-infinite optimization structure and stability of the feasible set. Journal of Optimization Theory and Applications 72(3), 529–552 (1992)
  * (17) Stein, O.: On constraint qualifications in nonsmooth optimization. Journal of Optimization Theory and Applications 121, 647–671 (2004)
  * (18) Stein, O., Steuermann, P.: The adaptive convexification algorithm for semi-infinite programming with arbitrary index sets. Mathematical Programming 136(1), 183–207 (2012)
  * (19) Leyffer, S., Menickelly, M., Munson, T., Vanaret, C., Wild, S.M.: A survey of nonlinear robust optimization. INFOR: Information Systems and Operational Research 58(2), 342–373 (2020)
  * (20) Kelley, J.E. Jr.: The cutting-plane method for solving convex programs. Journal of the Society for Industrial and Applied Mathematics 8(4), 703–712 (1960)
  * (21) Bonami, P., Cornuéjols, G., Lodi, A., Margot, F.: A feasibility pump for mixed integer nonlinear programs. Mathematical Programming 119(2), 331–352 (2009)
  * (22) Cheney, E.W., Goldstein, A.A.: Tchebycheff approximation and related extremal problems. Journal of Mathematics and Mechanics 14(1), 87–98 (1992)
  * (23) Fletcher, R., Leyffer, S.: Nonlinear programming without a penalty function. Mathematical Programming 91(2), 239–269 (2002)
  * (24) ApS, M.: MOSEK Optimizer API for C 9.3.21. (2022). https://docs.mosek.com/latest/capi/index.html#project
  * (25) Lobo, M.S., Fazel, M., Boyd, S.P.: Portfolio optimization with linear and fixed transaction costs. Annals of Operations Research 152, 341–365 (2007)
  * (26) Boyd, S., Busseti, E., Diamond, S., Kahn, R.N., Koh, K., Nystrup, P., Speth, J., et al.: Multi-period trading via convex optimization. Foundations and Trends in Optimization 3(1), 1–76 (2017)
  * (27) Bellett, A.: Australian Historical Stock Prices. https://www.kaggle.com/code/ashbellett/portfolio-optimisation/data. Accessed: 2022-07-04
  * (28) Dolan, E.D., Moré, J.J.: Benchmarking optimization software with performance profiles. Mathematical Programming 91(2), 201–213 (2002)
  * (29) PJM: Seven-Day Load Forecast. https://dataminer2.pjm.com/feed/load_frcstd_7_day/definition. Accessed: 2022-07-04