Max-affine estimators for convex stochastic programming

arXiv:1609.06331v1 [math.OC] 18 Sep 2016

Gábor Balázs and Csaba Szepesvári
Department of Computing Science, University of Alberta,
Edmonton, Alberta, T6G 2E8, Canada
András György
Department of Electrical and Electronic Engineering,
Imperial College London, London, SW7 2BT, United Kingdom
January 11, 2018

Abstract
In this paper, we consider two sequential decision making problems with a convexity
structure, namely an energy storage optimization task and a multi-product assembly
example. We formulate these problems in the stochastic programming framework and
discuss an approximate dynamic programming technique for their solutions. As the
cost-to-go functions are convex in these cases, we use max-affine estimates for their approximations. To train such a max-affine estimate, we provide a new convex regression
algorithm, and evaluate it empirically for these planning scenarios.

1

Introdution

This paper considers multi-stage stochastic programming problems (e.g., Shapiro et al.,
2009; Birge and Louveaux, 2011) with restricting the model so that the cost-to-go functions remain convex (Section 2). To motivate this framework, we provide two realistic
benchmark planning problems (Section 5), namely solar energy production with storage,
and the operation of a beer brewery.
To address these problems, we consider an approximate dynamic programming approach
(e.g., Bertsekas, 2005; Powell, 2011). For this, we estimate the cost-to-go functions by
convex regression techniques using max-affine representations (formed by the maximum
of finitely many affine functions). To train a max-affine estimator, we propose a novel
algorithm (AMAP, Section 3), which combines ideas of the convex adaptive partitioning
method (CAP, Hannah and Dunson, 2013) and the least squares partition algorithm (LSPA,
Magnani and Boyd, 2009), while learns the model size by cross-validation.
We discuss a full approximate dynamic programming approach (Section 4) that estimates the cost-to-go functions globally. It uses a forward pass to generate a data set for
the cost-to-go estimation by uniformly sampling the reachable decision space. Then it performs a backward pass to recursively approximate the cost-to-go functions for all stages.
The technical details of this algorithm is provided for polyhedral decision sets and convex
piecewise-linear cost-to-go representations.
Finally, we evaluate max-affine estimators in the contexts of our benchmark planning
problems through the full approximate dynamic programming algorithm (Section 5).

1

2

Convex stochastic programming

Consider a T -stage stochastic programming (SP) problem (e.g., Ruszczyński and Shapiro,
2003; Shapiro et al., 2009; Birge and Louveaux, 2011), where the goal is to find a decision
x∗1 solving the following:
h
i
.
x∗1 ∈ argmin J1 (x1 ) , Jt (xt ) = E ct (xt , Z t ) +
min
Jt+1 (xt+1 ) ,
(1)
xt+1 ∈Xt+1 (xt ,Z t )

x1 ∈X1 (x0 ,z0 )

.
.
with t = 1, . . . , T , some fixed initial values x0 , z0 , XT +1 (xT , Z T ) = {⊥}, JT +1 (⊥) = 0, and
a sequence of independent random variables Z 1 , . . . , Z T .
Notice that (1) includes discrete-time finite-horizon Markov decision process formulations (e.g., Puterman, 1994; Sutton, 1998; Bertsekas, 2005; Szepesvári, 2010; Powell, 2011)
after the state and action variables are combined into a single decision variable xt , and
the environment dynamics along with the action constraints are described by the decision
constraint functions Xt .
In this text, we consider only a subset of SP problems (1) when the cost functions
c1 , . . . , cT are convex in xt , and graph(Xt (xt , Z t )) are convex sets for all t = 1, . . . , T and
all Z t realizations,
where the graph of a set-valued function C : X → 2Y is defined as
. 
graph(C) = (x, y) ∈ X × Y : y ∈ C(x) . In this case, the cost-to-go functions Jt (·) are
convex for all t = 1, . . . , T (e.g., see Lemma 1), hence we call these SP problems convex.
One approach to deal with such SP problems (1) is to use approximate dynamic programming (ADP) methods (e.g., Bertsekas, 2005; Powell, 2011; Birge and Louveaux, 2011;
Hannah et al., 2014) which construct nested approximations to the cost-to-go functions,


ˆ
ˆ
Jt+1 (xt+1 ) ≈ Jt (xt ) ,
(2)
Jt (xt ) ≈ E ct (xt , Z t ) +
min
xt+1 ∈Xt+1 (xt ,Z t )

backwards for t = T, T −1, . . . , 1. Notice that for convex SP problems with convex cost-to-go
functions Jt , the estimates Jˆt can be restricted to be convex without loss of generality, and
so the minimization tasks in (2) can be solved efficiently.

3

Adaptive max-affine partitioning algorithm

