Global and Robust Optimisation for Non-Convex Quadratic Programs

Abstract

This paper presents a novel algorithm integrating global and robust optimisation methods to solve continuous non-convex quadratic problems under convex uncertainty sets. The proposed Robust spatial branch-and-bound (RsBB) algorithm combines the principles of spatial branch-and-bound (sBB) with robust cutting planes. We apply the RsBB algorithm to quadratically constrained quadratic programming (QCQP) pooling problems, utilising McCormick envelopes to obtain convex lower bounds. The performance of the RsBB algorithm is compared with state-of-the-art methods that rely on global solvers. As computational test bed for our proposed approach we focus on pooling problems under different types and sizes of uncertainty sets. The findings of our work highlight the efficiency of the RsBB algorithm in terms of computational time and optimality convergence and provide insights to the advantages of combining robustness and optimality search.

1 Introduction

One of the fundamental assumptions when constructing mathematical programming models is the degree of uncertainty in the input parameter data. If uncertainty is neglected, then a deterministic optimisation problem is generated, while for the case where uncertainty is considered to some degree, there are various methods to model uncertainty. There are two key drivers for the selection of the appropriate optimisation method: statistical data availability and desired degree of robustness. Robust optimisation has emerged as a leading method for problems with limited uncertain parameter data and requirements for a risk-averse solution. In a robust optimisation context, solutions that lead to attainable decisions under all parameter realisations are referred to as robust feasible, while the solution that also leads to the best objective value is called robust optimal. A general form of the robust optimisation problem is given by ($P_{rob}$).

$$\underset{x\in\mathcal{X}}{\min}~{}f(x) P_{rob}$$ ()

$$\text{s.t.}~{}g(x,u)\leq 0,~{}\forall u\in\mathcal{U}$$

where $x\in\mathcal{X}\subseteq\mathbb{R}^{n_{x}}$ denotes continuous decisions, and $u$ is the vector of uncertain parameters that resides within an uncertainty set $\mathcal{U}$. Despite the presence of an embedded optimisation problem within ($P_{rob}$), robust optimisation problems can yield tractable reformulations under restrictive assumptions about the problem structure. That is, if Problem ($P_{rob}$) is convex with respect to the decisions $x$ and concave in $u$, and $\mathcal{U}\subset$ is a convex set, then using duality theory one can derive a tractable robust counterpart reformulation of Problem ($P_{rob}$). Because of this, the field of convex robust optimisation is well-established, whereas the field of non-convex robust optimisation is still at its infancy [1].

State-of-the-art robust optimisation algorithms for non-convex problems typically rely on global solvers to obtain robust optimal solutions. The prevailing approach involves deriving the dual reformulation and utilising a global solver, either directly or adaptively. However, the dual reformulation can increase the problem complexity, and in the case of challenging non-convex problems may only yield a robust feasible solution. An alternative approach is based on an iterative robust cutting plane algorithm. In this case, significant computational time can be spent searching for a global solution that may be deemed robust infeasible.

Contributions of this work: Our research hypothesis is that integrated exploration of global and robust optimality can yield computational benefits. To this end, in this article we introduce a novel algorithm that conducts a concurrent global optimality and robustness search for continuous non-convex problems. The proposed approach integrates spatial branch-and-bound with robust cutting plane notions. The key idea is that, while exploring the branch-and-bound tree, we assess the robustness of nodes entailing the best-found solutions. We illustrate the performance and benefits of our proposed approach through benchmark Quadratically Constrained Quadratic Programs (QCQPs) of pooling problems. At each node, the non-convex problem is solved via a local solver. If the computed solution is as good as the best-found so far, an infeasibility test is performed to evaluate the robustness of the obtained solution. If not, then the corresponding cutting planes are added both to the non-convex and convex problems. The algorithm proceeds to the next step once no more violations are detected. The convex problem is solved next, and the solution of this problem is used to decide the most promising variable for branching into two child nodes. With the use of appropriate fathoming criteria, the tree nodes are exhausted, and the robust optimal solution is obtained.

Section 2 provides a literature review on state-of-the-art nonlinear robust optimisation methods and their applications, along with a brief overview of the pooling problem which serves as our case study. The outline of the robust cutting planes algorithm and the details for proposed Robust spatial branch-and-bound algorithm are introduced in Section 3. In Section 4 we present the pooling problem under uncertainty, and derive the corresponding formulations for the dual reformulation and RsBB solution methods. In Section 5 we evaluate the performance of the RsBB algorithm compared to the state-of-the-art for benchmark pooling problems under box, ellipsoidal and polyhedral sets of varying sizes.

2 Literature review

Robust optimisation (RO) is developing into one of the prevailing risk-averse methods to study problems under uncertainty in the optimisation literature [2, 3]. It is largely preferred for applications where there is no or limited statistical data for the uncertain parameters, or when dealing with hard feasibility constraints. In RO, the modeller must decide an appropriate, in size and shape, uncertainty set to characterise the uncertain parameters. The RO algorithms provide an optimal solution that corresponds to the worst-case uncertainty realisation and at the same time guarantee that the solution is feasible given any parameter value from the selected uncertainty set. As such, robust optimisation problems belong to the general class of Semi-Infinite Programming (SIP)[4]. The main challenge of RO is to find the appropriate way to transform the uncertain optimisation problem into a tractable deterministic formulation. The most prevailing methods in that direction, are dual reformulation [5, 6, 7] and robust cutting sets [8, 9]. The first comparison of the two methods can be traced to [10] who applied both methods in linear and mixed integer problems. Based on the computational experiments, cutting planes were found to be more suitable for the linear programs, while the dual reformulation performed better for the mixed integer programs. [11] quantified the performance of the two methods for different variants of the cutting planes, as well as that of a hybrid approach. In contrast to previous results in the literature, for the ellipsoidal uncertainty sets, the dual reformulation dominated in the linear programs and cutting planes for mixed-integer instances. For polyhedral uncertainty sets neither method was deemed to be superior.

Nonlinear robust optimisation problems have gained increasing attention from the research community over the past two decades [12]. For nonlinear convex problems under convex uncertainty sets both of the aforementioned methods can be applied. [13] highlight the computational challenge in solving nonlinear robust optimisation problems under general norms and propose an approximate reformulation that employs a linearisation of the uncertainty set. For nonlinear convex problems convex analysis may also be applied to retrieve tractable reformulations for specific classes of uncertainty sets [14, 6]. [15] presented a rigorous method to choose safety margins for uncertain parameters in nonlinear and non-convex problems via local linearisation around a nominal uncertainty value. Instead of using the nominal uncertainty point, [16] proposed an algorithm based on outer approximation over sampled uncertainty points and tested their methodology in non-convex process systems engineering problems. [17] performed a neighbourhood search using a robust local descent method to avoid locally infeasible regions. Their methodology obtained local robust solutions for non-convex and simulation-based problems. [18] used a sequential bilevel approach for non-convex min-max problems resulting in less conservative solutions compared to the state-of-the-art.

