Disciplined Convex Stochastic Programming:
A New Framework for Stochastic Optimization

Alnur Ali
Machine Learning Dept.
Carnegie Mellon University
alnurali@cmu.edu

J. Zico Kolter
Steven Diamond
Stephen Boyd
School of Computer Science Dept. of Computer Science Dept. of Electrical Engineering
Carnegie Mellon University
Stanford University
Stanford University
zkolter@cs.cmu.edu stevend2@stanford.edu boyd@stanford.edu

Abstract
We introduce disciplined convex stochastic programming (DCSP), a modeling framework thatcan significantly lower the barrier for modelers
to specify and solve convex stochastic optimization problems, by allowing modelers to naturally
express a wide variety of convex stochastic programs in a manner that reflects their underlying mathematical representation. DCSP allows
modelers to express expectations of arbitrary expressions, partial optimizations, and chance constraints across a wide variety of convex optimization problem families (e.g., linear, quadratic, second order cone, and semidefinite programs). We
illustrate DCSP’s expressivity through a numberof sample implementations of problems drawn
from the operations research, finance, and machine learning literatures.

1

INTRODUCTION

We introduce disciplined convex stochastic programming (DCSP), a modeling framework for specifying andsolving convex stochastic programs: convex optimization
problems that include random variables. DCSP builds onprinciples from stochastic optimization and convex analysis to allow modelers to naturally express a wide variety ofstochastic programs in a manner that reflects their underlying mathematical representation. At a high level, DCSP enables modelers to specify — in a straightforward way —and solve convex optimization problems that include (1)expectations of arbitrary expressions, (2) partial optimizations, optimizations over (only) a subset of the optimization
variables, which additionally pave the way for the specification of multi-stage stochastic programs (Sec. 2.1), and(3) chance constraints, constraints that are required to holdwith high probability — these three building blocks can beused to express a wide variety of stochastic optimization
problems.

Concurrently with this paper, we also make available anopen source Python implementation of DCSP, which werefer to as cvxstoc1 , that allows modelers to write andsolve stochastic programs — we present a variety of examples of using cvxstoc to model stochastic optimization
problems, drawn from the operations research, finance, andmachine learning literatures, in Sec. 4.
Related work Although other frameworks for stochastic
programming do exist ([24, 20, 11], and in Python mainly[26]), they often require significant effort from the modeler
to manipulate the optimization problem into an amenable
form, support a limited number of stochastic programming
constructs (e.g., [11] only supports chance constraints withuncertainty sets), and cannot express certain families of
convex optimization problems; indeed, checking the convexity of and solving (convex) optimization problems ingeneral is challenging. DCSP builds on (and extends) disciplined convex programming (DCP) [10], a recently introduced framework that makes it natural for modelers to express convex optimization problems, and additionally automates the tasks of verifying the convexity of these problems
and translating them into conic form (see, e.g., [8]). Thismeans that DCSP can be used to express and solve a wide
variety of stochastic convex optimization problems, including linear, quadratic, second order cone, and semidefinite programs. Probabilistic programming languages (e.g.,[9, 16, 13]) offer an alternative approach, but tend to focus on inference problems, and may not contain the features to capture traditional stochastic programming problem formulations; in contrast, convex modeling can be attractive because local solutions are global solutions, efficient solvers exist, and guarantees can often be obtained on
the optimality of a solution obtained by a solver.
This paper is structured as follows. In Sec. 2, we review background on stochastic programming, DCP, andcvxpy (an open source Python implementation of DCP).
In Sec. 3, we describe DCSP and cvxstoc’s syntax. In
1
cvxstoc is available as an extension of the cvxpy Python
package [7]: see http://www.cvxpy.org.

Sec. 4, we present a number of examples that illustrate ourframework.

Prob (f (x, ω) ≤ 0) ≥ η,

2

BACKGROUND

2.1

STOCHASTIC PROGRAMMING

where f is convex in x for each value of ω, and η is typically a large probability (e.g., 0.95); a chance-constrained
problem is an optimization problem with one or morechance constraints. Chance constraints are typically nonconvex, although effective convex approximations exist(see Sec. 3.3).

A convex optimization problem has the form
minimize f0 (x)subject to fi (x) ≤ 0,
hi (x) = 0,i = 1, . . . , m,i = 1, . . . , p,

2.2

where x ∈ Rn is the optimization variable, f0 : Rn → R isa convex objective function, fi : Rn → R, i = 1, . . . , m
are convex inequality constraint functions, and hi : Rn →R, i = 1, . . . , p are affine equality constraint functions.
A convex stochastic program has the form
minimizex

E f0 (x, ω)

subject to E fi (x, ω) ≤ 0, i = 1, . . . , mhi (x) = 0, i = 1, . . . , p,
n

(1)

q

where fi : R × R → R, i = 0, . . . , m are convex functions in x for each value of a random variable ω ∈ Rq , andhi : Rn → R, i = 1, . . . , p are (deterministic) affine functions; since expectations preserve convexity, the objective
and inequality constraint functions in (1) are (also) convex
in x, making (1) a convex optimization problem.
Two-stage stochastic programs An important special
case of (1) is a so-called two-stage stochastic program (alsoreferred to as an optimization problem with recourse) [6]:minimizex

