Polytopic Superset Algorithm for Nonlinear Robust Optimization

TL;DR

This paper introduces a novel polytopic superset algorithm for nonlinear robust optimization that guarantees feasible iterates and bounds on the optimal value, improving over traditional methods especially with many constraints.

Contribution

The proposed algorithm generates feasible solutions and bounds for NRO problems using polytopic supersets and a cutting plane method, with proven convergence to the optimum.

Findings

The superset algorithm converges to the optimal solution.

It outperforms Polak's outer approximation method with many constraints.

Numerical results demonstrate improved efficiency and feasibility.

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.

Tables

Table 1. Table 1 : Solution pair ( x k , u k , λ k ) subscript 𝑥 𝑘 subscript 𝑢 𝑘 subscript 𝜆 𝑘 (x_{k},u_{k},\lambda_{k}) and cutting planes obtained in the first three iterations of Algorithm  2 (i.e., k = 0 , 1 , 2 𝑘 0 1 2 k=0,1,2 ) where a 1 = 4 ​ 6 + 12 − 5 − 2 ​ 10 7 subscript 𝑎 1 4 6 12 5 2 10 7 a_{1}=\frac{4\sqrt{6}+\sqrt{12}-\sqrt{5}-2\sqrt{10}}{7} and a 2 = ( 10 − 3 ) ​ ( 1 + 2 ​ 2 ) 7 subscript 𝑎 2 10 3 1 2 2 7 a_{2}=\frac{(\sqrt{10}-\sqrt{3})(1+2\sqrt{2})}{7} .

Iteration $k$	$x_k$	$u_k$	${\lambda }_k$	$k$th Cutting Plane
$0$	$(2,1)$	$(1,2)$	$1/4$	$u_1+2u_2\le 3$
$1$	$(\sqrt{3},\sqrt{3})$	$(1,1)$	$\sqrt{3}/6$	$u_1+u_2\le 1.5$
$2$	$(2,2)$	$(0.75,0.75)$	$1/3$	$u_1+u_2\le 443/300$