In various applications, finding a local robust solution may not suffice, necessitating thus the need for coupling global and robust optimisation. Polynomial max-min and min-max problems are represented as general SDP problems and their corresponding robust optimal values can be approximated solving a hierarchy of SDP relaxations [19, 20]. [21] studied the scheduling of crude oil operations under demand uncertainty. The authors derived the deterministic robust counterpart of the uncertain constraints using a tailored branch and bound (BB) algorithm. [22] addressed the pooling problem under uncertainty comparing dual reformulation and cutting plane methods using a global solver. For the cutting planes method they relied on using a global solver and gradually decreasing the solver tolerance. Cutting planes outperformed the dual reformulation approach, since the latter increased the problem complexity for ellipsoidal and polyhedral uncertainty sets. [9] extending the work of [8], and proposed a robust cutting set algorithm to be able to certify robust solutions to non-convex problems entailing a large number of equality constraints. The quality of the final solution, i.e. being robust feasible or robust optimal, is solver dependant. On the same year, [23] employed enhanced normalised multiparametric disaggregation technique (ENMDT) and optimality-based bound tightening to solve the problem of refinery-wide planning operations under uncertainties. [24] introduced a BB algorithm with interval arithmetic to solve benchmark unconstrained non-convex optimisation problems of low-dimensionality under robust uncertainty. Most robust optimisation methods require the convexity of the uncertainty set or the inner-maximisation problem, however [25] proposed an adaptive bundle method for MINLP problems under robust uncertainty that allows the inner-maximisation problem to be non-convex. In the SIP literature, [26] proposed an iterative global optimisation method that relies on continuity rather than convexity for both the outer and inner optimisation problems. Their approach involved constructing a decreasing sequence of right-hand side constraint relaxations for the outer- minimisation problem. The feasibility of the obtained solution was evaluated in an inner-maximisation problem. Both optimisation steps were conducted using a global solver. This method was later generalised to address general SIP problems [27].

We derive our QCQP case study from the process systems engineering literature as the planning of pooling problems. The pooling problem belongs to the general class of non-convex quadratic problems and is proven to be NP-hard even for a small number of variables [28]. Historically, the most common solution approach to address non-convex NLPs, such as the pooling problem, is spatial branch-and-bound (sBB) [29]. In sBB, in contrast to classic branch-and-bound, branching is performed on continuous variables. The original non-convex problem constitutes the upper bound and the equivalent relaxed linear programming (LP) problem the lower bound. [30] used convex relaxations to approximate the bilinear terms within a branch-and-bound framework. [31] integrated the McCormick envelopes [32] with RLT constraints [33] in an sBB framework selecting the branching variable as the one with the worst-approximation error. Using RLT constraints result in tightening the convex relaxations of the bilinear terms without affecting the feasible region. Later on, [28] used both convex and concave envelopes following a novel branching strategy. For high-dimensional problems, the number of convex relaxation constraints increases drastically increasing the complexity of the problem. To that end, [34] proposed a branch-and-cut algorithm where only the violating linearisation constraints are added as cuts, instead of all convex and concave envelopes. Every generated cut is applied to all tree nodes. [35], replaced the non-convex pooling problem with a piecewise linear relaxation resulting in an MILP problem. [36] studied the impact of using piecewise linearisation for the linear envelopes, they also introduced pooling benchmark problems of increased dimensionality and complexity. [37] introduced the open-source MIQCQP solver GALINI that is based on a branch-and-cut algorithm using convex linear relaxations and mixed-integer restrictions for the pooling problems. The performance of the MILP pooling formulation can be improved using appropriate tightening methods and reduce computational time [38].

3 Methods

The proposed Robust spatial branch-and-bound (RsBB) integrates the robust cutting set algorithm within a spatial branch-and-bound global optimisation algorithm in search of the optimal robust solution. The robust cutting set algorithm [8] will be outlined with a theoretical study on the convergence of the integrated robust cutting planes and BB algorithms. In the following subsection, the proposed RsBB algorithm is introduced. The assumptions for our methodology are the following.

1. The domain of x, $\mathcal{X}\subset\mathbb{R}^{n_{x}}$ is bounded;

2. The objective and constraint functions $f,~{}g:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}$ are uniformly Lipschitz continuous on $x$;

3. The uncertainty set $\mathcal{U}$ is convex and generates a non-empty robust feasible region;

4. The uncertain constraints $g$ are concave in $u$;

5. The problem of maximising $g$ over $u\in\mathcal{U}$ is tractable and the exact solution $u^{*}$ can be obtained.

3.1 Robust cutting set algorithm

Consider the deterministic problem ($P_{nom}$) for the nominal uncertainty value $u_{nom}$.

$$\underset{x\in\mathcal{X}}{\min}~{}f(x) P_{nom}$$ ()

$$\text{s.t. }~{}g(x,u_{nom})\leq 0$$

Following [8], there exists a subset $\hat{\mathcal{U}}=\{u_{0},\cdots,u_{K}\}\subseteq\mathcal{U}$ involving $K$ finite uncertainty samples or realizations and the corresponding problem $P_{samp}$.

$$\underset{x\in\mathcal{X}}{\min}~{}f(x) P_{samp}$$ ()

$$\text{s.t.}~{}g(x,u)\leq 0,~{}u\in\hat{\mathcal{U}}$$

Let $u_{nom}=u_{0}$, the following uncertainty sample can be determined via solving the outer-minimisation problem $P_{nom}$, and evaluating the robustness of the solution $x^{*}_{nom}$ in the inner-maximisation problem T.

$$\underset{u}{\max}~{}g(x^{*},u)$$ (T)

$$\text{s.t. }\,u\in\mathcal{U}$$

Given specific properties of $g$ and $\mathcal{U}$ the exactness of the solution of problem (T) can be guaranteed [8]. We focus on the case that $g$ is affine in $u$ and we are evaluating box, polyhedral and ellipsoidal sets, hence exactness of the solution is satisfied. If a solution $u^{*}$ is detected, such that $g(x^{*},u^{*})>0$, then it is added on subset $\hat{\mathcal{U}}$. Let $\hat{G}(x)$ be the worst-case constraint function

$$\hat{G}(x)=\max\{g(x,u_{0}),\cdots,g(x,u_{K})\}$$

it is true that

$$g(x,u_{0})\leq\hat{G}(x)\leq G(x)$$

where $G(x)=\underset{u\in\mathcal{U}}{\sup}\,g(x,u)$.

Proposition 1.