i = 1, . . . , m,i = 1, . . . , p,

DISCIPLINED CONVEX PROGRAMMING

Disciplined convex programming (DCP) is a recently introduced modeling framework for specifying and solving convex optimization problems [10]. In a nutshell, DCP consists of a library of convex atomic functions, and a convex
rule-set that prescribes how these atomic functions may becomposed to express (more complex) convex optimization
problems.
Convex rule-set Verifying the convexity of arbitrary expressions is challenging; DCP checks convexity usingThm. 2.1, which is equivalent to enforcing a set of rules.
Theorem 2.1 ([10]). Suppose f = h(g1 (x), . . . , gk (x)),
where h : Rk → R is convex and gi : Rn → R, i =1, . . . , k, and one of the following holds for each i =1, . . . , k:
• gi is convex and h is nondecreasing in argument i• gi is concave and h is nonincreasing in argument i• gi is affine.

f0 (x) + E Q(x, ω)

subject to fi (x) ≤ 0,
hi (x) = 0,

Chance-constrained problems A chance constraint [5]is a constraint on the variable x of the form(2)

where Q(x, ω) = inf y {φ0 (x, y, ω) : φi (x, y, ω) ≤ 0,
ψj (x, y) = 0, i = 1, . . . , s, j = 1, . . . , w}
is the second stage problem, y ∈ Rr is the second stage
optimization variable, φ0 : Rn × Rr × Rq → R is thesecond stage objective function, and is convex in (x, y) for
each value of ω, φi : Rn × Rr × Rq → R, i = 1, . . . , s
are the second stage inequality constraint functions, also
convex in (x, y) for each value of ω, and ψi : Rn × Rr →R, i = 1, . . . , w are the second stage equality constraintfunctions, and are affine in (x, y). That is, Q is itself the
optimal value of another convex optimization problem, andis convex in x for each value of ω.
Two-stage stochastic programs model the uncertain consequences (in the second stage) of here-and-now decisionmaking (in the first stage): e.g., in a finance application, we
may wish to decide which assets to purchase now, while(also) factoring in how the asset prices might fluctate later.

Then f is convex2 .
Thm. 2.1 permits a wide variety of convex expressions: forexample, the maximum eigenvalue of a symmetric matrix,
λmax (2X − 4I), where X ∈ Sn , is recognized as convex.
(On the other hand, as an example of a limitation of the

Pn
2 1/2
rule-set, the expression, where x ∈ Rn , isi=1 xinot recognized as convex, although it is recognized as convex once reformulated as kxk2 .)Library of atoms Atomic convex functions3 are specified in DCP in their epigraph form: for example, the (convex) function f (x) = kxk1 is specified asminimizet

1T t

subject to −t  x  t,2

(3)

A similar result holds for concave functions.
See
http://www.cvxpy.org/en/latest/
tutorial/functions/index.html for a list of theconvex atoms available in cvxpy.
3

where x, t ∈ Rn . Thus, whenever a modeler writes theatom f (x) = kxk1 , DCP internally replaces it with (3),introduces the variable t, and can subsequently optimize
over (x, t).

its implementation in the cvxstoc Python package. In
a nutshell, DCSP consists of the addition of three operations to the disciplined convex programming (DCP) framework, which can be used to express a wide variety of convex stochastic programs: the ability to (1) compute (approximations to) expectations of arbitrary expressions, (2)handle partial optimization, and (3) compute (approximations to) chance constraints.

Disciplined convex programming DCP certifies a problem’s convexity by constructing an abstract syntax tree forthe objective and constraint functions, with atoms as internal nodes, and variables and constants as leaves, and thenapplying Thm. 2.1 recursively [10, 22].

3.1

RANDOM VARIABLES AND EXPECTATIONS

Random variables The most fundamental operations instochastic programs, and hence in DCSP, are the ability to
specify random variables, and compute (approximations to)
cvxpy [7] is an open source Python DCP implementation;
expectations of arbitrary expressions containing these ranwe briefly describe its syntax next.
dom variables. As in Sec. 2, DCSP assumes that all expresVariables are declared simply in cvxpy as follows: sions in a stochastic program are convex in the optimization
x = Variable()variable(s) for each value of the random variable(s) — thus,
x = NonNegative()from the point of view of DCSP, random variables do not
X = Semidefinite(n) affect the convexity of their parent expressions and can be
regarded as equivalent to constants, thereby requiring no
The first line declares x to be a variable in R, the second
additions to the DCP convex rule-set. (Practically speakdeclares x to be a variable in the nonnegative orthant R+ ,ing, DCSP permits the specification of a variety of random
and the third declares X to be a (n × n matrix) variable invariables; see Sec. 3.4.)nthe positive semidefinite cone S+ .
2.3

cvxpy

Convex expressions are specified by composing
conExpectations DCSP computes (approximations to) exPm
vex atoms; for example, the log loss
pectations of arbitrary expressions using simple Monte
i=1 log(1 +exp(−yi (wT xi + b))) can be specified by using the simCarlo evaluation, i.e.,pler log sum exp atom as follows:NXf (x, ωi ),
E
f
(x,
ω)
≈
(1/N
)
expr = [log_sum_exp(vstack(0, -y[i]*(w.T*x[i]+b)))


for i in range(m)]

