Deciding before the data arrives

How much to stock

A shop orders a product once a day, at 3 per unit, before it knows how many customers will come. Demand varies from day to day, roughly normally, around 100 units with a standard deviation of 15. Every unit a customer asks for that the shop does not have costs it 10 in lost profit and goodwill. Order too little and the shop loses sales; order too much and it pays for stock nobody buys. How much should it order?

This is a classic problem in operations research, known as the newsvendor problem, after a newspaper seller deciding how many papers to buy each morning.

Random variables and scenarios

The model has one decision, the stock, and one uncertain quantity, the demand:

m <- model()
stock  <- num_var(m, "stock", lower = 0, upper = 200)   # decided now
demand <- rand_var(m, "demand", normal(100, 15))       # learned later
set_scenarios(m, 512, seed = 42)

rand_var() declares a random variable. It appears in expressions just like a decision variable, but the service does not choose its value: the value is drawn from the distribution. normal(100, 15) takes the mean and the standard deviation, in the same order as rnorm().

set_scenarios() says how many draws to make: here, 512 possible days, each with its own demand. From now on, every expression that contains demand has 512 values, one per scenario. The seed belongs to the model, so solving the model again uses the same 512 days. The draws are made by the service, which means R’s set.seed() has no effect on them.

Planning for the average day

On a given day, the shortfall is the demand the shop cannot meet: max(demand - stock, 0), which is zero whenever the stock covers demand. The day’s cost is 3 * stock + 10 * shortfall. That cost is a different number in each of the 512 scenarios, but an objective has to be a single number. So the 512 values are summarized first, here by their average, with expectation():

shortfall <- max(demand - stock, 0)          # units short, in each scenario
minimize(m, 3 * stock + 10 * expectation(shortfall))

res <- solve(m)
res$solution[["stock"]]
#> [1] 107.642

The cost of buying, 3 * stock, is the same in every scenario, so only the shortfall needs averaging.

The answer has a simple explanation. Stocking one more unit costs 3, and saves 10 on each day on which demand exceeds the stock. That is worth it as long as demand exceeds the stock on more than 3 days out of 10. So the best stock is the level that demand stays below on 70% of days, the 70th percentile of demand. For a normal distribution with mean 100 and standard deviation 15, that is qnorm(0.7, 100, 15), about 107.9. The service worked from 512 sampled days rather than from the distribution itself, so its answer is close to that value but not equal to it.

An expression that contains a random variable is called random: it has one value per scenario. The objective and the constraints of a model must not be random, so they always contain a summary such as expectation(). quicopt stops with an error if one is missing, and is_random() tells you which kind of expression you have:

is_random(shortfall)
#> [1] TRUE
is_random(expectation(shortfall))
#> [1] FALSE

As mentioned in vignette("quicopt"), a model with random variables is solved by a search, so its status is "heuristic" rather than "optimal":

res$status
#> [1] "heuristic"

Meeting demand on 90% of days

With a stock of about 108, the shop runs out on roughly 3 days in 10. Suppose it wants to do better and meet demand on at least 90% of days. That is a requirement on a probability:

add(m, prob(demand <= stock) >= 0.9)

res <- solve(m)
res$solution[["stock"]]
#> [1] 118.7283

prob(demand <= stock) is the share of scenarios in which demand is at most the stock. The share is a single number, so it can appear in a constraint. A requirement like this is called a chance constraint.

The line holds two comparisons, and they do different jobs. The inner one, demand <= stock, is the event that is checked in each scenario. The outer one, >= 0.9, is the requirement on how often that event happens.

add() changed the model from the previous section, which kept its objective, so solving it again answers the new question. The stock rises to 118.7, the smallest stock that covers demand on 90% of the 512 sampled days. For the normal distribution itself, the answer would be qnorm(0.9, 100, 15), about 119.2.

Checking the answer on new scenarios

The service chose the stock by looking at 512 particular days. On those days the constraint holds, as it should. evaluate() computes any summary at a given solution, so it can show the share of days on which demand is met:

met <- prob(demand <= stock)
evaluate(m, res, met)
#> [1] 0.9003906

But will it hold on days the search never saw? Given a seed, evaluate() draws new scenarios and computes the summary on those instead. A larger sample than the original gives a more precise check:

evaluate(m, res, met, seed = 7, n_scenarios = 1000)
#> [1] 0.896

On new days the stock meets demand only 89.6% of the time. In this example we can also compute the exact answer, because demand is normal:

pnorm(res$solution[["stock"]], 100, 15)
#> [1] 0.8940854

The true probability is 0.894, short of the 90% the stock was chosen to meet. Nothing went wrong in the solve. The service found the smallest stock that met the target on its 512 days, so that stock is just enough for those particular days, and on other days it falls short about half the time. It is the same effect as a statistical model that fits its training data better than new data.

resample() runs this check for every constraint at once. It keeps the solution fixed, draws new scenarios, and reports whether the model’s constraints still hold:

chk <- resample(m, res, seed = 7, n_scenarios = 1000)
chk$feasible
#> [1] FALSE
chk$solver_data$max_violation
#> [1] 0.004

