{'fraction_non_alphanumeric': 0.02661321483771252, 'fraction_numerical': 0.006343379701184956, 'mean_word_length': 4.99604208900473, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient Stochastic Programming in Julia

We present StochasticPrograms.jl, a user-friendly and powerful open-source framework for stochastic programming written in the Julia language. The framework includes both modeling tools and structure-exploiting optimization algorithms. Stochastic programming models can be efficiently formulated using expressive syntax and models can be instantiated, inspected, and analyzed interactively. The framework scales seamlessly to distributed environments. Small instances of a model can be run locally to ensure correctness, while larger instances are automatically distributed in a memory-efficient way onto supercomputers or clouds and solved using parallel optimization algorithms. These structure-exploiting solvers are based on variations of the classical L-shaped and progressive-hedging algorithms. We provide a concise mathematical background for the various tools and constructs available in the framework, along with code listings exemplifying their usage. Both software innovations related to the implementation of the framework and algorithmic innovations related to the structured solvers are highlighted. We conclude by demonstrating strong scaling properties of the distributed algorithms on numerical benchmarks in a multi-node setup.


Introduction
Stochastic programming is an effective mathematical framework for modeling multi-stage decision problems that involve uncertainty [1]. It has been used to model complex real-world problems in diverse fields such as power systems [2,3,4], finance [5,6], and transportation [7,8]. The classical setting is linear stochastic programs where an actor takes decisions in two stages: initial decision x → observation ω → recourse action y(x, ξ(ω)) The actor first takes a decision x. Then, the realization of a random event ω alters the state of the world. The actor can observe ω and take a recourse action y with respect to x and the output of some random variable ξ(ω). We are interested in finding the optimal decision x, accounting for the ability to make a recourse action once ω has been observed. The notion of an optimal decision is captured by letting x and y be optimization variables in linear programs, where ξ(ω) parameterizes the second-stage problem for each event ω.
In applications, a stochastic program models some real-world decision problem under a statistical model of ξ. We can then compute approximations of optimal decision policies by solving approximated instances of the stochastic program. In brief, this involves computing a first-stage decisionx that is optimal in expectation over a set of secondstage scenarios ξ(ω i ) sampled from the model of ξ. This technique is known as sampled average approximation (SAA). In the linear setting, one can in principle formulate sampled instances on a extensive form that considers all available scenarios at once. This mathematical program can be solved using standard linear programming solvers, including both open-source solvers such as GLPK [9] and commercial solvers such as Gurobi [10]. However, the size of the extensive form grows linearly in the number of scenarios, and industry-scale applications typically involve 10,000+ scenarios. For example, the 24-hour unit commitment problem studied in [4] has 16,384 scenarios and the resulting extensive form has 4 billion variables. Solving the extensive form in such applications becomes practically infeasible. Moreover, the memory requirement for storing the stochastic program instances will eventually exceed the capacity of a single machine. This clarifies the need for a distributed approach for modeling large-scale stochastic programs. Structure-exploiting decomposition methods [11,12] that operate in parallel on distributed data become essential to solve large-scale instances.

Contribution
In this work, we present a user-friendly open-source software framework for efficiently modeling and solving stochastic programs in a distributed-memory setting. The framework allows researchers to formulate complex stochastic models and quickly typeset and test novel optimization algorithms. Stochastic programming educators will benefit from the clean and expressive syntax. Also, the framework supports analysis tools and stochastic programming constructs from classical theory and leading textbooks. Industrial practitioners can make use of the framework to rapidly formulate complex models, analyze small instances locally, and then run large-scale instances in production on a supercomputer or a cloud cluster. We implemented the framework in the Julia [13] programming language. Henceforth, we refer to the framework as SPjl. The framework is freely available through the registered Julia package StochasticPrograms.jl.
The design philosophy adopted during implementation of SPjl is centered around flexibility and efficiency, with the aim to provide a feature-rich and user-friendly experience. Also, the framework should be scalable to support largescale problems. With this in mind, we adhered to the fundamental principle that the optimization modeling should be separated from the data modeling. This design principle results in two key software innovations: deferred model instantiation and data injection. Optimization models are formulated in stages using a straightforward syntax that simultaneously specifies the data dependencies between the stages. The data structures related to future scenarios, and their statistical properties, are defined separately. An essential consequence of this design is that we can efficiently distribute stochastic program instances in memory, reducing interprocess communication to a minimum. Many computations involving distributed stochastic programs can then natively be run in parallel. Moreover, when the sample space is infinite, it becomes possible to adequately distinguish between the abstract representation of a stochastic program and finite sampled instances. The design also enables swift implementation of various constructs from classical stochastic programming theory. Another design choice is that the solver suites included in the framework are developed using policy-based techniques. We have shown in prior work how policy-based design can be used to create customizable and efficient optimization algorithms [14]. In short, SPjl is a powerful, versatile, and extensible framework for stochastic programming. It provides both an educational setting for new practitioners and a research setting were experts can further the field of stochastic programming.
We developed SPjl in Julia, which has several distinct benefits. Through just-in-time compilation and type inference, Julia can achieve C-like performance while being as expressive as similar dynamic languages such as Python or Matlab. Using the high-level metaprogramming capabilities of Julia, it is possible to create domain-specific tools with expressive syntax and high performance. Another benefit is access to Julia's large and rapidly expanding ecosystem of libraries, many of which play a central role in SPjl. For example, the parallel capabilities of SPjl are implemented using the standard library module for distributed computing, while optimization models are formulated using the JuMP [15] ecosystem. JuMP is an algebraic modeling language implemented in Julia using similar metaprogramming tools. It has been shown to achieve similar performance to AMPL [15], with syntax that is both readable and expressive. Also, it is possible to mutate model instances at runtime, which we utilize in the structure-exploiting algorithms. Recently, the backend of JuMP was redesigned into the new MathOptInterface [16]. The redesign introduces automatic reformulation bridges, which are used frequently in the current implementation of the SPjl framework. JuMP implements interfaces to many third-party optimization solvers, both open-source and commercial. These can be hooked in to solve extensive forms of stochastic programs or subproblems that arise in decomposition methods.

Related work
We give a short survey of similar software packages and highlight distinguishing features of SPjl. The most similar approach is the PySP framework [17], implemented in the Python language. Optimization models in PySP are created using Pyomo [18]; an algebraic modeling language also implemented in Python. In contrast, SPjl is written in the Julia language and formulates optimization models in JuMP, which has been shown to outperform Pyomo in various benchmarks [15]. In PySP, stochastic programs are composed of multiple .dat files and .py files, and the models are solved by running different solver scripts. In SPjl, all models are described in pure Julia and can be created, analyzed and solved in a single interactive session. Moreover, all operations are natively distributed in memory and run in parallel if multiple Julia processes are available. The parallel capabilities of PySP extend to running parallelized versions of the solver scripts. The primary function of PySP is to formulate and solve stochastic programs, while SPjl also provides a large set of stochastic programming constructs and analysis tools. The expressiveness of the modeling syntax can be compared by observing how the well-known farmer problem [1] is modeled using PySP [17] and how it is modeled using SPjl, as shown in Appendix 4. In particular, the PySP definition requires about 100 lines of code spread out over four different files, while SPjl requires 30 lines of code with the added benefit of being more readable. In addition, the resulting model can be analyzed interactively in Julia in a user-friendly way.
A more extensive list of similar software approaches is provided in [17], along with comparisons to PySP. This allows for a transitive comparison to SPjl. Other notable examples include the commercial FortSP solvers [19] coupled with the AMPL extension SAMPL for modeling. Out of all these approaches, SPjl has the most user-friendly interface and is also freely available.
The StructJuMP package [20] provides a simple interface to create block-structured JuMP models. The primary reason for developing StructJuMP was to facilitate a parallel modeling interface to existing structured solvers [21,4] that operate in computer clusters. These parallel solvers are implemented in C++ and are parallelized using MPI. This led to StructJuMP also making use of MPI to distribute stochastic programs in blocks. Apart from formulating distributed stochastic programs in a cluster, StructJuMP does not offer any modeling tools nor any way to generate the extensive form of a stochastic program. In comparison, SPjl provides numerous analysis tools as well as a compatible suite of structured solvers. In addition, SPjl natively distributes and solves stochastic programs using Julia, without relying on external software such as MPI.

Preliminaries
We give a short mathematical introduction to linear stochastic programming. The purpose is to provide background for the code examples presented in the subsequent section and also to keep this work self-contained. A more thorough introduction to the field is given in the textbook by [1].

Stochastic programming
A linear two-stage recourse model enables a simple but powerful framework for making decisions under uncertainty. We formalize this procedure in the following brief review. The first-stage decision made by the actor is denoted by x ∈ R n . We associate x with a linear cost function c T x that the actor pays after making the decision. Moreover, x is constrained to the standard polyhedron in linear programming, i.e.
where A ∈ R p×n and b ∈ R p . The recourse actions are represented by y ∈ R m . To describe the uncertainty in the decision problem, we consider some probability space (Ω, F , π) where Ω is a sample space, F is a σ-algebra over Ω and π : F → [0, 1] is a probability measure. Let ξ(ω) : Ω → R N be some random variable on Ω and let E ξ denote expectation with respect to ξ. We can now let ω ∈ Ω denote a scenario observed after deciding x. The scenario affects both cost and the constraints of the recourse action. Specifically, after realization of ω, the following second-stage problem is formulated to determine y with respect to x and ξ(ω): In other words, the random variable takes on the form ξ(ω) = (q ω T ω h ω ) T in this linear setting. Note that q ω ∈ R m , T ω ∈ R q×n and h ω ∈ R q are scenario-dependent while W ∈ R q×m is fixed. This is a standard setting in literature, which covers a wide range of problems [1]. It is possible to define W as scenario-dependent in the framework, but standard algorithms are then no longer certain to converge. Now, we formulate the two-stage recourse problem as follows. minimize The optimal value of (2) is referred to as the value of the recourse problem (VRP).
Apart from solving (2), we can compute two classical measures of stochastic performance. The first measures the value of knowing the random outcome before making the decision. This is achieved by taking the expectation in (2) outside the minimization, to obtain the wait-and-see problem: Now, the first-and second-stage decisions are taken with knowledge about the uncertainty. The difference between the expected wait-and-see value and the value of the recourse problem is known as the expected value of perfect information: The EVPI measures the expected loss of not knowing the exact outcome beforehand. It quantifies the value of having access to an accurate forecast.
Finally, we introduce the concept of decision evaluation to quantify the performance of a candidate first-stage decision x in the stochastic program (2). The expected result of x is given by This concept is used to compute the second classical measure. If the expectation in (2) is instead taken inside the second-stage objective function Q, we obtain the expected-value-problem: The solution to the expected-value-problem is known as the expected value decision, and is denote byx. The expected result of taking the expected value decision is known as the expected result of the expected value decision: The difference between the value of the recourse problem and the expected result of the expected value decision is known as the value of the stochastic solution: The VSS measures the expected loss of ignoring the uncertainty in the problem. A large VSS indicates that the second stage is sensitive to the stochastic data.
The EVPI, VSS, and VRP are important tools when gauging the performance of a stochastic model. All of these introduced measures are readily computed in the SPjl framework, which allows for easy analysis of user-defined models. Next, we discuss how to calculate the VRP, EVPI, and VSS depending on the form of the sample space Ω.

The finite extensive form and sample average approximation
If Ω is finite, say with n scenarios of probability π 1 , . . . , π n , then we can represent (2) compactly as We refer to this problem as the finite extensive form. It is often recognized in literature as the deterministic equivalent problem (DEP). Similar closed forms can be determined for the EVPI and the VSS. For small n, it is viable to solve this problem with standard linear programming solvers. For large n, decomposition approaches are required. In SPjl, the user provides a description of the abstract stochastic model (2) and a separate description of the uncertainty model of ξ. These are then combined internally to generate instances of the finite form (9), which are stored and solved efficiently on a computer or a compute cluster.
If Ω is not finite, the stochastic program (2) is exactly computable only under certain assumptions [1]. However, it is possible to formulate computationally tractable approximations of (2) using the finite form (9). The most common approximation technique is the sample average approximation (SAA) [22]. Assume that we sample n scenarios ω s , s = 1, . . . , n independently from Ω with equal probability. These scenarios now constitute a finite sample spaceΩ and we can use them to create a sampled model in finite extensive form (9). An optimal solution to this sampled model approximates the optimal solution to (2) in the sense that the empirical average second-stage cost V n = 1 n n s=1 q T sŷs , whereŷ s = arg min y∈R m {Q(x, ξ(ω s ))}, converges pointwise with probability 1 toV = E ξ [Q(x, ξ(ω))] as n goes to infinity [23]. Further, under certain assumptions it can be shown that in distribution as n goes to infinity [24]. This result provides a basis for calculating confidence intervals around the VRP of (2) [22,25], as well as around the EVPI and the VSS.

Structure-exploiting solvers
Efficient methods for storing and solving finite stochastic programs on the form (9) are key for high-performance stochastic programming computations. Therefore, this has been a main focus in the development of the SPjl framework. An important insight is that the finite extensive form (9) lends itself to block-decomposition approaches, which allow the stochastic program to be efficiently distributed in memory. Moreover, structure-exploiting solvers can be employed to solve the decomposed models efficiently. These approaches also readily extend to parallel settings where the stochastic program is distributed over several compute nodes. A key idea in the SPjl framework is to let the storage of the stochastic program depend on the type of optimizer used to solve it. In this way, the memory structure is optimized for the solver operation, and there is no redundant storage for other operations such as decision evaluation. We say that the underlying structure of the stochastic program is induced by the solver. Henceforth, we will refer to the treatment of (9) as one large optimization problem as the deterministic structure. This is the default structure for standard third-party solvers. For block-decomposition approaches, we adopt the terminology introduced in [17] and divide such strategies into two classes. In short, "vertical strategies decompose a stochastic program by stages" while "horizontal strategies decompose a stochastic program by scenarios" [17]. In the following, we will introduce two different solver algorithms that fall into these two categories and highlight the stochastic program structures they induce.

The L-shaped algorithm
The L-shaped algorithm is an efficient cutting-plane method for solving the finite extensive form (9) by decomposing into a master problem and a set of subproblems. The master problem has the form whereQ(x) is a lower bound on Here, each Q s (x) is the optimal value to a subproblem of the form (1). The idea of the L-shaped algorithm is to generate increasingly tight piecewise linear lower bounds on Q. We refer to the memory structure inferred by the L-shaped algorithm henceforth as the vertical structure.
During the L-shaped procedure, solution candidates x k are generated by solving the master problem (11), which are then used to parameterize subproblems of the form (1). Optimal dual variables in these subproblems are then used to improve the bound of Q(x) before the next solution candidate x k+1 is computed. Specifically, it follows from duality theory that λ T s (h s − T s x), where λ s is the dual optimizer of (1), is a valid support function for Q s (x), and hence, is a valid support function for Q(x). In the original formulation of the L-shaped algorithm [12], the above result is used at each iteration k to construct optimality cuts by introducing π s λ T s h s , and add to the master problem as the constraint ∂Q k x + θ ≥ q k . Aggregating the results from all subproblems in this way is known as the single-cut approach. This was later extended to a multi-cut variant where separate cuts are constructed for each subproblem [26]. If the iterate x k is not second-stage feasible, some subproblems will be infeasible. We handle this by solving the auxilliary problem: If w s > 0, then subproblem s is infeasible for the current iterate x k . Further, it follows from duality theory that where σ s is the dual optimizer of (10), is necessary for x to be second-stage feasible. The above result can be used to both check for second-stage infeasibility and construct feasibility cuts by introducing for all infeasible subproblems 1, . . . , f . Because W has a finite number of bases, finitely many feasibility cuts are required to completely describe the set of feasible first-stage decisions [12]. The optimality cuts and the feasibility cuts enter the master problem as follows: The master problem is then re-solved to generate the next iterate x k+1 , θ k+1 . This is repeated until the gap between the upper bound Q(x k ) and lower bound θ k+1 becomes small, upon which the algorithm terminates. Many variations can be introduced to improve the performance of the L-shaped algorithm. We provide an overview of such improvements available in SPjl in Section 4.1.

The progressive-hedging algorithm
The progressive-hedging algorithm was first introduced in [11]. In contrast to the L-shaped algorithm, applying progressive-hedging to solve (9) yields a complete decomposition over the n scenarios. The method is a specialization of the proximal-point algorithm [27], and convergence in the linear case (9) is derived in [11]. The main idea behind this approach is to introduce individual first-stage decisions x s to each scenario but force them to be equal. We then relax (dualize) these consistency constraints and solve the corresponding augmented Lagrangian problem. In other words, we consider the following problem: The consistency constraints x s = ξ, s = 1, . . . , n are called non-anticipative because they make the x s independent of scenario and enforce the fact that the first-stage decision is known when the second-stage uncertainty is realized. We refer to the memory structure inferred by the progressive-hedging algorithm henceforth as the horizontal structure. Separability across the n scenarios is achieved by introducing the following regularized relaxation of each subproblem: The algorithm now proceeds by iteratively alternating between generating new admissible solutions x k s , s = 1, . . . , n, and an implementable solution ξ k . In the two-stage setting, an admissible solution is feasible in every scenario, and an implementable solution is consistent in the sense that x s = ξ for all s. We obtain the implementable solution through aggregation: π s x k s and the Lagrange multipliers are updated scenario-wise through . Hence, the non-anticipative constraints are enforced while the dual variables converge. Progressive-hedging is a primal dual algorithm that is run until both the primal gap ξ k − ξ k−1

StochasticPrograms.jl
In this section, we showcase the capabilities of SPjl. We first give a brief overview of the framework and introduce the main functionality through a set of simple examples. Accompanying code excerpts are included. Next, we exemplify the effectiveness of SPjl model creation by giving a compact definition of the farmer problem. Finally, we summarize the algorithmic improvements and variations included in the framework.
SPjl extends the well-known JuMP syntax to support the definition of stages, decision variables, and uncertain parameters. Models are defined using the @stochastic_model macro. This creates a lightweight model object that can be used to instantiate finite stochastic programs by supplying a description of the uncertain parameters. Specifically, the user provides a list of discrete scenarios, or a sampler object capable of generating scenarios, to the model object.
The object then combines the model definition with the supplied uncertainty data and generates a finite stochastic program instance. The instantiated stochastic program can then be inspected, analyzed and solved in an interactive Julia session. This is useful in educational settings, but also for reasoning about complex models on a small scale. SPjl also supports reading problems specified in the SMPS format. For large-scale instances, SPjl provides scalable block-structured instantiation and structure-exploiting solvers that can operate in parallel. In addition, operations such as EWS calculation and decision evaluation are embarrassingly parallel over the subproblems. In other words, the workload is readily decoupled into independent subtasks that can be executed in parallel. This is leveraged when instantiating vertical or horizontal structures in distributed environments.
SPjl can be installed directly from the command line through Julia's package manager (pkg> add StochasticPrograms). Provided that a basic linear quadratic solver, such as GLPK or Ipopt, is installed, all code examples in this paper can be repeated by copying the lines verbatim. A more extensive introduction to the framework is given by the "Quick start" section of the online documentation 1 .
In SPjl, we create the model (13) in two steps. First, we formulate the optimization models as shown in Listing 1.
This creates a stochastic model where the two stages are given by the mathematical programs (13) and (14), expressed using an enhanced JuMP syntax. The @decision and @recourse lines work as standard @variable definitions in JuMP, but behind the scenes they also specify internal data dependencies between the first and second stage; and the @uncertain line annotates the random parameters and defines a point of data injection. The code specifies how the optimization models should be defined, but the actual model instantiation is deferred until we add a stochastic model of the uncertainties. We will consider two different distributions of ξ and use the same model object simple_model from Listing 1 to instantiate stochastic programs. This is a key feature in SPjl. The underlying stochastic model (2) object can be re-used to generate different finite stochastic program instances. Regardless of the distribution of ξ, a stochastic program instance is always a finite program of the form (9). This allows us to evaluate the same problem under different uncertainty models and to automatically adapt the underlying memory structure to optimize solver performance.

Finite sample space
First, let ξ be a discrete distribution, taking on the values ξ 1 = (500 100 24 28) T , ξ 2 = (300 300 28 32) T with probability 0.4 and 0.6 respectively. In Listing 2, an instance of the stochastic program (13) is created for this distribution. This code uses the model recipe created in Listing 1 to create second-stage models for each of the supplied scenarios. Here, we have used the default scenario constructor @scenario, where data values are named in accordance with the @uncertain annotation. The deterministic structure (extensive form) is used by default. Because this is a small example, correctness of the generated problem is easily verified. We can now set an optimizer and solve the model, as shown in Listing 3. The underlying memory structure can be set explicitly by setting the instantiation keyword to any of the supported structures during model instantiation. Alternatively, if an optimizer is chosen during instantiation, an appropriate structure is chosen automatically. For example, if we instantiate the same problem with an L-shaped optimizer the vertical structure is used instead, as can be seen in Listing 4. The same stochastic program has now been decomposed into a first-stage master problem and two second-stage subproblems. For completeness we also exemplify how the same problem is instantiated and solved using the progressive-hedging algorithm in Listing 5.

Infinite sample space
To demonstrate how SPjl handles continuous distributions for uncertain parameters, we assume that the uncertainties in our simple example follow a multivariate normal distribution, ξ ∼ N (µ, Σ). In general, there is no closed form solution of (2) when ξ has a continuous distribution. However, by the law of large numbers, a viable discrete approximation

✝ ✆
can be obtained by sampling scenarios from the continuous distribution. In SPjl, we achieve this by creating a sampler object associated with the defined scenario structure. In Listing 6, a sampler object for a multivariate distribution with is created and used to generate an instance of (13) with 100 sampled scenarios. Note that the same stochastic model object defined in Listing 1 is used in Listing 6 to generate the sampled instance.
With the ability to instantiate sampled models with an arbitrary number of scenarios, we can adopt the SAA methodologies developed in [22] to calculate confidence intervals around the optimal value of (13) as well as around the EVPI and the VSS. This is exemplified in Listing 7. These methods require re-solving sampled stochastic programs multiple times and the accuracy of the solution is increased by increasing the number of scenarios in the sampled models. Consequently, the parallel capabilities of SPjl become significant as these subproblems can become too large for single-core approaches. If multiple Julia processes are available, either locally or remotely, then the code in Listing 6 would automatically distribute the stochastic program on the available nodes in either a vertical or a horizontal structure. Although not practically required for this small example, this leads to significant performance gains for large-scale industrial models. See for example the scaling results presented in Section 6.

The farmer problem
To exemplify functional correctness, and allow for comparisons with similar tools, we consider the instructive farmer problem by [1]. Listing 8 shows a suggested code excerpt for how the farmer problem can be defined in SPjl and Listing 9 shows how the problem can be instantiated, solved, and analyzed using various solvers. The correctness of the numerical values can be verified in [1]. For comparison, the same problem is defined in PySP as outlined in [17] in about 100 lines spread out in separate files. Again, we stress that only 30 lines of Julia code are required to define the farmer problem in SPjl. Moreover, the optimal value, as well as the EVPI and VSS, can be calculated interactively in the same Julia session. This feature distinguished SPjl from other similar tools such as PySP. The time required to solve the farmer problem using the L-shaped algorithm was 0.57 seconds for both SPjl and PySP, measured on the same master node as the numerical benchmarks presented in the paper. Hence, there is no performance decrease from using SPjl instead of PySP for this small problem, with the added benefit of SPjl being more user-friendly.

Advanced solver configurations in the SPjl framework
The SPjl framework includes a variety of customizable improvements to the L-shaped and progressive-hedging algorithms. The possible variations of the classical algorithms included in the framework range from efficient implementations of influential research papers [28,29,30] to novel variants developed by the framework authors [31] or others [32,33,34]. We provide a summary of the improvements available for both L-shaped and progressive-hedging. In brief, each algorithm has a set of options that can be varied through a simple interface. In all examples, it is assumed that a given stochastic program instance sp has been instantiated with an appropriate optimizer. We can then use set_optimizer_attribute(sp, option, value) to customize the optimizer algorithm used by sp.

L-shaped
The L-shaped solver suite of SPjl includes a large set of customizable options. These are summarized below.
Regularization: A Regularization procedure limits the candidate search to a neighborhood of the current best iterate in the master problem. It tends to result in more effective cutting planes and improved performance of the L-shaped algorithm. Moreover, regularization enables warm-starting the L-shaped procedure with initial decisions. We have previously covered the regularization procedures in SPjl more in depth in [35].
The SPjl framework includes the following regularizations: Trust-region regularization [29], Regularized decomposition [28], Level set regularization [30]. Since the two latter techniques involve solving problems with quadratic penalty terms, the SPjl framework also provide an option for replacing quadratic penalties with various linear approximations, if only a linear solver is available.
Aggregation: Cut aggregation can reduce communication overhead and load imbalance and yield major performance improvements in distributed settings. In the classical L-shaped algorithm [12], all cuts are aggregated every iteration. The authors of [26] suggested a multi-cut variant where cuts are added separately in a disaggregate form, which on average yields faster convergence. We recently explored a novel set of aggregation approaches [31], which are all included in SPjl.
Consolidation: Cut consolidation, as proposed by [33], is also implemented in SPjl to reduce load imbalance by removing stale cuts from the master.
Execution: In a distributed environment with multiple Julia processes, the execution policy of the L-shaped algorithm can be executed in a serial, synchronous or asynchronous mode. The synchronous variant runs the L-shaped algorithm in parallel using a map-reduce pattern each iteration. The asynchronous scheme is appropriate in a heterogeneous environment where some workers may finish slower than others. We show how these algorithm policies can be applied to increase performance on large-scale problems in Section 6.

Progressive-hedging
The progressive-hedging solver suite shares a few options with the L-shaped suite. First, as each subproblem in the progressive-hedging procedure includes a quadratic penalty term, the same linear approximations as for L-shaped regularizations can be applied. Second, just like the L-shaped solvers, the progressive-hedging algorithms can be run serially, synchronously or asynchronously.
Penalization: The convergence rate of the progressive-hedging algorithm is sensitive to the choice of the penalty parameter r. The SPjl framework supports both a fixed penalty parameter and the adaptive strategy introduced in [32].

Implementation details
In this section, we provide a summary of the main software innovations in SPjl. We also discuss the implementation of the framework's distributed capabilities. The inner workings of SPjl are primarily based on two ideas: deferred model instantiation and data injection. In brief, a model definition in SPjl is a recipe for how to use data structures when building optimization models, while the actual model creation is deferred until data is provided. When a specific model is instantiated, the provided data is injected where required to construct the model. The main effect of this approach is that the stochastic model formulation is separated from the design of stochastic data parameters, which makes the SPjl framework versatile and flexible to use. For instance, it is possible to test small instances of a model locally to ensure that it is properly defined, and then run the same model in a distributed environment with a large set of scenarios. Deferred model instances and data injection also play a large role when distributing stochastic program instances in memory.

Deferred model instantiation
The advantages of deferred model instantiation is a smaller memory footprint and the ability to create various structures that use the first-and second-stage recipes as building blocks in a clever way. Examples include the deterministic, vertical, and horizontal structures, as well as wait-and-see problems and expected-value problems. The technique is also a premise for implementing data injection. In contrast to standard JuMP models, SPjl models defined through the @stage macros are not necessarily instantiated immediately. Instead, the user-defined Julia code that constructs the optimization problems is stored in lambda functions as model recipes. In other words, instead of creating and storing a JuMP object, the lines of code required to create the JuMP object is stored. This is achievable since Julia code is itself a data structure defined within the Julia language.
Deferred model instantiation is made possible through metaprogamming and the automatic reformulation bridges introduced in MathOptInterface [16]. These techniques allow us to add linking constraints between the stages that adhere to the data dependencies defined by the user. During model creation, any @decision line in a @stage definition creates special JuMP variables whose behaviour depends on the context of the instantiation. Any variable defined in this way can be included in @constraint definitions in subsequent stages. See for example Listing 1, where the last two constraint definitions in the second stage include references to x 1 and x 2 which were defined with @decision in the first stage. Next, we will discuss in more detail how instantiation is implemented for the main underlying structures: deterministic, vertical, and horizontal. In addition, we explain how decision evaluation is implemented in the different structures.

Deterministic structure
We construct the extensive form of a finite model (9) in steps using the stored model recipes. First, we generate the first-stage model in full using the corresponding recipe. Next, we process all available scenarios iteratively. For each scenario, we apply the second-stage recipe and append the resulting subproblem to the extensive model. In this context, any variables defined with @decision in the first stage are treated as regular JuMP variables. Before generating the subsequent scenario problem, we internally annotate the variables and constraints to associate them with the scenario they originated from. This labeling is visible in the printout shown in Listing 2. During decision evaluation, all variables defined with @decision are fixed to their corresponding values. The deterministic equivalent problem is then solved as usual, giving exactly (5).

Vertical structure
The vertical structure, introduced in Section 2.3.1, is also instantiated in steps. First, the first-stage master problem (11) is created using the corresponding recipe. Here, the @decision variables are again treated as regular JuMP variables. Next, subproblem instances of the form (1) are created for each possible scenario using the second-stage recipe. During second-stage generation, first-stage variables annotated with @decision enter the model as so called known decisions. These are not optimization variables, but rather parameters with given values. This design reflects the fact that the firststage decisions have already been taken when the second stage is reached. The values of the first-stage decisions can be entered into the second-stage constraints in which they appear through automatic reformulation bridges. Internally, all decisions defined in the first-stage are made known to the second stage by the @stochastic_model macro. It is also possible to explicitly add @known annotations to the second-stage definition to mark variables that originate from previous stages.
The subproblems are either stored in vector format on the master node or distributed on remote nodes as described in Section 5.3. We distribute new scenarios and generated subproblems as evenly as possible on remote nodes to achieve load balance. During decision evaluation, all variables defined with @decision are fixed to their corresponding values in the first stage. Further, these values are communicated to all subproblems, that can then update their respective second-stage constraints. The first stage and second stage problems are then solved separately, in parallel if possible, and the results are map-reduced to form (5).

Horizontal structure
Instantiation of the horizontal structure introduced in Section 2.3.2 is similar to instantiation of the vertical structure. They differ in that there is no master problem and in that the subproblems have the structure given in (12) instead of (1). The process for generating subproblems of this wait-and-see form is equivalent to one iteration of the finite extensive form generation. In short, the first-stage recipe is applied, followed by applying the second-stage recipe on the scenario data corresponding to the subproblem. Now, the variables defined with @decision are again treated as standard JuMP variables. Note that generation of the expected-value-problem (6) is equivalent to generating a wait-and-see model on the expected scenario of all available scenarios. The implementable solution ξ that enter the Listing 10: Simple showcase of data injection in SPjl ✞ ☎ @stochastic_model begin @stage 1 begin @decision(model, x) end @stage 2 begin @parameters d @known x @uncertain ξ @recourse(model, y <= d) @constraint(model, x + y <= ξ) end end ✝ ✆ horizontal form (12) through the non-anticipative constraints is added as a known decision to the subproblems. During the progressive-hedging procedure, the value of ξ can then be updated efficiently through bridges. This design is also used to implement proximal terms in the regularized variants of the L-shaped algorithm. Decision evaluation is performed similar to the other structures. In each subproblem, the first-stage decisions are fixed to their corresponding values and the subproblem is solved as usual. The results are then map-reduced to form (5). Again, the decision evaluation process is embarrassingly parallel in a distributed environment.

Data injection
Data injection is the second software pattern used to separate model and data design in SPjl. The aim is to make an object independent of how its dependencies are created. In SPjl, the dependencies consist of the data required to construct the optimization problems as described by the model recipes. The data includes uncertain parameters, as well as first-stage decisions and deterministic parameters. By adopting this approach, users of SPjl can focus on the design of the optimization model and the uncertainty model separately, while the framework is responsible for combining these designs into actual stochastic program instances. In the following, we describe the data injection functionality in more detail.
When an SPjl model is formulated using @stochastic_model, special annotations are used inside the @stage blocks to specify points of data injection. These annotations inform the framework which parameters are necessary to construct the model according to the @stochastic_model definition. The @stage macro transforms the stage blocks into anonymous lambda functions that map supplied data into optimization problems. Internally, when the user wants to instantiate the defined SPjl model, the required data is passed to the stored lambda functions according to one of the instantiation procedures outlined in the previous section.
We give a short review of the different types of data dependencies that can be specified in an SPjl model. Consider the simple second-stage formulation in Listing 10, which includes several data injection annotations.
Deterministic data: The @parameters annotation specifies scenario-independent data, i.e., deterministic parameters that are the same across all scenarios. Default parameter values can be specified inside the @parameters block. Otherwise, the values must be supplied during instantiation.
Uncertain data: The @uncertain annotation specifies the scenario-specific data. The scenario-dependent values are either created and supplied directly by the user or by some user-defined sampler object that models the uncertainty.
Decisions: The @known annotation makes the first-stage decision x available in the second-stage. Note again that @known annotations are implicitly added by @stochastic_model because of the @decision x in the first stage. When the second-stage generator is run, the framework will have already created a decision variable x using the first-stage generator, either as a standard JuMP variable or as a fixed known decision. All such first-stage variables are injected into the second-stage generator. These can then be used as if they were ordinary JuMP variables. See for example the last @constraint definition in the second stage of Listing 10.

Models:
The model keyword is a placeholder for a JuMP object that stores the actual optimization problem. In a deterministic structure, the model object is the same in every generator call. In the block-decomposition structures, the generators are instead applied to multiple JuMP models that form subproblems.
The use of data injection adds versatility to the framework. The user is only restricted to use the model keyword in the JuMP macro calls. Otherwise, all JuMP features are supported in the stage blocks. Also, there is no restriction on the scenario data types. Hence, instead of a simple structure with fields, it is possible to define a more complex data type that for example performs calculations at runtime to determine optimization parameters. In addition, any Julia methods defined on the scenario type become available in the stage blocks. This allows the user to design complex models of the uncertainty orthogonally to the definition of the stochastic program.
Because the model definition is decoupled from the data, it is possible to send the model recipe to a remote process where the scenario data resides and create the model from there. This is the foundation of the distributed implementation described next.

Distributed computations
SPjl has distributed capabilities for both modeling, analysis and optimization. All implementations rely on the Distributed module in Julia. This allows us to develop SPjl using high-level abstractions that utilize the efficient low-level communication protocols in Julia. In this way, the same codebase can be used to distribute computations locally, using shared-memory, and remotely, in a cluster or in the cloud.
Distributed computing in Julia is centered around the concepts of remote references and remote calls. Remote references are used to administer which node particular data resides on and to provide the remaining processes access to the remote data. Remote calls are used to schedule tasks on the nodes. Any process can wait on a remote reference, which blocks until data can be fetched, and then fetch the result when it is ready. The RemoteChannel objects are special remote references where processes can also put! data. Besides, specialized channel objects can be designed for specific data types. This feature is used frequently in the implementation of the distributed structured solvers.

Distributed stochastic programs
The distributed capabilities of SPjl were designed with the aim to minimize the amount of data passing. This is mainly achieved through the deferred instantiation and data injection techniques outlined above. In principle, a stochastic program can be instantiated in a distributed environment by passing all necessary data to each worker node. However, the data injection technique is independent of the way data is created. Therefore, a far more efficient approach is to let the workers generate the necessary scenario data and the optimization models themselves, with minimal data passing. This is possible since SPjl has support for passing lightweight sampler objects capable of randomly generating scenario data, such as the one in defined in Listing 6, along with passing the lightweight model recipes created in the @stage blocks. Scenario data and subproblems can then be generated in parallel on the worker nodes. The master keeps track of the scenario distribution and ensures that new scenarios and subproblems are generated on available workers in a way that promotes load-balance.
If multiple Julia processes are available, then any instantiated stochastic program in SPjl is automatically distributed in memory according to either a vertical structure or a horizontal structure. In a vertical structure, the master node administers the first-stage problem and schedules tasks and data transfers. In a horizontal structure, the master node is only responsible for task scheduling and data transfers. Aside from distributing the models in memory, SPjl parallelizes as many computations as possible. In many cases, speedups stem from subtasks being embarrassingly parallel over the independent subproblems. For example, this occurs during decision evaluation and calculation of EVPI and VSS. In these instances, the master schedules the same computation tasks on all workers using remote calls and then initiates any necessary reductions after the workers have finished using a standard map-reduce pattern. The more involved parallelization strategies in SPjl relate mostly to the structure-exploiting distributed solvers described in more detail in Appendix 5.3.2.

Distributed structured optimization algorithms
The implementations of the distributed structured solvers are also centered around remote calls and channels. Here, remote calls are used to initiate running tasks on every worker node, and the algorithm logic is driven by having the master and worker tasks wait on and write/fetch to/from specialized queue channels.
In the case of the L-shaped method, whenever the master node re-solves the master problem (11), it writes the new decision vector to a specialized Decision channel. It then sends a corresponding index to a Work channel on every remote node. Every worker continuously fetches tasks from its Work channel and uses the acquired index to fetch the latest decision vector from the master. Every new decision candidate infers a batch of subproblems to solve for each worker. After a worker has solved a subproblem (1), it sends the computed cutting planes to a CutQueue channel on the master. The master continuously fetches cuts from the CutQueue and appends them to the master problem. In the synchronous variant, the master only updates after all workers have finished their work for the current iteration. In other words, the synchronous algorithm is driven by the master node initiating and waiting for worker tasks through remote calls. In the asynchronous version, the master updates after it has received κn cuts, where n is the total number of subproblems. Timestamps are communicated throughout to keep track of the algorithm history and allow synchronized convergence checks. All subproblems are solved to completion each iteration regardless of the value of κ, to be able to check convergence properly. When the master has received all cuts corresponding to a specific iteration, it performs a convergence check and terminates if appropriate. For clarity, the procedure is illustrated in Fig. 1. A similar design is used to implement synchronous and asynchronous variants of the progressive-hedging algorithm.

Numerical benchmarks
We now evaluate the distributed performance of SPjl by benchmarking the structure-exploiting solvers on a large-scale planning problem. The numerical experiments are performed in a multi-node setup where a laptop computer acts as the master node and a desktop compute server of up to 32 cores provides worker nodes.

The SSN problem
We evaluate the solvers on the telecommunications problem SSN, first introduced in [36]. This problem is often included in similar benchmarks [29,34]. The SSN problem is formulated to plan bandwidth capacity expansion in a network before customer demands are known. The problem is freely available in the SMPS format 2 . The problem has 89 decision variables in the first stage, and 706 variables and 175 constraints in the second stage. We first run an SAA procedure to gauge the number of scenarios required to obtain a stable solution. The results are shown in Figure 2.
There is no visible improvement after 6000 scenarios. Moreover, the confidence interval around the optimal value is considered relatively tight at this point and is consistent with similar experiments [25]. With 6000 scenarios, the extensive form of the SAA model has 4.2 million variables and 5.3 million constraints, and about 20 minutes is required to build and solve the extensive form using Gurobi [10]. From this baseline, we run the distributed benchmarks.

Benchmarks
We evaluate the structured solvers by solving distributed SSN instances of 6000 scenarios. Benchmarks are performed using the Julia package BenchmarkTools.jl, which schedules multiple solve procedures and reports median computation times. Every solver runs until convergence criteria are reached with a relative tolerance of 10 −2 . The master node is a laptop computer with a 2.6 GHz Intel Core i7 processor and 16 GB of RAM. We spawn workers on a remote multi-core machine with two 3.1 GHz Intel Xeon processors (total 32 cores) and 128 GB of RAM. The two machines were 30 kilometers apart at the time of the experiments. The time required to pass a single decision or optimality cut at this distance is about 0.01 seconds. Hence, the communication latency is small, but not negligible as will be apparent in the results. For single-core experiments we only run the procedures once because the time to convergence is long and the measurement variance becomes relatively small. Throughout, the Gurobi optimizer [10] is used to solve emerging subproblems.
We first benchmark a set of L-shaped solvers. The nominal method is the multi-cut L-shaped algorithm without any advanced configuration. On average, this algorithm requires 19 iterations and 92 000 optimality cuts to solve an SSN instance of 6000 scenarios. This takes just over 30 minutes on the master node under serial execution. We run a strong scaling test where the number of worker cores on the remote machine is doubled in size up to 32 cores. Apart from multi-cut L-shaped, we also evaluate two variants with advanced algorithm policies. Specifically, one solver is configured to use trust-region regularization and partial cut aggregation with 32 cuts in each bundle. This aggregation scheme is static; the cuts are partitioned into groups of 32 in the same order each iteration. The second solver is configured to use level-set regularization and K-medoids cluster aggregation. This is a dynamic aggregation scheme where the cuts are clustered using the K-medoids algorithm based on a generalized cut distance matrix each iteration. We fix the partitioning scheme of the dynamic method after the first five iterations, as outline in [31]. All solvers are configured to use synchronous execution. The results from these experiments are shown in Figure 3.
First, we just consider the multi-cut method. The initial scaling is very poor with almost no speedup. We then observe speedups up to eight cores upon which the scaling curve flattens. The primary sources of inefficiency in distributed L-shaped algorithms are communication latency and load imbalance. This is especially true for multi-cut L-shaped because all cuts are passed separately and the master increases in size by the maximum number of constraints possible each iteration. We re-ran the two-core experiment on the master node with local threads as workers. In other words, without communication overhead. The time to convergence was then about 18 minutes. With 0.01 seconds required to pass a single cut and 92 000 cuts passed in total, this accounts for the extra 15 minutes required to converge in the multi-node setup. Therefore, we can conclude that much of the inefficiency stems from communication latency. The   Strong scaling Multi-cut L-shaped L-shaped with trust-region regularization and partial cut aggregation L-shaped with level-set regularization and K-medoids cluster aggregation  fact that the scaling curve flattens stems mostly from load imbalance. In the final iterations, most of the time is spent solving the now large master problem or passing cuts, so the worker nodes are not utilized optimally.
Next, we consider the advanced methods. The distributed performance is significantly improved compared to the multi-cut method. The main reason for this is that cut aggregation reduces both communication latency and load imbalance. Because cuts are aggregated, less data is passed each iteration. Further, the master problem does not grow as fast. Hence, the workload is more evenly spread out between master and workers, which improves parallel performance. In this particular case, the more advanced aggregation scheme yields slightly better performance, but it could also hold that level-set regularization is more performant than trust-region regularization on the SSN problem. Even with cut aggregation, the size of the master eventually exceeds the size of the subproblems and data passing still becomes a bottle-neck as the number of cores increase. Therefore, the scaling curves still flatten for larger numbers of cores. We do not claim that these configurations are the best possible. We can for example note that they are not optimal for single-core execution where both variants are outperformed by the multi-cut method. Also, the parallel efficiency increase as workers are added is not uniform. This is because the aggregation schemes are more optimal for some work granularities. We could possibly improve the convergence times further by parameter tuning. For this particular configurations, we could also let a processor on the remote machine act as the master node and remove communication latency all together. However, we believe that our results are a strong encouragement for the distributed capability of the SPjl framework. With non-negligible communication latency we are able to solve a large-scale planning problem in just over a minute by employing some of the readily available algorithm policies in the framework. This can be seen as a proof of concept for running industrial planning problems in a modern cloud architecture.
We tested the algorithms with asynchronous execution as well, but saw no performance improvements. Even though there is communication latency between the master node and the remote node, worker performance is even. Moreover, the subproblems are equally difficult to solve. There is therefore no immediate gain from introducing asynchrony and the overhead from doing so decreases performance. The asynchronous variants are expected to yield better performance in a more heterogeneous environment with stalling workers.
Next, we evaluate the performance of the progressive-hedging methods. Using the nominal method, we did not observe convergence even after long waiting times. Using the adaptive penalty policy eventually yields convergence. We configure the solvers to use adaptive penalty and synchronous execution and run the same strong scaling experiment as for the L-shaped methods. The results are shown in Figure 4.
Although at much worse time-to-solution than the L-shaped methods initially, the distributed progressive-hedging algorithm displays great scaling and outperforms the multi-cut L-shaped method after 16 cores. The efficiency probably stems from the problem being load-balanced across the workers. Communication latency again becomes a bottle-neck at 32 cores from which we attribute the worsened scaling. Again, the subproblems appear equally difficult as there were no stalling workers. Consequently, we did not observe any speedups from running the asynchronous variant.
The time to convergence is notably large and the progressive-hedging method is consistently outperformed by the advanced L-shaped methods. This is not surprising as we have spent more time on L-shaped improvements. Future work includes further algorithmic improvements to the progressive-hedging algorithms.

Conclusion
In this work, we have presented an open-source framework, StochasticPrograms.jl, for large-scale stochastic programming. It is written entirely in Julia and includes both modeling tools and solver algorithms. The framework is designed for distributed computations and naturally scales to high-performance clusters or the cloud. By using the extensive form, which is efficiently generated using metaprogramming techniques, stochastic program instances can be solved using open-source or commercial solvers. Through deferred model instantiation, data injection, and clever algorithm policies, the framework can operate in distributed architectures with minimal data passing. In addition, several analysis tools and stochastic programming constructs are included with efficient implementations, many of which can run in parallel.
The framework also includes a solver suite of scalable algorithms that exploit the structure of the stochastic programs. The structured solvers are shown to perform well on large-scale planning problems. High parallel efficiency is achieved for distributed L-shaped methods using cut aggregation techniques and regularizations. Moreover, distributed progressive-hedging algorithms are accelerated using an adaptive penalty procedure. The solver suites are made modular through a policy-based design, so that future improvements can readily be added.
There are several directions for future additions to the framework. First, SPjl does not yet fully support multi-stage problems. We have finished an infrastructure for representing multi-stage problems in a way that leverages the twostage design. Ongoing work involves designing a suitable Julian syntax for encoding transitive probabilities in a multi-stage scenario tree. Second, we will consider further algorithmic improvements to the existing L-shaped and progressive-hedging solvers. We also want to explore alternative sample-based approaches to the SAA method where the sampling is instead performed inside the structure-exploiting algorithm procedure. Examples of such approaches include L-shaped with importance sampling [37] or stochastic decomposition [38].
The framework is well-tested through continuous integration and is freely available on Github 3 . A comprehensive documentation is included 4 . The modeling framework, StochasticPrograms.jl, exists as a registered Julia package, which can be installed and run in any interactive Julia session.