i=1where f is (again) assumed to be convex in the optimizaAn objective is specified by instantiating a sense (i.e.,tion variable x for each value of the random variable ω,
Minimize or Maximize) with an expression: and ωi , i = 1, . . . , N are samples of ω; this approximationis referred to as the sample average approximation (SAA)
obj = x.T*c
Minimize(obj) in the stochastic programming literature, and methods thatuse it are often referred to as scenario-based methods. By
Constraints are specified by forming a list of expressions:the DCP rule-set, the nonnegative weighted sum of convex

functions is a convex function; thus, the expectation operconstrs = [x >= 0, x.T*numpy.ones((n,1)) == 1, ...]
 ator applied to an expression that is convex in x returns an
expression that is (also) convex in x.
A convex optimization problem, then, is specified by instantiating a Problem with an objective and a list of conThe SAA is, of course, a very simple method for approxstraints:
imating an expectation, and much more involved methods
for solving stochastic programs exist, but the clear advanprob = Problem(Minimize(obj), constrs)
prob.solve()
tage of this method is its simplicity: any random variable
can be included in a stochastic program as long as we areable to draw samples of it. If ω is a discrete random variThe last line solves the optimization problem.
able, then DCSP calculates its expectation exactly; otherwise, DCSP draws samples using Markov chain Monte
3 DISCIPLINED CONVEX STOCHASTIC
Carlo (MCMC) methods4 [15].

PROGRAMMING

In this section we present the chief methodological contribution of the paper: the disciplined convex stochastic programming (DCSP) framework, along with an overview ofIn the case of unconstrained stochastic programs, the SAA
objective value is (naturally) an unbiased estimator of the4
We implement MCMC by leveraging the PyMC Python package [15].

true objective value, E f0 (x, ω), with variance ∝ 1/N ,and an asymptotically normal distribution [21, chap. 5].
Thm. 3.1 additionally tells us that (roughly) both the optimal value and optimal set of a SAA converge almost surely
to the optimal value and optimal set of the true problem.

with a convex upper bound derived as follows [3]. Suppose
φ : R → R+ is a nonnegative, increasing convex function
with φ(0) = 1; then φ(z) ≥ 1(z ≥ 0), where 1(z ≥ 0)
equals 1 if z ≥ 0 and 0 otherwise, and so φ(z/α) ≥ 1(z ≥0), for some variable α ∈ R++ . Thus

Theorem 3.1 ([21, Thm. 5.3]). Define p̂∗N , ŜN and p∗ , S as
the optimal value and optimal set of a SAA with N samples
and of the true problem, respectively, and let K ⊂ Rn be acompact set. Suppose (a) S ⊆ K is nonempty, (b) ŜN ⊆ K
is nonempty(c) f0 is finite and continuous on K, andPa.s.,a.s.
N(d) (1/N ) i=1 f0 (x, ωi ) −−→ f (x) (uniformly) for x ∈
a.s.
a.s.
K. Then p̂∗N −−→ p∗ and supx∈ŜN inf y∈S kx−yk2 −−→ 0.

E φ(f (x, ω)/α) ≥ Prob (f (x, ω) ≥ 0) ,p̂∗N is also a downward biased estimator of p∗ , although
its bias decreases with N [21, Prop. 5.6]. In Sec. 3.4, we
empirically investigate the quality of the SAA.
3.2

PARTIAL OPTIMIZATIONDCSP adds a new partial optimization atom to the DCP
atom library, allowing modelers to express partial optimizations, i.e., optimizations over (only) a subset of the optimization variables; this atom also forms the basis for specifying two-stage stochastic programs.
We start with the observation that partial optimization is aconvex operation (see, e.g., [4, page 87]): i.e., if f is convex
in (x, y) and C is a nonempty convex set, theng(x) := inf {f (x, y) : (x, y) ∈ C} ,yis convex in x.
Accordingly, DCSP specifies a new partial optimization
atom that takes as input a convex optimization problem andreturns (the epigraph form for) another convex atom, which
complies with the DCP prescription for specifying atoms— this means that modelers can use partial optimizations
in stochastic programs as they would other atoms. In particular, two-stage stochastic programs, i.e., (2), can be naturally expressed using this atom; furthermore, the second
stage optimization problem Q need not be in standard form
(as required by other frameworks).

and soα E φ(f (x, ω)/α) ≤ α(1 − η)(5)=⇒ Prob (f (x, ω) ≥ 0) ≤ 1 − η,i.e., (5) is a conservative approximation to (4). Note that (5)is convex in (x, α): it is the perspective of the expectationof a convex increasing function, φ, of a convex function,
f .6
In (5), α can be interpreted as modulating the “steepness”
of the approximation; several choices of φ are possible, andare analogous, e.g., to different approximations to the zeroone loss common in machine learning. DCSP uses φ(z) =max{0, z + 1}, which roughly corresponds to a Markovinequality type bound7 on (4), and can also be interpreted
as the conditional value-at-risk of f (x, ω) [19]. Prop. 3.2
also tells us that this is the tightest possible choice of φ.
Practically speaking, DCSP approximates (5) with its SAA(at which point all the benefits of DCP readily apply), then
optimizes over (x, α) to obtain the tightest possible bound;in Sec. 4.3, we empirically investigate the quality of theseapproximations.
Proposition 3.2 ([14]). Suppose φ : R → R+ is a nonnegative, increasing convex function and φ(z) ≥ 1(z ≥0); then ∃α ∈ R+ such that E (f (x, ω)/α + 1)+ ≤E φ(f (x, ω)).
3.4

