Robust optimization for the pooling problem

TL;DR

This paper explores robust optimization techniques for the non-linear, non-convex pooling problem, accounting for parametric uncertainty, and compares reformulation and cutting plane methods on multiple instances.

Contribution

It introduces the application of robust optimization approaches to the pooling problem with uncertain parameters, a novel extension beyond traditional deterministic methods.

Findings

Reformulation and cutting plane methods show different computational efficiencies.

Accounting for uncertainty alters the optimal solutions.

Robust optimization is feasible for non-convex, global optimization problems.

Abstract

The pooling problem has applications, e.g., in petrochemical refining, water networks, and supply chains and is widely studied in global optimization. To date, it has largely been treated deterministically, neglecting the influence of parametric uncertainty. This paper applies two robust optimization approaches, reformulation and cutting planes, to the non-linear, non-convex pooling problem. Most applications of robust optimization have been either convex or mixed-integer linear problems. We explore the suitability of robust optimization in the context of global optimization problems which are concave in the uncertain parameters by considering the pooling problem with uncertain inlet concentrations. We compare the computational efficiency of reformulation and cutting plane approaches for three commonly-used uncertainty set geometries on 14 pooling problem instances and demonstrate how accounting for uncertainty changes the optimal solution.

Tables

Table 1. Table 1: Robust reformulations for box, ellipsoidal, and polyhedral uncertainty sets.

Uncertainty set geometry	$p$	${\lambda }_{j,k}$
Box	$\infty$	$\sum _i{\hat{C}}_{i,k}\cdot x_{i,j}$
Ellipsoid	$2$	$\sqrt{\sum _i{\hat{C}}_{i,k}^2\cdot x_{i,j}^2}$
Ellipsoid (correlation)	$2$	$\sqrt{\sum _i\sum _{i^{\prime }}{\sigma }_{i,i^{\prime }}\cdot {\hat{C}}_{i,k}\cdot {\hat{C}}_{i^{\prime },k}\cdot x_{i,j}\cdot x_{i^{\prime },j}}$
Polyhedron	$1$	$\underset{i}{\max }{\hat{C}}_{i,k}{x}_{i,j}$

Table 2. Table 2: Percentage of instances solved to global optimality within 1 hr time limit and median time taken for different uncertainty set types and different robust optimization approaches. The median times were calculated for only those instances which could be solved to optimality within 1 hr.

Uncert. set	% instances solved			median time [s]
	reform.	CP - multiple	CP - one	reform.	CP - multiple	CP - one
Box	100	100	91	0.14	0.72	2.13
Ellipsoid	66	88	75	1.42	2.91	4.75
Polyhedron	83	96	79	0.15	1.13	1.77

Variables
$𝝀$	vector of dual variables
$q_{i,l}$	fraction of flow entering pool $l$ from source $i$
$v_j$	total flow to terminal $j$
$y_{l,j}$	flow from pool $l$ to terminal $j$
$x_{i,j}$	total flow from source $i$ to terminal $j$
$𝒙$	vector of decision variables
$𝝃$	vector of uncertain parameters
$z_{i,j}$	direct flow from source $i$ to terminal $j$
Parameters
$A_i^L$	lower availability limit feed $i$
$A_i^U$	upper availability limit feed $i$
$b$	right-hand-side constant
$c_i$	per unit cost of material from source $i$
${\tilde{C}}_{i,k}$	uncertain inlet concentration of component $k$ at source $i$
$C_{i,k}$	nominal inlet concentration of component $k$ at source $i$
${\hat{C}}_{i,k}$	maximum deviation of inlet concentration of component $k$ at source $i$
$\delta$	relative termination tolerance
$d_j$	per unit price of product at terminal $j$
$D_j^L$	minimum demand product $j$
$D_j^U$	maximum demand product $j$
${\ell }_i$	location vector for source $i$
$P_{j,k}^L$	minimum concentration of component $k$ at terminal $j$
$P_{j,k}^U$	maximum concentration of component $k$ at terminal $j$
$r$	uncertainty set size parameter
$S_l$	maximum capacity pool $l$
$s$	safety factor
${\sigma }_{i,i^{\prime }}$	covariance between input concentrations at sources $i$ and $i^{\prime }$
Sets
$\Lambda$	feasible set of dual problem
$T_X$	contains $(i,l)$ if source $i$ is connected to pool $l$
$T_Y$	contains $(l,j)$ if pool $l$ is connected to terminal $j$
$T_Z$	contains $(i,j)$ if source $i$ is conntected to terminal $j$
$𝒰$	uncertainty set
$\widehat{𝒰}$	set containing finite number of uncertainty scenarios
$𝒳$	certain feasible set

Equations

$$\underset{\mathbf{x}\in \mathcal{X}}{\min }$$

$$g(\mathbf{x},\mathbf{\xi })\le b$$

$$\underset{\mathbf{x}\in \mathcal{X}}{\min }$$

$$\underset{\mathbf{\xi }\in \mathcal{U}}{\max }g(\mathbf{x},\mathbf{\xi })\le b.$$

$$\underset{\mathbf{x}\in \mathcal{X}}{\min }$$

$$\underset{\mathbf{\lambda }\in \Lambda }{\min }{g}^{\ast }(\mathbf{x},\mathbf{\lambda })\le b,$$

$$\underset{\mathbf{x}\in \mathcal{X},\mathbf{\lambda }\in \Lambda }{\min }$$