The optimal solution of the sampled problem provides a lower bound for the robust problem.

Proof.

Let $\mathcal{F}_{rob},~{}\mathcal{F}_{samp},~{}\mathcal{F}_{nom}$ the feasible sets and $z^{*}_{rob},~{}z^{*}_{samp},~{}z^{*}_{nom}$ the optimal values for the robust, sampled and nominal problems. Then, $\mathcal{F}_{nom}:=\{x\in\mathcal{X},u\in\hat{\mathcal{U}}_{0}:g(x,u)\leq 0\}$, $\mathcal{F}_{samp}:=\{x\in\mathcal{X},u\in\hat{\mathcal{U}}:g(x,u)\leq 0\}$ and $\mathcal{F}_{nom}:=\{x\in\mathcal{X},u\in\mathcal{U}:g(x,u)\leq 0\}$. By definition $\hat{\mathcal{U}}_{0}\subseteq\hat{\mathcal{U}}\subseteq\mathcal{U}$. Hence, $\mathcal{F}_{rob}\subseteq\mathcal{F}_{samp}\subseteq\mathcal{F}_{nom}$ and $z^{*}_{rob}\geq z^{*}_{samp}\geq z^{*}_{nom}$. ∎

Proposition 2.

If the optimal solution of the sampled problem is feasible for the robust problem, then it is also optimal for the robust problem.

Proof.

Let $x^{*}_{samp}\in\mathcal{F}_{samp}$ be the optimal solution for the sampled problem. From Proposition 1 $z^{*}_{rob}\geq z^{*}_{samp}$. If $G(x^{*}_{samp})>0$ then $x^{*}_{samp}\in\mathcal{F}_{rob}$ as the sampled solution belongs to the robust feasible region. Hence $z^{*}_{rob}=z^{*}_{samp}$ at $x^{*}_{samp}$. ∎

[8] have proved that there exist a finite number $K$ of cutting plane iterations to find the robust optimal solution using the robust cutting set method, without depending on convexity assumptions of $f,~{}g$ on $x$. An outline of the robust cutting set algorithm is provided in Algorithm 1. We note that [9] have implemented a generalised cutting set method in the PyROS solver enabling the solution of non-convex problems as well.

Algorithm 1 Robust Cutting Set Algorithm

1:Solve ($P_{nom}$) $\rightarrow x^{*}$

2:Solve (T) $\rightarrow u^{*}$

3:if $g(u^{*},x^{*})>0$ then

4:     Append the robust set $\hat{\mathcal{U}}(u)=\hat{\mathcal{U}}\cap\{u^{*}\}$

5:else

6:     Go to Step 9

7:end if

8:Solve ($P_{samp}$) $\rightarrow x^{*}$

9:Go to Step 2

10:Terminate: Robust solution obtained

3.2 Robust spatial Branch and Bound

The aim of this work, is to solve problems with non-convex bilinear functions $f,~{}g$ in $x$ under robust uncertainty using an sBB algorithm via McCormick relaxations. Let $x_{i},~{}i=1,\cdots,B$ be the index for the complicating variables participating in bilinear formulations where $B\leq n_{x}$. Consider the convex nominal problem $\tilde{P}_{nom}$ and sampled problem $\tilde{P}_{samp}$ augmented by $y$ for $n_{b}$ bilinear terms of $x$ using the McCormick envelopes [32].

$$\min~{}f(x,y) \tilde{P}_{nom}$$ ()

$$\text{s.t.}~{}g(x,y,u_{nom})\leq 0$$

$$\qquad\bm{\mu}(x,y)\leq 0$$

$$x\in\mathcal{X},~{}y\in\mathbb{R}^{n_{b}}$$

$$\min~{}f(x,y) \tilde{P}_{samp}$$ ()

$$\text{s.t.}~{}g(x,y,u)\leq 0,~{}u\in\hat{\mathcal{U}}$$

$$\qquad\bm{\mu}(x,y)\leq 0$$

$$x\in\mathcal{X},~{}y\in\mathbb{R}^{n_{b}}$$

Proposition 3.

The relaxed problem provides a lower bound for the the original robust problem.

Proof.

Let $\tilde{\mathcal{F}}_{rob},~{}\tilde{\mathcal{F}}_{samp},~{}\tilde{\mathcal{F}}_{nom}$ be the robust, sampled and nominal feasible regions of the relaxed convex problem. The non-convex problem can be reformulated using the auxiliary variables without affecting the feasible region of the problem $\mathcal{F}_{samp}:=\{x\in\mathcal{X},~{}y\in\mathbb{R}^{n_{b}},u\in\hat{\mathcal{U}}:g(x,u)\leq 0,~{}y=x_{i}x_{j},~{}i,j=1,\cdots,B\}$. Hence, the relaxed feasible region is $\tilde{\mathcal{F}}_{samp}=\{x\in\mathcal{X},~{}y\in\mathbb{R}^{n_{b}},u\in\hat{\mathcal{U}}:g(x,u)\leq 0,~{}\bm{\mu}(x,y)\leq 0\}$ which is a superset of the original problem $\mathcal{F}_{samp}\subseteq\tilde{\mathcal{F}}_{samp}$. From Proposition 1 the relaxed problem provides also a lower bound for the robust problem $\mathcal{F}_{rob}\subseteq\mathcal{F}_{samp}\subseteq\tilde{\mathcal{F}}_{samp}$. ∎

In the Robust spatial Branch and Bound (RsBB) algorithm we consider a two-level problem approach as in the robust cutting plane algorithm. The outer-level is comprised of the deterministic problems ($P_{samp}$) and ($\tilde{P}_{samp}$). The inner-level (T) problems evaluate the robustness of the solution of the outer-level non-convex problem. The RsBB algorithm is initialised in the root node for the original variable bounds $\mathcal{X}_{0}=\mathcal{X}$. Note that for both outer-level problems the domain $\mathcal{X}$ and the subset $\hat{\mathcal{U}}$ are the same at each node. The domain $\mathcal{X}$ is partitioned into subsets based on the sBB steps of the algorithm and the subset $\hat{\mathcal{U}}$ is augmented by any uncertainty samples that violate the robustness of the obtained solution. For a given node, $n$ problems $P_{samp}^{n}$ and $\tilde{P}_{samp}^{n}$ are defined for $\mathcal{X}_{n}$ and the corresponding $\mathcal{F}_{n},~{}\tilde{\mathcal{F}}_{n}$. The proposed algorithm for the RsBB implementation is displayed in Algorithm 3 and the Infeasibility Test in Algorithm 2.

Lemma 1.

The Infeasibility Test Algorithm 2 terminates in a finite number of iterations.

Proof.