cvxstoc

Next, we briefly detail the syntax of cvxstoc;
cvxstoc builds on cvxpy, and thus much of the usage
is similar.
Random variables are specified simply in cvxstoc as follows:omega = RandomVariableFactory().create_normal_rv(0,1)

3.3CHANCE CONSTRAINTSDCSP computes conservative approximations to chanceconstraints, as they are typically nonconvex5 ; in particular,
DCSP replaces
Prob (f (x, ω) ≥ 0) ≤ 1 − η,5(4)

One notable exception is chance constraints involving affine
functions of normal random variables, which can be expressed asa second order cone constraint (see, e.g., [4, page 157]). However, we favor the approximate approach described in this section
because it is substantially more general, and applies to any classof random variables.

Here, omega is a standard normal random variable.
cvxstoc includes a RandomVariableFactory object to simplify the specification of common random variables; see Sec. 4 for examples of the specification of otherrandom variables.
Expectations are specified by applying the expectation
atom:
6

Alternatively, the modeler can fix α, in which case the boundis convex in (x, η) if desired.
7
Alternatively, one could take φ(z) = exp z, which would be
analogous to a Chernoff-type inequality.







result = expectation(exp(omega*x), m)

Here, exp(omega*x) is the expression we wish to compute the expectation of, and m is the number of Monte Carlosamples to use when constructing the SAA.
Partial optimizations are specified by applying thepartial optimize atom:
atom = partial_optimize(prob, [y], [x])The first argument here is a Problem, the second is a listof variables to optimize over, and the third is a list of variables to not optimize over.
Two-stage stochastic programs can, in turn, be specified asfollows:
Figure 1: The optimal values of the stochastic least squares problem (6) (red) and a SAA to (6) (blue) with 95% confidence inter vals (light blue) vs. the number of Monte Carlo samples; in thiscase, m = 100, n = 50, although similar results hold across avariety of problem sizes.

Q = partial_optimize(prob2, [y], [x])prob1 = Problem(Minimize(f0 + expectation(Q(x),m)),constrs)Here, prob2 is the second stage problem, y is the second
stage variable, x is the first stage variable, and prob1 isthe first stage problem. Multi-stage stochastic programs
4 EXAMPLES
can be specified by iterating this construction:
Q2 = partial_optimize(prob3, [z], [x,y])At this point, we switch gears slightly and present sevprob2 = Problem(Minimize(phi0 + expectation(Q2(y),m)),,→ constrs2)
eral examples of stochastic programs along with their
Q1 = partial_optimize(prob2, [y], [x])corresponding cvxstoc implementations8 ; the majority
prob1 = Problem(Minimize(f0 + expectation(Q1(x),m)),constrs1)of these applications are well established or previously
known, though we also include some formulations that arenovel, to the best of our knowledge (namely, the precise
Chance constraints are specified by instantiating the prob
formulation of the stochastic optimal power flow problem
class and chaining it with an inequality: in Sec. 4.4, and the budgeted learning of a classifier in aprob(constr >= 0, m) <= 1-eta
 cascade problem in Sec. 4.6).
Here, constr is a (stochastic) constraint, and m is thenumber of Monte Carlo samples to use.
3.4.1

Quality of the sample average approximation

Here, we investigate the quality of the SAA employed by
cvxstoc in the special case of a (unconstrained) leastsquares problem
minimize E kAx − bk22 ,x(6)where the entries of A ∈ Rm×n ∼ Normal(µ1 , σ12 ), b ∈Rm ∼ Normal(µ2 , σ22 ), and x ∈ Rn , in which case theobjective has the analytic form
xT E AT Ax − 2 E bT Ax + E bT b,

4.1

YIELD-CONSTRAINED COSTMINIMIZATION

We begin with a simple example from the operations research literature (see, e.g., [4, page 107]). Consider the(general) problem of choosing the parameters x ∈ Rn
governing a manufacturing process so that our cost cT x,
where c ∈ Rn , is minimized, while the parameters lie
in a set of allowable values S; we can model noise in
the manufacturing process by expressing this constraint as
Prob(x + ω ∈ S) ≥ η, where ω ∈ Rn is a random vectorand η is a large probability (e.g., 0.95), which is referred
to as an η-yield constraint. Thus, we have the optimization
problemminimizex

assuming we know the second moments of (A, b). Fig. 1
plots the optimal value of the true problem (6) and a SAA
to (6): we see that the SAA obtains reasonable accuracy
after roughly 100 Monte Carlo samples.

cT xsubject to Prob(x + ω ∈ S) ≥ η.
8
Due to space constraints, we present one of these examples
in the supplementary material.

