Enhanced branch-and-bound algorithm for chance constrained programs with Gaussian mixture models

Abstract

We study a class of chance constrained programs (CCPs) where the underlying distribution is modeled by a Gaussian mixture model. As the original work, Hu et al. (IISE Trans 54(12):1117–1130, 2022. https://doi.org/10.1080/24725854.2021.2001608) developed a spatial branch-and-bound (B &B) algorithm to solve the problems. In this paper, we propose an enhanced procedure to speed up the computation of B &B algorithm. We design an enhanced pruning strategy that explores high-potential domains and an augmented branching strategy that prevents redundant computations. We integrate the new strategies into original framework to develop an enhanced B &B algorithm, and illustrate how the enhanced algorithm improves on the original approach. Furthermore, we extend the enhanced B &B framework to handle the CCPs with multiple chance constraints, which is not considered in the previous work. We evaluate the performance of our new algorithm through extensive numerical experiments and apply it to solve a real-world portfolio selection problem.

Appendices

Appendix A: The B &B algorithm

Appendix B: Proof of Proposition 3

The proof of Proposition 3 is adapted from the proof of Proposition 3 in Pang et al. (2023).

For any ${\varvec{x}}\in {\mathcal {X}}$ and $y_{pj}\in \Delta ^0$,

$$\begin{aligned} F_{pj}({\varvec{x}},y_{pj})-F_{pj}({\varvec{x}},{{\underline{y}}}^0_{pj})&=\left( \Phi ^{-1}( y_{pj})-\Phi ^{-1}({\underline{y}}^0_{pj})\right) \sqrt{f^p({\varvec{x}})^\top \Sigma _j f^p({\varvec{x}})}\\&=\left( y_{pj}-{\underline{y}}^0_{pj}\right) \frac{d}{dy}\Phi ^{-1}(y)\vert _{y={\hat{y}}_{pj}}\sqrt{f^p({\varvec{x}})^\top \Sigma _j f^p({\varvec{x}})}\\&\le ({{\overline{y}}}^0_{pj}-{\underline{y}}^0_{pj})d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}, \end{aligned}$$

where ${\hat{y}}_{pj}\in l _{y^0_{pj}}$. The last inequality holds since the derivative of $\Phi ^{-1}(y)$ is $\frac{1}{\phi \left( \Phi ^{-1}(y)\right) }$ and $\lambda ^{\star } I-\Sigma _j$ is positive semi-definite and thus $\lambda ^{\star } f({\varvec{x}})_{\max }^2 \ge \lambda ^{\star }\Vert f^p({\varvec{x}})\Vert ^2 \ge f^p({\varvec{x}})^{\top }\Sigma _j f^p({\varvec{x}})$.

Suppose that a domain $\Delta =\prod _{p=1}^{P}\prod _{j=1}^{K} l _{pj}$ satisfies

$$\begin{aligned} {{\overline{y}}}_{pj}-{\underline{y}}_{pj}\le \frac{\epsilon }{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}}, \forall j\in {\mathbb {K}}, \quad \forall p\in {\mathbb {P}}. \end{aligned}$$

If Problem (13) formulated over $\Delta$ is feasible, denote the optimal solution as $({{\varvec{x}}}^c,{{\varvec{y}}}^c)$. Since $F_{pj}({{\varvec{x}}}^c,{{{\overline{y}}}}_{pj})-F_{pj}({{\varvec{x}}}^c,{{\underline{y}}}_{pj})\le \epsilon$ for all $j\in {\mathbb {K}}$ and $p\in {\mathbb {P}}$, the solution is an $\epsilon$-optimal solution. If Problem (13) formulated over $\Delta$ is infeasible, Problem (2) formulated over $\Delta$ must be infeasible. In both cases, the domain $\Delta$ will not be divided into two smaller subdomain.