Algorithm 2 is applied to nodes for which a solution to $(P_{samp}^{n})$ is obtained. Since $(P_{samp}^{n})$ is assumed to be tractable, the algorithm remains well-defined at every iteration. Moreover, Algorithm 2 follows the same structure as Algorithm 1, which corresponds to the robust cutting set algorithm. By Assumptions 2 and 5, along with the tractability of $(P_{samp}^{n})$, the conditions required for finite convergence of Algorithm 1 are satisfied. Thus, by the results in [8] (Section 5.2), Algorithm 2 terminates in a finite number of iterations. ∎

The effectiveness of the sBB algorithm heavily depends on the choice of bounding, branching, and node update strategies. Our bounding strategy relies on McCormick relaxations [39, 40], which generate degenerate perfect sequences of under-estimators [41][Theorem 5]. For branching, we adopt the bisection method, which ensures an exhaustive subdivision of the original variable domain [42][Definition IV.10]. The node update strategy follows a best-first approach, selecting the node with the lowest lower bound. This method is node-improving [42](Definition IV.6), meaning it prioritises promising candidates for early exploration [43]. An additional decision lies in the choice of branching variable. In this work we consider the maximum violation strategy [44, 45, 46] and a pseudoscore strategy [47, 46]. The use the maximum approximation error strategy is prioritised and if the errors are below the selected tolerance, then the branching variable is selected based on the pseudoscore. We define the branching variable selected via the maximum approximation error as:

$$x_{br}=\operatorname*{arg\,max}_{{i,j}}(\max|y_{ij}-x_{i}x_{j}|),~{}i,~{}j\in\{1,\cdots,n_{x_{b}}\}$$ (1)

If the maximum approximation error is below the selected tolerance, the branching variable is selected based on a pseudoscore. For the pseudoscore calculation, strong branching is taking place on the evaluating node $n$ and the variable with the most promising lower bound improvement is selected. Using bisection to define the branching point $x_{i}^{b}$ at node $n$, the intervals for child nodes can be defined as $\delta_{i}^{n,+}=\overline{x}_{i}-x_{i}^{b}$ and $\delta_{i}^{n,-}=x_{i}^{b}-\underline{x}_{i}-$. Let $\phi_{i}^{n,+},~{}\phi_{i}^{n,-}$ lower bound values of the child nodes for branching in variable $i$ and $\phi^{n}$ the lower bound of the parent node. If any of the child nodes is infeasible, the corresponding lower bound is set to 0. The lower bound improvement of the pseudonodes is evaluated as $\Delta_{i}^{+}=\phi_{i}^{n,+}-\phi^{n}$ and $\Delta_{i}^{+}=\phi_{i}^{n,-}-\phi^{n}$. Let $\varsigma_{i}^{+}=\Delta_{i}^{+}/\delta_{i}^{+}$ be the objective gain for branching on variable $i$. Let $x_{br}$ be the variable selected for branching via pseudoscores:

$$x_{br}=\operatorname*{arg\,max}(\varsigma_{i}^{+},\varsigma_{i}^{-}),~{}i\in\{1,\cdots,n_{x_{b}}\}$$ (2)

Proposition 4.

The RsBB algorithm 3 terminates after finitely many steps.

Proof.

A sBB algorithm terminates after finitely many steps if the bounding operation is finitely consistent as stated in [42][Theorem IV. I]. Let $X_{n}$ denote the subset of $x$ corresponding to node $n$. Then $\beta(\mathcal{X}_{n})$ be the lower bound and $\alpha(\mathcal{X}_{n})$ the upper bound of $f(x,y),x\in\mathcal{X}_{n}$. A bounding operation is finitely consistent if $\underset{n\rightarrow\infty}{\lim}(\alpha(\mathcal{X}_{n})-\beta(\mathcal{X}_{n}))$. The McCormick relaxations are known to be degenerate perfect as established by [41][Theorem 5], which implies that $\alpha(\mathcal{X}_{n})\rightarrow f(x,y),x\in\mathcal{X}_{n}$ and $\beta(\mathcal{X}_{n})\rightarrow f(x,y),x\in\mathcal{X}_{n}$. Since both bounds converge to $f(x,y)$ the condition for finite consistency is satisfied. ∎

Proposition 5.

The optimal solution obtained after the exhaustion of the RsBB tree is also optimal for the robust problem.

Proof.

Let $z_{n}$ denote the objective value of the sampled problem at node $n$, and let $z^{BF}$ be the best feasible objective value found by the RsBB algorithm up to this node. Let $x^{*}_{n}$ be the corresponding optimal solution for the sampled problem at node $n$. By construction, RsBB performs a robustness evaluation for all nodes where $z_{n}\leq z^{BF}$. This ensures that any solution $x^{*}_{n}:z_{n}\leq z^{BF}$ is also feasible for the robust problem. Hence, the final optimal solution $x^{*}$ obtained by RsBB is necessarily feasible for the robust problem. By Proposition 2, the RsBB solution $x^{*}$ is also optimal for the robust problem. ∎

Algorithm 2 Infeasibility Test Algorithm

1:Solution of $(P_{samp}^{n})$ at node $n$, i.e., $x^{*}$

2:Subset of sampled uncertain parameters $\hat{\mathcal{U}}$

3:Updated set of sampled uncertain parameters $\hat{\mathcal{U}}^{\prime}$

4:Solve (T) $\rightarrow u^{*}$

5:if $g(x^{*},u^{*})>0$ then

6:     Add new sample to $\hat{\mathcal{U}}^{\prime}=\hat{\mathcal{U}}\cup\{u^{*}\}$

7:     Solve $(P_{samp}^{n})\rightarrow x^{*}$

8:     Go to Step 1

9:end if

10:Terminate: Robust node, continue to BB

Algorithm 3 The RsBB Algorithm

1:Non-convex problem $P_{samp}^{n}$ at node $n$ with feasible region $\mathcal{F}^{n}$ and optimal value $z_{n}$

2:Relaxed convex problem $\tilde{P}_{samp}^{n}$ at node $n$ with feasible region $\mathcal{F}^{n}$ and optimal value $\tilde{z}_{n}$

3:Set of waiting nodes $\mathcal{W}:=\{0\}$, set of current nodes $\mathcal{C}:=\emptyset$

4:Set of uncertainty samples $\hat{\mathcal{U}}:=\{u_{0}\}$

5:Initialize complicating variable bounds $\underline{x}_{i},\overline{x}_{i}$ for $i=1,\dots,B$

6:Initialize best possible solution $z^{BP}\leftarrow\infty$ and best found solution $z^{BF}\leftarrow\infty$

7:while $\mathcal{W}\neq\emptyset$ do

8:     $z^{BP}\leftarrow\underset{n\in\mathcal{W}}{\min}{\tilde{z}_{n}}$

9:     Reset $\mathcal{C}\leftarrow\emptyset$