Note that if the distribution over ω is log-concave and S
d = RandomVariableFactory().create_categorical_rv(,→ d_vals, d_probs)
is a convex set, then this constraint is convex in x. We# Create optimization variables
can directly express the yield-constrained cost minimizax = NonNegative()tion problem using cvxstoc; an implementation is given
y1, y2 = NonNegative(), NonNegative()in Listing 1 (S is taken to be an ellipsoid).
 # Create second stage problem
obj = -s*y1 - r*y2
constrs = [y1+y2<=x, y1<=d]p2 = Problem(Minimize(obj), constrs)
Q = partial_optimize(p2, [y1, y2], [x])# Create problem datan = 10c = numpy.random.randn(n)P, q, r = numpy.eye(n), numpy.random.randn(n), numpy.
,→ random.randn()
mu, Sigma = numpy.zeros(n), 0.1*numpy.eye(n)omega = RandomVariableFactory().create_normal_rv(mu,,→ Sigma)
m, eta = 100, 0.95
# Create and solve optimization problem
x = Variable(n)yield_constr = prob(quad_form(x+omega,P)+ (x+omega).T*q + r >= 0, m) <= 1-eta
p = Problem(Minimize(x.T*c), [yield_constr])p.solve()

Listing 1: A cvxstoc implementation of the yield-constrained
cost minimization problem.

4.2

THE NEWS VENDOR PROBLEM

The news vendor problem is a classic problem in thestochastic programming literature (see, e.g., [2, page 15]);in this problem, a vendor must decide how much newspaper to stock, so that profit is maximized while backorderand return fees (due to excess or insufficient demand, respectively) are minimized, in the face of uncertain demand.
Our optimization variables are the number of units ofstocked newspaper x ∈ R+ , the number of units purchasedby customers y1 ∈ R+ , and the number of unpurchased(surplus) units that must be returned by the vendor y2 ∈R+ . Our problem data are b, s, r ∈ R+ , which denote theprice to stock, sell, and return a unit of newspaper, respectively. Lastly, we let the random variable d ∼ Categoricalmodel the uncertain (newspaper) demand.
We can pose the news vendor problem as the following
two-stage stochastic program:bx + E Q(x)minimizex

subject to 0 ≤ x ≤ u,where Q(x) =min−(sy1 + ry2 )s.t.

y1 + y2 ≤ x
0 ≤ y1 ≤ dy2 ≥ 0.

y1 ,y2# Create and solve first stage problem
p1 = Problem(Minimize(b*x + expectation(Q(x), want_de=,→ True)), [x<=u])p1.solve()Listing 2: A cvxstoc implementation of the news vendor
problem.



We can also represent a stochastic program by means of an
influence diagram, a directed acyclic graph, where circular
nodes correspond to random variables, square nodes correspond to decision variables, diamond nodes correspond tocosts, and edges flow from node x to node y iff the valueof node y depends in some way on the value of node x;
Fig. 2 presents the influence diagram for the news vendor
problem.
x

d

y1

y2

−s

−r

b

Figure 2: The influence diagram for the news vendor problem.

4.3

PORTFOLIO OPTIMIZATIONIn portfolio optimization, we wish to maximize wealthwhile meeting certain restrictions on risk, in the face of
uncertain asset prices; we can pose a standard portfolio optimization problem [12], subject to two kinds of risk constraints, as a stochastic program.

The risk constraints we consider here are the value-at-risk(VaR) (see, e.g., [25, chap. 29]) and conditional value-atrisk (CVaR) (e.g., [19], [23, page 286]); intuitively, VaR
A cvxstoc implementation of the news vendor problem
allows the modeler to control the probability of a loss (onis given in Listing 2; in contrast, a PySP [26] implemenasset sales) beyond a (modeler-defined) threshold, and istation (see the supplementary material) required 111 lines
often nonconvex, while CVaR allows the modeler to control
spanning 6 files.
 the expected value of such a loss, and is convex.
# Create problem datab, s, r, u = 10, 25, 5, 150
d_probs = [0.3, 0.6, 0.1]
d_vals = [55, 139, 141]

Our optimization variables are the allocation vector (acrossa set of n assets) x ∈ Rn , and the CVaR β ∈ R; the prob-lem data is the loss threshold u ∈ R+ , and the vector of
returns p ∼ Normal(p̄, Σ).
We can pose a CVaR-constrained portfolio optimization
problem as

SAA also obtains reasonable accuracy after roughly 100
Monte Carlo samples. Fig. 4 plots the probability of a SAA
to (7) vs. the number of Monte Carlo samples, and has asimilar interpretation.

E −pT x

minimizex,βsubject to β + 1/(1 − η) E(−pT x − β)+ ≤ u
1T x = 1, x  0,

(7)

where (z)+ := max{0, z}.
A cvxstoc implementation of the CVaR-constrained
portfolio optimization problem is given in Listing 3.
# Create problem datan = 10pbar, Sigma = numpy.random.randn(n), numpy.eye(n)p = RandomVariableFactory().create_normal_rv(pbar,,→ Sigma)u, eta, m = numpy.random.rand(), 0.95, 100# Create optimization variables
x, beta = NonNegative(n), Variable()# Create and solve optimization problem
cvar = expectation(pos(-x.T*p - beta), m)
cvar = beta + 1/(1-eta)*cvar
prob = Problem(Minimize(expectation(-x.T*p,m)),[x.T*numpy.ones((n,1)) == 1, cvar<=u])prob.solve()

