Lecture 13: Ranking and Selection - Simulation and Optimization: Lecture Notes

Comparing Scenarios chose between two solutions. This notebook covers the second of the four types of problems: a small, discrete set of solutions. We discuss two ways to spend the simulation budget, an equal split over all solutions and a smarter one that discards bad solutions early, and apply both to the newsvendor problem.

Learning outcomes

On completion of this notebook, you will be able to:

formulate the newsvendor problem as a simulation optimization problem
apply ranking and selection with an equal split of the simulation budget, and explain its drawbacks
formulate the one-sided tests that discard solutions, and explain the Šidák correction for multiple testing
apply ranking and selection with discarding, and read and adapt the Python code that does so

The Setting¶

We solve the simulation optimization problem

\max_{\pi \in S} \ E[r(X, \pi)],

(1)

where $E[r(X, \pi)]$ can only be evaluated by simulation. In this notebook $S$ is discrete and relatively small, with $k = |S|$ solutions: the simulation budget of $m$ runs is much larger than the number of solutions ( $m \gg k$ ), so it suffices to simulate every solution at least a couple of times. The task is to select the best solution, and on the way we rank the solutions by their estimated performance, hence the name ranking and selection.

The Newsvendor Problem¶

The newsvendor problem is a typical ranking and selection problem. A newsvendor buys $\pi$ newspapers in the morning for a price $c$ each, and sells them during the day for a price $p > c$ . Demand $X$ is random, and newspapers that are left at the end of the day are worthless. The newsvendor sells $\min(X, \pi)$ papers, so the profit is

r(X, \pi) = p \min(X, \pi) - c\pi,

(2)

and the newsvendor chooses $\pi$ from $S = \{0, 1, \dots, K\}$ . The upper bound $K$ comes from the newsvendor’s stand: it has room for at most $K$ newspapers, so ordering more is not possible. The example in Comparing Scenarios was a newsvendor problem with only two order sizes, 9 and 10.

In the code below, simulate simulates the profit of one order size n times:

import numpy as np

price, cost = 1, 0.75
mean_demand = 10
orders = np.arange(21)  # S = {0, 1, ..., 20}


def profit(order, demand):
    return price * np.minimum(demand, order) - cost * order


def simulate(order, n, rng):
    """Simulated profits of n days with the given order size."""
    return profit(order, rng.poisson(mean_demand, n))

For this simple example, the expected profit of every order size can also be calculated exactly. How to do this is outside the scope of this course, and for real-life simulation optimization problems it is in general not possible. We use the exact values, stored in true_value by the hidden cell below, only to check how well the methods do.

# exact expected profit, only to check the methods (outside the scope of
# the course)
from scipy import stats

expected_sales = [
    stats.poisson.sf(np.arange(o), mean_demand).sum() for o in orders
]
true_value = price * np.array(expected_sales) - cost * orders
best_order = orders[true_value.argmax()]
print(f"best order size: {best_order}")
print(f"expected profit of order size {best_order}: {true_value.max():.3f}")
print(f"expected profit of order size 10: {true_value[10]:.3f}")

best order size: 8
expected profit of order size 8: 1.540
expected profit of order size 10: 1.249

Ordering 8 papers is best, not 10: with leftovers worthless and a small margin of €0.25 per sold paper, the newsvendor should order less than the mean demand.

Option 1: Split the Budget Equally¶

The natural first approach is to simulate every solution $m/k$ times, estimate $E[r(X, \pi)]$ for each $\pi$ by the average of its runs, and return the solution with the highest average. Let $y(\pi)$ and $s^2(\pi)$ denote the sample average and sample variance computed from the runs of $\pi$ (the $\bar X$ and $S_X^2$ of Comparing Scenarios). With $s^2(\pi)$ we can also give a 95% CI for each $E[r(X, \pi)]$ .

For the newsvendor, a budget of $m = 4200$ runs gives 200 runs per order size:

rng = np.random.default_rng(1)
budget = 4200
runs_per_order = budget // len(orders)

means = np.zeros(len(orders))
sds = np.zeros(len(orders))
for idx, order in enumerate(orders):
    profits = simulate(order, runs_per_order, rng)
    means[idx] = profits.mean()
    sds[idx] = profits.std(ddof=1)
half_widths = 1.96 * sds / np.sqrt(runs_per_order)

best = means.argmax()
print(f"selected order size: {orders[best]}")
print(f"estimated expected profit: {means[best]:.3f}")

selected order size: 8
estimated expected profit: 1.600

import plotly.graph_objects as go