10:     for $n\in\mathcal{W}$ do

11:         if $\tilde{z}_{n}=z^{BP}$ then

12:              Add current node to $\mathcal{C}\leftarrow\mathcal{C}\cup\{n\}$

13:         end if

14:     end for

15:     for $n\in\mathcal{C}$ do

16:         Remove $n$ from waiting nodes $\mathcal{W}\leftarrow\mathcal{W}\setminus\{n\}$

17:         Solve $P_{samp}^{n}\rightarrow z_{n}$

18:         if $z_{n}\leq z^{BF}$ then

19:              Perform Infeasibility Test (Algorithm 2) $\rightarrow\hat{\mathcal{U}}^{\prime}$

20:              if New violations detected: $|\hat{\mathcal{U}}^{\prime}|>|\hat{\mathcal{U}}|$ then

21:                  Update $\hat{\mathcal{U}}\leftarrow\hat{\mathcal{U}}^{\prime}$

22:                  Solve $\tilde{P}_{samp}^{n},~{}\forall n\in\mathcal{W}\rightarrow\tilde{z}_{n}$

23:              end if

24:              Set $z^{BF}\leftarrow\min\{z_{n},z^{BF}\}$

25:              Solve $\tilde{P}^{n}_{samp}\rightarrow\tilde{z}_{n}$

26:              Select branching variable

27:              Generate child nodes $n^{\prime},n^{\prime\prime}$

28:              Append $\mathcal{W}\leftarrow\mathcal{W}\cup\{n^{\prime},n^{\prime\prime}\}$

29:              for $w\in\mathcal{W}$ do

30:                  if $\tilde{z}_{w}(1-\epsilon)\geq z^{BF}$ then

31:                       Fathom node $\mathcal{W}\leftarrow\mathcal{W}\setminus\{w\}$

32:                  end if

33:              end for

34:         end if

35:     end for

36:end while

37:Terminate: Robust optimal solution obtained at $z^{BF}$

3.3 Motivating example

The implementation of the RsBB algorithm is demonstrated in a non-convex QCQP toy problem $(P_{toy})$ defined by Eqs (3)-(7). We consider box uncertainty for $u\in\mathcal{U}=[2,6]$ around the nominal value $u_{0}=4$. The corresponding relaxed problem using McCormick envelopes will be denoted as $(\tilde{P}_{toy})$. The feasible regions for the toy nominal, sampled and robust problems are depicted in Fig. 1. The robust feasible region corresponds to a set of equations when no more feasibility violations are detected, while the sampled region corresponds to intermediate robust cutting plane iterations. In essence, at the termination of our algorithm the sampled and robust regions overlap.

$$\min-2x_{1}x_{2}$$ (3)

$$\text{s.t.}~{}u_{0}x_{1}x_{2}+2x_{1}+2x_{2}\leq 3$$ (4)

$$0.1-(x_{1}-0.5)^{2}-(x_{2}-0.5)^{2}\leq 0$$ (5)

$$x_{2}-0.09x_{2}-0.5\leq 0$$ (6)

$$\textbf{x}\in\mathcal{X}=[0,1]^{2}$$ (7)

The contour plots for the objective value in the feasible regions is displayed in Figure 2 which also entails the solutions of the examined nodes in rhombuses. For the root node and node 1 a trivial solution was found by the local solver hence no robust cuts were added to the problem. Node 2 obtained the optimal solution for the nominal problem resulting in an infeasibility cut. The robust optimal solution was obtained on the same node for a tighter sampled feasible region. The descendand nodes are not depicted on the figure but can be found in the detailed RsBB tree representation of Figure 3. The implementation of the RsBB algorithm is the following:
Step 0: Initialise the sampled uncertainty set $\hat{\mathcal{U}}={4}$ with the nominal uncertainty value. Solve $(P_{toy})$ with local solver. The obtained solution at root node is $z_{0}=0$ which is set as the best found $z^{*}=0$. The solution trivially satisfies the uncertain constraint in Eq.(4) hence no cuts are added.
Step 1: Solve the relaxed problem $(\tilde{P}_{toy})$ and generate two new (1,2) child nodes based on maximum approximation error.
Step 2: Select new node as $n=1$ based on lowest-lower bound $\tilde{z}_{1}<\tilde{z}_{2}$.
Step 3: Solve $(P_{toy})$ with local solver. The obtained optimal value is $z_{1}=0$. The Infeasibility Test is trivial for this node as well.
Step 4: Solve the relaxed problem $(\tilde{P}_{toy})$ and generate two new (3, 4) child nodes based on maximum approximation error. None of the waiting nodes is fathomed.
Step 5: Select new node as $n=2$ based on lowest-lower bound $\tilde{z}_{2}<\tilde{z}_{3}<\tilde{z}_{4}$.
Step 6: Solve $(P_{toy})$ with local solver. The obtained optimal value is $z_{2}=-0.45$.
Step 7: Perform the Infeasibility Test. The sampled point $u^{*}=6$ was found to violate Eq. (4) and was appended in $\hat{\mathcal{U}}\cup\{6\}$.
Step 8: Solve $(P_{toy})$ with local solver and the new feasible region. The obtained optimal value is $z_{2}=-0.36$.
Step 9: No more violations are detected by the infeasibility test, hence we set the new best found solution as $z^{*}_{2}=-0.36$.
Step 10: Solve the relaxed problem $(\tilde{P}_{toy})$ and generate two new (5,6) child nodes based on maximum approximation error. Node 6 is fathomed due to infeasibility of the relaxed problem.
Step 11: Fathom nodes for which $\tilde{z}_{n}(1-\epsilon)>p^{*}$. Node 4 is fathomed.
Step 12: Select new node based on lowest-lower bound.
Steps 8 to 12 are repeated for the remaining nodes, the corresponding solutions for each node can be found in Figure 3. The RsBB terminates in node 14 and the robust optimal solution $z^{*}=-0.36$.

Figure 1: Feasible regions for nominal sampled and robust toy problems

Figure 2: Contour plot for toy problem over the nominal and sampled feasible regions. Rhombuses correspond to objective values found at denoted nodes.

Figure 3: RsBB tree for the toy problem. We denote as $z_{n}$ and $\tilde{z}_{n}$the solutions of $P_{toy}$ and $\tilde{P}_{toy}$ at node $n$. We note as $infeas$ the nodes for which $\tilde{P}_{toy}$ was found infeasible. Nodes with dashed lines correspond to fathomed nodes.

4 Problem description

For the purpose of this study, we evaluate the performance of RsBB for the pooling problems. We use the pq-formulation of the pooling problem, resulting in a quadratically constrained quadratic program (QCQP). bilinear terms are manifested both in the objective function and the quality constraints.

4.1 The pooling problem