To represent convex cost-to-go functions, we use max-affine maps formed by the maximum
of finitely many affine functions (hyperplanes). To train the parameters of such a maxaffine estimate, we present a new convex regression algorithm called Adaptive Max-Affine
Partitioning (AMAP), which combines the partitioning technique of the convex adaptive
partitioning algorithm (CAP, Hannah and Dunson, 2013), the least squared partition algorithm (LSPA, Magnani and Boyd, 2009), and learns the model size (number of hyperplanes)
by a cross-validation scheme.
Just as LSPA and CAP, AMAP also P
aims to reduce the empirical risk with respect to
. 1 n
the squared loss defined as Rn (f ) = n i=1 |f (xi ) − yi |2 for some function f and data
.
set Dn = {(xi , yi ) : i = 1, . . . , n} ⊆ (Rd × R)n with n samples and dimension d. For the
discussion of AMAP, denote the max-affine function of model M = {(ak , bk ) : k = 1, . . . , K}
.
d
by fM (x) = maxk=1,...,K a⊤
k x+bk for x ∈ R . Notice that each max-affine model M induces
a partition over the data set Dn as
.
P = {C1 , . . . , CK } ,

. 
Ck = i ∈ {1, . . . , n} x⊤
i ak + bk = fM (xi ) ,
2

(3)

for some K ∈ N and all k = 1, . . . , K, where ties are broken arbitrarily so that Ck 6= ∅ and
Ck ∩ Cl 6= ∅ ⇐⇒ k = l for all k, l = 1, . . . , K. Furthermore, each partition P induces a
max-affine model by fitting each cell of P using the linear least squares algorithm:
1 X
.
xi ,
∆ik = xi −
|Ck |
i∈Ck

−1 X
X
1
. X
.
ak =
+
βI
∆ik ∆⊤
∆
y
,
b
=
yi − a⊤
i
d
ik
k
ik
k xi ,
|Ck |
.
M = {(ak , bk ) : k = 1, . . . , K} ,

i∈Ck

i∈Ck

(4)

i∈Ck

where Id is the d × d identity matrix, and β > 0 is set to some small value for stability (we
.
use β = 10−6 ).
Similar to CAP, AMAP builds the model incrementally by cell splitting, and improves
the partition using LSPA by alternating steps (3) and (4). The AMAP model improvement
step is given by Algorithm 1. Notice that AMAP performs coordinate-wise cell splitting
1.

2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.

input: training set Dn , model M = {(ak , bk ) : k = 1, . . . , K},
partition P = {C1 , . . . , CK }, empirical risk E, minimum cell size s∗
{cell splitting}
M∗′ ← M, E∗′ ← E, P∗′ ← P
for all k = 1, . . . , K do
if |Ck | ≥ 2s∗ then
for all j = 1, . . . , d do
mj ← median({Xij : i ∈ Ck })
Cle ← {i ∈ Ck : Xij ≤ mj }, Cgt ← {i ∈ Ck : Xij ≥ mj }
(ale , ble ) ← least squares on {(X i , Yi ) : i ∈ Cle }
(agt , bgt ) ← least squares
 on {(X i , Yi ) : i ∈ Cgt }
M ′ ← M \ {(ak, bk )} ∪ {(ale , ble ), (agt , bgt )}
E ′ ← Rn fM ′ (·)
if E ′ < E∗′ then
M∗′ ← M ′ , E∗′ ← E ′ , P∗′ ← (P \ {Ck }) ∪ {Cle , Cgt }
end if
end for
end if
end for
{running LSPA}
repeat
M∗ ← M∗′ , E∗ ← E∗′ , P∗ ← P∗′
P∗′ ← induced partition of M∗ by (3)
M∗′ ← fitting of partition
P∗′ by (4)

′
E∗ ← Rn fM∗′ (·)
until minC∈P∗′ |C| ≥ s∗ and E∗′ < E∗
output: model M∗ , partition P∗ , empirical risk E∗
Algorithm 1: AMAP model improvement step.

(steps 5 to 15), just as CAP, however, AMAP makes the split always at the median (steps 6
and 7) instead of checking multiple cut points. This saves computation time, but can also
create worse splits. To compensate for this loss in quality, AMAP runs a restricted version
of LSPA (steps 18 to 23) not just for a single step as CAP, but until the candidate model
3

improves the empirical risk and its induced partition satisfies the minimum cell requirement
(step 23). We also mention that indices {i ∈ Ck : Xij = mj } are assigned to Cle and Cgt
(step 7) in order to preserve the minimum cell requirement.
Notice that the difference between M ′ and M is only two hyperplanes (step 10), so the
number of arithmetic operations for computing E ′ (step 11) can be improved from O(nKd)
to O(nd). Further, the cost of least squares regressions (steps 8 and 9) is O(|Ck |d2 ). Hence,
the computational cost of the entire cell splitting process (steps 2 to 17) is bounded by
O(max{K, d}d2 n). For the LSPA part, the partition fitting (step 21) is O(nd2 ) and the
error calculation (step 22) is O(nKd). So, the cost of a single LSPA iteration (steps 19
to 22) is bounded by O(max{K, d}dn),
implying that the cost of Algorithm 1 is bounded