Figure 3: The optimal values (higher means more wealth) of the
VaR-constrained portfolio optimization problem (8) (green), theCVaR-constrained portfolio optimization problem (7) (red), and a SAA to (7) (blue) with 95% confidence intervals (light blue)
Listing 3: A cvxstoc implementation of the CVaR-constrainedvs. the number of Monte Carlo samples; the problem size n = 50,
portfolio optimization problem.
although similar results hold across a variety of problem sizes.

We can also pose a VaR-constrained portfolio optimization
problem as
minimizex

E −pT xT

subject to Prob(p x ≤ 0) ≤ 1 − η
1T x = 1, x  0.

(8)

As per Sec. 3.3, DCSP replaces the chance constraint in(8) with a sample average approximation (SAA) to a (moreconservative) CVaR constraint (making (8) equivalent to(7)). We investigate the quality of this approximation in thespecial case where p ∼ Normal(p̄, Σ), in which case boththe VaR and CVaR constraints have analytic forms [18]: theVaR constraint can be expressed asp̄T x ≥ Φ−1 (η)kΣ1/2 xk2 ,(9)where Φ−1 (·) is the inverse standard normal cumulative
distribution function, while the CVaR constraint can be expressed as √
2π(1 − η)kΣ1/2 xk2 .
p̄T x ≥ exp −(Φ−1 (η))2 /2 /(10)
Fig. 3 plots the optimal value (i.e., wealth) of the VaRconstrained portfolio optimization problem (8), the CVaRconstrained portfolio optimization problem (7), and a SAA
to (7): we see that the wealth obtained by constraining VaR
is indeed less conservative than by constraining CVaR. The

4.4

OPTIMAL POWER FLOW

Consider a network G = (V, E), with a set of vertices V
and a set of edges E, that models an electrical grid: i.e., asubset of the vertices G ⊆ V are generators, which produce
power, the remaining vertices L = V \ G are loads, which
consume power, and an edge is drawn between a generator
and a load if and only if there is a (physical) transmission
line between them.
In the standard optimal power flow problem, we wish tominimize the total cost of generating power, while satisfying demand, subject to the topology of the network andper-generator capacity constraints. We often do not havecomplete control over all the generators in the grid, so wedenote the subset of generators that we do have control overas G1 , and also define G2 = G \ G1 ; we also define G = |G|,
G1 = |G1 |, G2 = |G2 |, and L = |L|. We write the pergenerator costs as cG1 ∈ RG1 and cG2 ∈ RG2 , and the
per-generator lower and upper (respectively) limits as l andu ∈ RG .
The topology/demand constraints can be expressed as
Aplin = (pG1 , pG2 , pL ), where A ∈ Rn×E is the incidence
matrix for the (directed) graph G, pG1 ∈ RG1 and pG2 ∈RG2 are variables denoting (nonnegative) power generation, pL ∈ RL are constants denoting the (non-positive)
power consumption at the loads, and plin ∈ RE are vari-# Create optimization variables
p_g1, p_g2 = NonNegative(), NonNegative()z = NonNegative(num_winds)
p_lines = Variable(E)p_w = RandomVariable(pymc.Lognormal(name="p_w", mu=1,
tau=1, size=num_winds))# Create second stage problem
p_g = vstack(p_g1, p_g2)
p = vstack(p_g1,
p_g2,
p[load_idxes[:-1]],
p_w-z,
p[load_idxes[-1]])
p2 = Problem(Minimize(p_g.T*c_g + z.T*c_w),[A*p_lines == p, p_g<=u_gens, z<=p_w,
abs(p_lines)<=u_lines])
Q = partial_optimize(p2, [p_g2, z, p_lines], [p_g1])# Create and solve first stage problem
p1 = Problem(Minimize(expectation(Q(p_g1), m)))p1.solve()

Figure 4: The probability of a SAA to (7) (blue) with 95% confidence intervals (light blue) and β (red) vs. the number of Monte
Carlo samples; the problem size n = 50, although similar results
hold across a variety of problem sizes.

Listing 4: A cvxstoc implementation of the optimal power
flow problem.

ables denoting the power flowing through each edge.
Now, additionally consider the presence of a set of renewable generators (e.g., wind farms), which we denote W,
whose (intermittent) generation an operator can either sellon the spot market [17], or use to power loads. We let
W = |W|, and model this situation with a random vectorpW ∈ RW , (pW )i ∼ LogNormal(µi , σi2 ), i = 1, . . . , W .
We can cast this as the following optimization problem:
minimize E Q(pG1 )
pG1

where

Q(pG1 ) =min

pG2 ,z,plin

s.t.

T 


pG1
+ cTW z
 pG2
pG1
 pG2 

Aplin = pL
pW − z
0  z  pW
|plin |  ulin pG1
lG  uG ,
pG2
cG1
cG2