In the standard pooling formulation, given an existing infrastructure, the optimal flows between input streams, mixing pools and output product streams need to be determined. The objective is to minimise the cost of the process under storage capacity, demand satisfaction and product quality constraints. [48] proposed the q-formulation of the pooling problem in order to derive smaller in size dual problems. Later on, [49] proposed an alternative pq-formulation via the use of two RLT constraints which improved the performance of the state-of-the-art solvers at that time, while not changing the problem structure. We introduce the auxiliary variables $v_{ilj}$ to replace the bilinear terms. Let $\Pi_{nom}$ be the original non-convex pooling problem of cost minimisation:

$${\text{min}}~{}\sum_{ilj}c_{i}v_{ilj}-\sum_{j}d_{j}\sum_{l}y_{lj}+\sum_{j}d_{j}\sum_{i}z_{ij}$$ (8)

$$\text{s.t.}~{}\sum_{lj}v_{ilj}+\sum_{j}z_{ij}\leq A_{i},\quad\forall i$$ (9)

$$\sum_{j}y_{lj}\leq S_{l},\quad\forall l$$ (10)

$$\sum_{l}y_{lj}+\sum_{i}z_{ij}\leq D_{j},\quad\forall j$$ (11)

$$\sum_{i}v_{ilj}=y_{lj},\quad\forall l$$ (12)

$$\sum_{l}y_{lj}+\sum_{i}z_{ij}=f_{j},\quad\forall j$$ (13)

$$P_{jk}^{L}f_{j}\leq\sum_{il}C_{ik}v_{ilj}+\sum_{i}C_{ik}z_{ij}\leq P_{jk}^{U}f_{j},\quad\forall j,k$$ (14)

$$v_{ilj}=q_{il}y_{l,j}~{}\forall i,l,j$$ (15)

To derive the convex relaxation $\tilde{\Pi}_{nom}$ of the pooling problem $\Pi_{nom}$, Eq. (15) is replaced by the equivalent McCormick envelopes Eqs.(16)-(19) [32].

$$v_{ilj}\geq\underline{q}_{il}y_{lj}+q_{il}\underline{y}_{lj}-\underline{q}_{il}\underline{y}_{lj},\quad\forall i,l,j$$ (16)

$$v_{ilj}\leq\overline{q}_{il}y_{lj}+q_{il}\underline{y}_{lj}-\overline{q}_{il}\underline{y}_{lj},\quad\forall i,l,j$$ (17)

$$v_{ilj}\geq\overline{q}_{il}y_{lj}+q_{il}\overline{y}_{lj}-\overline{q}_{il}\overline{y}_{lj},\quad\forall i,l,j$$ (18)

$$v_{ilj}\leq\underline{q}_{il}y_{lj}+q_{il}\overline{y}_{lj}-\underline{q}_{il}\overline{y}_{lj},\quad\forall i,l,j$$ (19)

In the pooling problem, uncertainty can be manifested in any of the model parameters, such as the feed component cost, product price and feed component quality. For this study we consider uncertainty in the feed component quality, i.e. $C_{ik}$ in Eq. (14). Following the notation in [7], the true value of the uncertain level quality can be formulated as follows:

$$\tilde{C}_{ik}=C_{ik}+\xi_{ik}\hat{C}_{ik}$$ (20)

where $C_{ik}$ represents the nominal value, $\hat{C}_{ik}$ the constant perturbation (which is positive) and $\xi_{ik}$ is a vector of random variables which are subject to uncertainty. Given a selected uncertainty set $\mathcal{U}_{p}$, the deterministic quality constraints in (14) are replaced by their robust formulations in Eqs.(21) and(22) resulting in the semi-infinite programming problems $\Pi_{rob},~{}\tilde{\Pi}_{rob}$ [50].

$$\sum_{il}\tilde{C}_{ik}v_{ilj}+\sum_{i}\tilde{C}_{ik}z_{ij}\leq P_{jk}^{U}f_{j},\quad\forall\bm{\xi}\in\mathcal{U}_{p},~{}\forall j,k$$ (21)

$$P_{jk}^{L}f_{j}\leq\sum_{il}\tilde{C}_{ik}v_{ilj}+\sum_{i}\tilde{C}_{ik}z_{ij},\quad\forall\bm{\xi}\in\mathcal{U}_{p},~{}\forall j,k$$ (22)

4.2 Uncertainty set modelling

In this study we consider the case where $\mathcal{U}_{p}$ is a convex uncertainty set with $p\in\{\infty,2,1\}$, i.e. box, ellipsoidal and polyhedral uncertainty. [51] was the first to study linear problems under box uncertainty using the dual reformulation method, resulting in the most conservative solution, i.e. worst-case solution. Later on, [52] aiming at a less conservative approach, proposed the ellipsoidal uncertainty set defined by the $2$-norm. Despite the fact that the ellipsoidal set reduces the conservatism of the problem, it introduces extra nonlinear terms. To this end, [53] proposed an alternative formulation assuming that it it unlikely that all uncertain parameters would be affected by the maximum perturbation. The corresponding polyhedral set is defined by the $1$-norm. Given that Eqs. (21) and (22) are linear in $\xi_{ik}$ both dual reformulation and robust cutting planes can be used to address the robust problem $\Pi_{rob}$. The deterministic dual reformulations of Eqs. (21) and(22) derived in Eqs. (23) and(24) using the $\rho_{jk}$ according to the selected uncertainty set as displayed in Table 1.

Uncertainty set	$\mathcal{U}$	$\bm{\rho}_{jk}$
Box	$\{\bm{\xi}|~{}|\xi_{i}|\leq\Psi\,~{}\forall i\}$	$\Psi\big{(}\sum_{il}\hat{C}_{ik}q_{il}y_{lj}+\sum_{i}\hat{C}_{ik}z_{ij}\big{)}$
Ellipsoidal	$\{\bm{\xi}|~{}\sqrt{\sum_{i}\xi_{i}^{2}}\leq\Omega\}$	$\Omega\bigg{(}\sqrt{\sum_{i}\hat{C}^{2}_{ik}\big{(}\sum_{l}q_{il}y_{lj}+z_{ij}\big{)}^{2}}\bigg{)}$
Polyhedral	$\{\bm{\xi}|~{}\sum_{i}|\xi_{i}|\leq\Gamma_{k}\}$	$\Gamma_{k}\underset{i}{\max}\hat{C}_{ik}(\sum_{l}q_{il}y_{lj}+z_{ij})$

Table 1: Uncertainty set definitions and dual reformulations for the examined uncertainty set types.

$$\sum_{il}C_{ik}q_{il}y_{lj}+\sum_{i}C_{ik}z_{ij}-P_{jk}^{U}v_{j}+\rho_{jk}\leq 0$$ (23)

