A parallel hub-and-spoke system for large-scale scenario-based optimization under uncertainty

Abstract

Practical solution of stochastic programming problems generally requires the use of parallel computing resources. Here, we describe the open source package mpi-sppy, in which efficient and scalable parallelization is a central feature. We report computational experiments that demonstrate the ability to solve very large stochastic programming problems—including mixed-integer variants—in minutes of wall clock time, efficiently leveraging significant parallel computing resources. We report results for the largest publicly available instances of stochastic mixed-integer unit commitment problems, solving to provably tight optimality gaps. In addition, we introduce a novel software architecture that facilitates combinations of methods for accelerating convergence that can be combined in plug-and-play manner. The mpi-sppy package is written in Python, leverages the widely used Pyomo (http://www.pyomo.org) library for modeling mathematical programs, builds on existing MPI implementations to ensure efficiency and scalability, and is available via http://github.com/Pyomo/mpi-sppy.

Appendices

A Details on the NetDes instances

We use the following network design problem (termed “NetDes”) for computational testing. Given a directed graph $G=(V,E)$, the decision maker selects a subset of the arcs E to build. Once the arcs have been built, two sets of uncertain parameters are realized: a demand value for each node (which may be positive or negative, i.e., a demand or a supply), and a capacity for each arc. The decision maker next selects the amount of flow to be sent along each arc that was built in the first stage, in order to satisfy the realized demand with the minimum expected cost, while respecting the capacity limit on each edge. This problem can be formulated as the following mixed-integer program:

$$\begin{aligned} \min \;&\sum _{e\in {E}}c_ex_e+\sum _{\xi \in \varXi }\Pr (\xi )\sum _{e\in {E}}f_e(\xi ){y_e(\xi )} \end{aligned}$$

(11a)

$$\begin{aligned} \mathrm {s.t.}\;&\sum _{e\in {N^+_v}}y_e(\xi )-\sum _{e\in {N^-_v}}y_e(\xi )=d_v(\xi ){} & {} \forall \;v\in {V},\xi \in \varXi \end{aligned}$$

(11b)

$$\begin{aligned}&y_e(\xi )\le {u_e(\xi )}x_e{} & {} \forall \;e\in {E},\xi \in \varXi \end{aligned}$$

(11c)

$$\begin{aligned}&x_e\in \{0,1\}{} & {} \forall \;e\in {E}\end{aligned}$$

(11d)

$$\begin{aligned}&y_e(\xi )\ge 0{} & {} \forall \;e\in {E},\xi \in \varXi \end{aligned}$$

(11e)

Here, $x_e=1$ if and only if edge $e\in {E}$ is built in the first stage, and $y_e(\xi )$ represents the amount of flow sent on arc $e\in {E}$ under scenario $\xi \in \varXi$, which occurs with probability $\Pr (\xi )$. The parameter $c_e$ is the cost to build edge $e\in {E}$, and $f_e(\xi )$ is the cost to send one unit of flow along arc $e\in {E}$ under scenario $\xi \in \varXi$. The set $N^+_v$ (resp. $N^-_v$) is the set of directed edges entering (resp. exiting) node $v\in {V}$. The constraints (11b) enforce conservation of flow for each node $v\in {V}$ and scenario $\xi \in \varXi$, where $d_v(\xi )$ is the “demand” (or, possibly, supply) of the node v under scenario $\xi \in \varXi$. The constraints (11c) ensure that we do not send any flow on an arc which we did not build, and further that the flow on each arc is less than the capacity $u_e(\xi )$.

The network design instances used in this paper were generated using the NETDES algorithm [20]. All instances are feasible under all scenarios—that is, there always exists a subset of edges which may be selected in the first stage that allows for a flow which satisfies the demand in all scenarios $\xi \in \varXi$.

B Scalability and overhead for small examples

To study the use of the library on shared-memory computers, we make use of the farmer instance from [3] modified using two instance creation parameters cropsmult and numscens. The original problem has three crops and three scenarios. The scalable instances have cropsmult sets of the original three crops with the same non-stochastic characteristics as in the original problem. Scenarios are also created in groups of three with a pseudo random number that is uniformly distributed added to the values for the original three scenarios. The problems have no practical interest, but they are easy to describe, scalable in number of scenarios, and size of scenario sub-problems and non-trivial to solve.

Table 5 gives results for $cropsmult=100$ and $numsens=2048$. Rather than solving individual scenarios as sub-problems, the hub PH algorithm solves bundles of 32 scenarios (see, e.g., [7, 29] for more discussion of bundling with PH). The spokes are Lagrangian and xhat-shuffle-looper. We limited the gurobi solver to one thread on each rank because we are not interested in solving the problem, just in looking at scalability. The PH algorithm is allowed to run for 10 and 100 iterations. We conducted our experiments on a workstation with 48 dual threaded Intel Xeon workstation operating at 2.1GHz, and 1TB of RAM. The column "Hub Ranks" gives the number ranks assigned to the hub, so the total number of ranks is three times that because there are two spokes.

Table 5 Fixed number of sub-problems, using bundling and varying the number of ranks for the scalable farmer instance with $cropsmult=100$ and $numsens=2048$

We now look at results of runs on a compute node with 24 dual threaded Intel CPUs operating at 2.3 GHz. Table 6 shows experiments with $cropsmult=1000$ and the number of processors set at three times the number of scenarios so that there is one scenario per rank in each of the three cylinders: hub, lower bound, upper bound. Our interest is not in solving the problem but in checking overhead and scalability on small systems.

Table 6 One sub-problem per stratum (i.e., one hub rank per sub-problem) for the scalable farmer instance with $cropsmult=1000$

The results in Tables 5 and 6 give information about overhead and scalability with respect to the number of sub-problems. There is also a practical question related to scalability with respect to the number of spokes. The effort to exchange data to spokes increases linearly with the number of spokes, so for a large enough number of spokes, the hub would spend most of its time writing and reading from spoke buffers and only some of its time working. In such cases, it would be necessary to implement a helper for the hub whose job would be to take data from the hub and distribute it to spokes.

The current implementation does not have such a helper and Table 7 illustrates that it is probably not needed for modest numbers of spokes. These experiments were done on a compute node with 24 dual threaded Intel CPUs operating at 2.3 GHz using scalable farmer with $cropsmult=1000$ and the number of processors set at 1 plus the number spokes so there are multiple scenarios per rank. For every experiment we have a PH hub cylinder and a Lagrangian lower bound spoke. We vary the number of xhat-shuffle-looper upper bound spokes. Our interest is not in solving the problem but in checking overhead and scalability on small systems. All experiments ran for 100 PH iterations. We see that increasing the number of spokes from 2 to 8 increases the runtime by less than 10%.

Table 7 Time in seconds for 100 iterations for the scalable farmer problem with $cropsmult=1000$

We note the coefficient of variation in run times for replicated experiments is less than 2%, so there is not much value in replicating these experiments.

C Details concerning aircond experiments

The aircond model is a multi-stage production planning problem with inventory, backorders (with linear and quadratic penalties), and production costs used for testing and experimentation. Full details are provided at https://github.com/DLWoodruff/aircond. This example makes use of a very new feature of mpi-sppy referred to as proper bundles. These bundles are created before the cylinder system runs, pickled (i.e., serialized), and stored for use during the execution by the cylinder system. Bundle pickling is parallelized and can use all ranks (e.g. 600 ranks for the 1 M experiments) and takes about 40 s for the experiments reported in this paper.

In the interest of repeatability, a slightly condensed form of the slurm script for 100k scenarios is given in Fig. 4.

Fig. 4

slurm script for the 100k scenario aircond example