$$g^{\ast }(\mathbf{x},\mathbf{\lambda })\le b.$$

$${\mathbf{x}}_n=\underset{\mathbf{x}\in \mathcal{X}}{\mathrm{argmax}}$$

$$g(\mathbf{x},\mathbf{\xi })\le b,\forall \mathbf{\xi }\in {\hat{\mathcal{U}}}_n,$$

$${\mathbf{\xi }}_{n+1}=\underset{\mathbf{\xi }\in \mathcal{U}}{\mathrm{argmax}}g({\mathbf{x}}_n,\mathbf{\xi }),$$

$$\underset{q_{i,l},y_{l,j},z_{i,j}}{\min }$$

$$\displaystyle\begin{array}[]{l}\text{Feed}\\ \text{Availab.}\end{array}$$

$$\displaystyle\begin{array}[]{l}\text{Pool}\\ \text{Capacity}\end{array}$$

$$\displaystyle\begin{array}[]{l}\text{Product}\\ \text{Demand}\end{array}$$

$$\displaystyle\begin{array}[]{l}\text{Simplex}\end{array}$$

$$\displaystyle\begin{array}[]{l}\text{Product}\\ \text{Quality}\end{array}$$

$$v_j=\sum _{l:(l,j)\in T_Y}{y}_{l,j}+\sum _{i:(i,j)\in T_Z}{z}_{i,j}\forall j.$$