$$P_{jk}^{L}v_{j}-\sum_{il}C_{ik}q_{il}y_{lj}-\sum_{i}C_{ik}z_{ij}+\rho_{jk}\leq 0$$ (24)

For the robust cutting planes algorithm, consider the sampled quality constraints Eqs. (25), (26). The robustness of each $\Pi_{samp}$ solution is evaluated by the infeasibility test problem $T_{pool}$, $\forall j,k$ and both upper and lower quality constraints. For each sampled $\bm{\xi}$ that violates a quality constraint, the sampled uncertainty set $\hat{\mathcal{U}}_{p}$ is augmented by the corresponding $\bm{\xi}$ value.

$$\sum_{il}\tilde{C}_{ik}v_{ilj}+\sum_{i}\tilde{C}_{ik}z_{ij}\leq P_{jk}^{U}f_{j},\quad\bm{\xi}\in\hat{\mathcal{U}}_{p},~{}\forall j,k$$ (25)

$$P_{jk}^{L}f_{j}\leq\sum_{il}\tilde{C}_{ik}v_{ilj}+\sum_{i}\tilde{C}_{ik}z_{ij},\quad\bm{\xi}\in\hat{\mathcal{U}}_{p},~{}\forall j,k$$ (26)

$$\begin{rcases}\max~{}\sum_{il}\tilde{C}_{ik}v^{*}_{ilj}+\sum_{i}\tilde{C}_{ik}z^{*}_{ij}-P_{jk}^{U}f^{*}_{j}\\ \text{s.t.}\quad\tilde{C}_{ik}=C_{ik}+\xi_{ik}\hat{C}_{ik}\\ ~{}\bm{\xi}\in\mathcal{U}_{p}\end{rcases}{\forall j,k} T_{pool}$$ ()

Table 2 summarises the notation for the different problems and the equations comprising them providing a correspondence to the general problem notation used for the RsBB methodology.

Optimisation level	Problem	Notation	Equations
Outer-level	Nominal pooling	$\Pi_{nom}$	(8)-(15)
	Robust pooling	$\Pi_{rob}$	(8)-(13), (15), (21), (22)
	Sampled pooling	$\Pi_{samp}$	(8)-(13), (15), (25), (26)
	Relaxed nominal pooling	$\tilde{\Pi}_{nom}$	(8)-(14), (16)-(19)
	Relaxed robust pooling	$\tilde{\Pi}_{rob}$	(8)-(13), (16)-(19), (21), (22)
	Relaxed sampled pooling	$\tilde{\Pi}_{samp}$	(8)-(13), (16)-(19), (25), (26)
Inner-level	Quality robustness	$T_{pool}$	($T_{pool}$)

Table 2: Overview of the examined pooling problem categorised by optimisation level. Each problem is defined by its notation and the corresponding set of equations that describe its feasible region..

5 Computational experiments

5.1 Case studies

The performance of the RsBB algorithm is evaluated for 10 benchmark pooling instances, details of the problems can be found in Table 3. For the uncertain inlet quality, we set $\hat{C}_{ik}=C_{ik}$. The uncertain parameters are defined by box, ellipsoidal and polyhedral uncertainty sets for six different uncertainty set sizes, i.e. $\Psi,~{}\Omega,~{}\Gamma:=\{0.05,0.1,0.15,0.2,0.25,0.3\}$. The selected modelling environment is Pyomo v6.7 [54]. For the RsBB algorithm HiGHS v.1.5.3 [55] is used as the LP solver and conopt4 as the local NLP solver. The proposed algorithm is compared with PyROS v1.2.8 using Gurobi v11.1 [56] as the global solver. The deterministic robust counterparts are evaluated as well using Gurobi v11.1 as global solver. We set a time limit of 1 hour for all of the examined methods.

Problem	Feeds	Pools	Products	Qualities	Variables	Equations QP	Equations LP
haverly1	3	1	2	1	10	17	29
haverly2	3	1	2	1	13	20	38
haverly3	3	1	2	1	10	17	29
bental4	4	1	2	1	13	21	39
bental5	13	3	5	2	92	121	301
foulds2	6	2	4	1	36	46	94
adhya1	5	2	4	4	33	62	122
adhya2	5	2	4	6	33	70	130
adhya3	8	3	4	6	52	94	190
adhya4	8	2	5	4	58	95	215

Table 3: Model statistics for examined pooling problems

5.2 Results

A comparison of the RsBB algorithm with state-of-the-art methods in terms of % of problems solved within 1 hour of CPU time can be found in Figure 4. Dual-Gurobi corresponds to solving the dual problem, using the dual constraints for the different uncertainty types, with global solver Gurobi. PyROS-Gurobi corresponds to using PyROS solver with Gurobi as a global solver both for the upper and lower level problems. The RsBB algorithm is using conopt4 for the upper and lower level problems for box and polyhedral sets and BARON for the lower level problems under ellipsoidal uncertainty. For each of the 10 benchmark problems of Table 3 6 different uncertainty sizes are evaluated, resulting in 60 problems in total for each uncertainty type. The RsBB and Dual algorithms solve the same number of problems for all uncertainty types, i.e. 100% for box and polyhedral and 92% for ellipsoidal. In contrast, PyROS solves all problems only for box, 90% for ellipsoidal and 40% for polyhedral set. Dual reformulation method retrieves robust optimal solutions in less than 1 sec for all problems with box set and 80% of problems with ellipsoidal and polyhedral set. RsBB algorithm for box and polyhedral sets requires few seconds to less than 5 minutes to solve all instances. For the ellipsoidal set, RsBB relies on using BARON to solve the lower level problems resulting in an increase of the total computational time, rendering PyROS marginally faster. The results of Figure 4 suggest that integrating the robust cutting planes in a sBB framework can improve the performance of solely using robust cutting planes. Further optimising the RsBB algorithm via more rigorous BB techniques could further improve the computational efficiency of the proposed algorithm.

Figure 4: Performance profiles for Dual reformulation, PyROS and RsBB methods for different uncertainty types

Statistics for RsBB algorithm can be evaluated via the variability of the number of tree nodes explored in Figure 5. As a general trend, for all set types the range of the nodes explored reduces with the increase of the uncertainty set size. As the uncertainty size increases, the corresponding robust cuts become more restraining reducing the feasible region more rigorously. For uncertainty set sizes greater than 0.15 the number of nodes explored does not depend on the set type and all problems are solved evaluating less than 15 nodes. However, for smaller set sizes, the variability of the nodes explored increases with the decrease of the conservativism of the uncertainty set (box$>$ellipsoidal$>$polyhedral). The box uncertainty set corresponds to the worst-case uncertainty realisation, were all uncertain parameters can obtain the maximum allowed perturbation. For ellipsoidal and polyhedral sets more RsBB iterations are required to find all combinations of maximum parameter violations.