cW is the (nonpositive) revenue obtained by selling renewable power, z ∈ RW is the decision vector for the renewable generators, and ulin are the limits on the power flowing
through each edge.
A cvxstoc implementation of the stochastic optimal
power flow problem is given in Listing 4. We solved thisproblem on the IEEE 14 Bus Test Case, i.e., with n = 14,
G1 = 1, G2 = 1, W = 1, and L = 10: Fig. 5 presents theresults.

Figure 5: The electrical grid and (optimal) power generation forthe optimal power flow problem on the IEEE 14 Bus Test Case.
Red vertices are generators: a positive number indicates the optimal power generation, while “sec. stg.” denotes an uncontrolled
generator. The blue vertex is a (stochastic) renewable generator:
its mean available (wind) power is shown above it. Other vertices
are loads: their (nonpositive) demanded powers are shown above
them.

4.5

ROBUST SUPPORT VECTOR MACHINE

Consider the problem of learning a support vector machine(SVM) from a set of m data points {(xi , yi )}m
i=1 . Suppose we would like to (additionally) model the fact thatour data collection process is noisy (in order to gain robustness in our solution), by incorporating the belief that(say) xi ∼ Normal(µ1 , Σ1 ) for all i where yi = 1 andxi ∼ Normal(µ2 , Σ2 ) for all i where yi = −1 into
the learning process. We can thereby pose the following
chance-constrained variant of the canonical (soft-margin)SVM optimization problem [1]Pm
minimize kwk22 + C i=1 ξiConcretely, we can write this optimization problem as
minimize L1 (a, b; {xi , yi }mi=1 ) + E Q(a, b),a,bw,b,ξi
subject to Prob yi (wT xi + b) ≥ 1 − ξi ≥ η,
ξi ≥ 0, i = 1, . . . , m,where Q(a, b) =min L2 (c, d; {zi , wi }pi=1 )
c,d

s.t.

n

where w ∈ R , b ∈ R, ξi ∈ R+ for i = 1, . . . , m, C
is the regularization trade-off parameter, and η is a large
probability (e.g., 0.95)9 .

kak1 + kck1 ≤ u,{(xi , yi )}m
i=1
n

is the (fixed) training set of m points in
R for the first stage classifier, and {(zi , wi )}pi=1 is the(stochastic) training set of p points in Rq for the second
stage classifier.

A cvxstoc implementation of the robust SVM problem
is given in Listing 5.
 A cvxstoc implementation, where L and L are (both)