by O max{tLSPA , d} max{K, d}dn , where tLSPA denotes the number of LSPA iterations.
Undesirably, coordinate-wise cell splitting does not provide rotational invariance. To fix
this, we run AMAP after a pre-processing step, which uses thin singular value decomposition
(thin-SVD) to drop redundant coordinates and align the data along a rotation invariant
basis. Formally, let the raw (but already centered) data be organized into X ∈ Rn×d and
.
y ∈ Rn . Then, we scale the values [y1 . . . yn ]⊤ = y/ max{1, kyk∞ }, and decompose X
by thin-SVD as X = U SV ⊤ , where U ∈ Rn×d is semi-orthogonal, S ∈ Rd×d is diagonal
with singular values in decreasing order, and V ∈ Rd×d is orthogonal. Coordinates related
.
to zero singular values are dropped1 and the points are scaled by S as [x1 . . . xn ]⊤ =
U S/ max{1, S11 }, where S11 is the largest singular value. Now observe that rotating the
raw points X as XQ with some orthogonal matrix Q ∈ Rd×d only transforms V to Q⊤ V
and does not affect the pre-processed points x1 , . . . , xn . Finally, we note that thin-SVD
can be computed using O(nd2 ) arithmetic operations (with n ≥ d),2 which is even less than
the asymptotic cost of Algorithm 1.
AMAP is presented as Algorithm 2 and run using uniformly shuffled (and pre-processed)
data Dn , and a partition {F1 , . . . , Fu∗ } of {1, . . . , n} with equally sized cells (the last one
might be smaller) defining the cross-validation folds.
For model selection, AMAP uses a u∗ -fold cross-validation scheme (steps 8 to 19) to
find an appropriate model size that terminates when the cross-validation error (step 14) of
the best model set M∗ (steps 7 and 16) cannot be further improved for twait iterations. At
the end, the final model is chosen from the model set M∗ having the best cross-validation
error, and selected so to minimize the empirical risk on the entire data (step 20). For this
.
.
scheme, we use the parameters u∗ = 10 and twait = 5.
AMAP starts with models having a single hyperplane (steps 2 to 6) and increments each
model by at most one hyperplane in every iteration (step 11). Notice that if AMAP cannot
find a split for a model Mu to improve the empirical risk Eu , the update for model Mu
(steps 11 and 12) can be skipped in the subsequent iterations as Algorithm 1 is deterministic.
.
We also mention that for the minimum cell size, we use s∗ = max{2(d + 1), ⌈log 2 (n)⌉}
allowing model sizes up to O nd/(d+4) / ln(n) , which is enough for near-optimal worst-case
performance (Balázs et al., 2015, Theorems 4.1 and 4.2).
1

By removing columns of U and V , and columns and rows of S.
First decompose X by a thin-QR algorithm in O(nd2 ) time (Golub and Loan, 1996, Section 5.2.8) as
X = QR, where Q ∈ Rn×d has orthogonal columns and R ∈ Rd×d is upper triangular. Then apply SVD for
R in O(d3 ) time (Golub and Loan, 1996, Section 5.4.5).
2

4

1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.