Figure 5: Variability of tree nodes explored for set size and set type for problems solved to robust optimality within 1 hour

Figure 6 illustrates the optimality gap closure of the RsBB algorithm for problems of varying dimensionality under a box uncertainty set with $\Psi=0.05$. The computational time for the depicted problems and total number of cuts is summarised in Table 4. An initial observation suggests that the number of nodes explored is not proportional to the problem size (haverly1$<$ foulds2$<$ adhya4$<$ bental5). Pooling problems have multiple local minima which are independent of the problem dimensionality, hence larger problems may be solved in less CPU time. This observation is further supported by the comparison of the corresponding CPU times in Table 4. Both in RsBB and Dual methods the higher CPU time is required for adhya4 followed by bental5, haverly1 and foulds2. While for PyROS, the problem size and CPU time are positively correlated, the CPU time entails the time required for Gurobi to find an optimal solution and the total robust cut iterations. The delay of PyROS can also be attributed to the smaller number of total cuts added compared to RsBB. In the proposed method, all uncertain parameter violations are included in the sampled uncertainty set, whereas PyROS considers only the uncertain parameters that cause the maximum violation. Selecting only a subset of violating cuts, rather than all, can be advantageous for larger problems where multiple constraints are impacted by uncertainty. This approach helps mitigate the exponential growth of constraints in the sampled problem. Another noteworthy aspect of Figure 6 is the significant variability in the optimality gap at the root node across the examined problems, which directly influences the algorithm’s convergence.

Figure 6: RsBB convergence for different pooling problems under box uncertainty set

Problem	CPU Dual [s]	CPU PyROS [s]	CPU RsBB [s]	Cuts PyROS	Cuts RsBB
haverly1	0.4	5.5	4.4	2	2
foulds2	0.3	20.0	4.2	5	8
adhya4	0.7	164.2	73.9	9	21
bental5	0.3	205.7	30.2	17	37

Table 4: Computational time and total cuts used to solve different pooling problems under box uncertainty set of size $\Psi=0.05$

Focusing on a specific problem, i.e. foulds2, we can investigate the dependence of the RsBB convergence on the uncertainty set size and set type as in Figure 7. Increasing the uncertainty size results in a decrease of the RsBB iterations and at the same time a decrease in the optimality gap at the root node. Fixing the same uncertainty size to 0.15, faster convergence is achieved for the box uncertainty set followed by the ellipsoidal and the polyhedral, analogous to the decrease of the set conservativism. For foulds2 with $\Psi=0.15$, RsBB terminates at the robust optimal solution on the root node, providing grounds for further investigation. Figure 8 depicts the increase of the objective value in robust cut iterations at the root node of foulds2 from Figure 7. For $\Psi=0.15$, RsBB performs 11 robust cut iterations to obtain the robust optimal solution at the root node. In contrast, for $\Psi=0.1$, fewer cut iterations lead to a higher optimality gap at the root node. For the polyhedral set and $\Gamma=0.15$ only 3 cut iterations take place in the root node, resulting in an initial gap of 9.5%. It is important to highlight that the objective values retrieved at the root node in 8b correspond to the robust optimal solution and the extra tree nodes are used to tighten the solution of the linear relaxation. For the box and polyhedral uncertainty sets solving the lower level problem has a finite number of solutions, which are the vertices of the polytopes. Finiteness has also been proved in the case of ellipsoidal set by [57], however as displayed in Figure 8b there can be multiple degenerate lower level violating solutions that do not affect the feasibility region of the upper level problem. The optimal objective value of foulds2 for the nominal uncertain parameter is at -1100. The right axis of Figure 8 depicts the % increase of the cost as different uncertainty set and sizes and considered. Allowing a 15% perturbation with the polyhedral set will result in an increase of 38% of the objective value, 50% for the ellipsoidal set and 80% for the box set. Throughout the examined results is has been highlighted that increasing the conservativism of the problem either via set size or set type, facilitates and accelerates the performance of the RsBB algorithm. Notably, robust cuts serve as a powerful feasibility bounds tightening method, significantly reducing the feasible region for the sBB tree nodes, which enhances the overall efficiency and accuracy of the algorithm.

Figure 7: foulds2 convergence a) for different box sizes, b) and different set types

Figure 8: foulds2 objective value in robust cut iterations at root node a) for different box sizes, b) and different set types

6 Conclusion

In this work, we introduced a novel algorithm that integrates global and robust optimisation methods. The proposed RsBB algorithm combines the principles of spatial branch-and-bound (sBB) and robust cutting planes to solve continuous non-convex problems under convex uncertainty sets. We applied our approach to QCQP pooling problems utilising McCormick envelopes to obtain convex lower bounds. The inlet quality parameter was treated as uncertain, with its range defined by box, ellipsoidal, and polyhedral uncertainty sets of varying sizes. The performance of the RsBB algorithm was compared to state-of-the-art methods that rely on global solvers. Our computational results demonstrate that the proposed method can solve the 97% of the tested problems to robust optimality while being computationally competitive with the state-of-the-art. We observed that using robust cutting planes in the sBB algorithm serves as an effective feasibility bounds tightening method. Finally, the integration of robust cuts is proven to have a beneficial impact on the sBB algorithm facilitating the optimality search both in terms of problem tractability and CPU time. Ongoing work in our work is looking into extending this approach to general non-convex problems.

Nomenclature

Sets

$i$	Feed components
$l$	Pools
$j$	Products
$k$	Qualities

Parameters

$\xi_{k}$	Random uncertain variable of quality
$\Psi$	Adjustable parameter controlling the size of the uncertainty set
$A_{i}$	Upper bound for component availability
$D_{j}$	Upper bound for product demand
$S_{l}$	Upper bound for pool size
$P_{jk}^{U}$	Upper bound for product quality
$P_{jk}^{L}$	Lower bound for product quality
$C_{ik}$	Nominal level of quality for component
$\tilde{C}_{ik}$	True level of quality for component
$\hat{C}_{ik}$	Perturbed level of quality for component
$c_{i}$	Unit price for component
$d_{j}$	Unit price for product

Variables

$q_{il}$	Fraction of flow of component $i$ entering pool $l$
$y_{lj}$	Flow from pool $l$ to product $j$
$z_{ij}$	Flow from component $i$ to product $j$
$f_{j}$	Total flow to product $j$
$v_{ilj}$	Auxiliary variable for flow from component $i$ passing pool $l$ to product $j$

Data availability statement Data and software will be made available by the authors upon reasonable request. Acknowledgements Financial support under the EPSRC grants ADOPT (EP/W003317/1) and RiFtMaP (EP/V034723/1) is gratefully acknowledged by the authors.