Sunday, November 6, 2022

A Bicriterion IP Model for a Game Strategy

A question on Mathematics Stack Exchange asks how to solve an optimization problem related to a game apparently called "speedrunner". The problem boils down to selecting integer values for variables $n_a,$ $n_c$ and $n_{b,i}$ ($i=1,\dots,n_c$) so as to minimize elapsed time $$t_{\textrm{total}} =t_{a}n_{a}+\sum_{i=1}^{n_{c}}t_{b}n_{b,i}+t_{c}n_{c},$$

where  $t_a$, $t_b$ and $t_c$ are parameters with specified values. The sole constraint specified is that winnings, defined by

$$x_{\textrm{total}} = x_{0}n_{a}\prod_{i=1}^{n_{c}}n_{b,i}$$(where $x_0$ is another parameter), must be at least $\$1$billion. We can infer from the requirement of positive earnings that$n_a \ge 1,n_c \ge 1,$and$n_{b,i} \ge 1$for$i=1,\dots,n_c.$For modeling purposes, we will allow the index$i$to exceed$n_c$and require that$n_{b,i}=0$for$i \gt n_c.$An integer linear programming model would be a no-brainer were it not for two things: the product form of the expression for winnings is nonlinear; and it involves the product of a variable number of variables. The approach I chose was to express the original variables in terms of binary variables. To do that, we first need upper bounds on the original variables. We get those by guessing an upper bound$T$on$t_{\textrm{total}}.$(If we guess too low, we will either get a solution that exceeds that upper bound or the model will become infeasible, either of which will tip us off to try a larger guess.) Once we have that upper bound, we get an upper bound for each variable by setting the other two variables as small as possible. That leads to the following upper bounds:$$N_{a} =\left\lfloor \frac{T-t_{b}-t_{c}}{t_{a}}\right\rfloor$$ $$N_{b} =\left\lfloor \frac{T-t_{a}-t_{c}}{t_{b}}\right\rfloor$$ and $$N_{c} =\left\lfloor \frac{T-t_{a}}{t_{c}+t_{b}}\right\rfloor .$$ The denominator in the third equation arises from the fact that adding 1 to$n_c$not only directly costs$t_c$units of time but also adds another$n_{b,i}$variable with minimum value 1, costing$t_b$units of time. Armed with those upper bounds, we can introduce binary variables$y_j(j=1,\dots,N_a),w_j(j=1,\dots, N_c)$and$z_{i,j}(i=1,\dots,N_c$and$j=0,\dots,N_b)$along with constraints making them type 1 special ordered sets (i.e., each set contains exactly one variable with value 1) and expanding the original variables in terms of them. Specifically: •$n_{a}=\sum_{j=1}^{N_{a}}j\cdot y_{j}$with constraint$\sum_{j=1}^{N_{a}}y_{j}=1;$•$n_{c}=\sum_{j=1}^{N_{c}}j\cdot w_{j}$with constraint$\sum_{j=1}^{N_{c}}w_{j}=1;$and •$n_{b,i}=\sum_{j=0}^{N_{b}}j\cdot z_{i,j}$for$i=1,\dots,N_{c}$with the following constraints for all$i\in\left\{ 1,\dots,N_{c}\right\}:$•$\sum_{j=0}^{N_{b}}z_{i,j}=1$(the value of$n_{b,i}$is uniquely defined); and •$z_{i,0}+\sum_{k=i}^{N_{c}}w_{k}=1$($n_{b,i}=0$if and only if$n_{c}<i$). Armed with all that, the original objective function can be expressed as minimizing$$t_{a}n_{a}+t_{b}\sum_{i=1}^{N_{c}}n_{b,i}+t_{c}n_{c},$$where the only tweak is that we now sum over a constant number$N_c$of terms rather than a variable number$n_c,$relying on the fact that$n_{b,i}=0$for$i>n_c.$That brings us to the nonlinear formula for winnings. The original constraint is$ x_{\textrm{total}} \ge 10^9,$which we can rewrite as$\log(x_{\textrm{total}}) \ge \log(10^9)$using whatever flavor logarithm you like. Expanding the logs of$n_a$and$n_{b,i}$in terms of the binary variables, we get$$\log(x_{\textrm{total}}) = \log(x_{0})+\sum_{j=1}^{N_{a}}\log(j)\cdot y_{j} +\sum_{i=1}^{N_{c}}\sum_{j=1}^{N_{b}}\log(j)\cdot z_{i,j}.$$So the constraint that log winnings equal or exceed the log of$10^9$is linear in the binary variables. The solution in the Math Stack Exchange post$(n_a = 7,n_c = 3$and$n_{b,1}=n_{b,2}=n_{b,3}=18)$turns out to be optimal given the parameter values in the question, and the IP model here will correctly achieve the minimum time (approximately 52.3 minutes) ... but it will not necessarily reproduce the best winnings within that time (approximately$\$1.02$ billion). That is because there are multiple feasible solutions that hit the correct time value, with $n_a = 7$ and with $n_c = 3$ but with the 54 cumulative $n_b$ units allocated differently (for instance, $\{19, 19, 16\}$ rather than $\{18, 18, 18\}.$ To get the best possible solution, we can use a bicriterion model with lexicographically ordered objectives, first minimizing time and then maximizing winnings (by minimizing the negative of the log of winnings).

I tested that model using the Java API to CPLEX 22.1. My code is available from my university repository.

Tuesday, November 1, 2022

Holt Double Exponential Smoothing

Given a time series $x_t,\, t=1,2,\dots,$ presumed to contain linear trend, Holt's double exponential smoothing method computes estimates $s_t$ of the level (mean) at time $t$ and $b_t$ of the (presumed constant) slope using the following formulas: $$s_t = \alpha x_t + (1-\alpha)(s_{t-1} + b_{t-1})$$ and $$b_t = \beta(s_t - s_{t-1}) + (1-\beta)b_{t-1}$$ where $\alpha,\beta \in (0,1)$ are smoothing weights. A recent question on OR Stack Exchange asked why the second formula is based on the level estimate and not the observed value. In other words, the proposed alternative to the trend update was $$b_t = \beta(x_t - x_{t-1}) + (1-\beta)b_{t-1}.$$

The intuition for doing it Holt's way is fairly simple. If exponential smoothing is working as intended (meaning smoothing things), then the difference in level estimates $s_t - s_{t-1}$ should be less variable than the difference in observed values $x_t - x_{t-1}.$ A formal proof probably involves induction arguments, requiring more functioning brain cells than I had available at the time, so I was a bit loose mathematically in my answer on OR SE. Just to confirm the intuition, I did some Monte Carlo simulations in R. The notebook containing the experimental setup, including code, is available here. It requires the dplyr and ggplot2 library packages.

The following plots show confidence intervals over time for the errors in the level and trend estimates using both Holt's formula and what I called the "variant" method. They are from a single experiment (100 independent time series with identical slope and intercept, smoothed both ways), but other trials with different random number seeds and changes to the variability of the noise produced similar results.

In both cases, the estimates start out a bit wobbly (and the Holt estimates may actually be a bit noisier), but over time both stabilize. There does not seem to be much difference between the two approaches in how noisy the level estimates are, at least in this run. The Holt estimates may have slightly narrower confidence intervals, but that is not clear, and the difference if any seems pretty small. The Holt trend estimates, however, are considerably less noisy than those of the variant method, supporting the intuitive argument.