TEXT_COLOR = "#111827"
PLOT_LAYOUT = {
    "paper_bgcolor": "rgba(0,0,0,0)",
    "plot_bgcolor": "rgba(0,0,0,0)",
    "font": {"color": TEXT_COLOR, "size": 13},
    "margin": {"l": 50, "r": 20, "t": 20, "b": 50},
    "height": 360,
    "legend": {
        "bgcolor": "rgba(0,0,0,0)",
        "orientation": "h",
        "yanchor": "top",
        "y": -0.2,
    },
}
AXIS_STYLE = {"gridcolor": "rgba(128,128,128,0.25)"}
# static pictures: no hover, zoom or drag
PLOT_CONFIG = {"displayModeBar": False, "staticPlot": True}


def add_cis(fig, x, y, half_width, color, name):
    fig.add_trace(
        go.Scatter(
            x=x,
            y=y,
            mode="markers",
            marker={"color": color, "size": 7},
            error_y={
                "type": "data",
                "array": half_width,
                "color": color,
                "thickness": 2,
                "width": 4,
            },
            name=name,
        )
    )


def finish_ci_plot(fig):
    fig.update_layout(**PLOT_LAYOUT)
    fig.update_xaxes(title="order size π", dtick=2, **AXIS_STYLE)
    fig.update_yaxes(title="expected profit (€)", **AXIS_STYLE)
    fig.show(config=PLOT_CONFIG)

option1_fig = go.Figure()
add_cis(option1_fig, orders, means, half_widths, "#1f77b4", "95% CI")
option1_fig.add_trace(
    go.Scatter(
        x=orders,
        y=true_value,
        mode="lines",
        line={"color": "grey", "dash": "dash", "width": 2},
        name="true expected profit",
    )
)
finish_ci_plot(option1_fig)

Loading...

Figure 1:95% CIs for the expected profit of every order size, each based on 200 runs, with the true expected profit.

Figure 1 shows the CIs. The CIs for small order sizes are very narrow (see Exercise 1). Around the optimum, the CIs of order sizes 7, 8 and 9 overlap: the budget is not enough to tell them apart. At the same time, a large part of the budget went to order sizes 14 to 20, whose CIs lie far below the others after a few runs already.

Because of the noise, option 1 shows both problems mentioned in About Simulation Optimization: it may select the wrong order size, and the estimated profit of the selected order size tends to be too high. To see how often this happens, we repeat option 1 1000 times:

def option_1(budget, rng):
    """Equal split: return the selected order size and its estimate."""
    n = budget // len(orders)
    estimates = [simulate(order, n, rng).mean() for order in orders]
    best = int(np.argmax(estimates))
    return orders[best], estimates[best]


rng = np.random.default_rng(2)
results_1 = [option_1(budget, rng) for _ in range(1000)]
selected_1 = np.array([order for order, _ in results_1])
overestimate = [est - true_value[order] for order, est in results_1]
print(f"fraction best order size selected: {np.mean(selected_1 == 8):.2f}")
print(f"average overestimation of the profit: {np.mean(overestimate):.3f}")

fraction best order size selected: 0.54
average overestimation of the profit: 0.056

Option 1 selects the best order size 8 only about half of the time; otherwise it mostly selects 7 or 9. And the estimate of the selected order size’s expected profit is too high on average, because we select the order size whose estimate happened to come out high.

Common Random Numbers¶

The code above draws new demands for every order size. As in Comparing Scenarios, we can also use common random numbers: simulate one set of demands and evaluate every order size on the same demands. This does not make the individual CIs narrower, but it reduces the variance of the differences between order sizes, and these differences decide which order size is selected. In a discrete-event simulation it is the same: reuse the same customer arrival moments for every candidate solution.

Option 2: Discard Bad Solutions Early¶

Option 1 spends as many runs on clearly bad solutions as on promising ones. Option 2 first simulates all solutions a limited number of times, discards the solutions that are significantly worse than another solution, and spends the rest of the budget on the remaining ones. In general terms, with $k = |S|$ solutions and a budget of $m$ runs:

Simulate every solution $m_0$ times, where $m_0 < \lfloor m/k \rfloor$ is the number of runs per solution in this first round ( $\lfloor \cdot \rfloor$ rounds down to an integer). Compute the average $y(\pi)$ and the sample variance $s^2(\pi)$ of each solution.
For every pair of solutions $(\pi, \pi')$ , test whether $\pi$ is significantly worse than $\pi'$ .
Discard every solution that is significantly worse than at least one other solution. The solutions that remain form the set $I$ .
Divide the remaining budget $m - k m_0$ evenly over the solutions in $I$ , and return the solution in $I$ with the highest average.

When to Discard a Solution¶

Consider first one pair of solutions $\pi$ and $\pi'$ , each simulated $m_0$ times, independently. We maximize, so $\pi$ is worse than $\pi'$ if $E[r(X, \pi)] < E[r(X, \pi')]$ . We only discard $\pi$ if the data give strong evidence for this, so this is the alternative hypothesis of a one-sided test:

\begin{aligned} H_0(\pi, \pi') &: E[r(X, \pi)] \ge E[r(X, \pi')] \quad (\pi \text{ is not worse than } \pi'), \\ H_1(\pi, \pi') &: E[r(X, \pi)] < E[r(X, \pi')] \quad (\pi \text{ is worse than } \pi'). \end{aligned}

(4)

The difference of the averages $y(\pi') - y(\pi)$ has variance $(\sigma^2(\pi) + \sigma^2(\pi'))/m_0$ , which we estimate with the sample variances. So the test statistic

T = \frac{y(\pi') - y(\pi)}{\sqrt{s^2(\pi) + s^2(\pi')}/\sqrt{m_0}}

(5)

has approximately a $t$ -distribution with $m_0 - 1$ degrees of freedom if the two expectations are equal, and a large value of $T$ is evidence for $H_1$ . As on the slides, let $\beta_q$ be the inverse of this $t$ -distribution at $q$ : the number $b$ with $P(T \le b) = q$ . In Python, $\beta_q$ is stats.t.ppf(q, m0 - 1). A one-sided test at significance level $\alpha$ rejects $H_0(\pi, \pi')$ if $T \ge \beta_{1-\alpha}$ , that is, if

y(\pi) \le y(\pi') - \beta_{1-\alpha}\, \frac{\sqrt{s^2(\pi) + s^2(\pi')}}{\sqrt{m_0}}.

(6)

For $\alpha = 0.05$ and $m_0 = 100$ , $\beta_{0.95} = 1.66$ . For large $m_0$ , the $t$ -distribution is close to the standard normal distribution, so the normal critical values are a good approximation (1.64 here); the exercises use them for calculations by hand.

The Multiple-Testing Problem and the Šidák Correction¶