$$\sum _{\begin{array}{c}l:(l,j)\forall T_Y\\ i:(i,l)\forall T_X\end{array}}{\tilde{C}}_{i,k}\cdot q_{i,l}\cdot y_{l,j}+\sum _{i:(i,j)\in T_Z}{\tilde{C}}_{i,k}\cdot z_{i,j}\{\begin{array}{l}\ge P_{j,k}^L\cdot v_j\\ \le P_{j,k}^U\cdot v_j\end{array}\forall {\tilde{C}}_{i,k}\in {\mathcal{U}}_p,\forall j,k.$$

$${\mathcal{U}}_p=\left\{{\tilde{C}}_{i,k}=C_{i,k}+{\hat{C}}_{i,k}\cdot {\xi }_{i,k}|\Vert {\mathbf{\xi }}_k{\Vert }_p\le r\right\},$$

$$k(\mathbf{\ell },{\mathbf{\ell }}^{\prime })={\sigma }^2\exp \left(\frac{\Vert \mathbf{\ell }-{\mathbf{\ell }}^{\prime }{\Vert }_2^2}{2l^2}\right)$$

$${\sigma }_{i,i^{\prime }}=k({\mathbf{\ell }}_i,{\mathbf{\ell }}_{i^{\prime }}),$$

$$\left\{{\tilde{C}}_{i,k}=C_{i,k}+{\hat{C}}_{i,k}\cdot {\xi }_{i,k}|{\mathbf{\xi }}_k^{⊺}{\Sigma }^{-1}{\mathbf{\xi }}_k\le r^2\right\}.$$

$$\text{mathbbm}1(s\in S)=\{\begin{array}{ll}1& s\in S\\ 0& s\notin S,\end{array}$$

$$\begin{array}[]{ll}\sum\limits_{i}{\tilde{C}_{i,k}\cdot x_{i,j}}\begin{cases}\geq P^{L}_{j,k}\cdot v_{j}\\ \leq P^{U}_{j,k}\cdot v_{j},\end{cases}&\forall\tilde{C}_{i,k}\in\mathcal{U}_{p},\forall j,k,\end{array}$$

$$x_{i,j}=\sum _{l:(l,j)\in T_Y}\left[\text{mathbbm}1\left((i,l)\in T_X\right)\cdot y_{l,j}\cdot q_{i,l}\right]+\text{mathbbm}1\left((i,j)\in T_Z\right)\cdot z_{i,j}.$$

$$\sum _i{C}_{i,k}\cdot x_{i,j}\{\begin{array}{l}\le P_{j,k}^U\cdot v_j-r\cdot {\lambda }_{j,k}\\ \ge P_{j,k}^L\cdot v_j+r\cdot {\lambda }_{j,k}\end{array}$$

$$\displaystyle\begin{array}[]{l}\sum\limits_{\begin{subarray}{c}l:(l,j)\in T_{Y}\\ i:(i,l)\in T_{X}\end{subarray}}{C}_{i,k}\cdot q_{i,l}\cdot y_{l,j}\\ +\sum\limits_{i:(i,j)\in T_{Z}}{C}_{i,k}\cdot z_{i,j}\end{array}\begin{cases}\leq P^{U}_{j,k}\cdot v_{j}-r\left[\sum\limits_{\begin{subarray}{c}l:(l,j)\in T_{Y}\\ i:(i,l)\in T_{X}\end{subarray}}\hat{C}_{i,k}\cdot q_{i,l}\cdot y_{l,j}+\sum\limits_{i:(i,j)\in T_{Z}}\hat{C}_{i,k}\cdot z_{i,j}\right]\\ \geq P^{L}_{j,k}\cdot v_{j}+r\left[\sum\limits_{\begin{subarray}{c}l:(l,j)\in T_{Y}\\ i:(i,l)\in T_{X}\end{subarray}}\hat{C}_{i,k}\cdot q_{i,l}\cdot y_{l,j}+\sum\limits_{i:(i,j)\in T_{Z}}\hat{C}_{i,k}\cdot z_{i,j}\right]\end{cases}$$

$$\sum _{\begin{array}{c}l:(l,j)\in T_Y\\ i:(i,l)\in T_X\end{array}}\left(C_{i,k}+r\cdot {\hat{C}}_{i,k}\right)\cdot q_{i,l}\cdot y_{l,j}+\sum _{i:(i,j)\in T_Z}\left(C_{i,k}+r\cdot {\hat{C}}_{i,k}\right)\cdot z_{i,j}\le P_{j,k}^U\cdot v_j\forall j,k,$$

$$\sum _{\begin{array}{c}l:(l,j)\in T_Y\\ i:(i,l)\in T_X\end{array}}\left(C_{i,k}-r\cdot {\hat{C}}_{i,k}\right)\cdot q_{i,l}\cdot y_{l,j}+\sum _{i:(i,j)\in T_Z}\left(C_{i,k}-r\cdot {\hat{C}}_{i,k}\right)\cdot z_{i,j}\ge P_{j,k}^L\cdot v_j\forall j,k,$$

$$r^2\sum _i{\hat{C}}_{i,k}^2{x}_{i,j}^2\{\begin{array}{l}\le (P_{j,k}^U\cdot v_j-\sum _i{C}_{i,k}\cdot x_{i,j}{)}^2\\ \le (\sum _i{C}_{i,k}\cdot x_{i,j}-P_{j,k}^L\cdot v_j{)}^2,\end{array}$$

$$\sum _{\begin{array}{c}l:(l,j)\in T_Y\\ i:(i,l)\in T_X\end{array}}{C}_{i,k}\cdot q_{i,l}\cdot y_{l,j}+\sum _{i:(i,j)\in T_Z}{C}_{i,k}\cdot z_{i,j}$$

$${\lambda }_{j,k}$$

$$\underset{q_{i,l},y_{l,j},z_{i,j}}{\min }$$

$$\sum _i{\tilde{C}}_{i,k}\cdot x_{i,j}\{\begin{array}{l}\ge P_{j,k}^L\cdot v_j\\ \le P_{j,k}^U\cdot v_j\end{array}$$

$$\underset{j,k}{\max }\underset{{\tilde{C}}_{i,k}\in {\mathcal{U}}_p}{\max }\max \left\{\sum _i{\tilde{C}}_{i,k}\cdot x_{i,j}-P_{j,k}^U\cdot v_j,P_{j,k}^L\cdot v_j-\sum _i{\tilde{C}}_{i,k}\cdot x_{i,j}\right\},$$

Full text

1 Introduction

Robust optimization has become increasingly popular as a tool for planning and scheduling of chemical and manufacturing processes under parametric uncertainty 1, 2, 3, 4, 5, 6, 7, 8, 9. The key assumption in robust optimization is that the uncertain parameters are within some known, bounded uncertainty set. Constraints containing the uncertain parameters are robustified by requiring that they hold for all uncertainty set values. Benefits of the robust approach are that (i), for a rich class of problems and suitable uncertainty sets, considering uncertainty induces only moderate increases in the number of variables and constraints3, 10, 11 and (ii) it does not require probabilistic information 12 , even though such information may be used to construct uncertainty sets with probabilistic guarantees 3, 13, 14.

The robust optimization literature has traditionally studied polynomially-solveable problems. While starting with a focus on linear programming problems (LPs)15, 16, 17, robust optimization has since evolved to encompass many non-linear problem classes18, 19, 20. In particular, a rich class of problems which are convex in the decision variables and concave in the uncertain parameters are often tractable given a suitable uncertainty set21. In addition, the process systems engineering literature has used robust optimization for many applications in mixed-integer linear programming (MILP)2, 13.

This work considers the connection of robust optimization and global optimization in the context of the non-convex pooling problem. The pooling problem is relevant in oil and gas refining, water systems, supply chains, and more22, 23. It is non-linear and non-convex due to bilinear mixing terms and is strongly $\mathcal{NP}$-hard24. The pooling problem is frequently used as a case study in global optimization 25, 26, 27, 28, 29, 30, 31, 24, 32, 33. Parametric uncertainty in the pooling problem has been previously considered: Li et al.34 apply robust optimization to demand uncertainty in the linear product demand constraints while Barton, Li, and co-workers develop a two-stage stochastic programming approach for the extended pooling problem with uncertain inlet concentrations and apply their methodology more generally to several applications 35, 36, 37, 38, 39, 40. A similar approach has also been proposed in the context of water network synthesis41.

While multi-stage approaches can lead to less conservative solutions because they can incorporate information which becomes available during operation, they tend to assume perfect knowledge of the uncertain parameters after the uncertainty has been revealed. In practice, real-time measurements may be able to reduce parametric uncertainty but not entirely eliminate it. In the context of pooling, the uncertain inlet concentrations may also change over time. Rigorously addressing these temporal changes would require a multi-stage, e.g., approximate dynamic programming, approach to managing the pooling network (due to accumulation in the pools). The single-stage robust approach offers an alternative to this which is especially relevant in applications, e.g., oil/natural gas networks or water networks, where some parametric uncertainty remains unresolved until the final blending step. This conservatism is relevant (i) when quality constraint violations may lead to unusable products or (ii) in health and safety-critical domains.

The pooling problem is an interesting case study for global robust optimization, because it is non-linear, non-convex in the decision variables but linear (and therefore concave) in all potential uncertain parameters. Therefore many robust optimization techiques, such as duality-based reformulations and robust cutting plane approaches, are applicable. The resulting global optimization problem is not significantly more difficult than the nominal pooling problem without consideration of uncertainty. While any of the problem’s parameters could be considered uncertain, we focus on the interesting case of uncertain component inlet concentrations, which occur on the left hand side of the problems quality constraint in combination with bilinear terms in the decision variables. This work extends a previous conference paper42 which developed robust reformulations for the pooling problem:

• •

In addition to the duality-based robust reformulations from Wiebe et al. 42, we apply robust cutting planes to the non-convex pooling problem and develop a termination tolerance-based strategy for their integration with global optimization.

• •

We develop a simple but practical, safety factor inspired approach as a benchmark and as a computationally effective, conservative approximation to the robust problem.

• •

We assess and compare the suitability of these techniques in the context of global optimization and analyze their computational efficiency for several commonly used uncertainty set geometries (box, ellipsoidal, and polyhedral sets).

• •

We show how the optimal solution to the pooling problem is affected by varying degrees of uncertainty in the source component concentrations.

The first part of this paper reviews the different methods that have been developed for the robust optimization problem and discusses the situations in which each is applicable. It furthermore briefly discusses the integration of global and robust optimization in the context of robust cutting planes. The second part introduces the pooling problem, reviews its robust reformulations, and develops the safety factor-inspired approach. The final part compares the efficiency of the different methods for solving the robust pooling problem and highlights the effect of parametric uncertainty on the optimal flow configuration.

2 Robust optimization

Consider a generic robust optimization problem:

[TABLE]

where $\boldsymbol{x}\in\mathbbm{R}^{n}$ is the vector of decision variables, $\boldsymbol{\xi}\in\mathbbm{R}^{m}$ the vector of uncertain parameters, $f(\cdot)$ the objective function, $g(\cdot,\cdot)$ an uncertain constraint, $\mathcal{U}$ a non-empty, bounded set containing possible realizations of $\boldsymbol{\xi}$, and $\mathcal{X}$ is the certain feasible region (which does not depend on $\boldsymbol{\xi}$). Note that the objective does not depend on $\boldsymbol{\xi}$ and that we are only considering a single constraint. This formulation is without loss of generality, because an uncertain objective can always be transformed into an uncertain constraint using an epigraph formulation and because robust optimization usually considers each uncertain constraint separately (although there are settings where constraints can be modeled jointly18, 43). In robust optimization, it is often useful to rewrite Constraint (1b) as an inner maximization problem, leading to an equivalent bilevel formulation44:

[TABLE]

The complexity of Problem (1) and the applicable solution techniques largely depend on the convexity/concavity of $f(\boldsymbol{x})$, $g(\boldsymbol{x},\boldsymbol{\xi})$, $\mathcal{X}$, and $\mathcal{U}$ with respect to $\boldsymbol{x}$ and $\boldsymbol{\xi}$ 21. Four categories of problems exist:

$g(\boldsymbol{x},\boldsymbol{\xi})$ and $f(\boldsymbol{x})$ are convex in $\boldsymbol{x}$ and concave in $\boldsymbol{\xi}$ and $\mathcal{U}$ and $\mathcal{X}$ are convex. 2. 2.

$g(\boldsymbol{x},\boldsymbol{\xi})$ is concave in $\boldsymbol{\xi}$ and $\mathcal{U}$ is convex, but $g(\boldsymbol{x},\boldsymbol{\xi})$ or $f(\boldsymbol{x})$ is non-convex in $\boldsymbol{x}$ or $\mathcal{X}$ is non-convex. 3. 3.

$g(\boldsymbol{x},\boldsymbol{\xi})$ and $f(\boldsymbol{x})$ are convex in $\boldsymbol{x}$ and $\mathcal{X}$ is convex, but $g(\boldsymbol{x},\boldsymbol{\xi})$ is non-concave in $\boldsymbol{\xi}$ or $\mathcal{U}$ is non-convex. 4. 4.

$g(\boldsymbol{x},\boldsymbol{\xi})$ is non-concave in $\boldsymbol{\xi}$ or $\mathcal{U}$ is non-convex and $g(\boldsymbol{x},\boldsymbol{\xi})$ or $f(\boldsymbol{x})$ is non-convex in $\boldsymbol{x}$ or $\mathcal{X}$ is non-convex.

Fig. (1) summarizes which methods are applicable to each category. Most work in the robust optimization literature falls into Category 1, also called convex robust optimization. Category 1 includes linear programming problems with different types of convex uncertainty sets, e.g., box sets 15, ellipsoidal sets 16, and polyhedral sets 17. It also includes non-linear convex problems, e.g., convex-quadratic problems 18, 20, semi-definite programs 18, or general convex constraints 45. Both robust reformulation approaches and robust cutting planes have been applied extensively Category 1 and solutions can often be obtained in polynomial time.

Categories 3 and 4 are generally much harder, because the inner maximization Problem (2b) is a non-convex optimization problem. Furthermore, even finding a feasible solution to Problem (1) requires solving this non-convex problem to global optimality. Problems in these categories can generally only be solved using semi-infinite programming techniques 46, 47, 48, which are limited to small scale problems, or approximate robust optimization schemes 49, 50, 51, 52.

Category 2 has received limited attention in the robust optimization and semi-infinite programming communities53, 54, 55. A particular subclass, (non-convex) MILP problems, has been studied extensively in the process systems engineering literature, e.g., several robust planning and scheduling applications 2, 13, 56, 57. We focus instead on the subclass of continuous but non-linear, non-convex problems which are concave in the uncertain parameters, which includes the pooling problem. The same methods (both reformulation and cutting plane approaches) which are used to solve problems in Category 1 are generally also applicable to Category 2, but solutions can rarely be found in polynomial time. Solving such problems effectively therefore requires a tighter integration between robust optimization and global optimization techniques.

2.1 Reformulation approaches

The robust reformulation approach was first proposed by Soyster15 for LPs with a box uncertainty set and has since been developed for many classes of problems18, 16, 17, 45, 20. These approaches generally utilize strong duality by replacing the inner maximization in Problem (2) with its dual minimization problem:

[TABLE]

where $g^{*}(\cdot,\cdot)$ is the objective of the dual problem, $\boldsymbol{\lambda}$ are the dual variables, and $\Lambda$ is the feasible space. Problem (3) can clearly be written as a single level optimization problem which can often be solved using off-the-shelf solver software:

[TABLE]

2.2 Cutting plane approaches

Robust cutting plane approaches outer-approximating the feasible space by replacing the the uncertainty set $\mathcal{U}$ with a finite number of uncertainty realizations $\hat{\mathcal{U}}_{n}=\left\{\boldsymbol{\xi}_{1},\ldots,\boldsymbol{\xi}_{n}\right\}$. A robustly feasible solution is found by iteratively adding uncertainty realizations which violate the constraint. The key idea is to alternately solve a master problem:

[TABLE]

and a separation problem:

[TABLE]

which determines whether the current solution $\boldsymbol{x}_{n}$ is robustly feasible ($\max\limits_{\boldsymbol{\xi}\in U}g(\boldsymbol{x}_{n},\boldsymbol{\xi})\leq 0$) and, if it is not, generates a new cut $\boldsymbol{\xi}_{n+1}$.

Note that the master and separation problem are both convex optimization problems when $f(\boldsymbol{x})$ and $g(\boldsymbol{x},\boldsymbol{\xi})$ are convex in $\boldsymbol{x}$ and concave in $\boldsymbol{\xi}$ and $\mathcal{U}$ and $\mathcal{X}$ are convex. When $f(\boldsymbol{x})$ or $g(\boldsymbol{x},\boldsymbol{\xi})$ is non-convex in $\boldsymbol{x}$, e.g., in the pooling problem, the master problem becomes a global optimization problem which has to be solved at every iteration. This could be disadvantageous because the global optimization solver may spend a lot of time proving optimality of a solution which is not robustly feasible. One way to avoid this would be integrating the robust cutting planes generation into the global optimization procedure. The separation problem could be solved to generate a new cut whenever an improved incumbent solution is found by the global optimization solver, e.g., using lazy constraints. These cuts have been proposed for MILP 58, but could also be useful in the global optimization setting.

Unlike mature MILP solver software, however, global optimization solvers do not routinely offer callbacks for adding lazy or incumbent constraints. In light of this, we propose a simpler strategy for integration of global optimization and robust cutting planes by adaptively changing the termination tolerance for the master problem. The general procedure is outlined in Algorithm (1).

The master problem is first solved to a large tolerance $\delta_{0}$ Whenever a robustly feasible solution is found, the tolerance $\delta$ is reduced by a factor $\lambda\leq 1$ until it reaches the final tolerance $\delta^{*}$.

Robust cutting planes are similar to the concept of backoff in the control literature59, 60, 61, 62. The backoff approach generally also starts by solving the nominal problem and then subsequently “backing off” to suboptimal but more robust solutions which remain feasible in light of uncertainty. While the robust optimization approach assumes that the uncertain parameter is within a known uncertainty set, uncertainty in the backoff approach is usually due to process dynamics.

3 Robust pooling

Consider the $q$-formulation of the standard pooling problem 63, 64:

[TABLE]

The decision variables are the fraction of flow $q_{i,l}$ from source $i$ to pool $l$, the flow $y_{l,j}$ from pool $l$ to terminal $j$, and the direct flow $z_{i,j}$ from source $i$ to terminal $j$. The sets $T_{X},T_{Y},T_{Z}$ describe the connections between sources, pools and terminals. $A^{L}_{i}/A^{U}_{i}$, $D^{L}_{j}/D^{U}_{j}$ and $P^{L}_{j,k}/P^{U}_{j,k}$ are the lower and upper limits for feed availability, product demand, and product quality respectively. The pool capacity is limited by $S_{l}$, the inlet concentration of component $k$ at source $i$ is ${C}_{i,k}$, product $j$ is sold at price $d_{j}$ and the cost of material from source $i$ is $c_{i}$. The objective minimizes negative profit. For notational brevity we introduce the total flow $v_{j}$ at terminal $j$ in Constraint (7w):

[TABLE]

While non-linearities (bilinear mixing terms) in the decision variables occur in the objective function, feed availability, and product quality constraints, the problem is linear in all (potentially uncertain) parameters. The problem therefore belongs to Category 2 outlined above, independent of which parameter is considered to be uncertain. Furthermore, unlike the original $p$-formulation of the pooling problem, there are no potentially uncertain equality constraints. In fact, the only equality constraint is the simplex which does not contain parameters. This makes the $q$-formulation particularly suitable for robust optimization, because uncertain equality constraints generally either have to be eliminated or treated using adjustable robust optimization 65, 66.

While we could consider any or all of the problems parameters to be uncertain, we focus on robustifying the product quality Constraints (7w) considering uncertain input component concentrations $\tilde{C}_{i,k}$:

[TABLE]

This constraint is particularly interesting, because the uncertain parameters occur on the left hand side and in combination with the bilinear terms $q_{i,l}\cdot y_{l,j}$.

3.1 Uncertainty set

We consider a commonly used class of uncertainty sets $\mathcal{U}_{p}$ defined by the $p$-norm13:

[TABLE]

where $\boldsymbol{\xi}_{k}=[\ldots,\xi_{i,k},\ldots]^{T}$, ${C}_{i,k}$ is the nominal value of $\tilde{C}_{i,k}$, $\hat{C}_{i,k}$ its maximum deviation, and $r$ a size parameter which determines how much of the uncertainty is considered (For $r=0$ the constraint is equivalent to the nominal constraint and for $r=1$ all possible realizations of $\tilde{C}_{i,k}$ are in the uncertainty set. The set $\mathcal{U}_{p}$ corresponds to a box uncertainty set for $p=\infty$, an ellipsoid for $p=2$, and a polyhedron for $p=1$. All of these uncertainty set geometries have been discussed extensively in the robust optimization literature 13.

While these uncertainty sets can capture different types of uncertainty scenarios, they do not utilize any known correlation between the random variables, e.g., the geographic structure of the pooling network. For example, inlet concentrations at geographically neighboring sources may be correlated more strongly than geographically distant sources. To leverage such dependence structure, we constructe ellipsoidal uncertainty sets based on distance between sources using kernel functions. Kernel functions, e.g., the squared exponential kernel:

[TABLE]

where $\sigma^{2}>0$ is the signal variance and $l>0$ is the length scale, can be used to model correlation between two geographic locations $\boldsymbol{\ell}$ and $\boldsymbol{\ell}^{\prime}$67. We construct a covariance matrix $\Sigma$ with elements $\sigma_{i,i^{\prime}}$ describing the correlation between each pair of network sources $i$ and $i^{\prime}$:

[TABLE]

where $\ell_{i}$ is the location of source $i$. Note that this approach could be combined with Gaussian process regression (Kriging) to construct $\Sigma$ based on available concentration data from different locations67. Using $\Sigma$, Eqn. (11) constructs an ellipsoidal uncertainty set exploiting the correlation between sources:

[TABLE]

We note that, if $\sigma=1$ in Eqn. (10) and the sources are very far from each other, i.e., there is no correlation, $\Sigma$ becomes the identity matrix and the uncertainty set in Eqn. (11) is equivalent to $\mathcal{U}_{2}$ (Eqn. 9).

3.2 Reformulation

Because Constraint (8) is linear in the uncertain parameters $\tilde{C}_{i,k}$ and $\mathcal{U}_{p}$ is convex, duality-based robust reformulation techniques are applicable. To this end, we introduce the indicator function:

[TABLE]

allowing us to rewrite Constraint (8) in standard linear form:

[TABLE]

where $x_{i,j}$ is the total flow from source $i$ to terminal $j$:

[TABLE]

Using robust optimization results 13, the semi-infinite Constraint (12) can be rewritten as an equivalent deterministic constraint:

[TABLE]

where $\lambda_{j,k}$ depends on the selected uncertainty set geometry $p$ as indicated in Table (1).

Note that no absolute values are required because the flow $x_{i,j}$ is always positive.

For the box set substituting Eqn. 13 for $x_{i,j}$ leads to:

[TABLE]

which simplifies to:

[TABLE]

and

[TABLE]

respectively. In the simple case of the box uncertainty set, the worst case uncertainty scenario is clearly always ${C}_{i,k}+r\cdot\hat{C}_{i,k}$ for the upper quality constraint and ${C}_{i,k}-r\cdot\hat{C}_{i,k}$ for the lower quality constraint. The reformulation of the ellipsoidal uncertainty set can potentially be improved by squaring each side of the quality constraints:

[TABLE]

This constraint eliminates the square root which can be numerically challenging, but is only valid when the nominal Constraint (7w) is also added to the model, increasing the number of equations. The same applies to the ellipsoidal uncertainty set with correlation. For the polyhedral set, Constraint (8) can be written as:

[TABLE]

While the box uncertainty set does not change the problem complexity with respect to the nominal case, the ellipsoidal and polyhedral uncertainty sets both add (potentially many) convex constraints to the model. These constraints do not change the complexity class, but they may introduce practical difficulties. While the ellipsoidal set adds only $\mathcal{O}(|i|\cdot|j|)$ constraints, the polyhedral set adds $\mathcal{O}(|i|\cdot|j|\cdot|k|)$ constraints and variables. An advantage of the polyhedral set, however, is that the problem remains bilinear, while the ellipsoidal set introduces higher order terms.

3.3 Cutting planes

The cutting plane approach outlined above is directly applicable to the pooling problem. The master problem is:

[TABLE]

where $\hat{\mathcal{U}}_{0}=\left\{{C}_{i,k}\right\}$ contains only the nominal values. The separation problem is:

[TABLE]

which can be solved by solving $2\cdot|J|\cdot|K|$ convex optimization problems. For the uncertainty sets considered in this work, it can also be reformulated as a MILP (for box and polyhedral sets) or a mixed-integer quadratically constrained program (MIQCP) using a big-M reformulation which can be solved effectively with existing MIP solvers.

We explore two different ways of adding cuts: adding a single cut in every round, consisting of the $j^{*}$-$k^{*}$-quality constraint which is most violated, and adding multiple cuts in every round, consisting of all $2\cdot|J|\cdot|K|$ quality constraints for the worst case scenario $\tilde{C}_{i,k}^{*}$, where ($j^{*},k^{*},\tilde{C}_{i,k}^{*}=\operatorname*{arg\,max}$ Eqn. (20)). Note that an intermediate approach which adds the first $n$ most violated cuts would also be possible. Furthermore, $\hat{\mathcal{U}}_{0}$ could also be initialized with other values. For example, for the box uncertainty set the maximum and minimum concentrations could be used. This would make the master problem essentially equivalent to the reformulation (Eqn. 17) and the cutting plane algorithm would always converge after the first iteration.

While there are examples for which Algorithm 1 does not converge 68, it often performs similar to the reformulation approach in practice 69, 58. Furthermore, the algorithm is guaranteed to converge for the robust pooling problem with the Table 1 uncertainty sets. Since Eqn. 2b is linear in $\tilde{C}_{i,k}$, the separation problem generally has a finite number of solutions for the box and polyhedral sets, the vertices of the polytope. For the ellipsoidal set, Mínguez and Casero–Alonso70 show that convergence is guaranteed when the constraint obtained with the reformulation approach is convex by showing that the cutting plane approach is equivalent to solving the reformulation using outer-approximation. Since all of the non-convexities are handled by the global optimization solver and the robust reformulation (Eq. 14) is convex in $x_{i,j}$, the Mínguez and Casero–Alonso70 result also applies to the pooling problem.

3.4 Safety factors approach

A very simple but practically relevant approach to design under uncertainty is the concept of safety factors 71, 72, 73. In this approach, the space of feasible designs is restricted by scaling bounds with a fixed safety factor $s>0$. In the case of the pooling problem with uncertain inlet concentration, the upper (and lower) bound $P^{U}_{j,k}$ (and $P^{L}_{j,k}$) may be scaled to discard solutions which may easily become infeasible, while $\tilde{C}_{i,k}$ is replaced by its nominal value ${C}_{i,k}$:

[TABLE]

A benefit of safety factors is that they do not increase the complexity of the problem at all. A disadvantage is, however, that one generally has to rely on experience when choosing a suitable value for $s$. Note that Problem (21) is equivalent to a robust version of Problem (7) with uncertain bounds $P^{U}_{j,k}$ and $P^{L}_{j,k}$ and a box uncertainty set, e.g.:

[TABLE]

For $s=1$, Problem (21) is equivalent to the nominal problem and for $s\to\infty$ the optimal solution is $y_{l,j}=z_{i,j}=0$ for all $(l,j)$, $(i,j)$, i.e., the safest solution is producing nothing.

To compare the safety factor approach with the robust optimization approach, we find the smallest safety factor $s^{\text{min}}$ which gives a feasible solution for the robust problem with a given uncertainty set $\mathcal{U}_{p}$. This is outlined in Algorithm (2). We employ a combination of bisection and secant search to find the safety factor for which Eqn. (20) is just feasible.

4 Results

The robust reformulation, robust cutting plane, and safety-factor approach outlined above were applied to 14 literature pooling instances 74, 75, 64, 63, 76. The model was implemented in GAMS 26.1.0 and solved on an i7-6700 CPU with $8\times 3.4$GHz and 16GB RAM. NLP problems were solved using ANTIGONE 1.1 77 while MILP problems were solved using CPLEX 12.8. Each instance was solved for 30 values of $r$ for the box, ellipsoidal, and polyhedral uncertainty set. A time limit of 1 hr, a relative termination tolerance of $10^{-6}$, and $\hat{C}_{i,k}={C}_{i,k}$ were used. For the cutting plane approach, a maximum of 200 cuts were added.

Fig. 3 shows the optimal solution to instance Adhya 1 for three uncertainty set sizes $r=0.0,0.1,0.3$ using the polyhedral uncertainty set. The width of the arrows between nodes indicates the fraction of flow from/to each pool. In the nominal case, i.e., no uncertainty is considered ($r=0$), Products 2 and 4 are produced using Sources 1, 2, 4, and 5. As $r$ is increased, hedging against more uncertainty in the inlet concentrations, at first the total amount of Products 2 and 4 produced remains constant, but Source 1 is gradually replaced by Source 2 which is more expensive but of better quality. When $r$ reaches $0.1$, Source 1 stops being useful altogether. Past $r=0.14$ the production of Product 2 ends, because its quality cannot be guaranteed anymore for all parameter values in the uncertainty. Stopping Product 2 production causes a sharp decrease in the expected profit from $446.2$ to $65.9$. Fig. (3) shows the objective value of instance Adhya 1 (scaled by the nominal objective value) as a function of the uncertainty set size $r$. The three scenarios shown in Fig. (3) and discussed above are indicated by vertical dashed lines. Results are shown for all three uncertainty set geometries. For a given uncertainty set geometry, as $r$ increases, the fraction of the nominal objective value achieved always decreases, as would be expected. Notice that, independent of the geometry, the profit first always decreases slowly as cheap sources are substituted by higher quality sources and then suddenly drops as the production of Product 2 becomes infeasible.

Fig. (4) shows similar trends for a number of other instances.

Figs. (3) and (4) show a clear ordering between uncertainty sets geometries. By construction the polyhedral uncertainty set is always smallest and the box uncertainty set largest for a given $r$. Note that this does not necessarily mean that the polyhedral set is always superior. While it is always less conservative, i.e., achieves a better worst case objective value, it may also be less robust. Black points in Fig. (4) indicate instances which could not be solved to optimality within 1 hour.

Fig. (5) shows a comparison for instance Foulds 3 between the reformulation and the multi-cut version of the robust cutting plane approach. As would be expected, the two approaches lead to identical objective values, except for some small deviations with the ellipsoidal sets in regions where neither approach converges to the optimal solution. The cutting plane approach, however, converges for many more instances than the reformulation approach. For this particular instance with the polyhedral uncertainty set, it converges for all but one value of $r$ while the reformulation approach only converges for very large values of $r$. A similar trend, albeit not as strong, can be seen for the ellipsoidal uncertainty set.

Table 2 shows this trend generally holds across instances. The multi-cut cutting plane strategy solves a larger percentage of instances than the reformulation approach for all uncertainty set geometries. Note that the single-cut cutting plane approach, which adds only one cut in every round, performs poorly in comparison to both alternatives. The uncertainty set geometry also clearly has an effect on the tractability of the robust problem. The box uncertainty set-constrained problem is generally easy to solve. This is expected, especially for the reformulation approach, because the problem difficulty is equivalent to the nominal case. The reformulation approach performs particularly poorly for the ellipsoidal set. This is most likely due to the higher order non-linear terms in this formulation. Table 2 also shows the median time taken across instances for those instances which could be solved to optimality.

While the cutting plane approach could potentially add a large number of constraints, in practice it often adds a small number of constraints. Across instances, the multi-cut strategy needs an average of 13.3 iterations to converge. This leads to an average final problem size of 335 constraints which is only slightly larger than the average number of constraints in the reformulation approach: 328 for the polyhedral, 166 for the ellipsoidal, and 154 for the box uncertainty set. The cutting plane approach does tends to take longer to solve than the reformulation approach. This may, however, be at least partially because this approach solves a larger fraction of the more difficult instances to optimality than the reformulation approach.

While these results suggest that, at least for the pooling problem with the uncertainty sets considered in this work, the cutting plane approach is the preferable strategy for global robust optimization, the reformulation approach has one important advantage: any feasible solution to the robust reformulation is a valid feasible (albeit potentially non-optimal) solution to the robust problem. In contrast, intermediate solutions in the cutting plane algorithm are often not robustly feasible.

Fig. (6) shows a comparison between the multi-cut robust cutting plane approach and the safety factor approach outlined above. For instance Haverly 2, the safety factor approach (with optimized safety factors) is almost identical to the robust approach. This means that, for this instance, even the ellipsoidal and polyhedral uncertainty set on the inlet concentration $\tilde{C}_{i,k}$ can be approximated by a box uncertainty set on $P^{U}_{j,k}$ and $P^{L}_{j,k}$. However, e.g., for instances Adhya 4 and Bental 5, the safety-factor approach is not equivalent anymore and leads to a conservative solution for both the ellipsoidal and polyhedral set. For instance Foulds 3, the safety factor approach is equivalent to the box uncertainty set. Overall, while the safety factor is very effective in the sense that it does not increase the complexity of the problem, it is generally a conservative approximation to the robust problem.

Fig. (7) demonstrates correlation between sources. The figure considers three different geographic layouts for instance Adhya 1 with varying degrees of proximity between sources. The for the “far” scenario are very similar to those of the ellipsoidal uncertainty set without correlation. With decreasing distance between sources, i.e., increasing correlation, the worst case profit generally decreases because correlation between sources will make it more likely that, e.g., inlet concentrations at multiple sources are all large simultaneously.

5 Conclusion

The pooling problem is an interesting case study for robust optimization because it is non-convex in the decision variables but concave in all potentially uncertain parameters. Applying robust reformulation or cutting plane approaches to this problem leads to a global optimization problem with increased complexity compared to the nominal case. The tractability of the resulting optimization problem with global optimization solvers is highly dependent on the selected uncertainty set. While simple box uncertainty sets hardly increase computational time, more advanced ellipsoidal and polyhedral sets are more difficult to solve. For the pooling problem, the cutting plane approach generally solves a larger percentage of instances to optimality while the reformulation approach has the advantage that it produces guaranteed feasible intermediate solutions. The relevance of optimization under uncertainty is supported by the significant changes in the solution observed when different degrees of uncertainty are taken into account.

6 Acknowledgements

This invited contribution is part of the I&EC Research special issue for the 2019 Class of Influential Researchers. This work was funded by the Engineering & Physical Sciences Research Council (EPSRC) Center for Doctoral Training in High Performance Embedded and Distributed Systems (EP/L016796/1), an EPSRC/Schlumberger CASE studentship to J.W. (EP/R511961/1, voucher 17000145), and an EPSRC Research Fellowship to R.M. (EP/P016871/1).

7 Nomenclature