1
2
w, b, xi = Variable(n), Variable(), NonNegative(m)taken to be the `2 -regularized log loss, is given in Listing
constr = []6.
Sigma = 0.1*numpy.eye(n)for i in range(m):mu = numpy.array(X[i])[0]x = RandomVariableFactory().create_normal_rv(mu,,→ Sigma)chance = prob(-y[i]*(w.T*x+b) >= (xi[i]-1), ns)
constr += [chance <= eta]

p = Problem(Minimize(norm(w,2) + C*sum_entries(xi)),
constr)p.solve()

Listing 5: A cvxstoc implementation of the robust SVM
problem.

4.6

BUDGETED LEARNING OF A CLASSIFIERIN A CASCADE

Suppose we are interested in learning a (single) classifierthat is part of a system (cascade) of classifiers; i.e., we areinterested in estimating the parameters a ∈ Rn and b ∈ R
of a first stage classifier, whose output is to be (somehow)
combined with the output of a second stage classifier, before presenting the combined output to a user10 .
If we knew the second stage classifier’s parameters, thenour learning task would be trivial. Instead, we choose tomodel our uncertainty as follows: we assume that we doknow the second stage classifier’s loss function, but remain
uncertain of its feature representation. We can pose this asa two-stage stochastic program, where the expectation in
the second stage is taken over all possible feature representations for the second stage classifier; for instance, if thecascade is being used for document classification, then wemight posit that each possible feature representation in thesecond stage is a function of a sample of a word from agenerative model (e.g., latent Dirichlet allocation).
We additionally assume that there is some overall test timebudget on the cascade, which we express as an upper boundu ∈ R+ on the quantity kak1 + kck1 , where c ∈ Rq are theparameters of the second stage classifier [27].
9

We note that this formulation is quite fine-grained, in thesense that per-data point noise models/distributions, as well as
mistake probabilities, may be specified.
10
Such scenarios are common in web search: see, e.g., [27].

# Create optimization variables
a, b = Variable(n), Variable()c, d = Variable(q), Variable()# Create second stage problem
obj2 = [log_sum_exp(vstack(0, -w[i]*(c.T*z[i]+d)))
for i in range(p)]budget = norm1(a) + norm1(c)p2 = Problem(Minimize(sum(obj2) + C*norm(c,2)),[budget<=u])Q = partial_optimize(p2, [c,d], [a,b])# Create and solve first stage problem
obj1 = [log_sum_exp(vstack(0, -y[i]*(x[i]*a+b)))
for i in range(m)]p1 = Problem(Minimize(sum(obj1) + C*norm(a,2) +expectation(Q(a,b), ns)), [])p1.solve()Listing 6: A cvxstoc implementation of the budgeted learningof a classifier in a cascade problem.

5

CONCLUSION

We described disciplined convex stochastic programming (DCSP), a modeling framework that can significantly
lower the barrier for modelers to specify and solve convex stochastic programs. We presented a number of sample implementations of stochastic programs that illustrated
DCSP’s expressivity; in constrast, other frameworks often
require significantly more effort from the modeler to express the problem and/or manipulate it into standard form,
support a limited number of stochastic programming constructs, and cannot express certain families of convex optimization problems.
Acknowledgements
We thank John Duchi and the reviewers for helpful discussions. This work was supported by a Dept. of Energy
Computational Science Graduate Fellowship under grant
number DE-FG02-97ER25308, and by the National Science Foundation under grant number IIS-1320402.

References[1] A. Ben-Tal, S. Bhadra, C. Bhattacharyya, and J. Nath.
Chance-constrained uncertain classification via robust optimization.
Mathematical Programming,
127(1):145–173, 2011.
[2] J. Birge and F. Louveaux. Introduction to Stochastic Programming. Springer Series in Operations Research and Financial Engineering. Springer, 1997.
[3] S. Boyd.
Chance-constrained optimization.
http://stanford.edu/class/ee364a/
lectures/chance_constr.pdf,
January
2015.
[4] S. Boyd and L. Vandenberghe. Convex Optimization.
Cambridge University Press, 2004.
[5] A. Charnes and W. Cooper. Chance-constrained programming. Management Science, 6(1):73–79, 1959.

[15] A. Patil, D. Huard, and C. Fonnesbeck. PyMC:
Bayesian stochastic modelling in Python. Journal of
Statistical Software, 35(4):1–81, 7 2010.
[16] A. Pfeffer. Figaro: An object-oriented probabilistic programming language. Charles River Analytics
Technical Report, page 137, 2009.
[17] D. Phan and S. Ghosh. Two-stage stochastic optimization for optimal power flow under renewablegeneration uncertainty. ACM Transactions on Modeling and Computer Simulation, 24(1):2:1–2:22, January 2014.
[18] A. Prékopa and R. Wets. Stochastic Programming,volume 27. North-Holland, 1986.
[19] R. Rockafellar and S. Uryasev. Conditional value-atrisk for general loss distributions. Journal of Banking
and Finance, pages 1443–1471, 2002.
[20] Richard E Rosenthal. GAMS — A user’s guide. 2004.

[6] G. Dantzig. Linear programming under uncertainty.
Management Science, 50(12 Supplement):1764–
1769, December 2004.
[7] S. Diamond, E. Chu, and S. Boyd. CVXPY: A
Python-embedded modeling language for convex optimization, version 0.2. http://www.cvxpy.
org, May 2014.
[8] A. Domahidi, E. Chu, and S. Boyd. ECOS: An SOCP
solver for embedded systems. In Proceedings of theEuropean Control Conference, 2013.
[9] N. Goodman, V. Mansinghka, D. Roy, K. Bonawitz,
and J. Tenenbaum. Church: A language for generative models with non-parametric memoization andapproximate inference. In Proceedings of the 24th
Conference on Uncertainty in Artificial Intelligence,
2008.
[10] M. Grant. Disciplined Convex Programming. PhD
thesis, Stanford University, 2004.
[11] J. Löfberg. YALMIP: A toolbox for modelingand optimization in MATLAB. In Proceedings of
CCA/ISIC/CACSD, September 2004.

[21] A. Shapiro, D. Dentcheva, and A. Ruszczynski. Lectures on Stochastic Programming: Modeling andTheory. MOS-SIAM Series on Optimization. Society for Industrial and Applied Mathematics, 2009.
[22] M. Udell, K. Mohan, D. Zeng, J. Hong, S. Diamond,and S. Boyd. Convex optimization in Julia. In Proceedings of the 1st First Workshop for High Performance Technical Computing in Dynamic Languages,
2014.
[23] S. Uryasev and P. Pardalos. Stochastic Optimization.
Applied Optimization. Springer, 2001.
[24] C. Valente, G. Mitra, M. Sadki, and R. Fourer. Extending algebraic modelling languages for stochastic programming. INFORMS Journal on Computing,21(1):107–122, 2009.
[25] S. Wallace and W. Ziemba. Applications of Stochastic Programming. MPS-SIAM Series on Optimization. Society for Industrial and Applied Mathematics,
2005.

[12] H. Markowitz. Portfolio selection. The Journal of
Finance, 7(1):77–91, 1952.

[26] J. Watson, D. Woodruff, and W. Hart. PySP: Modeling and solving stochastic programs in Python. Mathematical Programming Computation, 4(2):109–149,
2012.

[13] T. Minka, J. Winn, J. Guiver, S. Webster, Y. Zaykov, B. Yangel, A. Spengler, and J. Bronskill. Infer.NET 2.6, 2014. Microsoft Research Cambridge.
http://research.microsoft.com/infernet.

[27] Z. Wu, M. Kusner, K. Weinberger, M. Chen, andO. Chapelle. Classifier cascades and trees for minimizing feature evaluation cost. Journal of Machine
Learning Research, 15:2113–2144, 2014.

[14] A. Nemirovski and A. Shapiro. Convex approximations of chance-constrained programs. SIAM Journal
on Optimization, 17(4):969–996, 2006.