The constraint fails, and max_violation says by how much: the share of days on which demand is met falls short of 0.9 by 0.004.

A safety margin

The share computed from 512 scenarios is an estimate, and its uncertainty can be quantified. Its standard error is sqrt(0.9 * 0.1 / 512), about 0.013. add() takes a margin argument for this: with margin = 2, the service must meet the target plus two standard errors, about 0.927, on its sample. Then the true probability is at least 0.9 with a confidence of about 98%.

To compare the two versions, here is the model as a function with the margin as its argument:

newsvendor <- function(margin = 0) {
  m <- model()
  stock  <- num_var(m, "stock", lower = 0, upper = 200)
  demand <- rand_var(m, "demand", normal(100, 15))
  set_scenarios(m, 512, seed = 42)
  minimize(m, 3 * stock + 10 * expectation(max(demand - stock, 0)))
  add(m, prob(demand <= stock) >= 0.9, margin = margin)
  m
}

res_safe <- solve(newsvendor(margin = 2))
res_safe$solution[["stock"]]
#> [1] 121.036

The safer stock is 121.0, about 2.3 units more. To check it, evaluate it in the original model m, since the question is whether it meets the original target of 90%, on the same 1000 new days as before:

evaluate(m, res_safe, met, seed = 7, n_scenarios = 1000)
#> [1] 0.924
resample(m, res_safe, seed = 7, n_scenarios = 1000)$feasible
#> [1] TRUE
pnorm(res_safe$solution[["stock"]], 100, 15)
#> [1] 0.9196016

This time the target holds on new days, and the exact probability is 0.920.

Several uncertain quantities

In vignette("quicopt"), a contractor chose which of six jobs to accept for a 12-hour day, knowing how long each job takes. In practice the lengths are estimates. Suppose each job takes its estimated time give or take an hour, normally distributed with a standard deviation of 1, and the contractor wants the day to end on time in at least 90% of scenarios:

pay      <- c(9, 6, 7, 4, 3, 2)     # what each job pays
estimate <- c(6, 5, 4, 3, 3, 2)     # how long each job is expected to take

m <- model()
take  <- bin_var(m, "take", length = 6)
hours <- rand_var(m, "hours", normal(estimate, 1))
set_scenarios(m, 500, seed = 1)

on_time <- prob(sum(hours * take) <= 12)
maximize(m, sum(pay * take))
add(m, on_time >= 0.9)

res <- solve(m)
which(res$solution == 1)
#> take[1] take[3] 
#>       1       3
res$objective
#> [1] 16

Because estimate is a vector of six means, normal(estimate, 1) declares six random variables at once, one per job, each drawn independently. sum(hours * take) is the length of the day, which differs from scenario to scenario, and on_time is the share of scenarios in which it fits into 12 hours.

The answer is not the plan from vignette("quicopt"), jobs 1, 3 and 6. That plan was expected to take exactly 12 hours, so it runs late whenever the jobs together take longer than estimated, which happens on about half the days. evaluate() accepts a solution written out by hand too, as a named vector:

old_plan <- c("take[1]" = 1, "take[2]" = 0, "take[3]" = 1,
              "take[4]" = 0, "take[5]" = 0, "take[6]" = 1)
evaluate(m, old_plan, on_time)
#> [1] 0.522

Dropping job 6 gives up a pay of 2, and the day then ends on time in more than 90% of the scenarios, both those of the model and new ones:

evaluate(m, res, on_time)
#> [1] 0.922
evaluate(m, res, on_time, seed = 2)
#> [1] 0.92
pnorm(12, mean = 6 + 4, sd = sqrt(2))     # the exact probability
#> [1] 0.9213504

Unlike the stock, which could be set to any number, a set of jobs either fits or does not. There is usually no plan that meets 90% exactly, and the best one keeps some room to spare: here its true probability is 0.921.

Planning for bad days

An average says little about the bad days. A shop that can absorb an ordinary day but not a disastrous one may care more about its worst scenarios. cvar(cost, 0.95) is the average cost over the worst 5% of scenarios, a summary known in finance as the conditional value at risk. Minimizing it asks for the stock that keeps those bad days as cheap as possible:

m <- model()
stock  <- num_var(m, "stock", lower = 0, upper = 200)
demand <- rand_var(m, "demand", normal(100, 15))
set_scenarios(m, 512, seed = 42)
cost <- 3 * stock + 10 * max(demand - stock, 0)     # the day's cost, in each scenario

minimize(m, cvar(cost, 0.95))
solve(m)$solution[["stock"]]
#> [1] 132.5795

The stock is higher than the 107.6 that minimizes the average cost, because the worst days are the days of high demand, and more stock is what makes them cheaper.

The most cautious summary is the single worst scenario. scenario_max(cost) is the highest cost among the 512 scenarios. Calling minimize() again replaces the objective:

minimize(m, scenario_max(cost))
solve(m)$solution[["stock"]]
#> [1] 143.4227

This stocks for the highest demand among the 512 sampled days. Be careful with it: a single scenario decides the answer, so it depends on the sample much more than an average does. A larger sample will usually contain a higher demand, and the stock will grow with it.

