Advanced Topics in Stochastic Programming - Optimization for Decision Science

This notebook contains advanced topics for optimization under uncertainty.

Sample Average Approximation¶

Sample average approximation (SAA) method is an approach for solving stochastic optimization problems by Monte Carlo simulation. It approximates the expected objective function of the stochastic problem by a sample average estimate derived from a random sample. The resulting sample average approximation problem is then solved by deterministic optimization techniques.

Consider a stochastic program in the following form:

$z = \inf \int_{\Xi} g(x,\xi) P(d\xi)$

The optimal solution $x^*$ will have optimal value $z^*$ . SAA consider a sample ${\xi^i}$ of independent observations of $\xi$ :

$z^{\nu} = \inf \frac{1}{\nu} \sum_{i=1}^{\nu} g(x,\xi^i)$

Where $x^{\nu}$ is the random vector of solutions with independent random samples, $\xi^i$ , $i=1,...,\nu$ .

Sparse Grids (Worst Case)¶

Instead of Monte Carlo sampling, sparse grids can be used to obtain efficient characterizations of the integrals in stochastic programs. Many computational problems are solved on full grids. While this is feasible if the dimensionality $d$ of the problem is low, 2 or 3, full grids become very expensive when facing a higher dimensional problem $d>4$ . This is due to the curse of dimensionality, which states that the complexity of full grids grows exponentially with $d$ .

Sparse grids defeat the curse of dimensionality and allow the solution of the problems with much smaller effort. It is constructed by extending one-dimensional quadrature rules to higher dimensions. Gaussian quadrature rule is one of these rules to approximate the definite integral of a function. It is usually a weighted sum of function values at specified points within the domain of integration, stated as:

$\int f(x)dx \approx \sum_{i=1}^{n} \omega_i f(x_i)$

Comparison between 2-point Gaussian and trapezoidal quadrature. For the polynomial:

$y(x) = 7x^3-8x^2-3x+3$

whose integral in [-1,1] is $\frac{2}{3}$ .

The trapezoidal rule returns the integral of the orange dashed line:

$y(-1)+y(1) = -10$

The 2-point Gaussian quadrature returns the integral of the black dashed curve, equal to:

$y(-\sqrt{\frac{1}{3}}) + y(\sqrt{\frac{1}{3}}) = \frac{2}{3}$

Exact cubic together with the two-point trapezoidal and Gauss--Legendre interpolants.

The figure is generated from source in this repository. The two-node Gauss--Legendre rule is exact because the integrand is cubic.

import Tasmanian
import numpy as np
import matplotlib.pyplot as plt

# define bounds for integration
UB = 1
LB = -1

# function
y = lambda x: 7 * x**3 - 8 * x**2 - 3 * x + 3

# integrated function
int_y = lambda x: 1.75 * x**4 - 8 / 3 * x**3 - 1.5 * x**2 + 3 * x

# real area
real_S = int_y(UB) - int_y(LB)
print("The area:", real_S)

# area calculated by trapezoidal rule
trap_S = (y(UB) + y(LB)) * (UB - LB) / 2
print("The area calculated by trapezoidal rule:", trap_S)

# generate sparse grids
range_p = np.array([[LB, UB]])
grid_p = Tasmanian.SparseGrid()
grid_p.makeGlobalGrid(1, 0, 1, "level", "gauss-legendre")
grid_p.setDomainTransform(range_p)
points_p = grid_p.getPoints()
weights_p = grid_p.getQuadratureWeights()
# area calculated by sparse grids
gauss_S = sum(y(point) * weights_p[i] for i, point in enumerate(points_p))
print("The area calculated by Gauss rule:", gauss_S)

How to choose quadrature rules for sparse grids depends on how the problem is formed. A guidance on how to choose the one dimensional quadrature rule can be found in: https://tasmanian.ornl.gov/documentation/group__SGEnumerates.html#ga145e27d5ae92acdd5f74149c6d4f2ca2

The following example shows several example quadrature rules:

def show_sparse_grids(range_p, dim=2, output=0, depth=5, rule="gauss-legendre"):
    """
    This function shows the sparse grids generated with different rules
    Arguments:
        range_p: dimension ranges
        dim: sparse grids dimension
        output: output level
        depth: depth level
        rule: quadrature rules
    Return:
        None
    Other:
        A figure shows 2D sparse grids
    """
    grid_p = Tasmanian.SparseGrid()
    grid_p.makeGlobalGrid(dim, output, depth, "level", rule)
    grid_p.setDomainTransform(range_p)
    points_p = grid_p.getPoints()
    weights_p = grid_p.getQuadratureWeights()

    # One scatter call, not one per point: the loop drew every point in a
    # different colour from the cycle, which carried no information and did
    # not survive greyscale printing.
    plt.figure()
    plt.scatter(points_p[:, 0], points_p[:, 1], color="black", marker="o")
    plt.xlabel("Dimension 1")
    plt.ylabel("Dimension 2")
    plt.title("Sparse grids of " + rule)
    plt.show()


range_fix = np.array([[50, 220], [300, 600]])
show_sparse_grids(range_fix)
show_sparse_grids(range_fix, rule="gauss-hermite")
show_sparse_grids(range_fix, rule="chebyshev")

The following conceptual comparison uses the same number of points in two dimensions. The sparse grid is a union of tensor products of nested Clenshaw--Curtis nodes; the Monte Carlo points use a fixed seed. Sparse-grid convergence can be faster for sufficiently smooth integrands, but that conclusion depends on dimension, regularity, and the grid construction rather than on visual coverage alone.

Nested two-dimensional sparse-grid points beside the same number of seeded Monte Carlo samples.

The figure is generated from source in this repository. Reference [2] remains a motivating application of sparse grids to stochastic controller tuning.

Reference¶

Biegler, L.T., 2010. Nonlinear programming: concepts, algorithms, and applications to chemical processes. Society for Industrial and Applied Mathematics.

Birge, J.R. and Louveaux, F., 2011. Introduction to stochastic programming. Springer Science & Business Media.

Related sources:

[1] By Paolostar - Own work, CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=46820806

[2] Renteria, J.A., Cao, Y., Dowling, A.W. and Zavala, V.M., 2018. Optimal pid controller tuning using stochastic programming techniques. AIChE Journal, 64(8), pp.2997-3010.