According to the original B &B algorithm, the initial domain $\Delta ^0$ will not be divided into two smaller rectangles from the middle point of the longest edge if this edge is not larger than $\frac{\epsilon }{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}}$. Now, we consider the p-jth edge of a subdomain with the length of $\vert l _{y_{pj}}\vert ={{\overline{y}}}_{pj}-{\underline{y}}_{pj}> \frac{\epsilon }{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}}$ and then the partitioning continues. At the beginning of the B &B algorithm, we note that there is only one domain with $\vert l _{y^0_{pj}}\vert$ edge, i.e., the initial domain $\Delta ^0$. According to the original branching rule, partitioning in the middle to generate all domains with the length of this edge equal to $\left( \frac{1}{2}\right) ^{1}\vert l _{y^0_{pj}}\vert$, it is only required $2^0=1$ time of division. There are two subdomains generated with $\left( \frac{1}{2}\right) ^{1}\vert l _{y^0_{pj}}\vert$. If all subdomains on this edge are equal to $\left( \frac{1}{2}\right) ^2\vert l _{y^0_{pj}}\vert$, it is required $2^0+2^1=3$ times of division. It is not difficult to conclude that $2^n-1$ times of division can obtain all subdomains that the length of this edge is equal to $\frac{1}{2}^{n}\vert l _{y^0_{pj}}\vert$.

Let $n^{\star }=\log _2\left\lfloor \frac{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}\left| { l _{y^0_{pj}}}\right| }{\epsilon }+1\right\rfloor$. Then the subdomians with the length of the p-jth edge are equal to

$$\begin{aligned} \left| l _{y^0_{pj}}\right| \frac{1}{2}^{\log _2\left\lfloor \frac{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}\left| { l _{y^0_{pj}}}\right| }{\epsilon }+1\right\rfloor }&=\left| l _{y^0_{pj}}\right| \frac{1}{\left\lfloor \frac{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}\left| l _{y^0_{pj}}\right| }{\epsilon }+1\right\rfloor }\\&<\left| l _{y^0_{pj}}\right| \frac{1}{ \frac{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}\left| l _{y^0_{pj}}\right| }{\epsilon }}\\&=\frac{\epsilon }{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}}. \end{aligned}$$

Therefore this edge is no longer divided and the corresponding number of the partition is

$$\begin{aligned} 2^{n^\star }-1=\left\lfloor \frac{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}\vert l _{y^0_{pj}}\vert }{\epsilon }+1\right\rfloor -1 =\left\lfloor \frac{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}\vert l _{y^0_{pj}}\vert }{\epsilon }\right\rfloor \end{aligned}$$

Note that the each division corresponds to one computation of an edge of $\Delta ^0$ for Problem (13). This implies that the number of division is equal to the number of edges. Therefore, the number of Problem (13) solved is at most

$$\begin{aligned} \left\lfloor \frac{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}\vert l _{y^0_{11}}\vert }{\epsilon }\right\rfloor \times \cdots \times \left\lfloor \frac{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}\vert l _{y^0_{PK}}\vert }{\epsilon }\right\rfloor \le \left\lfloor \frac{d^{\star }f({\varvec{x}})_{\max }\sqrt{\lambda ^{\star }}\vert l ^\star _{y}\vert }{\epsilon }\right\rfloor ^{P\times K}. \end{aligned}$$

The proof is completed.

Appendix C: Comparing the performance of B &B and enhanced B &B algorithms

In Fig. 8, the y-axis represents the computation time, with the blue dotted line and orange line marked with stars stand for the B &B algorithm and enhanced B &B algorithm, respectively. The computation times of the B &B algorithm and the enhanced B &B algorithm under $K\in \{2, 5\}$, $d=50$ are $\{0.2\ \textrm{s}, 0.5\ \textrm{s}\}$ and $\{30.1\ \textrm{s}, 171.4\ \textrm{s}\}$ separately. The two lines almost coincide, indicating that each iteration of the two algorithms take roughly the same amount of time. However, the enhanced B &B algorithm avoids solving for redundant solutions, thus the number of iterations is largely reduced, resulting in significant advantage in computation time.

Fig. 8

The computation times of the B &B and the eB &B for Problem (6) with quadratic constraints (nonconvex subproblem) for $\alpha =0.1$ and $d=50$: a $K=2$ and b $K=5$

Appendix D: The statistics information of 48 Industry Portfolios

In Table 12, the abbreviations “M”, “S.D.”, “SK” and “K” refer to “Mean”, “Standard deviation”, “Skewness” and “Kurtosis”. We employ the normal distribution to model the distribution of the average monthly return rates of the 48 Industry Portfolio, as depicted in Fig. 9. It is noteworthy that the probability density function (pdf) manifests a long-tail characteristic, a feature not entirely captured by a normal distribution fit.

Fig. 9

Modeling the average monthly return rates of the 48 Industry Portfolios using normal distribution

Table 12 Descriptive statistics of monthly return rate of 48 Industry Portfolios