Table 2. (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\le 3$
$1$	$(\sqrt{3},\sqrt{3})$	$(1,1)$	$\sqrt{3}/6$	$u_1+u_2\le 1.5$
$2$	$(2,2)$	$(0.75,0.75)$	$1/3$	$u_1+u_2\le 443/300$

Table 3. (b) Euclidean Projection Cut.

Iteration $k$	$x_k$	$u_k$
0	$(2,1)$	$(1,2)$
1	$(\frac{\sqrt{30\sqrt{5}}}{5},\frac{2\sqrt{15\sqrt{5}}}{5})$	$(\sqrt{5}\left(\sqrt{2}-1\right),\frac{\sqrt{10}\left(\sqrt{2}-1\right)}{2})$
2	$(\sqrt{\frac{6a_2}{a_1{a}_2+a_1^2}},\sqrt{\frac{6a_1}{a_1{a}_2+a_2^2}})$	$(a_1,a_2)$

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

	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

Table 5. (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

Table 6. (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

Table 7. 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

Table 8. 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

Equations

$$\underset{x}{\min }f(x)\text{mbox}s.t.u_i^{\top }{h}_i(x)\le b_i(x),\forall u_i\in {\mathcal{U}}_i,i=1,2,...,I,$$

$$\underset{x}{\min }f(x)\text{mbox}s.t.u^{\top }h(x)\le b(x),\forall u\in \mathcal{U},$$

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

$$\nabla f(x^{\ast })+\sum _{u\in A^{\ast }}{\lambda }_u^{\ast }\left(\nabla h(x^{\ast })u-\nabla b(x^{\ast })\right)=0.$$

$$errm(x,{\lambda }_u,A)$$

$$+{\left[\underset{u\in \mathcal{U}}{\max }{u}^{\top }h(x)-b(x)\right]}^{+}$$

$$+\sum _{u\in A}[-{\lambda }_u{]}^{+}$$

$$+\sum _{u\in A}|u^{\top }h(x)-b(x)|+\sum _{u\in A}\Vert u-{\text{proj}}_{\mathcal{U}}(u){\Vert }_2.$$

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

$$\underset{x}{\min }f(x)\text{mbox}s.t.u^{\top }h(x)\le b(x),\forall u\in U_k.$$

$$\underset{t,u\in \mathcal{U}}{\max }t\text{mbox}s.t.u^{\top }h(x_k)-b(x_k)\ge t.$$

$$\underset{x}{\min }f(x)\text{mbox}s.t.u^{\top }h(x)\le b(x),\forall u\in {\mathcal{S}}_k.$$

$$\displaystyle\min_{x,\gamma_{k}}\$$

$${\gamma }_k^{\top }{d}_k\le b(x),$$

$$(B_k{)}^{\top }{\gamma }_k=h(x),$$

$${\gamma }_k\ge 0.$$

$$\underset{u\in {\mathcal{S}}_k}{\max }{u}^{\top }h(x_k)\le b(x_k)\iff u^{\top }h(x_k)\le b(x_k),\forall u\in {\mathcal{S}}_k.$$

$$\underset{\gamma \ge 0,(B_k{)}^{\top }\gamma =h(x_k)}{\min }{\gamma }^{\top }{d}_k=\underset{u\in {\mathcal{S}}_k}{\max }{u}^{\top }h(x_k).$$

$$\nabla f(x_k)+\nabla h(x_k)v_k-{\mu }_k\nabla b(x_k)=0,$$

$$-B_k{v}_k+{\mu }_k{d}_k-{\nu }_k=0,$$

$${\mu }_k\ge 0,{\nu }_k\ge 0,{\mu }_k\left({\gamma }_k^{\top }{d}_k-b(x_k)\right)=0,{\nu }_k^{\top }{\gamma }_k=0,$$

$${\gamma }_k^{\top }{d}_k\le b(x_k),{\gamma }_k\ge 0,B_k^{\top }{\gamma }_k=h(x_k).$$

$$\text{eqref}stt3$$

$$\text{eqref}stt2\text{mbox}and\text{eqref}stt4$$

$$\text{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),$$

$$\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$},$$

$$\displaystyle\nabla f(x_{k})=0,\quad\mbox{if $\mu_{k}=0$}.$$

$${\lambda }_k=\{\begin{array}{ll}{\mu }_k>0,& \text{if }{\mu }_k>0\\ 0,& \text{otherwise}\end{array},u_k=\{\begin{array}{ll}\frac{v_k}{{\mu }_k},& \text{if }{\mu }_k>0\\ \text{any }u\in \mathcal{U},& \text{otherwise.}\end{array}$$

$$\nabla g(u_k{)}^{\top }u\le \nabla g(u_k{)}^{\top }{u}_k-g(u_k).$$

$$\nabla g(z_k{)}^{\top }u\le \nabla g(z_k{)}^{\top }{z}_k-g(z_k).$$

$$(u_k-z_k{)}^{\top }u\le (u_k-z_k{)}^{\top }{z}_k.$$

$$\underset{x_1,x_2}{\min }$$

$$x_1^2{u}_1+x_2^2{u}_2\le 6,\forall (u_1,u_2{)}^{\top }\in \mathcal{U},$$

$$errm(\bar{x},\bar{\lambda },\{\bar{u}\})=\Vert \bar{u}-{\text{proj}}_{\mathcal{U}}(\bar{u}){\Vert }_2.$$

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

$$\underset{x,p\ge 0}{\min }p\text{mbox}s.t.u^{\top }h(x)\le b(x)+p,\forall u\in \mathcal{U}.$$

Full text

1 Introduction

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

[TABLE]

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:

[TABLE]

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):

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

[TABLE]

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

[TABLE]

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$:

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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.

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:

[TABLE]

Based on (8), we propose the following 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:

[TABLE]

**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

[TABLE]

Next, by linear programming duality, we have

[TABLE]

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:

[TABLE]

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

[TABLE]

Further, we have

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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)$.

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).

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:

[TABLE]

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:

[TABLE]

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):

[TABLE]

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

[TABLE]

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.

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})$.

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:

Kelley’s cutting plane (17):

[TABLE]

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

[TABLE] 2. 2.

Euclidean projection cut (18):

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE] 3. 3.

Gradient-free cut (19):

[TABLE]

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

[TABLE]

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

[TABLE]

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:

The sequence of $\{\|z_{{k_{p}}}-u_{{k_{p}}}\|_{2}\}$ converges to [math]. 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

[TABLE]

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

[TABLE]

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.

[TABLE]

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.

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.

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:

[TABLE]

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.

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.

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.

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.

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.

\bmhead

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.