Between these two extremes are scenario_quantile(cost, 0.95), the cost that 95% of scenarios stay at or below, and variance() and std_dev(), which measure how much the cost varies from scenario to scenario. Each of them can be used in the objective or in a constraint. The table at the end of this article lists all the summaries.

Putting a price on an event

prob() counts how often something happens. Sometimes you want to use the event inside the cost instead. holds(demand > stock) is 1 in the scenarios in which the shop runs out and 0 in the others. Suppose that running out costs an extra 200 each time it happens, for an emergency delivery, on top of the 10 per missing unit:

m <- model()
stock  <- num_var(m, "stock", lower = 0, upper = 200)
demand <- rand_var(m, "demand", normal(100, 15))
set_scenarios(m, 512, seed = 42)

shortfall <- max(demand - stock, 0)
stockout  <- holds(demand > stock)          # 1 when the shop runs out, else 0
minimize(m, 3 * stock + 10 * expectation(shortfall) + 200 * expectation(stockout))
solve(m)$solution[["stock"]]
#> [1] 121.0856

The fixed cost of running out does what the 90% requirement did, but through the price instead of a rule: the stock rises until running out is rare enough not to be worth avoiding further. expectation(holds(x)) is the same number as prob(x); holds() is useful because you can multiply the event by a cost first.

Uncertainty from data

Often the uncertainty is not a distribution you assume but data you have observed. set_empirical() turns the columns of a data frame into random variables: each column becomes one, named after the column, and each row becomes one scenario. Because all the columns are read row by row, values observed on the same day stay together, and any correlation between the columns carries over into the model.

Here are 200 days of demand and price, simulated in R so that the example is self-contained (so set.seed() matters here: it fixes this data, not the service’s draws):

set.seed(1)
z <- rnorm(200)
history <- data.frame(demand = round(100 + 15 * z + 5 * rnorm(200)),
                      price  = round(12 - 0.03 * (15 * z) + rnorm(200), 2))
cor(history$demand, history$price)
#> [1] -0.2968849

On days of high demand the price tends to be lower. The shop sells min(demand, stock) units at the day’s price and pays 3 per unit stocked:

m <- model()
stock <- num_var(m, "stock", lower = 0, upper = 200)
set_empirical(m, history)
maximize(m, expectation(m$price * min(m$demand, stock)) - 3 * stock)

res <- solve(m)
res$solution[["stock"]]
#> [1] 109

set_empirical() declares the random variables without returning them, so they are retrieved from the model by name, as m$price and m$demand.

No distribution had to be chosen or fitted: the 200 observed days are the scenarios.

Other distributions

Besides normal(), the service draws from three distributions, each with the same parameters as R’s corresponding random number function: uniform(min, max) as runif(), exponential(rate) as rexp(), and bernoulli(prob), which is 1 with probability prob and 0 otherwise. Here is the shop again, with a demand that is exponential with mean 100:

m <- model()
stock  <- num_var(m, "stock", lower = 0, upper = 600)
demand <- rand_var(m, "demand", exponential(1 / 100))
set_scenarios(m, 512, seed = 42)
minimize(m, 3 * stock + 10 * expectation(max(demand - stock, 0)))
solve(m)$solution[["stock"]]
#> [1] 120.7402

The reasoning from “Planning for the average day” applies to any distribution: the best stock is the 70th percentile of demand, which here is qexp(0.7, 1 / 100), about 120.4.

A parameter can also depend on a decision. bernoulli(0.2 - 0.015 * spend) is a machine breakdown that becomes less likely the more is spent on maintenance, so the model can weigh the cost of maintenance against the cost of breakdowns.

Any other distribution can be used by drawing a sample in R and passing it to empirical(). Here, a demand with a lognormal distribution:

set.seed(2)
m <- model()
stock  <- num_var(m, "stock", lower = 0, upper = 400)
demand <- rand_var(m, "demand", empirical(rlnorm(512, log(100), 0.3)))
minimize(m, 3 * stock + 10 * expectation(max(demand - stock, 0)))
solve(m)$solution[["stock"]]
#> [1] 121.5136

The sample’s length sets the number of scenarios, so set_scenarios() is not needed. Such a sample is fixed data, which has three consequences. set.seed() determines it, not the model’s seed. resample() cannot draw new values for it, so to check a solution on new data, draw a second sample in R. And its parameters cannot depend on a decision, so a distribution like the maintenance example above needs one of the four built-in distributions.

Solving again from an earlier answer

A model is often solved several times, with a constraint tightened, a cost updated, or more scenarios. set_start() gives the previous answer to the next solve as its starting point, so the search does not start from scratch:

m <- newsvendor()
res <- solve(m)                          # 512 scenarios
set_scenarios(m, 1000, seed = 42)        # the same model, more scenarios
set_start(m, res)                        # start where the last solve ended
solve(m)$solution[["stock"]]
#> [1] 117.6806

Summaries at a glance

Each of these turns a random expression x into a single number that can be used in an objective or a constraint.

The reported objective of a model with random variables is computed on the scenarios the solution was found on, and so tends to look better than the solution will do on new data. evaluate() and resample() give the figures on new scenarios.