--- title: "Deciding before the data arrives" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Deciding before the data arrives} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- Many decisions have to be made before all the data is in. A shop orders stock before it knows how much will sell. A contractor plans a day of jobs before knowing how long each will really take. This article shows how to model such a decision with quicopt, solve it, and check whether the answer can be trusted. It assumes you have read `vignette("quicopt")`. The idea in brief: you describe each uncertain quantity by a probability distribution. The service draws many possible outcomes from it, called **scenarios**, and looks for the decision that does best across all of them. You decide what "best across the scenarios" means: best on average, best on the bad days, or good enough on at least 90% of them. ## 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: ``` r 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()`: ``` r 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: ``` r 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"`: ``` r 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: ``` r 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: ``` r 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: ``` r 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: ``` r 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: ``` r 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: ``` r 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: ``` r 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: ``` r 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: ``` r 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: ``` r 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: ``` r 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: ``` r 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: ``` r 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): ``` r 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: ``` r 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: ``` r 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: ``` r 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: ``` r 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. | Summary | What it measures | |---------------------------|-----------------------------------------------------| | `expectation(x)` | the average over the scenarios | | `prob(a <= b)` | the share of scenarios in which the comparison holds | | `cvar(x, 0.95)` | the average over the worst 5% of scenarios | | `scenario_quantile(x, 0.95)` | the value that 95% of scenarios stay at or below | | `scenario_max(x)`, `scenario_min(x)` | the largest and the smallest value | | `variance(x)`, `std_dev(x)` | how much `x` varies from scenario to scenario | 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.