With $k$ solutions, we discard $\pi$ as soon as one of the $k - 1$ tests against the other solutions rejects $H_0(\pi, \pi')$ . Even if $\pi$ is the best solution, each of these tests can reject by chance, and the more tests, the larger the probability that at least one of them does. Compare it with a multiple-choice exam that you take together with $k - 1$ monkeys that answer at random: the more monkeys, the larger the chance that one of them beats you by luck, and we should not conclude from that that this lucky monkey knows more about the exam’s topic than you do. In the same way, a solution that beats $\pi$ in one of many tests by luck is not necessarily better than $\pi$ . So if every test is done at level $\alpha = 0.05$ , the probability of discarding the best solution is much larger than 5%.

We want the probability of discarding the best solution to be at most $\alpha$ . Therefore we do each test at a stricter, “conservative” level $\alpha^*$ . If $\pi$ is the best solution, all $k - 1$ hypotheses $H_0(\pi, \pi')$ are true, and each test does not reject with probability (at least) $1 - \alpha^*$ . The $k - 1$ tests are not independent: they all use the same average $y(\pi)$ , so if $y(\pi)$ happens to be low, $\pi$ looks bad in all comparisons at once. This dependence is positive (the tests tend to reject together), and as a result the probability that none of the $k - 1$ tests rejects is at least what it would be for independent tests:

P(\text{the best solution is not discarded}) \ge (1 - \alpha^*)^{k-1}.

(7)

The Šidák correction chooses $\alpha^*$ such that this lower bound equals $1 - \alpha$ :

(1 - \alpha^*)^{k-1} = 1 - \alpha \quad\Longrightarrow\quad \alpha^* = 1 - \sqrt[k-1]{1-\alpha}.

(8)

Then the best solution survives with probability at least $1 - \alpha$ .

You may know the Bonferroni correction from the statistics part of the course, which takes $\alpha^* = \alpha/(k-1)$ . It also guarantees that the best solution survives with probability at least $1 - \alpha$ , but the Šidák level is always a bit larger: $1 - \sqrt[k-1]{1-\alpha} > \alpha/(k-1)$ . So with the Šidák correction each test is slightly less strict and discards bad solutions a bit more often, with the same guarantee for the best solution. Here the difference is small: for $k = 21$ and $\alpha = 0.05$ , Bonferroni gives $\alpha^* = 0.0025$ and Šidák $\alpha^* \approx 0.0026$ .

For the newsvendor, with $k = 21$ and $\alpha = 0.05$ , this gives $\alpha^* \approx 0.0026$ and, with $m_0 = 100$ , $\beta_{1-\alpha^*} \approx 2.86$ instead of $\beta_{0.95} = 1.66$ : a much stricter bar for declaring a solution worse, because it is tested against many other solutions at once. The set of solutions that survive the first round is

I = \left\{ \pi \in S \;\middle|\; y(\pi) > y(\pi') - \beta_{1-\alpha^*} \frac{\sqrt{s^2(\pi)+s^2(\pi')}}{\sqrt{m_0}} \text{ for all } \pi' \ne \pi \right\}.

(9)

Option 2 for the Newsvendor¶

We use the same budget $m = 4200$ as for option 1, with $m_0 = 100$ runs per order size in the first round. That leaves $4200 - 21 \times 100 = 2100$ runs for the second round. Steps 1 to 3:

from scipy import stats

rng = np.random.default_rng(3)
budget = 4200
m0 = 100
k = len(orders)
alpha = 0.05
alpha_star = 1 - (1 - alpha) ** (1 / (k - 1))
beta = stats.t.ppf(1 - alpha_star, m0 - 1)
print(f"alpha* = {alpha_star:.4f}, beta = {beta:.2f}")

# step 1: simulate every order size m0 times
first_round = {order: simulate(order, m0, rng) for order in orders}
y = {order: runs.mean() for order, runs in first_round.items()}
s2 = {order: runs.var(ddof=1) for order, runs in first_round.items()}

# steps 2 and 3: discard every order size that is worse than another one
survivors = []
for order in orders:
    discard = False
    for other in orders:
        if other != order:
            margin = beta * np.sqrt(s2[order] + s2[other]) / np.sqrt(m0)
            if y[order] <= y[other] - margin:
                discard = True
    if not discard:
        survivors.append(order)
print("I =", [int(order) for order in survivors])

alpha* = 0.0026, beta = 2.86
I = [7, 8, 9, 10, 11]

first_half = {o: 1.96 * np.sqrt(s2[o] / m0) for o in orders}
option2_fig = go.Figure()
for kept, color, name in [
    (True, "#1f77b4", "kept (in I)"),
    (False, "grey", "discarded"),
]:
    shown = [o for o in orders if (o in survivors) == kept]
    add_cis(
        option2_fig,
        shown,
        [y[o] for o in shown],
        [first_half[o] for o in shown],
        color,
        name,
    )
finish_ci_plot(option2_fig)

Loading...

Figure 2:95% CIs for the expected profit of every order size after the first round of option 2 (100 runs each). The grey order sizes are discarded.

Figure 2 shows which order sizes survive. Small order sizes are discarded: because their CIs are so narrow, the data show clearly that they are worse than the best order sizes. In step 4 we divide the remaining 2100 runs over the survivors, add them to the runs of the first round, and select the order size with the highest average:

n_each = (budget - k * m0) // len(survivors)
final = []
for order in survivors:
    extra_runs = simulate(order, n_each, rng)
    all_runs = np.concatenate([first_round[order], extra_runs])
    final.append(all_runs.mean())
best = int(np.argmax(final))
print(f"{n_each} extra runs per surviving order size")
print(f"selected order size: {survivors[best]}, estimate: {final[best]:.3f}")

420 extra runs per surviving order size
selected order size: 8, estimate: 1.498

The five survivors now get 520 runs each instead of 200. To compare the two options fairly, we repeat both 1000 times with the same budget, with option 2 wrapped in a function option_2 that follows the code above line by line:

def option_2(budget, m0, alpha, rng):
    """Discard bad order sizes first: return the selected order size."""
    k = len(orders)
    alpha_star = 1 - (1 - alpha) ** (1 / (k - 1))
    beta = stats.t.ppf(1 - alpha_star, m0 - 1)
    first_round = {order: simulate(order, m0, rng) for order in orders}
    y = {order: runs.mean() for order, runs in first_round.items()}
    s2 = {order: runs.var(ddof=1) for order, runs in first_round.items()}
    survivors = []
    for order in orders:
        discard = False
        for other in orders:
            if other != order:
                margin = beta * np.sqrt(s2[order] + s2[other]) / np.sqrt(m0)
                if y[order] <= y[other] - margin:
                    discard = True
        if not discard:
            survivors.append(order)
    n_each = (budget - k * m0) // len(survivors)
    final = []
    for order in survivors:
        extra_runs = simulate(order, n_each, rng)
        all_runs = np.concatenate([first_round[order], extra_runs])
        final.append(all_runs.mean())
    return survivors[int(np.argmax(final))]

rng = np.random.default_rng(4)
selected_2 = np.array([option_2(budget, m0, alpha, rng) for _ in range(1000)])
print(f"option 1: best order size selected {np.mean(selected_1 == 8):.2f}")
print(f"option 2: best order size selected {np.mean(selected_2 == 8):.2f}")

option 1: best order size selected 0.54
option 2: best order size selected 0.66

With the same budget, option 2 finds the best order size more often than option 1, because it spends more runs where it matters: on the order sizes that are hard to tell apart. It does not remove the problems of noise completely: the expected profits of order sizes 7 and 9 are only €0.03 and €0.08 below that of order size 8, and telling them apart reliably needs many more runs.

References¶

Koole, G. (2019). An Introduction to Business Analytics. §8.2.
Šidák, Z. (1967). “Rectangular confidence regions for the means of multivariate normal distributions.” Journal of the American Statistical Association, 62:626–633.