input: training set Dn , minimum cell size s∗ , folds F1 , . . . , Fu∗ , iterations twait
{initialization}
for u = 1, . . . , u∗ do
.
Pu ← {Fu } with Fu = {1, . . . , n} \ Fu
u )} ← least squares on {(x , y ) : i ∈ F }
Mu = {(au1 , bP
i i
1
Pu
Eu ← |Fu |−1 i∈Fu |fMu (xi ) − yi |2 , Eu ← |Fu |−1 i∈Fu |fMu (xi ) − yi |2
end for
P ∗
Eu
M∗ ← {M1 , . . . , Mu∗ }, E∗ ← u1∗ uu=1
{cross-validation training}
t ← 1, tmax ← twait
while t ≤ min{tmax , ⌈nd/(d+4) ⌉} do
for u = 1, . . . , u∗ do
(Mu , Pu , Eu ) P
← update by Algorithm 1 using {(xi , yi ) : i ∈ Fu }, Mu , Pu , Eu , s∗
2
−1
Eu ← |Fu |
i∈Fu |fMu (xi ) − yi |
end for P
∗
E ← u1∗ uu=1
Eu
if E < E∗ then
M∗ ← {M1 , . . . , Mu∗ }, E∗ ← E, tmax ← t + twait
end if
t←t+1
end while
{choosing the final model} 
M∗ ← argminM ∈M∗ Rn fM (·)
output: model M∗
Algorithm 2: Adaptive max-affine partitioning (AMAP).

4

Approximate dynamic programming

Here we use max-affine estimators to approximate the cost-to-go functions of convex SP
problems. First, notice that solving (1) is equivalent to the computation of π1 (x0 , z0 ),
where
.
πt (xt−1 , zt−1 ) =
argmin
Jt (xt ) ,
xt ∈Xt (xt−1 ,zt−1 )
(5)


Jt (xt ) = E ct (xt , Z t ) + Jt+1 (πt+1 (xt , Z t )) ,

.
.
for all t = 1, . . . , T , and πT +1 (·, ·) = ⊥ with JT +1 (⊥) = 0. The sequence π = (π1 , . . . , πT )
represents an optimal policy.
We only consider SP problems with convex polyhedral decision constraints written as
Xt+1 (xt , Z t ) = {xt+1 : Qt+1 xt+1 + Wt+1 (Z t )xt ≤ ct+1 (Z t )} which are non-empty for all
possible realizations of xt and Z t . As the coefficient Qt+1 of the decision variable xt+1 is
independent of random disturbances Z t and the constraint xt ∈ Xt (xt−1 , zt−1 ) for policy
πt is feasible for any xt−1 and zt−1 , these SP problems are said to have a fixed, relatively
complete recourse (Shapiro et al., 2009, Section 2.1.3). We will exploit the fixed recourse
property for sampling (6), while relatively complete recourse allows us not to deal with
infeasibility issues which could make these problems very difficult to solve.3
3

Infeasible constraints can be equivalently modeled by cost-to-go functions assigning infinite value for
points outside of the feasible region. Then notice that the estimation of functions with infinite magnitude
and slope can be arbitrarily hard even for convex functions (Balázs et al., 2015, Section 4.3).

5

In order to construct an approximation Jˆt to the cost-to-go function Jt , we need “realizable” decision samples xt,i at stage t. We generate these incrementally during a forward pass
for t = 1, 2, . . . , T , where given n decisions xt,1 , . . . , xt,n and m disturbances zt,1 , . . . , zt,m
at stage t, we uniformly sample new decisions for stage t + 1 from the set

 [
.
Xt+1 (xt,i , zt,j )
X̂t+1 = conv
i=1,...,n,
j=1,...,m

n

= xt+1 : Qt+1 xt+1

(6)

o

≤ max ct+1 (zt,j ) − Wt+1 (zt,j )xt,i ,
i=1,...,n,
j=1,...,m

where conv(·) refers to the convex hull of a set, and the maximum is taken component-wise.
To uniformly sample the convex polytope X̂t+1 , we use the Hit-and-run Markov-Chain
Monte-Carlo algorithm (Smith, 1984, or see Vempala, 2005) by generating 100 chains (to
reduce sample correlation) each started from the average of 10 randomly generated border
points, and discarding d2t+1 samples on each chain during the burn-in phase,4 where dt+1 is
the dimension of xt+1 .
Then, during a backward pass for t = T, T − 1, . . . , 1, we recursively use the cost-to-go
estimate of the previous stage Jˆt+1 (·) to approximate the values of the cost-to-go function Jt
at the decision samples xt,i generated during the forward pass, that is
m

. 1 X
Jˆt+1 (xt+1 )
ct (xt,i , zt,j ) +
min
Jt (xt,i ) ≈ yt,i =
m
xt+1 ∈Xt+1 (xt,i ,zt,j )

(7)

j=1

for all t = 1, . . . , T , and JˆT +1 (·) ≡ 0. This allows us to set up a regression problem with
data {(xt,i , yt,i ) : i = 1, . . . , n}, and to construct an estimate Jˆt (·) of the cost-to-go Jt (·).
We call the resulting method, shown as Algorithm 3, full approximate dynamic programming (fADP) because it constructs global approximations to the cost-to-go functions.
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.

input: SP problem, trajectory count n, evaluation count m, estimator REG
x0,i ← x0 for all i = 1, . . . , n
z0,j ← z0 for all j = 1, . . . , m
{forward pass}
for all t = 0, 1, . . . , T − 1 do
sample xt+1,1 , . . . , xt+1,n from X̂t+1 by (6) using xt,1 , . . . , xt,n , and zt,1 , . . . , zt,m
sample zt+1,1 , . . . , zt+1,m from the distribution of Z t+1
end for
{backward pass}
JˆT +1 (·) ← 0
for all t = T, T − 1, . . . , 1 do
compute yt,1 , . . . , yt,n by (7) using Jˆt+1 (·), xt,1 , . . . , xt,n , and zt,1 , . . . , zt,m
Jˆt (·) ← REG({(xt,i , yt,i ) : i = 1, . . . , n})
end for
output: cost-to-go functions Jˆt (·), t = 1, . . . , T
Algorithm 3: Full approximate dynamic programming (fADP).

When the cost-to-go functions Jt are approximated by “reasonably sized” convex piecewise linear representations Jˆt , the minimization problem in (7) can be solved efficiently by
4

The choice was inspired by the O(d2t+1 ) mixing result of the Hit-and-run algorithm (Lovász, 1999).

6

linear programming (LP). In the following sections, we exploit the speed of LP solvers for
fADP using either AMAP or CAP as the regression procedure. Then, the computational
cost to run Algorithm 3 is mostly realized by solving T nm LP tasks for (7) and training T
estimators using the regression algorithm REG.
Here we mention that max-affine estimators using as many hyperplanes as the sample
size n can be computed by convex optimization techniques efficiently up to a few thousand
samples (Mazumder et al., 2015). These estimators provide worst-case generalization error
guarantees, however, the provided max-affine models are large enough to turn the LP tasks
together too costly to solve (at least using our hardware and implementation).
The situation is even worse for nonconvex estimators REG for which LP has to be replaced for (7) by a much slower nonlinear constrained optimization method using perhaps
randomized restarts to minimize the chance of being trapped in a local minima. Furthermore, when the gradient information with respect to the input is not available (as
for many implementations), the minimization over these representations require an even
slower gradient-free nonlinear optimization technique. Hence, multivariate adaptive regression splines, support vector regression, and neural networks were impractical to use for
fADP on our test problems.

5

Experiments

To test fADP using max-affine estimators, we consider two SP planning problems, namely
solar energy production with storage management (Section 5.1), and the operation of a
beer brewery (Section 5.2).
For our numerical experiments, the hardware has a Dual-Core AMD Opteron(tm) Processor 250 (2.4GHz, 1KB L1 Cache, 1MB L2 Cache) with 8GB RAM. The software uses
MATLAB (R2010b), and the MOSEK Optimization Toolbox (v7.139).
To measure the performance of the fADP algorithm, we evaluate the greedy policy
with respect to the learned cost-to-go functions {Jˆt : t = 1, . . . , T }. More precisely, we
.
run π̂ = (π̂1 , . . . , π̂T ) with π̂t (xt−1 , zt−1 ) ∈ argminxt ∈Xt (xt−1 ,zt−1 ) Jˆt (xt ) on 1000 episodes,
.
(e) (e) 
1 P1000 PT
and record the average revenue (negative cost) as REV = − 1000
e=1
t=1 ct xt , zt
 (e) (e) 
: t = 1, . . . , T , e = 1, . . . , 1000. We repeat
over the episodes’ trajectories
xt , zt
this experiment 100 times for each regression algorithm REG,5 and show the mean and
standard deviation of the resulting sample.

5.1

Energy storage optimization

Inspired by a similar example of Jiang and Powell (2015, Section 7.3), we consider an
energy storage optimization problem where a renewable energy company makes a decision
every hour and plans for two days (T = 48). The company owns an energy storage with
state s which can be charged with maximum rate rc , using the company’s renewable energy
source (E) or the electrical grid that the company can buy electricity from while paying the
retail price (p). The goal is to maximize profit by selling electricity to local clients on retail
price (p) according to their stochastic demand (D) or selling it back to the electrical grid
on wholesale price (w). Electricity can be sold directly from the renewable energy source
5

The random seeds are kept synchronized, so every algorithm is evaluated on the same set of trajectories.
Furthermore, fADP algorithms with the same n and m parameters use the same training data xt,1 , . . . , xt,n
and zt,1 , . . . , zt,m for all t = 0, . . . , T .

7

or from the battery with maximum discharge rate rd . The energy flow control variables,
fes , fed , feg , fsd , fsg , fgs , are depicted on Figure 1.

Figure 1: Flow diagram of a convex energy storage problem.
The SP model (1) of the energy storage problem can be formulated by using the decision
.
.
variable xt = [st fes,t fed,t feg,t fsd,t fsg,t fgs,t ]⊤ ∈ [0, ∞)7 and setting Z t = [Et+1 Dt+1 ]⊤ .
Furthermore, the cost function is defined as
.
ct (xt , Z t ) = pt (fgs,t − fed,t − fsd,t ) − wt (feg,t + fsg,t ) ,
for all t = 1, . . . , T , and the dynamics and control constraints are described by



s










f
s
=
s
+
f
−
f
−
f
+
f
,


es
t
es,t
sg,t
gs,t
sd,t










fes , fed , feg , fsd , fsg , fgs ≥ 0,

 fed 
. 

,
0 ≤ s + fes − fsd − fsg + fgs ≤ smax ,
Xt+1 (xt , Z t ) =  feg 






f
f
+
f
≤
r
,
f
+
f
≤
r
,


es
gs
c
sg
sd
d
 sd 





 fsg  fes + fed + feg ≤ Et+1 , fed + fsd ≤ Dt+1 






fgs

.
.
for all t = 0, . . . , T − 1. To initialize the system, define x0 = [s0 0 . . . 0]⊤ and z0 = [d1 e1 ]⊤ ,
where s1 = s0 ∈ [0, smax ] is the current storage level and d1 , e1 ≥ 0 are the currently
observed demand and energy production, respectively.
To set the parameters of this energy storage optimization problem, we consider a solar
energy source, a discounted nightly electricity pricing model (Economy 7 tariff), and plan.
ning for a two days horizon on hourly basis (T = 48). Retail and expected wholesale price
curves, along with the electricity demand and energy production distributions of this model
are shown on Figure 2 for the two consecutive sunny days. The distributions are truncated
normal with support D ∈ [0, 15] and E ∈ [0, 12] for demand and energy production, respec.
.
tively and the storage has capacity smax = 20 with charge and discharge rates rc = 4 and
.
.
.
.
rd = 10, respectively. The model is initialized by s0 = 0, d1 = E[D1 ], and e1 = E[E1 ].
To evaluate fADP on this problem, we use the CAP6 and AMAP convex regression
algorithms, and multiple configurations determined by the number of trajectories n generated for training (which is the sample size for the regression tasks as well), and the number
of evaluations m used to approximate the cost-to-go functions Jt at a single point (7). The
result is presented on Figure 3, which also includes a “heuristic” algorithm to provide a
6

As the implementation of CAP is not too reliable for highly-correlated features, we combined it with
the data preprocessing step of AMAP (see Section 3) which improved its stability.

8

retail price (p)
wholesale price (w)

demand (D)
energy (E)

25
20
15
10
5
0
1

12

24

36

48

hours
Figure 2: Parameters of the energy storage optimization problem. Retail (p) and expected
wholesale (w) price curves, energy demand (D) and production (E) distributions with mean
and standard deviation are shown for two-day long period.
fADP(CAP,m = 10)
fADP(AMAP,m = 10)
30

3700

minutes

revenue

3800

fADP(CAP,m = 100)
fADP(AMAP,m = 100)

heuristic

3600
3500

15

0
25

50

100

25

150

50

100

150

training samples per stage (nm, in 1000s)
Figure 3: Energy storage optimization results for the fADP algorithm using AMAP or CAP
as the inner convex regression procedure. Results show the total revenue and the training
time in minutes with mean and standard deviation for the number of trajectories n and
cost-to-go evaluations m.
baseline. The heuristic uses a fixed policy of immediately selling the solar energy preferrably
for demand (fed → max, feg ≥ 0, fes = 0), selling from the battery during the day when
demand still allows (fgs = 0, fsd ≥ 0), charging the battery overnight (fgs → max, fsd = 0),
and selling everything close to the last stage (fgs = 0, fsd → max, fsg → max). This policy
is much better than the optimal policy without storage 7 which scores 3227 ± 6.
The results of Figure 3 show that fADP using either CAP or AMAP significantly outperforms the heuristic baseline algorithm when the sample size is large enough. The regression
algorithms prefer larger sample sizes n to better sample quality m, although this significantly increases the computation time for CAP and provides only an insignificant revenue
increase compared to AMAP.
.
Because p ≥ w, the optimal policy for smax = 0 minimizes the immediate cost by fed = min{E, D},
.
.
.
.
.
= max{0, E − fed }, and fes = fsd = fgs = fsg = 0.

7

feg

9

5.2

Beer brewery operation

Inspired by Bisschop (2016, Chapter 17), we consider the multi-product assembly problem
of operating a beer brewery which makes a decision in every two weeks and plans for about
one year (48 weeks, T = 24). The factory has to order ingredients (stratch source, yeast,
hops) to produce two types of beers (ale and lager) which have to be fermented (for at least
2 weeks for ale and 6 weeks for lager) before selling. The states and actions of this process
are illustrated on Figure 4.

Figure 4: Diagram of the beer brewery problem. The state includes three types of ingredients and two types of bottled beer in the warehouse, and the four fermentation tanks.
Actions are formed by ingredient orders ur , brewing amounts ub , and beer sales us .
The decision variable xt is a 16 dimensional vector with the following components:


stratch source in storage


yeast in storage




hops in storage


 ale beer fermented for less than 2 weeks 




produced ale beer


 lager beer fermented for less than 2 weeks 




lager beer fermented for 2 to 4 weeks




lager
beer
fermented
for
4
to
6
weeks


. 
 ∈ [0, ∞)16 .
xt = 

produced lager beer




stratch source order






yeast order




hops order




ale beer brewing




lager beer brewing




ale beer sales
lager beer sales
The first 9 coordinates are state variables, while the last 7 coordinates represent actions.
.
⊤
⊤
The cost functions (which may take negative values) are ct (xt , Z t ) = [h⊤
t ct −rt ]xt , where
9
5
ht ∈ [0, ∞) is the storage cost, ct ∈ [0, ∞) is the market price of the ingredients with
the brewing costs (adjusted by the water price), and rt ∈ [0, ∞)2 is the selling price of the

10

beers, for stages t = 1, . . . , T . The constraint on the dynamics of the system is given by



F xt + Rur + Bub − Sus



ur , ub ≥ 0, us ∈ [0, Dt+1 ], 



u
.
r

Xt+1 (xt , Z t ) = 
F
x
+
Bu
−
Su
≥
0,
,
t
s
b


ub




F xt + Rur + Bub ≤ kt+1 

us

for all t = 1, . . . , T −1, where Z t = Dt+1 ∈ [0, ∞)2 is the beer demand, kt+1 ∈ [0, ∞]9 is the
capacity bound, 0 denotes a zero vector of appropriate size, and the fermentation matrix
F ∈ {0, 1}9×16 , the brewing matrix B ∈ R9×2 , the storage loading matrix R ∈ {0, 1}9×3
and the selling matrix S ∈ {0, 1}9×2 are defined as






−ba −bl
03×2
03×13
I3
 1 0 


 00 

00
. 




0 0 
B=
,
02×4



 10 

1
1
 0 1 


. 
. 
,
S=
F =
0000
00 


 06×3
,
06×7 
03×2
 00 

1
0
0
0




04×2

 00 



0100
I3
.
,
R=
0011
01
06×3
where 0m×n is a zero matrix of size m × n, and ba , bl ∈ [0, ∞)3 are the required ingredients
.
.
for brewing ales and lagers, respectively. To initialize the system, we use x0 = 0 and z0 = 0.
To test the fADP algorithm on this problem, we set the horizon to 48 weeks horizon
.
on a fortnight basis (T = 24). The demand distributions Z t for lager and ale beers are
shown on Figure 5 with mean and standard deviation. Both distributions are truncated
lager demand

ale demand

12
10
8
6
4
2
0
1

6

12

18

24

fortnights
Figure 5: Lager and ale beer demand distributions for the beer brewery optimization problem with mean and standard deviation are shown for a 48 weeks horizon.
normal with support [0.1, 12]. The cost vectors are set to fixed values for all t = 1, . . . , T
.
.
.
as ht = [1 0.2 0.2 1 2 1 1 1 2]⊤ , ct = [20 10 5 1 1]⊤ , and rt = [90 50]⊤ . Furthermore, the
.
.
ingredient requirement vectors for brewing are ba = [1 1 1]⊤ and bl = [0.5 0.9 0.8]⊤ for
.
ale and lager, respectively, and the capacity vector is k = [10 10 10 10 ∞ 10 ∞ ∞ ∞]⊤ to
ensure the feasibility (relatively complete recourse) requirements, as discussed in Section 4.
11

Similar to the energy optimization case, we use the CAP and AMAP estimators for
fADP with various trajectory set sizes n and cost-to-go evaluation numbers m. As a baseline, we use the optimal solution for the deterministic version of (1) where Z t is replaced by
its expectation EZ t for all t = 1, . . . , T . The results are presented on Figure 6, showing that
fADP(CAP,m = 10)
fADP(AMAP,m = 10)

5000
4000
3000
2000
1000
0

3
hours

revenue

fADP(CAP,m = 100)
fADP(AMAP,m = 100)

2
1
0

25 50

100

150

200

250

25 50

100

150

200

250

training samples per stage (nm, in 1000s)
Figure 6: Beer brewery optimization results for fADP algorithm using AMAP and CAP
convex regression to approximate the convex cost-to-go functions. Results show the revenue
with mean and standard deviation, and the training time in minutes for trajectories n and
cost-to-go evaluations m.
AMAP improves the performance significantly by collecting revenue over 4100 compared
to CAP which stays around 3600.
However, the result also shows that the running time of AMAP become significantly
larger than CAP. This indicates that larger max-affine models are trained by AMAP compared to CAP to improve the accuracy of the cost-to-go approximations, that increases the
computational cost of the LP tasks of (7), and eventually slows down the fADP algorithm
significantly.
Finally, notice that using larger trajectory sets n for AMAP provide better performance
for low sample levels nm, but the improved sample quality with m = 100 eventually achieves
the same using significantly less computational resources. So it seems to be still possible
for max-affine estimators to find better tradeoffs between accuracy and model size.

6

Conclusions and future work

In this paper, we considered solving convex SP problems by the combination of an approximate dynamic programming method (fADP) and convex regression techniques. For the
latter, we proposed a new state-of-the-art max-affine estimator (AMAP), and combined it
with an approximate dynamic programming algorithm to address two benchmark convex
SP problems of moderate size.
Clearly, scaling up the fADP algorithm for larger problems would require further improvements. One of these could be using more expressive convex piecewise-linear representations (e.g., sum-max-affine maps), which might compress the LP tasks enough for the
current solvers. For this, Hannah and Dunson (2012) considered various ensemble tehniques
(bagging, smearing, and random forests) to enhance the performance of the CAP estimator.
However, these techniques still seem to construct too large models to increase the accuracy
significantly, that makes the vast amount of LP tasks impractical to solve. Maybe, LP algo12

rithms which can more efficiently solve large number of similar LP problems with different
right hand sides (e.g., Gassmann and Wallace, 1996) could help with this issue.
But eventually, it would become inevitable to localize cost-to-go approximations to a
fraction of the decision space, perhaps by running fADP iteratively alternating between
sampling and estimation phases, and exploring at the boundary of the accurately approximated region of cost-to-go functions in order to find and avoid delayed rewards and costs.
However, this is left for future work.

A

Auxiliary tools

The following result is a slight generalization of Theorem 5.3 in Rockafellar (1972).
Lemma 1. Let X, Y be two convex sets and f : X × Y be a jointly–convex function in its
arguments. Additionally, let C : X → 2Y be a set–valued function for which graph(C) is
.
convex. Then g(x) = inf y∈C(x) f (x, y) is a convex function.
Proof. Let x1 , x2 ∈ X, y1 , y2 ∈ Y and λ ∈ (0, 1). As graph(C) is convex,
y1 ∈ C(x1 ), y2 ∈ C(x2 )

⇒


λy1 + (1 − λ)y2 ∈ C λx1 + (1 − λ)x2 .

Using this with the fact that the infimum on a subset becomes larger, and the joint–
convexity of f , we get


inf
f λx1 + (1 − λ)x2 , z
g λx1 + (1 − λ)x2 =
z∈C(λx1 +(1−λ)x2 )

≤ inf
inf f λx1 + (1 − λ)x2 , λy1 + (1 − λ)y2
y1 ∈C(x1 ) y2 ∈C(x2 )

≤

inf

inf

y1 ∈C(x1 ) y2 ∈C(x2 )

λf (x1 , y1 ) + (1 − λ)f (x2 , y2 )

= λg(x1 ) + (1 − λ)g(x2 ) ,
which proves the convexity of g.



References
G. Balázs, A. György, and C. Szepesvári. Near-optimal max-affine estimators for convex
regression. In G. Lebanon and S. Vishwanathan, editors, The 18th International Conference on Artificial Intelligence and Statistics (AISTATS), volume 38 of JMLR W&CP,
pages 56–64, 2015.
D. P. Bertsekas. Dynamic Programming and Optimal Control, Volume I. Athena Scientific,
3rd edition, 2005.
J. R. Birge and F. Louveaux. Introduction to Stochastic Programming. Springer, 2011.
J. Bisschop. AIMMS Optimization Modeling. 2016. http://www.aimms.com.
H. I. Gassmann and S. W. Wallace. Solving linear programs with multiple right-hand sides:
Pricing and ordering schemes. Annals of Operation Research, 64:237–259, 1996.
G. H. Golub and C. F. V. Loan. Matrix Computations. The Johns Hopkins University
Press, 3rd edition, 1996.
13

L. A. Hannah and D. B. Dunson. Ensemble methods for convex regression with applications
to geometric programming based circuit design. In J. Langford and J. Pineau, editors,
The 29th International Conference on Machine Learning (ICML), pages 369–376, 2012.
L. A. Hannah and D. B. Dunson. Multivariate convex regression with adaptive partitioning.
Journal of Machine Learning Research, 14:3261–3294, 2013.
L. A. Hannah, W. B. Powell, and D. B. Dunson. Semiconvex regression for metamodelingbased optimization. SIAM Journal on Optimization, 24(2):573–597, 2014.
D. R. Jiang and W. B. Powell. An approximate dynamic programming algorithm for
monotone value functions. CoRR, 2015. http://arxiv.org/abs/1401.1590v6.
L. Lovász. Hit-and-run mixes fast. Mathematical Programming, 86(3):443–461, 1999.
A. Magnani and S. P. Boyd. Convex piecewise-linear fitting. Optimization and Engineering,
10(1):1–17, 2009.
R. Mazumder, A. Choudhury, G. Iyengar, and B. Sen. A computational framework for multivariate convex regression and its variants. 2015. http://arxiv.org/abs/1509.08165.
W. B. Powell. Approximate Dynamic Programming: Solving the Curses of Dimensionality.
John Wiley & Sons, 2nd edition, 2011.
M. L. Puterman. Markov Decision Processes: Discrete Stochastic Dynamic Programming.
John Wiley & Sons, 1994.
R. T. Rockafellar. Convex Analysis. Princeton University Press, 1972.
A. P. Ruszczyński and A. Shapiro. Stochastic Programming. Elsevier, 2003.
A. Shapiro, D. Dentcheva, and A. Ruszczyński. Lectures on Stochastic Programming, Modeling and Theory. Society for Industrial and Applied Mathematics and the Mathematical
Programming Society, 2009.
R. L. Smith. Efficient monte carlo procedures for generating points uniformly distributed
over bounded regions. Operations Research, 32(6):1296–1308, 1984.
R. S. Sutton. Reinforcement Learning: An Introduction. MIT Press, 1998.
C. Szepesvári. Algorithms for Reinforcement Learning: Synthesis Lectures on Artificial
Intelligence and Machine Learning. Morgan & Claypool Publishers, 2010.
S. Vempala. Geometric random walk: A survey. Combinatorial and Computational Geometry, 52:573–612, 2005.

14

