Backward-Recurrent Capacity Provisioning vs. Lazy Capacity Provisioning. LCP stays within the bounds “lazily” moving forwards in time. The optimal solution stays within the bounds moving backwards in time.

<cit.> If there exists an index $t \in [1, \tau-1]$ such that $X_{\tau,t+1}^U < X_{\tau,t}^U$ or $X_{\tau,t+1}^L > X_{\tau,t}^L$, then $(\hat{X}_{\tau,1},\dots,\hat{X}_{\tau,t}) := (X_{\tau,1}^L,\dots,X_{\tau,t}^L) = (X_{\tau,1}^U,\dots,X_{\tau,t}^U)$, and no matter what the future arrival is, solving the optimization in $[1,\tau']$ for $\tau' > \tau$ is equivalent to solving two optimizations: one over $[1,t]$ with initial condition $X_0$ and final condition $\hat{X}_{\tau,t}$ and the second over $[t+1,\tau']$ with initial condition $\hat{X}_{\tau,t}$.

While not changing the worst-case complexity, this significantly improves the practical complexity in the application of right-sizing data centers as diurnal load patterns typically ensure that less than a day needs to be considered <cit.>. We denote by $X_{\tau}^{L,(t,x_0)}$ and $X_{\tau}^{U,(t,x_0)}$ the bounds resulting from optimizations beginning at time slot $t$ with initial condition $x_0$.  <cit.> showed that lazy capacity provisioning is $3$-competitive and also proved that this result is tight.

Lazy Capacity Provisioning <cit.>
$\mathcal{I}_{\text{SSCO}} = (\tau \in \mathbb{N}, m \in \mathbb{N}, \beta \in \mathbb{R}_{>0}, (f_1, \dots, f_{\tau}) \in (\mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0})^{\tau})$
$t_0 \gets 0$
$x_0 \gets 0$
$t \gets \tau-1$ $2$
$X_{t,t}^U < X_{t,t-1}^U \lor X_{t,t}^L > X_{t,t-1}^L$
$t_0 \gets t$
$x_0 \gets X_{t,t-1}^U$

find $X_{\tau,\tau}^{L,(t_0,x_0)}$ using the optimization described by <ref>
find $X_{\tau,\tau}^{U,(t_0,x_0)}$ using the optimization described by <ref>

§.§.§ Integral Algorithm

<cit.> applied lazy capacity provisioning to the integral variant Int-SSCO using their graph-based offline algorithm discussed in <ref> to compute the integral lower and upper bounds. It is apparent that this immediately yields a deterministic online algorithm for Int-SSCO.  <cit.> showed that similar to lazy capacity provisioning, their algorithm is $3$-competitive. Due to the changed method of determining the bounds, its runtime is $\mathcal{O}(\tau^2 C \log_2 m)$. Note that it is impossible to cache the intermediate results of the dynamic program (see <ref>) as the binary search over possible configurations considers different vertices depending on the obtained schedule, which changes over time. Thus, for large $\tau$, it may be beneficial to use caching instead of binary search resulting in a worst-case runtime of $\mathcal{O}(\tau C m)$. By using the same method of shortening the used history that was proposed by  <cit.>, we can reduce this time complexity drastically in practice (for large $\tau$). Thus, the adopted algorithm is still described by <ref>. We simply need to slightly modify the graph-based algorithm computing optimal offline solutions to allow for initial conditions other than $0$.

§.§ Memoryless Algorithm

<cit.> showed that for SSCO, a competitive ratio of $3$ can also be attained by a memoryless algorithm. In a memoryless online algorithm for smoothed convex optimization, the configuration $X_{\tau}$ only depends on the preceding configuration $X_{\tau-1}$ and the current hitting cost $f_{\tau}$. This generally allows for a more space and time-efficient algorithm, which is important when choosing a small time slot length $\delta$ to be more responsive to changes in load.

The algorithm proposed by  <cit.> works as follows. Let $\hat{x}$ be the minimizer of $f_{\tau}(x)$, i.e. $\hat{x} = \argmin_{x \in \mathcal{X}} f_{\tau}(x)$. The algorithm moves into the direction of the minimizer until it either reaches $\hat{x}$, or it reaches a configuration $x$ where its switching cost equals twice the hitting cost of $x$. <Ref> gives an example of a step of the algorithm. We observe that this is equivalent to the following convex optimization:
\begin{align}\label{eq:ud:memoryless}\begin{aligned}
&\min_{x \in \mathcal{X}} &&f_{\tau}(x) \\
&\text{subject to}        &&\beta |x - X_{\tau-1}| \leq \frac{f_{\tau}(x)}{2}.
\end{aligned}\end{align}
Originally,  <cit.> proposed this algorithm for a restricted variant of uni-dimensional SSCO where the decision space is unbounded, i.e., $\mathcal{X} = \mathbb{R}$, and the switching costs are given by the $\ell_1$ norm. In particular, they choose $\beta = 1$. First, it is easy to see that we can adapt the algorithm for a bounded decision space by bounding the feasible region of the optimization problem in <ref>. Second, we observe that $\beta$ can simply be interpreted as the weight that we associate with smoothing (i.e., minimizing movement) instead of minimizing hitting costs. This is shown by the following equation, which is obtained by dividing the cost of <ref> by $\beta$:
\begin{align}\label{eq:simplified_smoothed_convex_optimization:without_beta}
\sum_{t=1}^T \frac{1}{\beta} f_t(X_t) + \sum_{k=1}^d (X_{t,k} - X_{t-1,k})^+.
\end{align}
The cost associated with this equation is the cost of <ref> linearly scaled by $1 / \beta$. Especially, this argument shows the following lemma.

A schedule is optimal with respect to <ref> if and only if it is optimal with respect to <ref>.

Therefore, without loss of optimality, we can incorporate the weight of the switching cost $\beta$ into the hitting costs. Further, note that in their model,  <cit.> consider the absolute movement, i.e. $|X_{t,k} - X_{t-1,k}|$, rather than only positive movements, i.e. $(X_{t,k} - X_{t-1,k})^+$. However, with <ref> in <ref>, we have shown that these switching costs only differ by a constant factor (namely $1/2$).

The resulting algorithm is simply given by determining $\hat{x}$ based on the convex optimization in <ref>, see <ref>. Thus, the time (and space) complexity of this memoryless algorithm is $\mathcal{O}(C O_{\epsilon}^1)$ for finding $\epsilon$-optimal solutions.

The memoryless algorithm moves towards the minimizer of the hitting cost, balancing hitting and movement costs.

Memoryless algorithm <cit.>
$\mathcal{I}_{\text{SSCO}} = (\tau \in \mathbb{N}, m \in \mathbb{N}, \beta \in \mathbb{R}_{>0}, (f_1, \dots, f_{\tau}) \in (\mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0})^{\tau})$
$\hat{x}$ such that that $\hat{x}$ is the result of the optimization in <ref>

§.§ Probabilistic Algorithm

Next, we discuss a $2$-competitive algorithm developed by  <cit.>, which works by maintaining a probability distribution over configurations. Using this probability distribution, they describe a randomized algorithm which they subsequently translate into a deterministic algorithm. We discuss how to gather a randomized and then a deterministic algorithm from a probability distribution. Then, we describe how  <cit.> determine the probability distribution and how it can be computed in practice.

§.§.§ From Probability Distribution to Deterministic Algorithm

Let's suppose we have given a probability distribution $p$ over configurations $x \in \mathcal{X}$. A randomized algorithm is then described by initially picking a number $\gamma \in [0,1]$ uniformly at random and then maintaining the invariant that at time $\tau$ the chosen configuration $x_{\tau}$ has the property that the probability mass to the left of $x_{\tau}$ with respect to $p$ is exactly $\gamma$ <cit.>. Crucially, this approach only works in the fractional setting. Also, note that $\gamma$ is chosen only once prior to running the algorithm. This describes how we obtain a randomized algorithm from a probability distribution over configurations.

Next,  <cit.> show the following theorem, which describes how we can obtain a deterministic algorithm from a randomized algorithm.
<cit.> For the problem of (fractional) online convex optimization, if there exists a $\rho$-competitive randomized algorithm $\mathcal{R}$ then there exists a $\rho$-competitive deterministic algorithm $\mathcal{D}$.

<cit.> prove this theorem using Jensen's inequality. In the setting of a probability space, Jensen's inequalityJensen's inequality claims that given a convex function $\varphi$ and a random variable $X$ we have
\begin{align}
E(\varphi(X)) \geq \varphi(E X)
\end{align}
provided both expectations exist, i.e. $E |X|$ and $E |\varphi(X)| < \infty$ <cit.>.

Let $X_{\tau}$ be a random variable denoting the configuration of the randomized algorithm $\mathcal{R}$ at time $\tau$. Then, the deterministic algorithm $\mathcal{D}$ of  <cit.> sets their configuration to $x_{\tau} = E X_{\tau}$. The cost of $\mathcal{D}$ is thus given by $f_{\tau}(x_{\tau}) + (x_{\tau} - x_{\tau-1})^+$ and the cost of $\mathcal{R}$ is given by $E(f_{\tau}(X_{\tau})) + E((X_{\tau} - X_{\tau-1})^+)$. We observe that both $f_{\tau}$ and $(\cdot)^+$ are convex functions, implying that the cost of $\mathcal{R}$ is at least $f_{\tau}(E X_{\tau}) + (E X_{\tau} - E X_{\tau-1})^+$ which equals the cost of $\mathcal{D}$. Summing over all $t$ completes the proof.

Hence, we have seen that a deterministic algorithm can be obtained from a randomized algorithm by, in each time slot, choosing the expected configuration of the randomized algorithm.

§.§.§ Assumptions

In the description of their algorithm,  <cit.> consider a restricted variant of uni-dimensional SSCO. Similar to their memoryless algorithm, which we discussed in <ref>, they consider an unbounded decision space, i.e., $\mathcal{X} = \mathbb{R}$, and the $\ell_1$ norm as switching costs. Further, for their description of this probabilistic algorithm, they assume that the minimizer $\hat{x}$ of $f_{\tau}$ is unique and bounded and that the hitting costs $f_{\tau}$ are continuous and smooth, i.e., are infinitely many times continuously differentiable. In particular, they assume the first-order and second-order derivatives of $f_{\tau}$ are well-defined and continuous. In <ref>, we discussed the assumption of differentiability and how it relates to our data center model.

Our implementation generalizes their algorithm to instances of SSCO with a bounded decision space $\mathcal{X}$, variable switching costs $\beta$, and piecewise linear functions. The second assumption, namely that the minimizer of the hitting cost is bounded, is natural in a data center setting as revenue loss increases for small configurations, whereas energy costs increase for large configurations.

In summary, the final algorithm is $2$-competitive for arbitrary instances of uni-dimensional SSCO with the restriction that hitting costs must either be piecewise linear or smooth.

§.§.§ The Probability Distribution

For any time $\tau$, the algorithm maintains a probability distribution $p_{\tau}$ over configurations $x \in \mathcal{X}$. So $\int_a^b p_{\tau}(x) \,dx$ represents the probability that $X_{\tau} \in [a,b]$ for any two $a, b \in \mathcal{X}$. At each time step $\tau$ we first find the minimizer of $f_{\tau}$, $\hat{x} = \argmin_{x \in \mathcal{X}} f_{\tau}(x)$. Then, we find a point $x_r \geq \hat{x}$ such that
\begin{align}\label{eq:ud:probabilistic:right}
\frac{1}{2} \int_{\hat{x}}^{x_r} \diff[2]{f_{\tau}}{y}(y) \,dy = \beta \int_{x_r}^{\infty} p_{\tau-1}(y) \,dy
\end{align}
and a point $x_l \leq \hat{x}$ such that
\begin{align}\label{eq:ud:probabilistic:left}
\frac{1}{2} \int_{x_l}^{\hat{x}} \diff[2]{f_{\tau}}{y}(y) \,dy = \beta \int_{-\infty}^{x_l} p_{\tau-1}(y) \,dy.
\end{align}
Note that we use <ref> to linearly scale the hitting cost $f_{\tau}$ by $1 / \beta$ to allow for $\beta \neq 1$. We then simply moved the constant factor outside of the derivative and integral. The probability distribution is updated as follows:
\begin{align}\label{eq:ud:probabilistic:update}
p_{\tau}(x) = \begin{cases}
p_{\tau-1}(x) + \frac{1}{2 \beta} \diff[2]{f_{\tau}}{x}(x) & x \in [x_l,x_r] \\
0 & \text{otherwise}
\end{cases}
\end{align}
where $p_0$ is a discrete distribution concentrating all probability mass in the point $0$. Note that  <cit.> do not assume any particular initial distribution, yet in our original problem statement we assumed $X_0 = \mathbf{0}$. The continuous extension of this distribution can be approximated as $p_0 \sim \text{Unif}(0, \epsilon)$ for a suitably small $\epsilon > 0$. In our implementation we choose $\epsilon = 10^{-5}$.

§.§.§ The Algorithm

To begin with, recall that the algorithm developed by  <cit.> operates on an unbounded decision space. To translate the algorithm to a setting with a bounded decision space, it is easy to see that we need to ensure that the underlying probability distribution does not assign positive probability to $x \not\in \mathcal{X}$. This can be achieved by introducing the additional restrictions $0 \leq x_l$ and $x_r \leq m$ which requires the assumption $\hat{x} \in [0,m]$. This is not a restriction as we simply define $\hat{x}$ as the minimizer of the hitting cost $f_{\tau}$ with respect to the decision space $\mathcal{X}$.

We use a convex optimization (as described in <ref>) to find the minimizer of the hitting cost $\hat{x}$. To determine $x_l$ and $x_r$, we use Brent's method with suitably defined functions and intervals to find a root. We will define these intervals and functions and describe how they can be computed in the following. Before beginning their description, note that we can only use a bracketed root finding method as we assumed that our decision space is bounded. If the decision space was not bounded, $x_l$ could be determined using a search for a local optimum minimizing $x$ starting from $\hat{x}$ under the equality constraint given in <ref>. This works because the equality constraint reduces the dimensionality of the optimization to $0$, resulting in a single feasible point. The analogous approach can be used to determine $x_r$.

We begin by describing how $x_r$ can be determined by a local search for a root. First, note that $x_r \in [\hat{x},m]$. Next, we restate <ref> as a function of $x_r$:
\begin{align*}
&& \frac{1}{2} \int_{\hat{x}}^{x_r} \diff[2]{f_{\tau}}{y}(y) \,dy =&\ \beta \int_{x_r}^{\infty} p_{\tau-1}(y) \,dy \\
\iff&& \left(\diff{f_{\tau}}{x}(x_r) - \diff{f_{\tau}}{x}(\hat{x})\right) =&\ 2 \beta \int_{x_r}^{\infty} p_{\tau-1}(y) \,dy \\
\iff&& g(x_r) := \diff{f_{\tau}}{x}(x_r) - 2 \beta \int_{x_r}^{\infty} p_{\tau-1}(y) \,dy =&\ 0.
\end{align*}
Here, we used the fundamental theorem of calculus and that the first order derivative of $f_{\tau}$ at $\hat{x}$ is $0$ as $\hat{x}$ is the minimizer of $f_{\tau}$. We observe that $g(\hat{x}) \leq 0$ and $g(m) \geq 0$. As $g$ is continuous, we know that we can be sure to find a root $x_r$ of $g$ on the interval $[\hat{x},m]$.

We take the analogous approach to determine $x_l \in [0,\hat{x}]$. Using <ref>, we obtain:
\begin{align*}
&& \frac{1}{2} \int_{x_l}^{\hat{x}} \diff[2]{f_{\tau}}{y}(y) \,dy =&\ \beta \int_{-\infty}^{x_l} p_{\tau-1}(y) \,dy \\
\iff&& \left(\diff{f_{\tau}}{x}(\hat{x}) - \diff{f_{\tau}}{x}(x_l)\right) =&\ 2 \beta \int_{-\infty}^{x_l} p_{\tau-1}(y) \,dy \\
\iff&& h(x_l) := 2 \beta \int_{-\infty}^{x_l} p_{\tau-1}(y) \,dy - \diff{f_{\tau}}{x}(x_l) =&\ 0.
\end{align*}
Again, we observe that $h(0) \leq 0$, $h(\hat{x}) \geq 0$, and $h$ is continuous, implying that we can be sure to find a root $x_l$ of $h$ on the interval $[0,\hat{x}]$. <Ref> illustrates the choice of $x_l$ and $x_r$.

Visualization of the choice of $x_l$ and $x_r$ of the probabilistic algorithm.

Root Finding The previous arguments show that a bracketed root finding method can be used to find $x_l$ and $x_r$. We use Brent's methodBrent's method for root finding. Brent's method combines the bisection method with higher-order methods to guarantee convergence to the root, yet at a higher rate than if only bisection were used <cit.>.  <cit.> “recommend it as the method of choice for general one-dimensional
root finding where a function’s values only (and not its derivative or functional form)
are available”. We denote the convergence rate of approximating the root with tolerance $\epsilon$ by $\mathcal{O}(R_{\epsilon})$.

Numerical Differentiation We use the five-point stencilfive-point stencil
\begin{align*}
\diff{f_{\tau}}{x}(x) \approx \frac{-f_{\tau}(x - 2h) + 8 f_{\tau}(x + h) - 8 f_{\tau}(x - h) + f_{\tau}(x - 2h)}{12h}
\end{align*}
to find a finite difference approximation of order $\mathcal{O}(h)$ of the first order derivative of $f_{\tau}$ at configurations $x \in \mathcal{X}$ <cit.>. To match the accuracy of our convex optimizations we set $h := \epsilon / 10$. To approximate the second order derivative of $f_{\tau}$ at a configuration $x \in \mathcal{X}$ we use
\begin{align*}
\diff[2]{f_{\tau}}{x}(x) \approx \frac{-f_{\tau}(x + 2h) + 16 f_{\tau}(x+h) - 30 f_{\tau}(x) + 16 f_{\tau}(x-h) - f_{\tau}(x - 2h)}{12 h^2}
\end{align*}
which yields an approximation of order $\mathcal{O}(h^4)$ <cit.>. Thus, we set $h := (\epsilon / 10)^{-1/4}$. We are thus able to compute these approximations in $\mathcal{O}(C)$ time.

Next, we describe how the constraints from <ref> and <ref> can be computed numerically.

Numerical Integration In our implementation, we need to compute both finite and semi-infinite integrals over the probability distribution $p$.

We use the Tanh-sinh quadratureTanh-sinh quadrature (also known as the double exponential method) to compute finite integrals.  <cit.> describe the convergence and error of this method in more detail. They conclude that “overall, the tanh-sinh scheme appears to be the best
for integrands of the type most often encountered in experimental math research” and highlight that it has “excellent accuracy and runtime performance” <cit.>.

We use the Gauss-Laguerre quadratureGauss-Laguerre quadrature for semi-infinite integrals, which approximates values of integrals of the kind
\begin{align}\label{eq:gauss_laguerre}
\int_0^{\infty} e^{-x} f(x) \,dx
\end{align}
<cit.>. It is easy to see that we can eliminate the weights by multiplying the integrand $g$ with $e^x$. Let $\text{GL}(g)$ denote the approximation of <ref> obtained by the Gauss-Laguerre quadrature. We are then able to compute any right-open integral over the interval $[{a,\infty})$ with integrand $p$ by setting $g(x) := p(a+x)$ and any left-open integral over the interval $({-\infty,b}]$ with integrand $p$ by setting $g(x) := p(b-x)$.

Crucially, for numeric stability in both integration schemes, we need that the integrands are continuous. It is easy to see that, in general, this is not the case for our probability distribution $p$. We describe in the following paragraph how integrals can be suitably discretized to allow for stable numeric results. Moreover, as the integration schemes are not universal, probability distributions exist for which the used quadratures cannot find the integral. In such a case, one would have to resort to another integration scheme. We denote the convergence rate of approximating the integral with tolerance $\epsilon$ by $\mathcal{O}(I_{\epsilon})$.

Piecewise Linear Hitting Costs We have already discharged the assumptions that $\beta = 1$ and that the decision space is unbounded. It remains to extend this algorithm to allow for piecewise linear hitting costs.

In our description of the adaption to piecewise linear hitting costs $f_{\tau}$, we refer to the non-continuous or non-smooth points of $f_{\tau}$ as breakpointsbreakpoint. For a piecewise linear function, we can discretize the integral into a summation, replace the first-order derivative at a breakpoint by the difference of consecutive points, and replace the second-order derivative at a point by the difference in slopes of consecutive points <cit.>. Let $B_{f_{\tau}}$ denote the set of breakpoints of $f_{\tau}$. Further, we denote by $x_{\tau,l}$ and $x_{\tau,r}$ the values of $x_l$ and $x_r$ at time $\tau$, respectively. It is easy to see that the set of breakpoints of $p_{\tau}$ is then given by
\begin{align*}
B_{p_{\tau}} := \{0, m\} \cup \left(B_{f_1} \cup \{x_{1,l}, x_{1,r}\}\right) \cup \dots \cup \left(B_{f_{\tau}} \cup \{x_{\tau,l}, x_{\tau,r}\}\right).
\end{align*}

Let $B_{p_{\tau}}^I := B_{p_{\tau}} \cap I$. The integral of $p_{\tau}$ over $I \subseteq \mathbb{R} \cup \{-\infty, \infty\}$ can then be computed using the quadrature methods described previously by integrating piecewise:
\begin{align*}
\int_I p_{\tau}(y) \,dy = \int_{\min I}^{\min B_{p_{\tau}}^I} p_{\tau}(y) \,dy + \int_{\max B_{p_{\tau}}^I}^{\max I} p_{\tau}(y) \,dy + \sum_{(i, j) \in \text{sort}(B_{p_{\tau}}^I)} \int_i^j p_{\tau}(y) \,dy
\end{align*}
where $\text{sort}(A)$ denotes the pairs of consecutive elements of some set $A \subset \mathbb{R}$ in ascending order. We can continue to use finite difference approximations for the first-order and second-order derivatives.

Note that the computation of $p_{\tau}$ requires $\mathcal{O}(\tau)$ many approximations of the second-order derivative of $f_{t}$. Therefore, any evaluation of $p_{\tau}$ requires $\mathcal{O}(\tau C)$ time and we are thus able to compute the $\epsilon$-optimal integral in $\mathcal{O}(\tau C I_{\epsilon} |B_{p_{\tau}}|) = \mathcal{O}(\tau^2 C I_{\epsilon} |B_{f_0}|)$ time, assuming $B_{f_{\tau}}$ add a constant number of new breakpoints in each time step and with $B_{f_0}$ denoting the number of breakpoints that is shared between multiple $f_{\tau}$. Overall, the described convex optimizations can be solved $\epsilon$-optimally in $\mathcal{O}(\tau^2 C I_{\epsilon} |B_{f_0}| R_{\epsilon} O_{\epsilon}^1)$ time. This is also the time complexity of the algorithm. This shows that, similarly to lazy capacity provisioning, the computational complexity grows polynomially with time. However, unlike lazy capacity provisioning, we cannot regularly reset the history, rendering this algorithm computationally inefficient for short time slot lengths. The algorithm is described in <ref>.

Probabilistic algorithm <cit.>
$\mathcal{I}_{\text{SSCO}} = (\tau \in \mathbb{N}, m \in \mathbb{N}, \beta \in \mathbb{R}_{>0}, (f_1, \dots, f_{\tau}) \in (\mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0})^{\tau})$
$\hat{x} \gets \argmin_{x \in \mathcal{X}} f_{\tau}(x)$
find $x_r$ using the optimization described by <ref> subject to $x \in \mathcal{X}$
find $x_l$ using the optimization described by <ref> subject to $x \in \mathcal{X}$
set $p_{\tau}$ based on the update rule in <ref>
$\int_{x_l}^{x_r} y \cdot p_{\tau}(y) \,dy$

Updating the probability distribution $p$ can be done in constant time as this does not require any function evaluations. As discussed in the beginning of this subsection, given a uniformly picked $\gamma \in [0,1]$, the randomized algorithm chooses $X_{\tau}$ (randomly) such that $\int_{-\infty}^{X_{\tau}} p_{\tau}(y) \,dy = \gamma$. In other words, $P_{\tau}(X_{\tau}) \sim \text{Unif}(0,1)$ where $P_{\tau}$ is the cumulative distribution function of $p_{\tau}$. By the universality of the uniform, $X_{\tau}$ is $P_{\tau}$-distributed, i.e. $X_{\tau} \sim P_{\tau}$. Hence, $E X_{\tau} = \int_{x_l}^{x_r} y \cdot p_{\tau}(y) \,dy$ computes the configuration for time slot $\tau$ as $p_{\tau}(y) = 0$ for $y \not\in [x_l, x_r]$. We then return $E X_{\tau}$. Similarly to the previously discussed integrals, this integral can be computed $\epsilon$-optimally in $\mathcal{O}(\tau^2 C I_{\epsilon} |B_{f_0}|)$ time, not affecting the asymptotic time complexity of the algorithm.

§.§ Randomly Biased Greedy Algorithm

We have seen many constant-competitive online algorithms for uni-dimensional smoothed convex optimization. However, we have not yet paid much attention to minimizing regret. This is also partially because sublinear regret can be achieved easily using online gradient descent. As this approach generalizes to the multi-dimensional setting, we discuss this approach in <ref>, where we look more generally at descent methods.

Still, one important question regarding the uni-dimensional setting remains, namely, whether there is an algorithmic framework that achieves both a constant-competitive ratio and sublinear regret. While it is impossible to achieve both simultaneously, it is possible to develop an algorithmic framework that balances both performance metrics.

In their paper, where the incompatibility between the competitive ratio and regret was first introduced,  <cit.> already proposed an algorithmic framework for SCO balancing the two notions in the uni-dimensional setting. Their approach is to scale the norm used to penalize movement in the decision space with a parameter $\theta \geq 1$. If $\theta = 1$, i.e., the algorithm solves the original problem, their algorithm is $2$-competitive and has linear regret. In contrast, for $\theta > 1$, movement in the decision space is penalized more. This allows reducing the regret to an arbitrary amount (which still depends linearly on $T$) while maintaining a constant competitive ratio <cit.>.

Note that while regret is understood as introduced in <ref>, the competitive ratio is understood with respect to the modified problem with lookahead $1$. In general, with lookaheadlookahead $i$, the environment plays actions $i$ steps before the agent follows suit. In other words, a step at time $i$ is evaluated using the cost function from time $t-i$, and the initial step $X_i$ is $\mathbf{0}$. For lookahead $i$, we consider the modified overall cost
\begin{align*}
\sum_{t=1}^T f_t(X_{t+i}) + \norm{X_{t+i} - X_{t+i-1}}.
\end{align*}
Note that in this definition we assume that $f_t$ is known after $X_t$ is played, whereas in our original definition of SOCO we assumed that $f_t$ is known before $X_t$ is played. Hence, for $i=1$, this modified problem is equivalent to SOCO and similar to metrical task systems where, in both cases, the environment plays first. In contrast, with online convex optimization, the agent plays first, resulting in a lookahead of $i=0$ <cit.>. Note that online convex optimization is the more restricted setting as the agent has less knowledge than in metrical task systems, implying that given an algorithm for lookahead $i$, the corresponding algorithm for lookahead $0$ is as competitive. Given an algorithm with lookahead $i$, the corresponding algorithm with lookahead $0$ is obtained simply by shifting determined configurations $i$ time slots into the past.

<cit.> mention in their paper that the claims on this algorithm were withdrawn, however,  <cit.> clarified their proof, showing it is correct as stated in the original paper <cit.>.

The algorithm of  <cit.> is initialized with a random parameter $r$ which is uniformly sampled from $({-1,1})$, i.e. $r \sim \text{Unif}(-1, 1)$. They define the work function
\begin{align}\label{eq:randomly_biased_greedy:work_function}
w_{\tau}(x) = \min_{y \in \mathcal{X}} w_{\tau-1}(y) + f_{\tau}(y) + \theta \norm{x - y}
\end{align}
where $w_0(x) = \theta \norm{x}$. During each time slot $\tau$ the algorithm moves to the configuration $x$ minimizing $w_{\tau-1}(x) + r \theta \norm{x}$. The resulting algorithm is described in <ref>.

Randomly Biased Greedy <cit.>
$\mathcal{I}_{\text{SCO}} = (\tau \in \mathbb{N}, \mathcal{X} \subset \mathbb{R}, \norm{\cdot}, (f_1, \dots, f_{\tau}) \in (\mathcal{X} \to \mathbb{R}_{\geq 0})^{\tau}), \theta \geq 1, r \sim \text{Unif}(-1,1)$
$x \gets \argmin_{x \in \mathcal{X}} w_{\tau-1}(x) + r \theta \norm{x}$
set $w_{\tau}$ as described in <ref>

Observe that as expected, the algorithm returns $X_1 = \mathbf{0}$ for the initial step.

<cit.> show that given a $\theta \geq 1$ their algorithm attains the $\alpha$-unfair competitive ratio $(1+\theta) / \min \{\theta, \alpha\}$ and regret $\mathcal{O}(\max \{T / \theta, \theta\})$. Hence, for $\alpha = 1$ and $\theta = 1$ the algorithm is $2$-competitive. For any $\alpha > 0$, the optimal $\alpha$-unfair competitive ratio is $1 + 1 / \alpha$ and obtained by setting $\theta = \alpha$. In contrast, $\theta = 1 / \epsilon$ for some $\epsilon > 0$ yields the minimal regret $\mathcal{O}(\epsilon T)$ <cit.>. In general, when $T$ is known in advance,  <cit.> show that for $\theta \in \mathcal{O}(\sqrt{T})$, their algorithm achieves a $\mathcal{O}(\sqrt{T})$ $\alpha$-unfair competitive ratio and $\mathcal{O}(\sqrt{T})$-regret.

It is easy to see that the work function itself is convex as it can be interpreted as the inf-projection of the convex function $f(x,y) = w_{\tau-1}(y) + f_{\tau}(y) + \theta \norm{x - y}$. This fact is shown in proposition 2.22 of <cit.>. The evaluation of $w_{\tau}$ requires $\mathcal{O}(\tau)$ recursive evaluations of the work function each of which uses a convex optimization and evaluates the hitting costs. Thus, an $\epsilon$-optimal evaluation of $w_{\tau}$ can be obtained in $\mathcal{O}(C (O_{\epsilon}^1)^{\tau})$ time. Hence, the overall time complexity of the algorithm is in $\mathcal{O}(C (O_{\epsilon}^1)^{\tau+1})$. We improve the practical runtime by memoizing the work function.

§.§ Randomized Integral Relaxation

<cit.> use the probabilistic algorithm described in <ref> in their randomized algorithm for Int-SSCO achieving the optimal competitive ratio $2$ against an oblivious adversary. Although they used the probabilistic algorithm in their paper, their proof generalizes to any $2$-competitive fractional algorithm, so, in particular, the randomly biased greedy algorithm can be used as well. Roughly, the algorithm works by solving the relaxed problem using the probabilistic algorithm of  <cit.> and then randomly rounding the resulting fractional schedule.

Let $\bar{\mathcal{I}} = (T, m, \beta, \bar{F})$ with $\bar{F} = (\bar{f}_1, \dots, \bar{f}_T)$ be the fractional relaxation of the instance $\mathcal{I} = (T, m, \beta, F)$ of Int-SSCO and let $\bar{\mathcal{X}} = [0,m]$ denote the decision space of $\bar{\mathcal{I}}$.  <cit.> define the relaxed operating costs $\bar{f}_{\tau} : \bar{\mathcal{X}} \to \mathbb{R}_{\geq 0}$ as the linear interpolation of the integral operating costs $f_{\tau}$:
\begin{align*}
\bar{f}_{\tau} := \begin{cases}
f_{\tau}(x) & x \in [m]_0 \\
(\lceil x \rceil - x) f_{\tau}(\lfloor x \rfloor) + (x - \lfloor x \rfloor) f_{\tau}(\lceil x \rceil) & \text{otherwise}.
\end{cases}
\end{align*}
Note that $\bar{f_{\tau}}$ are continuous and piecewise linear with the set of breakpoints $[m]_0$. Hence, we are able to use <ref> to obtain the configuration $\bar{X}_{\tau} \in \bar{\mathcal{X}}$ at time $\tau$ for the relaxed problem instance  $\bar{\mathcal{I}}$. Further, let $\text{frac}(x) = x - \lfloor x \rfloor$ be the fractional part of $x$ and let $\bar{X}'_{\tau-1} = (\bar{X}_{\tau-1})_{\lfloor\bar{X}_{\tau}\rfloor}^{\lceil\bar{X}_{\tau}\rceil}$ be the projection of the preceding relaxed configuration onto the discrete interval of the current relaxed configuration.

The randomized algorithm distinguishes between time slots where the configuration is increased and time slots where the configuration is decreased. In the first case, i.e. $\bar{X}_{\tau-1} \leq \bar{X}_{\tau}$, if $X_{\tau-1} = \lceil\bar{X}_{\tau}\rceil$  the configuration remains unchanged. Otherwise, $X_{\tau}$ is set to $\lceil\bar{X}_{\tau}\rceil$ with probability
\begin{align*}
p_{\tau}^{\uparrow} := \frac{\bar{X}_{\tau} - \bar{X}'_{\tau-1}}{1 - \text{frac}(\bar{X}'_{\tau-1})}
\end{align*}
and to $\lfloor\bar{X}_{\tau}\rfloor$ with probability $1 - p_{\tau}^{\uparrow}$. Conversely, if $\bar{X}_{\tau-1} > \bar{X}_{\tau}$, the configuration remains unchanged if $X_{\tau-1} = \lfloor\bar{X}_{\tau}\rfloor$, and otherwise with probability
\begin{align*}
p_{\tau}^{\downarrow} := \frac{\bar{X}'_{\tau-1} - \bar{X}_{\tau}}{\text{frac}(\bar{X}'_{\tau-1})}
\end{align*}
the configuration is set to $\lfloor\bar{X}_{\tau}\rfloor$ and with probability $p_{\tau}^{\downarrow}$ the configuration is set to $\lceil\bar{X}_{\tau}\rceil$. The resulting algorithm is shown in <ref>.

Randomized integral relaxation <cit.>
$\mathcal{I}_{\text{Int-SSCO}} = (\tau \in \mathbb{N}, m \in \mathbb{N}, \beta \in \mathbb{R}_{>0}, (f_1, \dots, f_{\tau}) \in (\mathbb{N}_0 \to \mathbb{R}_{\geq 0})^{\tau})$
$\bar{X}_{\tau} \gets \text{\cref{alg:ud:probabilistic}}(\bar{\mathcal{I}}_{\text{Int-SSCO}})$
$\bar{X}_{\tau-1} \leq \bar{X}_{\tau}$
$X_{\tau-1} = \lceil\bar{X}_{\tau}\rceil$

$\gamma \sim \text{Unif}(0,1)$
$\gamma \leq p_{\tau}^{\uparrow}$

$X_{\tau-1} = \lfloor\bar{X}_{\tau}\rfloor$

$\gamma \sim \text{Unif}(0,1)$
$\gamma \leq p_{\tau}^{\downarrow}$

We use the universality of the uniform to simulate Bernoulli-distributed random variables with parameters $p_{\tau}^{\uparrow}$ and $p_{\tau}^{\downarrow}$, respectively. Any pseudo-random number generator can be used to produce the uniformly distributed $\gamma$. It is easy to see that the time complexity is given by the time complexity of <ref>, i.e., $\mathcal{O}(\tau^2 m C I_{\epsilon} R_{\epsilon} O_{\epsilon}^1)$ with $|B_{f_0}| \in \mathcal{O}(m)$, or the complexity of <ref>, i.e., $\mathcal{O}(C (O_{\epsilon}^1)^{\tau+1})$ depending on which algorithm is used for the relaxed problem.

§ MULTI-DIMENSIONAL

§.§ Lazy Budgeting

To begin with our discussion of the multi-dimensional setting, we examine two algorithms developed by  <cit.> for a restricted class of convex cost functions. Their first algorithm is $2d$-competitive for SLO, i.e., load and time-independent costs, and their second algorithm is $(2d+1)$-competitive for SBLO. Note that we only defined SBLO and SLO for the integral case. As no online algorithm for SLO can attain a competitive ratio smaller than $2d$, their first algorithm is optimal, and their second algorithm is nearly optimal <cit.>.

The idea behind both algorithms is to calculate optimal schedules up to the current time slot. Depending on this schedule, the algorithm decides if a server is powered up. To perform the smoothing, the algorithm remembers how long a server was idling and powers this server down if the idle duration surpasses a threshold.  <cit.> do not name their algorithms, yet, within this work, we refer to them as lazy budgeting methodslazy budgeting.

§.§.§ Lazy Budgeting for Smoothed Load Optimization

Let $\mathcal{I} = (d, \tau, m, \beta, \Lambda, c)$ be an instance of SLO.  <cit.> focus on a setting without inefficient server types. A server type $k$ is called inefficient if there exists another server type $k'$ where both the operating cost and the switching cost is lower, i.e. $c_k \geq c_{k'}$ and $\beta_k \geq \beta_{k'}$. In practice, this is not a restriction as a server of an inefficient server type is only ever powered up if all more efficient servers are already active as there is no trade-off between operating and switching costs. Furthermore, typically servers with a lower operating cost also have a higher switching cost. In addition,  assume that there are no duplicated server types, i.e., server types with equal operating and switching costs.

We assume that the server types are sorted in descending order by their operating costs, i.e. $c_1 > \dots > c_d$. As we excluded inefficient server types, switching costs are in ascending order, $\beta_1 < \dots < \beta_d$.

The algorithm of  <cit.> separates a problem instance into $m := \sum_{k=1}^d m_k$ lanes. Recall that with SLO, we assume that each active server can handle a single job during each time slot. The algorithm uses that at time slot $t$ there is a job in line $j$ if and only if $j \leq \lambda_{t}$. Thus, all servers represented by lines $j > \lambda_{t}$ are either inactive or idling.

Let $y_{t,j}$ denote the server type that handles the $j$-th lane during time slot $t$ given some underlying schedule $X$. We say $y_{t,j} = 0$ if there is no active server in lane $j$ during time slot $t$, which is only the case if $j > \lambda_{t}$ as we can assume that $\lambda_{t} \leq m$ holds for all time slots $t \in [T]$.  <cit.> give the following formal definition:
\begin{align*}
y_{t,j} := \begin{cases}
\max \{k \in [d] \mid \sum_{k' = k}^d X_{t,k'} \geq j\} & k \in \left[\sum_{k=1}^d X_{t,k}\right] \\
0 & \text{otherwise}.
\end{cases}
\end{align*}
By this definition, we prefer to use servers of the server type with the lowest operating cost (and largest switching cost). In other words, the server types handling each lane $y_{t,1}, \dots, y_{t,m}$ are sorted in descending order, i.e. $y_{t,j} \geq y_{t,j'}$ for $j < j'$. We denote by $\hat{y}_{t,j}^{\tau}$ the server type in lane $j$ during time slot $t$ induced by some optimal schedule $\hat{X}^{\tau}$ up to time $\tau$ and by $\widetilde{y}_{t,j}$ the server type in lane $j$ during time slot $t$ as assigned by the algorithm.

The algorithm begins by finding an optimal schedule $\hat{X}^{\tau}$ up to time slot $\tau$. This schedule is chosen such that the server type in a lane of $\hat{X}^{\tau}$ is never reduced compared to the previously used optimal schedule $\hat{X}^{\tau-1}$, i.e. $\hat{y}_{t,j}^{\tau} \geq \hat{y}_{t,j}^{\tau-1}$ for all time slots $t \in [\tau]$ and lanes $j \in [m]$. Moreover, we assume that $\hat{X}^{\tau}$ is a schedule that powers up servers as late as possible and powers down servers as early as possible. This is necessary in case $c_k = 0$ for some server type $k \in [d]$. We observe that these properties are fulfilled by all optimal schedules that <ref> obtains. We cache the results of the optimal graph-based algorithm such that in every iteration of the online algorithm, only one dynamic update needs to be performed. The asymptotic time complexity of this dynamic update is thus given as $\mathcal{O}(|\mathcal{M}| C d)$ where $|\mathcal{M}| \in \mathcal{O}(\prod_{k=1}^d m_k)$.

Now, the algorithm ensures that no server type is used for lane $j \in [m]$ that is smaller than the server type used by $\hat{X}^{\tau}$, i.e. $\widetilde{y}_{\tau,j} \geq \hat{y}_{\tau,j}^{\tau}$. If $\widetilde{y}_{\tau-1,j} < \hat{y}_{\tau,j}^{\tau}$, a server of type $\widetilde{y}_{\tau-1,j}$ is powered down and a server of type $\hat{y}_{\tau,j}^{\tau}$ is powered up. A server of type $k$ that is not replaced by a greater server type remains active for $\bar{t}_k := \lfloor \beta_k / c_k \rfloor$ time slots <cit.>. If $\hat{X}^{\tau}$ uses a smaller server type $k' \leq k$ in the meantime, then the server of type $k$ will run for at least $\bar{t}_{k'}$ further time slots.

The algorithm computes $\widetilde{y}_{\tau,j}$ directly. The corresponding number of active servers of type $k$ of the underlying schedule $\widetilde{X}$ can be obtained by $\widetilde{X}_{\tau,j} = |\{j \in [m] \mid \widetilde{y}_{\tau,j} = k\}|$. The resulting algorithm is shown in <ref>. Here, $h_j$ denotes the time until the server handling line $j \in [m]$ is powered down.

Lazy Budgeting for SLO <cit.>
$\mathcal{I}_{\text{SLO}} = (d \in \mathbb{N}, \tau \in \mathbb{N}, m \in \mathbb{N}^d, \beta \in \mathbb{R}_{>0}^d, \Lambda \in \mathbb{N}_0^{\tau}, c \in \mathbb{R}_{\geq 0}^d)$
Update the previously found optimal schedule $\hat{X}^{\tau-1}$ to $\hat{X}^{\tau}$ such that $\hat{y}_{t,j}^{\tau} \geq \hat{y}_{t,j}^{\tau-1}$ for all  $j \in [m]$
$j \gets 1$ $m$
$\widetilde{y}_{\tau-1,j} < \hat{y}_{\tau,j}^{\tau}$ $t \geq h_j$
$\widetilde{y}_{\tau,j} \gets \hat{y}_{\tau,j}^{\tau}$
$h_j \gets \tau + \bar{t}_{\hat{y}_{\tau,j}^{\tau}}$

$\widetilde{y}_{\tau,j} \gets \widetilde{y}_{\tau-1,j}$
$h_j \gets \max \{h_j, \tau + \bar{t}_{\hat{y}_{\tau,j}^{\tau}}\}$ where $\bar{t}_0 = 0$

$k \in [d]$
$X_{\tau,k} \gets |\{j \in [m] \mid \widetilde{y}_{\tau,j} = k\}|$

BuildLanes($x, d, m$)
$y \gets \mathbf{0}$
$j \gets 1$ $m$
$j \leq \sum_{k=1}^d x_k$
$k \gets 1$ $d$
$\sum_{k'=k}^d x_{k'} \geq j$
$y_j \gets k$

Schedules can be converted to lanes in $\mathcal{O}(m d^2)$ time as described by <ref> and lanes can be converted back to schedules in $\mathcal{O}(m)$ time by iterating over all lanes. Hence, the overall asymptotic time complexity of the algorithm is described by adding the time required to find the optimal schedule, i.e., $\mathcal{O}(m d^2 + C d \prod_{k=1}^d m_k)$.

§.§.§ Randomized Lazy Budgeting for Smoothed Load Optimization

<cit.> describe how the competitive ratio of the previously described <ref> can be improved to $\frac{e}{e-1}d \approx 1.582d$ against an oblivious adversary by randomizing the running time of a server. Before execution, the randomized algorithm chooses $\gamma \in [0,1]$ according to the probability density function
\begin{align*}
f_{\gamma}(x) = \begin{cases}
e^x / (e-1) & x \in [0,1] \\
0 & \text{otherwise}.
\end{cases}
\end{align*}
Then, the running time of a server of type $k \in [d]$, $\bar{t}_k$, is set to $\lfloor \gamma \cdot \beta_k / c_k \rfloor$.

To sample $\gamma$ we first seek to find the cumulative distribution function $F_{\gamma}$. For $x \in [0,1]$ we have
\begin{align*}
F_{\gamma}(x) &= \int_0^x f_{\gamma}(t) \,dt \\
&= \frac{1}{e-1} \int_0^x e^t \,dt \\
&= \frac{1}{e-1} (e^x - 1).
\end{align*}
By the universality of the uniform, realizations of $F_{\gamma}$ can be simulated given $U \sim \text{Unif}(0,1)$ as $F_{\gamma}^{-1}(U) \sim F_{\gamma}$. It is now easy to see that $F_{\gamma}^{-1}(x) = \ln (x (e - 1) + 1)$.

This completes the description of the implementation of the randomized variant of lazy budgeting for SLO. The asymptotic time complexity is given by the time complexity of <ref>.

§.§.§ Lazy Budgeting for Smoothed Balanced-Load Optimization

In a subsequent paper,  <cit.> modified their lazy budgeting method to SBLO, allowing for more complex cost functions. The method remains similar: First, an optimal schedule is found which ends at the current time slot. Then, the algorithm ensures for each server type that the number of active servers is at least as large as the number of active servers in the optimal schedule. Lastly, to make the algorithm competitive, if the accumulated idle operating cost of a server of type $k$, $g_{t,k}(0)$ exceeds the switching cost $\beta_k$, a server of type $k$ is powered down.

First,  <cit.> developed a $(2d+1)$-competitive online algorithm for a setting where the operating costs are time-independent. This allows to determine the runtime of a server in advance as both $g_{t,k}(0)$ and $\beta_k$ are known. Then, they extend their algorithm to a setting that allows for time-dependent operating costs where their algorithm attains a competitive ratio of $2d+1+\epsilon$ for any $\epsilon > 0$.

Time-Independent Operating Costs

Again, we denote by $\hat{X}^{\tau}$ the optimal schedule with information up to time slot $\tau$ and by $X$ the schedule obtained by the algorithm. If the optimal schedule has more active servers of some type $k \in [d]$ than the schedule obtained by the algorithm, i.e. $\hat{X}_{\tau,k}^{\tau} > X_{\tau-1,k}$, $(\hat{X}_{\tau,k}^{\tau} - X_{\tau-1,k})^+$ servers of type $k$ are powered up. After being active for $\bar{t}_k := \lfloor \beta_k / g_k(0) \rfloor$ time slots, a server of type $k$ is powered down again. The resulting algorithm is shown in <ref>. Here, $h_{t,k}$ denotes the number of servers of type $k$ that were powered up at time $t$. We assume that $h_{t,k}$ is initialized with $0$ for all $t \in \mathbb{Z}$ and $k \in [d]$.

Lazy Budgeting for SBLO (for time-independent operating costs) <cit.>
$\mathcal{I}_{\text{SBLO}} = (d \in \mathbb{N}, \tau \in \mathbb{N}, m \in \mathbb{N}^d, \beta \in \mathbb{R}_{>0}^d, \Lambda \in \mathbb{N}_0, G \in (\mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0}^d)^d)$
Update $\hat{X}^{\tau-1}$ to $\hat{X}^{\tau}$
$k \gets 1$ $d$
$X_{\tau,k} \gets X_{\tau-1,k} - h_{\tau - \bar{t}_k, k}$
$X_{\tau,k} < \hat{X}_{\tau,k}^{\tau}$
$X_{\tau,k} \gets \hat{X}_{\tau,k}^{\tau}$
$h_{t,k} \gets \hat{X}_{\tau,k}^{\tau} - X_{\tau,k}$

Similar to our implementation of Lazy Budgeting for SLO (<ref>), we cache intermediate results of the algorithm computing the optimal schedule up to time slot $\tau$ (<ref>). It is therefore enough to perform a single dynamic update to obtain $\hat{X}^{\tau}$ which is possible in $\mathcal{O}(|\mathcal{M}| C d)$ time where $|\mathcal{M}| \in \mathcal{O}(\prod_{k=1}^d m_k)$. It is easy to see that this is also the time complexity of <ref>.

Time-Dependent Operating Costs

Next,  <cit.> extend <ref> to an algorithm which supports time-dependent operating costs and achieves a competitive ratio of $2d + 1 + c(\mathcal{I})$ where $c(\mathcal{I}) := \sum_{k=1}^d \max_{t \in [T]} \frac{g_{t,k}(0)}{\beta_k}$. Now, the idle operating cost $g_{t,k}(0)$ is not constant over time anymore. Since at time $\tau$ we only know the operating costs up to time $\tau$, we cannot pre-determine the number of time slots that a server of type $k$ is active until powered down, yet at the time slot when a server is powered down we know all relevant cost functions.  <cit.> formally define the maximal number of time slots such that the sum of the idle operating costs beginning from the next time slot, $\tau+1$, is smaller than or equal to the switching cost as
\begin{align*}
\bar{t}_{t,k} := \max \left\{\bar{t} \in [T - t] \mid \sum_{t' = t+1}^{t+\bar{t}} g_{t,k}(0) \leq \beta_k\right\}
\end{align*}
In contrast to <ref> where servers of type $k$ that were powered up at time $\tau$ were active for $\bar{t}_{\tau,k}$ time slots, they are now active for $\bar{t}_{\tau,k} + 1$ time slots.

Let $W_{\tau,k}$ denote the set of all time slots $t$ such that servers of type $k$ that were powered up at time $t$ are powered down during time slot $\tau$, i.e., $t + \bar{t}_{t,k} + 1 = \tau$. Again, we denote by $h_{t,k}$ the number of servers of type $k$ that were powered up at time $t$. Using the same powering-up policy as <ref>, the updated algorithm is described in <ref>.

Lazy Budgeting for SBLO (for time-dependent operating costs) <cit.>
$\mathcal{I}_{\text{SBLO}} = (d \in \mathbb{N}, \tau \in \mathbb{N}, m \in \mathbb{N}^d, \beta \in \mathbb{R}_{>0}^d, \Lambda \in \mathbb{N}_0^{\tau}, G \in (\mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0}^d)^{d^{\tau}})$
Calculate $\hat{X}^{\tau}$
$k \gets 1$ $d$
$W_{\tau,k} \gets \left\{t \in [\tau-1] \mid \sum_{t'=t+1}^{\tau-1} g_{t',k}(0) \leq \beta_k < \sum_{t'=t+1}^{\tau} g_{t',k}(0)\right\}$
$X_{\tau,k} \gets X_{\tau-1,k} - \sum_{t \in W_{\tau,k}} h_{t, k}$
$X_{\tau,k} < \hat{X}_{\tau,k}^{\tau}$
$X_{\tau,k} \gets \hat{X}_{\tau,k}^{\tau}$
$h_{t,k} \gets \hat{X}_{\tau,k}^{\tau} - X_{\tau,k}$

$W_{\tau,k}$ can be computed in $\mathcal{O}(\tau^2 C)$ time. Thus, the overall time complexity of the algorithm is $\mathcal{O}(\tau^2 |\mathcal{M}| C d)$.

Reducing the Competitive Ratio

We now describe how  <cit.> improve the competitive ratio of their algorithm to $2d + 1 + \epsilon$ for any $\epsilon > 0$. The idea is to consider a modified problem instance $\widetilde{\mathcal{I}} = (d, \widetilde{\tau}, m, \beta, \widetilde{\Lambda}, \widetilde{G})$ which divides each time slot $t$ of the original problem instance $\mathcal{I}$ into $\widetilde{n}_t$ sub time slots, allowing for up to $\widetilde{n}_t$ intermediate state changes. Our goal is to choose $\widetilde{\mathcal{I}}$ such that $c(\widetilde{\mathcal{I}})$ becomes arbitrarily small.

In our description, we refer to time slots of $\mathcal{I}$ by $t$ and to sub time slots of $\widetilde{\mathcal{I}}$ by $u$. The total number of sub time slots is given by $\widetilde{\tau} := \sum_{t=1}^{\tau} \widetilde{n}_t$. We denote by $U(t) = [u+1 : u+\widetilde{n}_t] = \{u+1, u+2, \dots, u+\widetilde{n}_t\} \subseteq [\widetilde{\tau}]$ with $u = \sum_{t'=1}^{t-1} \widetilde{n}_t$ the set of sub time slots corresponding to time slot $t$. In contrast, let $U^{-1}(u)$ be the time slot $t \in [\tau]$ such that $u \in U(t)$.

The operating cost $g_{t,k}(l)$ of a server of type $k$ during time slot $t$ under load $l$ is divided equally among all sub time slots $\widetilde{n}_t$, i.e. $\widetilde{g}_{u,k}(l) := g_{U^{-1}(u),k}(l) / \widetilde{n}_{U^{-1}(u)}$. The job volume does not change, so $\widetilde{\lambda}_u := \lambda_{U^{-1}(u)}$.

<cit.> show that for
\begin{align*}
\widetilde{n}_t = \left\lceil\frac{d}{\epsilon} \cdot \max_{k \in [d]} \frac{g_{t,k}(0)}{\beta_k}\right\rceil
\end{align*}
the cost of the resulting schedule is at most $2d + 1 + \epsilon$ times larger than the cost of an optimal solution. For $\epsilon \to 0$, the competitive ratio converges to $2d + 1$.

Let $\widetilde{X}$ be the schedule obtained by <ref> for $\widetilde{\mathcal{I}}$. Then, the schedule $X$ for the original problem instance $\mathcal{I}$ can be obtained by setting $X_{\tau} = \widetilde{X}_{\mu(\tau)}$ where $\mu(\tau) = \argmin_{u \in U(\tau)} \widetilde{f}_{u}(\widetilde{X}_u)$ is the configuration that minimizes the operating cost during $U(t)$ <cit.>.

The resulting algorithm is shown in <ref>. First, the modified problem instance $\widetilde{\mathcal{I}}$ is created, and the next $\widetilde{n}_{\tau}$ time steps are simulated using <ref>. Then, the resulting schedule $X$ is constructed from $\widetilde{X}$.

Lazy Budgeting for SBLO <cit.>
$\mathcal{I}_{\text{SBLO}} = (d \in \mathbb{N}, \tau \in \mathbb{N}, m \in \mathbb{N}^d, \beta \in \mathbb{R}_{>0}^d, \Lambda \in \mathbb{N}_0^{\tau}, G \in (\mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0}^d)^{d^{\tau}})$
$\widetilde{n}_{\tau} \gets \left\lceil d / \epsilon \cdot \max_{k \in [d]} g_{\tau,k}(0) / \beta_k\right\rceil$
Extend the modified problem instance $\widetilde{\mathcal{I}}$ by $\widetilde{n}_{\tau}$ additional time slots
Update $\widetilde{X}$ by executing the next $\widetilde{n}_{\tau}$ time slots in <ref>
$X_{\tau} \gets \widetilde{X}_{\mu(\tau)}$ with $\mu(\tau) = \argmin_{u \in U(\tau)} \widetilde{f}_u(\widetilde{X}_u)$

The number of sub time slots $\widetilde{n}_t$ can be determined in $\mathcal{O}(C d)$ time. Creating the modified problem instance takes $\mathcal{O}(\widetilde{n}_{\tau})$ time. Simulating <ref> for $\widetilde{n}_{\tau}$ sub time slots takes $\mathcal{O}(\widetilde{n}_{\tau} \widetilde{\tau}^2 |\mathcal{M}| C d)$ time. The obtained schedule can be constructed in $\mathcal{O}(\widetilde{n}_{\tau} C)$ time. Hence, the overall asymptotic time complexity of <ref> is given by $\mathcal{O}(\widetilde{n}_{\tau} \widetilde{\tau}^2 |\mathcal{M}| C d)$, which is inversely proportional to $\epsilon^3$.

§.§ Descent Methods

When we seek fractional solutions, and the cost functions are differentiable, a promising method is to descent towards the minimizer of the current cost function. This is an approach that is commonly used in online convex optimization, yielding algorithms with sublinear regret <cit.>. As mentioned in <ref> on the Randomly Biased Greedy algorithm, movement costs are not considered in the online convex optimization setting, and the agent picks a point in the decision space before the cost function is revealed. A commonly used algorithm that achieves no-regret in the setting of online convex optimization is online gradient descent (OGD).

§.§.§ Online Gradient Descent

Online gradient descent works by selecting an arbitrary initial point $X_1 \in \mathcal{X}$ and then choosing, at time $\tau > 1$, $X_{\tau} = \Pi_{\mathcal{X}}(X_{\tau-1} - \eta_{\tau-1} \nabla f_{\tau-1}(X_{\tau-1}))$ where $\eta_t$ are the learning rates and $\Pi_{\mathcal{X}}(x)$ is the euclidean projection of $x$ onto $\mathcal{X}$ <cit.>. The algorithm is described in <ref>. Similar to RBG, OGD operates with lookahead $1$. In the multi-dimensional setting, the euclidean projectioneuclidean projection of a point $x \in \mathbb{R}^d$ onto a convex set $K$ is given as $\Pi_{K}(x) = \argmin_{y \in K} \norm{y - x}_2$.

Online Gradient Descent <cit.>
$\mathcal{I}_{\text{SCO}} = (\tau \in \mathbb{N}, \mathcal{X} \subset \mathbb{R}^d, \norm{\cdot}, (f_1, \dots, f_{\tau}) \in (\mathcal{X} \to \mathbb{R}_{\geq 0})^{\tau}), \eta > 0$
$X_{\tau} \gets \Pi_{\mathcal{X}}(X_{\tau-1} - \eta \nabla f_{\tau-1}(X_{\tau-1}))$

With appropriate learning rates OGD achieves no-regret for online convex otpimization. For example, OGD obtains $\mathcal{O}(\sqrt{T})$-regret for $\eta_t \in \Theta(1 / \sqrt{t})$ <cit.>.  <cit.> showed that an OGD algorithm with $\mathcal{O}(\rho_2(T))$-regret in the online convex optimization setting, and $\sum_{t=1}^T \eta_t \in \mathcal{O}(\rho_1(T))$, achieves $\mathcal{O}(\rho_1(T) + \rho_2(T))$-regret in the online smoothed convex otpimization setting. In particular, for learning rates $\eta_t \in \Theta(1 / \sqrt{t})$, OGD obtains $\mathcal{O}(\sqrt{T})$-regret when movement costs are considered and the agent picks a point after the hitting costs were revealed <cit.>.

Using finite difference methods, $\nabla f_{\tau-1}(X_{\tau-1})$ can be computed in $\mathcal{O}(d C)$ time. The euclidean projection can be computed $\epsilon$-optimally in $\mathcal{O}(O_{\epsilon}^d)$ time. Thus, the overall asymptotic time complexity of OGD is $\mathcal{O}(d C O_{\epsilon}^d)$.

§.§.§ Online Mirror Descent

Thus, concerning regret, where algorithms seek to minimize the hitting cost immediately, the smoothing property of SCO does not require the development of new algorithms. In contrast, competitive algorithms need to wait before moving to the minimizer of the hitting cost until the movement costs are amortized. In other words, in each step, a competitive algorithm has to decide how far to move into the direction of the minimizer to balance hitting cost and movement cost. Crucially, where to move depends on the geometry of the cost function. As is shown in <ref>, rather than moving towards the minimizer directly, it is advantageous to move to a projection onto some sub-level set of the cost function. This approach minimizes the movement costs that are required to reach a point with the same hitting cost. This approach can balance hitting costs and movement costs and is discussed in <ref>, where we introduce the online balanced descent framework.

Level sets of $f_{\tau}$ in two dimensions. The blue arrows show projections of $X_{\tau-1}$ onto some level set. The red line visualizes the projection of $X_{\tau-1}$ onto all level sets $\{x \in \mathbb{R}^2 \mid f_{\tau}(x) = l\}$ for $l \in [\hat{x}, f_{\tau}(X_{\tau-1})]$. The step in the direction of the minimizer of $f_{\tau}$ is shown in black. Note that it is not optimal to move directly in the direction of the minimizer as there likely exists a closer point on the same level set. Online balanced descent picks a point on the red line.

The mirror descent framework is an extension of gradient descent, allowing to adapt to the underlying “geometry” of a problem <cit.>. The original gradient descent algorithm uses Euclidean geometry. This is shown by a slightly modified form of its update rule:
\begin{align*}
X_{\tau} = \argmin_{x \in \mathcal{X}} \eta_{\tau-1} \langle\nabla f_{\tau-1}(X_{\tau-1}), x\rangle + \frac{1}{2} \norm{x - X_{\tau-1}}_2^2.
\end{align*}
Observe that OGD uses the squared Euclidean distance as a regularizer (i.e., a function ensuring that we remain close to the point $X_{\tau-1}$) which can be replaced by another distance to obtain different algorithms <cit.>.

Proximal Point View The Bregman divergence is a commonly used class of distance functions. Given a strictly convex distance-generating functiondistance-generating function $h$, the Bregman divergence measures the deviation of $h$ from its linear approximation.

Bregman divergence
<cit.> The Bregman divergence from a point $x$ to a point $y$ with respect to a strictly convex function $h$ is given as
\begin{align*}
D_h(x,y) = h(x) - h(y) - \langle\nabla h(y), x - y\rangle.
\end{align*}

Bregman Divergence $D_h(x,y)$ for a function $h : \mathbb{R} \to \mathbb{R}$ <cit.>. The Bregman divergence measures how much a function differs at $x$ from its linear approximation at $y$.

The definition of the Bregman divergence of a univariate function $h$ is visualized in <ref>.

The modified variant of OGD, which uses a Bregman divergence as a regularizer, is known as online mirror descentmirror descent (OMD) or online proximal gradient descent. Note that OMD is parametrized by $h$, which is used to describe the underlying geometry of the problem <cit.>.

For the function $h(x) = \frac{1}{2} \norm{x}_2^2$ from $\mathbb{R}^d$ to $\mathbb{R}$ (the squared $\ell_2$ normsquared $\ell_2$ norm), the Bregman divergence is the Euclidean distance, i.e. $D_h(x,y) = \frac{1}{2} \norm{x-y}_2^2$ <cit.>. Hence, OMD reduces to OGD if it is parametrized with the provided definition of $h$, i.e., a Euclidean geometry is used.

Mirror Map View The proximal point view yields just one perspective of mirror descent. Another perspective that is used frequently is the perspective of mirror maps.

Recall that in OGD, starting from some point $X_{\tau-1}$, we moved into the direction of the gradient of $f$. However, note that $\nabla f_{\tau-1}(X_{\tau-1})$ belongs to the dual space[The dual spacedual space of a vector space $V$ over some field $\mathbb{F}$ is the set of all linear maps from vectors in $V$ to scalars in $\mathbb{F}$ (which are called linear functionalslinear functionals) <cit.>] of $\mathbb{R}^d$. In the Euclidean space, this is not a problem as the dual space of the Euclidean space is the Euclidean space itself <cit.>. However, when working with normed spaces that are not self-dual, this is problematic.

Instead of adding elements from the dual space to elements from the primal space, mirror descent maps points from the primal space to the dual space, performs the gradient step in the dual space, and then maps the resulting point back to the primal space.

Visualization of a step of Mirror Descent. The previous point $X_{\tau-1}$ is first mapped to the dual space, $\theta_{\tau-1}$. Then, a step is taken into the direction of the gradient, $\theta_{\tau}$, and the resulting point is mapped back to the primal space. Finally, the resulting point $X'_{\tau}$ is projected back onto the feasible region $K$, resulting in the next point $X_{\tau}$ <cit.>.

mirror map
<cit.> Given a norm $\norm{\cdot}$ and a differentiable and $\alpha$-strongly convex function $h : \mathbb{R}^d \to \mathbb{R}$, the associated mirror map is $\nabla h : \mathbb{R}^d \to \mathbb{R}^d$ and the inverse mirror map is $(\nabla h)^{-1} : \mathbb{R}^d \to \mathbb{R}^d$.

For $h(x) = \frac{1}{2} \norm{x}_2^2$ the mirror map and its inverse are the identity map <cit.>. A complete description of the mirror descent framework is given in <ref>. Note that the choice of the mirror map is central as it describes the dual space (also called mirror image) where the gradient step is taken. <Ref> visualizes a step of mirror descent.

Online Mirror Descent <cit.>
$\tau \in \mathbb{N}, K \subset \mathbb{R}^d, \norm{\cdot}, h \in \mathbb{R}^d \to \mathbb{R}$
map to the dual space $\theta_{\tau-1} \gets \nabla h (X_{\tau-1})$
take a gradient step in the dual space $\theta_{\tau} \gets \theta_{\tau-1} - \eta_{\tau-1} \nabla f_{\tau-1}(X_{\tau-1})$
map back to the primal space $X'_{\tau} \gets (\nabla h)^{-1}(\theta_{\tau})$
project $X'_{\tau}$ onto a point $X_{\tau} \in K$ using the Bregman projection $X_{\tau} \gets \Pi_K^h(X'_{\tau})$

Bregman projection
<cit.> The Bregman projection of a point $x$ onto a convex set $K \subseteq \mathbb{R}^d$ given the distance-generating function $h$ is
\begin{align*}
\Pi_K^h(x) = \argmin_{y \in K} D_h(y,x).
\end{align*}

Note that when $h$ is the squared $\ell_2$ norm, the Bregman projection is equivalent to the Euclidean projection. Using finite difference methods, the Bregman projection can be computed similarly to the Euclidean projection, $\epsilon$-optimally in $\mathcal{O}(O_{\epsilon}^d)$ time, assuming the runtime of $h$ is constant.

Next, we describe how the ideas from mirror descent are adapted for the smoothed online convex optimization setting, where agents operate with lookahead $1$ and movement costs need to be considered.

§.§.§ Online Balanced Descent

The online balanced descent (OBD) algorithms that were developed by  <cit.> are a special case of online mirror descent (OMD) with lookahead $1$. In practice, to achieve the one-step lookahead, OBD moves to a point $X_{\tau}$ on a level set of $f_{\tau}(\cdot)$ such that the step is normal to the contour line of $f_{\tau-1}$. In contrast, OMD steps into a direction that is normal to the contour line of $f_{\tau-1}(X_{\tau-1})$. In other words, OMD takes a step with respect to its starting point on some level set, whereas OBD takes a step with respect to its destination on some level set. <Ref> shows how these updates compare.

Update of OMD. The contour lines represent level sets of $f_{\tau-1}$.

Update of OBD. The contour lines represent level sets of $f_{\tau}$.

Comparison of an update of OMD and OBD in two dimensions assuming the distance-generating function $h(x) = \frac{1}{2} \norm{x}_2^2$. OMD (red) takes a step in a direction normal to the contour line of $f_{\tau-1}$ at $X_{\tau-1}$. OBD (blue) takes a step in a direction normal to the contour line of $f_{\tau}$ at $X_{\tau}$ <cit.>. The step in the direction of the minimizer of $f_{\tau}$ is shown in black. Note that it is not optimal to move in the direction of the minimizer to a point on some level set as there likely exists a closer point on the same level set.

The algorithmic framework for OBD can roughly be divided into two parts. First, the projection of the previous point onto some level set of the cost function.  Second, the strategies to choose the specific level set and geometry to balance hitting costs and movement costs. The algorithm for (1) is also called the meta algorithm as it is parametrized with a concrete level set and geometry (i.e., mirror map).

§.§.§ Meta Algorithm

Similar to OMD, the meta algorithm of OBD chooses the next point $X_{\tau}$ in the dual space. However, whereas OMD takes an arbitrary step into the direction of the gradient, OBD takes the shortest step onto some sub-level set $K_l = \{x \in \mathcal{X} \mid f_{\tau}(x) \leq l\}$ of the revealed hitting cost $f_{\tau}$. In other words, we seek to find the Bregman projection of $X_{\tau-1}$ onto the sub-level set $K_l$. The first-order condition of the corresponding optimization in the dual space implies that
\begin{align}\label{eq:pbd:first_order_condition}
\nabla h(X_{\tau}) = \nabla h(X_{\tau-1}) - \eta_{\tau} \nabla f_{\tau}(X_{\tau})
\end{align}
must be satisfied by $X_{\tau}$ where $\eta_{\tau}$ is the optimal slack of the inequality constraint $f_{\tau}(x) \leq l$ <cit.>. Note that this corresponds to a variant of OMD with lookahead $1$.

OBD requires $h$ to be $\alpha$-strongly convex and $\beta$-Lipschitz smooth in the norm $\norm{\cdot}$ that is used to obtain the movement costs. The meta algorithm is described in <ref>.

Online Balanced Descent (meta algorithm) <cit.>
$\mathcal{I}_{\text{SCO}} = (\tau \in \mathbb{N}, \mathcal{X} \subset \mathbb{R}^d, \norm{\cdot}, (f_1, \dots, f_{\tau}) \in (\mathcal{X} \to \mathbb{R}_{\geq 0})^{\tau}), l \geq 0, \text{distance-generating function } h$
$X_{\tau} \gets \Pi_{K_l}^h(X_{\tau-1})$

As an example, we consider the Euclidean space with the $\ell_2$ norm. In this setting $h(x) = \frac{1}{2} \norm{x}_2^2$ is 1-strongly convex and 1-Lipschitz smooth <cit.>. As the corresponding mirror map, $\nabla h$ is the identity map, the first-order condition <ref> reduces to
\begin{align*}
X_{\tau} = X_{\tau-1} - \eta_{\tau} \nabla f_{\tau}(X_{\tau})
\end{align*}
corresponding to OGD with lookahead $1$.

We can generally choose $h$ to either perform well for the competitive ratio or regret. In their initial paper,  <cit.> propose two algorithms that balance hitting and movement costs in the primal and dual space and perform well concerning the competitive ratio and regret, respectively.

§.§.§ Primal Algorithm

Primal online balanced descent (P-OBD) balances hitting and movement costs in the primal space. Let $\hat{x} = \argmin_{x \in \mathcal{X}} f_{\tau}(x)$ and $x(l) = \text{Meta-OBD}(\mathcal{I}, l, h) = \Pi_{K_l}^h(X_{\tau-1})$. Given some $\beta > 0$, the balance parameter $l$ is chosen such that a balance condition $g(l) = \norm{x(l) - X_{\tau-1}} \leq \beta l$ is satisfied. More formally, $l$ is chosen such that either $x(l) = \hat{x}$ and $g(l) < \beta l$ or $g(l) = \beta l$ hold <cit.>.

<cit.> show that the balance function $g(l)$ is continuous in $l$. We observe that $l$ is lower bounded by $f_{\tau}(\hat{x})$. As we assume that $\hat{x}$ is unique, $x(l) = \hat{x}$ iff $l = f_{\tau}(\hat{x})$ and the first condition is fulfilled if and only if $g(f_{\tau}(\hat{x})) < \beta f_{\tau}(\hat{x})$. If the first condition is not satisfied, we can efficiently determine an $l$ fulfilling the second condition using a bracketed root finding method on $g(l) - \beta l$ within the interval $[f_{\tau}(\hat{x}), \gamma]$. Here, it suffices to choose $\gamma$ “large enough” to ensure that $g(\gamma) \leq \beta \gamma$. Observe that for $\gamma = f_{\tau}(X_{\tau-1})$, this is trivially satisfied, as $x(\gamma) = X_{\tau-1}$ and as such $g(\gamma) = 0$. The algorithm is described in <ref>.

Primal Online Balanced Descent <cit.>
$\mathcal{I}_{\text{SCO}} = (\tau \in \mathbb{N}, \mathcal{X} \subset \mathbb{R}^d, \norm{\cdot}, (f_1, \dots, f_{\tau}) \in (\mathcal{X} \to \mathbb{R}_{\geq 0})^{\tau}), \beta > 0, \text{distance-generating function } h$
$\hat{x} = \argmin_{x \in \mathcal{X}} f_{\tau}(x)$
$g(\hat{x}) \leq \beta f_{\tau}(\hat{x})$

$l \gets $ root of $g(l') - \beta l'$ for $l' \in [f_{\tau}(\hat{x}), f_{\tau}(X_{\tau-1})]$
$\text{Meta-OBD}(\mathcal{I}, l, h)$

The balancing is chosen such that the movement cost is upper bounded by the constant $\beta$ times the hitting cost <cit.>.

<cit.> show that P-OBD attains a competitive ratio of at most $3 + \mathcal{O}(1 / \alpha)$ for some $\beta > 0$, $\alpha$-locally polyhedral cost functions, and the $\ell_2$ norm as movement cost. They also show that local polyhedrality is useful when other norms like the $l_{\infty}$ norm are used. Note that the algorithm is memoryless and, therefore, nearly optimal, as  <cit.> showed in the uni-dimensional setting that no memoryless algorithm can attain a better competitive ratio than $3$ (which also holds for locally polyhedral cost functions) <cit.>. When used with the $\ell_1$ norm, which is relevant in the application of right-sizing data centers, P-OBD attains a competitive ratio of $\mathcal{O}(\sqrt{d})$ if $\alpha$ is fixed <cit.>. The memoryless algorithm (<ref>) of  <cit.> can be seen as a special case of P-OBD for $d = 1$ and $\beta = \frac{1}{2}$ <cit.>.

The minimizer can be found $\epsilon$-optimally in $\mathcal{O}(C O_{\epsilon}^d)$ time. Assuming the runtime of $h$ is constant, the balance parameter $l$ can be determined $\epsilon$-optimally in $\mathcal{O}(O_{\epsilon}^d R_{\epsilon})$ time. Overall, P-OBD runs in $\mathcal{O}(C O_{\epsilon}^d + O_{\epsilon}^d R_{\epsilon})$ time.

§.§.§ Dual Algorithm

<cit.> also developed an algorithm called dual online balanced descent (D-OBD), which balances the movement cost in the dual space with the gradient of the hitting cost (which is also in the dual space). Before describing the algorithm, we must first describe how the movement cost can be represented in the dual space. We thus introduce the notion of a dual norm.

dual norm
<cit.> Given some norm $\norm{\cdot}$ on $\mathbb{R}^d$, its dual norm $\norm{\cdot}_*$ is defined as
\begin{align*}
\norm{y}_* = \sup_{x \in \mathbb{R}^d} \{\langle x, y \rangle \mid \norm{x} \leq 1\}.
\end{align*}

The $\ell_2$ norm is self-dual <cit.>. In general, the dual norm for a concrete point $y \in \mathbb{R}^d$ can be computed with the following convex optimization:
\begin{align*}
&\max_{x \in \mathbb{R}^d} &&\langle x, y \rangle \\
&\text{subject to}         &&\norm{x} \leq 1.
\end{align*}

Returning to the description of D-OBD, for some fixed learning rate $\eta$, $l$ is now chosen such that
\begin{align*}
\norm{\nabla h(x(l)) - \nabla h(X_{\tau-1})}_* = \eta \norm{\nabla f_{\tau}(x(l))}_*
\end{align*}
holds. Let $g_1(l) = \norm{\nabla h(x(l)) - \nabla h(X_{\tau-1})}_*$ and $g_2(l) = \norm{\nabla f_{\tau}(x(l))}_*$. Again,  <cit.> show that the balance function $\frac{g_1(l)}{g_2(l)}$ is continuous in $l$ under the assumption that $h$ and $f_{\tau}$ are continuously differentiable on $\mathcal{X}$. Similar to our analysis of P-OBD, we observe  that $l$ is lower bounded by $f_{\tau}(\hat{x})$. We can determine $l$ using a bracketed root finding method on $g_1(l) - \eta g_2(l)$ within the interval $[f_{\tau}(\hat{x}), \gamma]$. We observe that for $l = f_{\tau}(\hat{x})$, $g_1(l) \geq \eta g_2(l) = 0$. Therefore, we need to choose $\gamma$ such that $g_1(l) \leq \eta g_2(l)$ is satisfied. Similar to our argument for P-OBD, it suffices to choose $\gamma = f_{\tau}(X_{\tau-1})$, resulting in $x(\gamma) = X_{\tau-1}$, and implying $g_1(\gamma) = 0$. The resulting algorithm is described in <ref>.

Dual Online Balanced Descent <cit.>
$\mathcal{I}_{\text{SCO}} = (\tau \in \mathbb{N}, \mathcal{X} \subset \mathbb{R}^d, \norm{\cdot}, (f_1, \dots, f_{\tau}) \in (\mathcal{X} \to \mathbb{R}_{\geq 0})^{\tau}), \eta > 0, \text{distance-generating function } h$
$\hat{x} = \argmin_{x \in \mathcal{X}} f_{\tau}(x)$
$l \gets $ root of $g_1(l') - \eta g_2(l')$ for $l' \in [f_{\tau}(\hat{x}), f_{\tau}(X_{\tau-1})]$
$\text{Meta-OBD}(\mathcal{I}, l, h)$

<cit.> show that the $L$-constrained dynamic regret of D-OBD is upper bounded by $\frac{G L}{\eta} + \frac{T \eta}{2 \alpha}$ where $h$ is $\alpha$-strongly convex in $\norm{\cdot}$, $\norm{\nabla h(x)}_*$ is upper bounded by $G$, and $\nabla h(0) = 0$. When $G$, $L$, and $T$ are known, $\eta$ can be chosen optimally as $\eta = \sqrt{\frac{2 G L \alpha}{T}}$, resulting in an $L$-constrained dynamic regret that is upper bounded by $\sqrt{\frac{2 G L T}{\alpha}}$ <cit.>. Further, in this setting, D-OBD achieves static regret $\mathcal{O}(\sqrt{T})$ <cit.>.

An evaluation of the dual norm can be computed $\epsilon$-optimally in $\mathcal{O}(O_{\epsilon}^d)$ time. Thus, assuming the runtime of $h$ is constant, $l$ can be found in $\mathcal{O}((O_{\epsilon}^d)^2 R_{\epsilon})$ time. Overall, the asymptotic time complexity of D-OBD is given as $\mathcal{O}(C O_{\epsilon}^d + (O_{\epsilon}^d)^2 R_{\epsilon})$.

§.§.§ Greedy and Regularized Algorithms

Later,  <cit.> proposed two additional algorithms using the OBD framework, greedy online balanced descent (G-OBD) and regularized online balanced descent (R-OBD). Both algorithms yield strong guarantees for the competitive ratio in the setting of squared $\ell_2$ norm movement costs and $\alpha$-strongly convex hitting costs where the optimal competitive ratio is $\mathcal{O}(1 / \sqrt{\alpha})$ as $\alpha$ approaches zero. G-OBD achieves this competitive ratio for quasiconvex[A function is quasiconvex iff it has a unique global minimum.] hitting costs that are $\alpha$-strongly convex around their minimizer and squared $\ell_2$ norm movement costs. R-OBD achieves this competitive ratio for $\alpha$-strongly convex hitting costs and arbitrary Bregman divergences as movement costs.

Both algorithms take an additional step of size $\mathcal{O}(\sqrt{\alpha})$ towards the minimizer of the hitting cost. G-OBD works by first taking a regular P-OBD step to some level set of the hitting cost. Then, it takes an additional step towards the minimizer of the hitting cost with a step size based on the convexity parameter $\alpha$. In contrast, R-OBD picks the next point by minimizing a weighted sum of hitting and movement costs. It uses an additional regularization term that encourages the algorithm to pick a point closer to the minimizer of the hitting cost <cit.>. We do not discuss G-OBD and R-OBD in more detail, as their theoretical guarantees do not cover the application of right-sizing data centers, but we provide implementations of them.

§ PREDICTING

In practice, we can attempt to use predicted hitting costs to improve the performance of online algorithms. Predicting future incoming loads to a high degree of accuracy in the data-center setting is often possible. Using predicted hitting costs and their uncertainty distributions, online algorithms can make more informed decisions in practice. In this section, we begin by discussing prediction windows. Then, we describe approaches for time-series predictions and end with discussing algorithms that use such predictions.

§.§ Prediction Window

A natural model to allow incorporating predictions is the use of a finite prediction window $w$. A prediction window bridges the gap between offline and online algorithms. Whereas an online algorithm only knows the hitting costs $f_t$ for $t \in [\tau]$ and an offline algorithm knows the hitting costs $f_t$ for all $t \in [T]$, an online algorithm with prediction windowprediction window of length $w$ knows all hitting costs $f_t$ up to $\tau + w$, i.e. $t \in [\tau + w]$. In other words, the prediction window $w$ represents the number of upcoming time slots at which the algorithm is assumed to have perfect knowledge of the future.

§.§.§ Lazy Capacity Provisioning with Prediction Window

<cit.> extend their algorithm lazy capacity provisioning, which we discussed in <ref> to support the prediction window by changing the update rule to
\begin{align*}
X_{\tau} = \begin{cases}
0 & \tau \leq 0 \\
(X_{\tau-1})_{X_{\tau+w,\tau}^L}^{X_{\tau+w,\tau}^U} & \tau \geq 1
\end{cases}
\end{align*}

The optimal schedules now need to obtained for $\tau + w$ rather than $\tau$ time slots. Thus, the time complexity changes to $\mathcal{O}((\tau + w) C O_{\epsilon}^{\tau + w})$ and $\mathcal{O}((\tau + w)^2 C \log_2 m)$ in the fractional and integral case, respectively.

The assumption of perfect knowledge of the future is sure to be violated when an online algorithm is used in practice. Still,  <cit.> show that lazy capacity provisioning with a prediction window is robust to this assumption in practice.  <cit.> and  <cit.> showed that using a finite prediction window does not improve the worst-case performance of the online algorithm for the fractional and integral case, respectively. In other words, the competitive ratio of lazy capacity provisioning is $3$ regardless of whether it uses a finite prediction window. In practice, however,  <cit.> show that a prediction window significantly improves the algorithm's performance.

There are two main drawbacks to using a finite prediction window. First, predictions windows are finite and typically constrained to a short period as they are assumed to be perfect. In contrast, predictions can be made for much longer time horizons, albeit with decreasing accuracy. Second, it completely disregards any knowledge or assumptions of the certainty and noise of the predictions by assuming the predictions to be perfect.

§.§ Making Predictions

There exist multiple paradigms for making time-series predictions. Due to much recent engagement in the field of deep learning generally and time-series predictions specifically, multiple approaches perform well in practical settings. Most algorithms separately tune parameters of individual models for short-term and long-term trends as well as seasonality <cit.>.

A fundamental difference between models is Bayesianness, i.e., whether they use an underlying uncertainty distribution within the model. Facebook's Prophet algorithm is Bayesian, whereas LinkedIn's Greykite algorithm is not <cit.>.

For Bayesian models, online algorithms can use the uncertainty distribution to consider outliers appropriately. For non-Bayesian models, additive white Gaussian noise can be added to the prediction to achieve a similar effect. In general, many strategies can be used to obtain a single representative prediction of the underlying distribution. In our experiments, we use the mean prediction to ensure appropriate consideration of outliers. The median or 90th percentile predictions are alternatives that are more robust to outliers.

Note that, in principle, predictions can be made arbitrarily far into the future. However, at some point, they become too uncertain to be valuable. For example, infeasible load profiles may be assigned a positive probability, which would result an infinite cost if we use the mean to obtain a representative prediction, even if all servers are active. Thus, we also use a prediction window, which needs to be set appropriately to account for the uncertainty distribution of the predicted loads.

§.§ Receding Horizon Control

Receding horizon controlreceding horizon control (RHC) (or model predictive control) is a methodology for making decisions based on predictions of the future that is commonly used to control data centers <cit.>. In RHC, an agent predicts their action up to some fixed point in time, referred to as the prediction window. Based on this prediction, the agent adjusts their action for the current time slot. In the next time slot, this process repeats <cit.>.

<cit.> previously investigated the performance of RHC in the context of right-sizing data centers. RHC works by solving a convex optimization from time $\tau$ to time $\tau + w$ starting from the initial configuration $X_{\tau-1}$. We set $X_0 = \mathbf{0}$. Similar to our analysis of capacity provisioning in <ref>, we describe by $X^{\tau}(X_{\tau-1}) \in \mathcal{X}^{w+1}$ the optimal schedule for times $\tau$ through $\tau+w$. This schedule is obtained by minimizing
\begin{align}\label{eq:rhc}
\sum_{t=\tau}^{\tau+w} f_t(X_t) + \norm{X_t - X_{t-1}}
\end{align}
over configurations $X_{\tau}, \dots, X_{\tau+w} \in \mathcal{X}$. This optimization has $\mathcal{O}(d w)$ dimensions and thus can be computed $\epsilon$-optimally in $\mathcal{O}(C O_{\epsilon}^{dw})$ time. Now, RHC simply picks the first predicted action. RHC is described in <ref>.

Receding Horizon Control <cit.>
$\mathcal{I}_{\text{SSCO}} = (\tau \in \mathbb{N}, m \in \mathbb{N}, \beta \in \mathbb{R}_{>0}, (f_1, \dots, f_{\tau}) \in (\mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0})^{\tau})$
$X_{\tau} = X_{\tau}^{\tau}(X_{\tau-1})$

<cit.> prove the competitive ratio of RHC in the application of right-sizing data centers. They show that in the uni-dimensional setting, RHC attains a competitive ratio of $1 + \mathcal{O}(1/w)$ which is strictly better than the optimal competitive ratio (for deterministic algorithms without predictions) of $2$ and $3$ for memoryless algorithms for $w > 1$ and $w > \frac{1}{2}$, respectively. However, in a multi-dimensional setting, RHC is $(1 + \max_{k \in [d]} \beta_k / e_k(0))$-competitive where we defined $\beta_k$ as the switching cost of a server of type $k$ and $e_k(0)$ as the average energy cost of an idling server of type $k$. Importantly, this competitive ratio does not depend on the size of the prediction window $w$.

§.§ Averaging Fixed Horizon Control

In their paper,  <cit.> present another algorithm, averaging fixed horizon control (AFHC), which attains a competitive ratio of $1 + \max_{k \in [d]} \frac{\beta_k}{(w+1) e_k(0)}$. In particular, AFHC is $(1 + \mathcal{O}(1/w))$-competitive. However,  <cit.> find that in many realistic settings, RHC performs better than AFHC.

At time $\tau$, AFHC works by performing $w + 1$ individual RHC steps starting from $t_0 = \tau-w$ up to $t_0 = \tau$ and averaging the results. Each individual step is also referred to as an iteration of fixed horizon control (FHC).

We describe the sub-iterations of AFHC using $k \in [w+1]$. We set $t_0 = \tau+k-(w+1)$, ensuring that $t_0 \in [\tau-w,\tau]$. We denote by $X^{t_0}(X_{t_0-1}^{(k)})$ the optimal schedule for times $t_0$ through $t_0+w$ which is obtained analogously to <ref>. We also set $X_t = \mathbf{0}$ and $X_t^{(k)} = \mathbf{0}$ for all $t \leq 0$ and $k \in [w+1]$. AFHC is described in <ref>.

Averaging Fixed Horizon Control <cit.>
$\mathcal{I}_{\text{SSCO}} = (\tau \in \mathbb{N}, m \in \mathbb{N}, \beta \in \mathbb{R}_{>0}, (f_1, \dots, f_{\tau}) \in (\mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0})^{\tau})$
$k \in [w+1]$
$t_0 \gets \tau+k-(w+1)$
$X^{(k)} \gets X^{t_0}(X_{t_0-1}^{(k)})$

$X_{\tau} = \frac{1}{w+1} \sum_{k=1}^{w+1} X_{\tau}^{(k)}$

Intuitively, AFHC can be interpreted as performing $w+1$ FHC-steps in parallel, where each FHC-step starts from a different $t_0 \in [\tau-w,\tau]$, and then averaging all configurations for time $\tau$. Note that RHC is equivalent to the last FHC-step with initial time $t_0 = \tau$, i.e., $k = w+1$. The asymptotic time complexity of AFHC is given as $\mathcal{O}(w C O_{\epsilon}^{dw})$.

<cit.> show that AFHC achieves sublinear regret and a constant competitive ratio using a prediction window of constant length.  <cit.> introduce a class of “forward-looking” algorithms that can consider cost functions within some prediction window but are only allowed to use a constant limited number of past cost functions. They show that among these algorithms, AFHC achieves optimal regret.

In <cit.>,  generalize RHC and AFHC to a class of algorithms called committed horizon control (CHC), which consist of $v \in [w+1]$ sub-iterations of FHC. Note that RHC corresponds to CHC with parameter $v = 1$, whereas AFHC corresponds to CHC with parameter $v = w+1$. They investigate how $v$ can be chosen optimally based on the noise distribution of predictions.

<cit.> extend AFHC to a new algorithm called synchronized fixed horizon control (SFHC), which is $(1 + \mathcal{O}(1/w))$-competitive for both convex and non-convex cost functions.  <cit.> propose two new gradient-based online algorithms, receding horizon gradient descent (RHGD) and receding horizon accelerated gradient (RHAG), and show that the dynamic regret of RHAG is near-optimal when compared to a class of online algorithms that includes CHC.

CHAPTER: IMPLEMENTATION

As part of our work, we implemented the algorithms described in <ref> and <ref>. Our implementation is written in Rust and has Python bindings to interact with some components <cit.>. Detailed documentation of the implementation is available <cit.>. In this chapter, we discuss the general architecture and the points of focus of our implementation.

§ ARCHITECTURE

Our implementation is separated into four main components. First, a component that encompasses the implementation of the general problems described in <ref>. Second, a related component containing the implementations of the offline and online algorithms from <ref> and <ref>, respectively. Third, a component including abstractions to generate and update problem instances over time (referred to as modelsmodels). This includes the models of data centers described in <ref>. Fourth, a component that implements the algorithms' practical use. Among other things, this component is used to execute online algorithms in real-time by sequentially updating the underlying problem instance as new information becomes available (referred to as streamingstreaming an online algorithm).

§.§ Problems and Reductions

This component includes data structures representing instances of the problems introduced in <ref>. Note-worthy aspects are the definition of norms for SCO and hitting costs for SCO, SSCO, and SBLO. We also implemented the reductions from SLO to SBLO, SBLO to SSCO, and SSCO to SCO.

Norms Our implementation includes the Manhattan norm, Euclidean norm, Mahalanobis norm, and a dimension-dependently scaled variant of the Manhattan norm, which is used in the reduction from SSCO to SCO. We also include implementations of the square of a norm as well as the dual of a norm.

Hitting Costs The primary data structure used across SCO, SSCO, and SBLO is the definition of hitting costs. Here, our implementation has to allow for sequentially arriving subsets of the hitting cost. Consider the example of an online algorithm with a prediction window $w$. Then, at each time slot $\tau$, the algorithm expects to receive a hitting cost with domain $[\tau : \tau + w]$ which conflicts with previous and subsequent hitting costs. To complicate the matter further, offline algorithms only require a single hitting cost covering the entire domain $[T]$.

In our implementation, we refer to these hitting costs, which may only cover definitions of the hitting cost for a subset of all time slots as single hitting costssingle hitting cost. In the online setting, the domain of a single hitting cost may cover future time slots, in which case the hitting cost is uncertain. From our discussion of predictions in <ref>, it naturally follows to model the distribution of hitting costs for a future time slot as a vector of samples that can then be used to estimate the density. Thus, a single cost function arriving at time slot $\tau$ is formally described as a function $\chi_{\tau} : [T] \times \mathcal{S} \to \bigcup_{n=1}^{\infty} \mathbb{R}_{\geq 0}^n$ mapping a time slot $t \in [T]$ and some point $x$ in the space $\mathcal{S}$ to a prediction with a varying sample size $n$. In the case of a certain prediction, the sample size is $1$. For SCO and SSCO we set $\mathcal{S} = \mathcal{X}$ while for SBLO we set $\mathcal{S} = \mathbb{R}$.

We implement the hitting cost $f_t(x)$ using a B-Tree-Map[A B-Tree-Map is a map based on a B-Tree balancing cache-efficiency and search performance <cit.>] mapping the time slot of their arrival to single hitting costs. Crucially, there need not be a single hitting cost for every time slot, as the example of an offline algorithm illustrates. The concrete single hitting cost is determined by choosing the last single hitting cost arriving during a time slot in $[t]$. In an online setting, this ensures that if $t$ is in the future, the current single hitting cost is used, and if $t$ is in the past, the single hitting cost from time slot $t$ (or the closest previous single hitting cost) is used.

§.§ Algorithms

The algorithms component encompasses the implementation of all offline and online algorithms described in <ref> and <ref>. For both classes of algorithms, we define a common interface.

Offline An offline algorithm receives as input a problem instance and some algorithm-specific options. It returns an arbitrary data structure that can be used to obtain the determined schedule. In our implementation, we generalized the algorithms slightly so as to support the following uses:

First, we adapted the graph-based algorithms for SSCO to support inverted movement costs, i.e., paying the switching cost for shutting down a server rather than powering up a server. Second, we adapted all algorithms to allow for computing the $\alpha$-unfair optimal offline solution where movement costs are scaled by a factor $\alpha$. Third, we adapted the general algorithm for the multi-dimensional fractional case to support computing the $L$-constrained optimal offline solution where movement costs are upper-bounded by the constant $L$. The $L$-constrained optimal offline solution cannot easily be computed using the graph-based approaches for the integral case as this additional constraint prevents the application of Bellman's optimality principle to find a shortest path using dynamic programming. Fourth, we extended the implementations to allow starting from an arbitrary initial time slot. Fifth, we adapted the uni-dimensional optimal graph search (<ref>) to support arbitrary initial configurations in the first time slot.

Some of these extensions are required to implement some online algorithms, and others were merely added for analysis purposes. We also provide functions to compute the static fractional and integral optima.

Online An online algorithm receives as input a problem instance, the current time slot $\tau$, the schedule for all previous time slots, some algorithm-specific memory, and some algorithm-specific options. The algorithm returns the configuration for time slot $\tau$ as well as the updated memory. The practical use of online algorithms is described in greater detail in <ref>.

§.§ Models

In practice, it is not very useful to generate and update the problem instances directly. In the application of right-sizing data centers, we have discussed models of operating costs and switching costs in detail in <ref>. These models can be used to generate instances of the problems discussed in <ref>. This approach is not limited to the application of right-sizing data centers. A model is a data structure with associated functions to produce an associated problem instance and update an existing problem instance online. As seen in the example of right-sizing data centers, these generators may require additional inputs, which are provided as the load profiles for each time slot. We refer to the inputs required for the initial generation of a problem instance as offline inputsoffline inputs and the inputs required to update a problem instance online as online inputsonline inputs. Offline inputs encapsulate information for all time slots with respect to some time horizon $[T]$ whereas online inputs for some time slot $\tau$ encapsulate information for the current time slot $\tau$ as well as all time slots in the prediction window $[\tau + 1 : \tau + w]$.

In the application of right-sizing data centers, an offline input is a vector of load profiles $\mathcal{I} = (\lambda_1, \dots, \lambda_T)$ with $\lambda_t \in \mathbb{N}_0^e$ for all $t \in [T]$. In contrast, an online input is a vector of predicted load profiles $\mathcal{I} = (\lambda_{\tau}, \mathcal{P}_{\tau + 1}, \dots, \mathcal{P}_{\tau + w})$ with $\lambda_{\tau} \in \mathbb{N}_0^e$ and where $\mathcal{P}_t \in \left(\bigcup_{n=1}^{\infty} \mathbb{N}_0^n\right)^e$ is a predicted load profilepredicted load profile, i.e., a vector of predicted loads for each job type. Similar to our implementation of single hitting costs, we use a varying number of $n$ samples to describe the distribution of predicted loads. The load profile for the current time slot $\tau$ is certain.

We observe that given a predicted load profile $\mathcal{P}_t$ with $n_i$ samples of loads with type $i \in [e]$, the sampled load profiles are all combinations of the samples of individual loads, resulting in $\prod_{i=1}^e n_i$ sampled load profiles. To reduce the number of considered samples, we randomly select $n$ sampled loads of each type to produce $n$ load profile samples. This approach is reasonable as the samples of the individual loads were only an approximate prediction, to begin with.

Another core element of a model are model outputsmodel outputs. These are results that are returned from the model, which reach beyond simply some cost. For our data center model, we implemented outputs to return the energy cost instead of the revenue loss and the determined optimal assignment of jobs to server types.

The implementation of our data center model consists of separate models for energy consumption, energy cost, delay, revenue loss, and switching cost, which can be combined in various ways to obtain a concrete model of a data center or a network of data centers. For the generation of instances of SBLO, the model must be limited to a single location, source, and job type. For the generation of instances of SLO, the model further assumes full utilization of every active server during each time slot and averages the energy cost over the time horizon.

Crucially, we assume in our model that job arrivals fall precisely onto the beginning of new time slots. Without this approximation, we cannot ensure that all arriving jobs can be completed by the end of the time slot. In practice, this requires matching jobs to time slots either by delaying them to the next time slots (in which case they are available from the beginning) or simply assuming they are available from the beginning of the current time slot.

§.§ Streaming and Practical Use

Using an offline algorithm is relatively straightforward. Given a model and offline inputs, we can generate a problem instance that the offline algorithm can then solve. With online algorithms, this becomes more involved. For example, some uses require the immediate streaming of an online algorithm from an initial to a final time slot. This is required by the Lazy Budgeting algorithm (<ref>), which streams an intermediate online algorithm for the sub time slots of the current time slot.

However, in general, the challenge in streaming an online algorithm is to remember the algorithm's state with infrequent incoming iterations (once per time slot). To this end, we implemented a simple client-server architecture. Here, the server runs in the background and remembers the problem instance and state of the algorithm. Whenever new information arrives, the client requests the next iteration from the server. On initialization of the server, the online algorithm can be streamed up to the current time slot given previous offline inputs. This architecture has two main benefits: First, the responsibility for memorizing the problem instance and state of the algorithm is offloaded entirely to the server. Second, the initialization of the server only requires a model and (optionally) an offline input. Requesting the next iteration only requires the relevant online input. In particular, the problem instance remains opaque, drastically simplifying the interaction with the interface and reducing the amount of data that needs to be remembered, and serialized, sent, and deserialized from the client to the server. This simplification also makes it feasible to stream an online algorithm relying entirely on Python bindings.

§ FOCUSES

We now turn to describe more general aspects of our implementation.

§.§ Adaptability

A central point of focus was to maximize the adaptability and general applicability of our implementation. More concretely, we focused on minimizing the number of cross-dependencies between the components described in <ref> to maximize the reusability across all components of the architecture. For example, an entirely new application can be supported solely by implementing a new model. Once a new algorithm is implemented, it can immediately be substituted for any other algorithm using any previously implemented model while benefiting from the infrastructure built around streaming online algorithms. This allows for faster prototyping of new algorithms while assessing their practical performance, particularly when using predictions. Also, this allows determining which algorithm performs best in a concrete application instance.

Moreover, using the client-server architecture, algorithms for smoothed online convex optimization can be utilized without much effort in any practical application. As the implementation is in Rust, it can be interfaced with any other programming language. For our case studies in <ref>, we already implemented Python bindings <cit.>.

§.§ Efficiency

Next to the strong safety guarantees of Rust, the main benefit of an implementation in Rust is its memory efficiency and speed. Rust achieves similar performance to C or C++ <cit.>. At the same time, Rust is readable and allows for high-level abstractions and polymorphism, which we use heavily to maximize adaptability.

In our implementation, we heavily parallelize tasks to achieve optimal performance <cit.>. For example, we parallelize calculating hitting costs for multiple predicted samples or determining the optimal predecessor of a vertice in iterations of the dynamic program finding a shortest path in a graph.

§.§ Numeric Computations and Rounding

In implementations involving numeric computations to some precision $\epsilon$ and frequent flooring or ceiling operations, it is essential to round numeric results to precision; otherwise, the results are not numerically stable. In our experimens, we use a precision of $\epsilon = 10^{-2}$, which performs well when results are rounded to integral solutions. For example, if the result of some numeric computation with precision epsilon is $10^{-3}$, which is then ceiled, we would falsely obtain $1$ as the result if we did not apply the precision to the result before proceeding with the algorithm.

For convex programs, we additionally interpret $\epsilon$ as relative to the absolute value of the optimization to maintain good performance.

CHAPTER: CASE STUDIES

This chapter examines the performance of the previously described models and algorithms using real server traces. Thereby we focus on two aspects: First, we are interested in how well the discussed algorithms compare in absolute terms and relative to each other. Second, we are interested in the general promise of dynamically right-sizing data centers, which we study by conservatively estimating cost savings and relating them to previous research.

§ METHOD

First, we describe our experimental setup. We begin with a detailed discussion of the characteristics of the server traces, which we use as a basis for our analysis. Then, we examine the underlying assumptions of our analysis. This is followed by a discussion of alternative approaches to right-sizing data centers, which we use as a foundation for estimating the cost savings resulting from dynamic right-sizing of data centers. Next, we describe the general model parameters we use in our analysis and relate them to previous research. Lastly, we introduce the precise performance metrics used in the subsequent sections.

Throughout our experiments, we seek to determine conservative approximations for the resulting performance and cost savings. Our experimental results were obtained on a machine with 16 GB memory and an Intel Core i7-8550U CPU with a base clock rate of 1.80GHz.

§.§ Traces

We use several traces with varying characteristics for our experiments. Some traces are from clusters rather than individual data centers. However, to simplify our analysis, we assume traces apply to a single data center without restricting the considered server architectures.

<cit.> showed that the characteristics of traces vary drastically even within a single trace when different subsets are considered individually. Their observation shows that it is crucial to examine long server traces and various server traces from different sources to gain a proper perspective of real-world performance.

§.§.§ Characteristics

To understand the varying effectiveness of dynamic right-sizing for the considered traces, we first analyze the properties of the given traces.

The most immediate and fundamental properties of a trace are its duration, the number of appearing jobs, the number of job types, and the underlying server infrastructure – especially whether this infrastructure is homogeneous or heterogeneous.

Then, we also consider several more specific characteristics. The interarrival timeinterarrival time (or submission rate) of jobs is the distribution of times between job arrivals. This distribution indicates the average system load as well as load uniformity. The peak-to-mean ratio (PMR)peak-to-mean ratio is defined as the ratio of the maximum load and the mean load. It is a good indicator of the uniformity of loads. We refer to time slots as peakspeak load when their load is greater than the 0.9-quantile of loads. We call the ratio of the 0.9-quantile of loads and the mean load true peak-to-mean-ratio (TPMR)true peak-to-mean ratio as it is less sensitive to outliers than the PMR. We refer to periods between peaks as valleysvalley. More concretely, we refer to the time between two consecutive peaks as peak distancepeak distance and the number of consecutive time slots up to a time slot with a smaller load as valley lengthvalley length. Further, we say that a trace follows a diurnal patterndiurnal pattern if during every 24 hours, excluding the final day, there is at least one valley spanning 12 hours or more. We exclude the final day as the final valley might be shortened by the end of the trace.

We also consider some additional information included in some traces, such as the measured scheduling rate (or queuing delay), an indicator for utilization.

§.§.§ Overview

We now give an overview of all used traces. For our initial analysis, we use a time slot length of 10 minutes.

MapReduce[MapReduce is a programming model for processing and generating large data sets in a functional style <cit.>] Workload from a Hadoop[Apache Hadoop is an open-source software for managing clusters] Cluster at Facebook <cit.> This trace encompasses three day-long traces from 2009 and 2010, extracted from a 6-month and a 1.5-month-long trace containing 1 million homogeneous jobs each. The traces are visualized in <ref> and summarized in <ref>. The cluster consists of 600 machines which we assume to be homogeneous. For the trace from 2010, we adjust the maximum number of servers to 1000 as otherwise the trace is infeasible under our models. <Ref> visualizes the corresponding dynamic and static offline optimal schedules under our second model (which is described in <ref>). The trace was published by  <cit.> as part of the SWIM project at UC Berkeley.

Facebook MapReduce workloads.

Optimal dynamic and static offline schedules for the last day of the Facebook workloads. The left y axis shows the number of servers of the static and dynamic offline optima at a given time (black). The right y axis shows the number of jobs (i.e., the load) at a given time (red).

characteristic     2009-0     2009-1     2010
duration     1 day     1 day     1 day

number of jobs     6 thousand     7 thousand     24 thousand

median interarrival time     7 seconds     7 seconds     2 seconds

PMR     3.91     2.97     2.2

TPMR     2.04     1.93     1.69

mean peak distance     95 minutes     106 minutes     87 minutes

mean valley length     44 minutes     36 minutes     35 minutes

Characteristics of Facebook's MapReduce workloads.

Los Alamos National Lab HPC Traces <cit.> This trace comprises two separate traces from high-performance computing clusters from Los Alamos National Lab (LANL). The traces were published by  <cit.> as part of the Atlas project at Carnegie Mellon University.

The first trace is from the Mustang cluster, a general-purpose cluster consisting of 1600 homogeneous servers. Jobs were assigned to entire servers. The dataset covers 61 months from October 2011 to November 2016 and is shown in <ref>. Note that the PMR is large at $622$ due to some outliers in the data. The median job duration is roughly 7 minutes, although the trace includes some extremely long-running outliers, resulting in a mean job duration of over 2.5 hours. In the trace, jobs were assigned to one or multiple servers. To normalize the trace, we consider each job once for each server it was processed on. <Ref> shows the dynamic and static offline optimal schedules under our second model.

The second trace is from the Trinity supercomputer. This trace is very similar to the Mustang trace but includes an even more significant number of long-running jobs. We, therefore, do not consider this trace in our analysis.

Microsoft Fiddle Trace <cit.> This trace consists of deep neural network training workloads on internal servers from Microsoft. The trace was published as part of the Fiddle project from Microsoft Research. The jobs are run on a heterogeneous set of servers which we group based on the number of GPUs of each server. There are 321 servers with two GPUs and 231 servers with eight GPUs. The median job duration is just below 15 minutes. The load profiles are visualized in <ref>.

The CPU utilization of the trace is extremely low, with more than 80% of servers running with utilization 30% or less <cit.>. However, memory utilization is high, with an average of more than 80% indicating that overall server utilization is already very high <cit.>. Again, the PMR is rather large at 89.43 due to outliers.

In our model, we adjust the runtime of jobs relative to the number of available GPUs in the respective server, i.e., the average runtime of jobs on a 2-GPU-server is four times as long as the average runtime of jobs on an 8-GPU-server. We adjust for the increased energy consumption of a server with eight GPUs by increasing the energy consumption of servers with two GPUs by a factor of 4.2. We also associate a fifteen times higher switching cost with servers with eight GPUs.

The dynamic and static offline optimal schedules under our second model are shown in <ref>. Note that under the given load servers with two GPUs are preferred to servers with eight GPUs when they are only needed for a short period due to their lower switching costs. This might seem counterintuitive at first, as 2-GPU-servers seem to be strictly better than 8-GPU-servers as the operating and switching cost of 8-GPU-servers is worse by a factor greater than four than the respective cost of 2-GPU-servers. However, we assume an average job runtime of 7.5 minutes on 8-GPU-servers as opposed to an average job runtime of 30 minutes on 2-GPU-servers (a factor of four), implying that 8-GPU-servers can process more than four jobs in an hour without a significant increase in delay, whereas 2-GPU-servers are limited to one job per time slot.

Alibaba Trace <cit.> This trace consists of a mixture of long-running applications and batch jobs. We are using their trace from 2018, covering eight days. The trace is visualized in <ref>, the dynamic and static offline optimal schedules under our second model are shown in <ref>. The jobs are processed on 4000 homogeneous servers. In our models, we assume a total of 10,000 servers to ensure that the number of servers is not a bottleneck. Jobs themselves are grouped into 11 types which we further simplify to 4 types based on their average runtime. We consider short, medium, long, and very long jobs. Their average runtime in the trace is shown in <ref>. The mean job duration is just below 15 minutes. The median job duration is 8 seconds, and the mean job duration is just over 1.5 minutes.

job type     mean runtime
short     68 seconds

medium     196 seconds

long     534 seconds

very long     1180 seconds

Characterization of the job types of the Alibaba trace.

Data from a previous trace indicates that mean CPU utilization varies between 10% and 40% while mean memory utilization varies between 40% and 65% <cit.>. This indicates that the overall server utilization is not optimal.

In our model, we scale job runtimes by a factor of 2.5 from short to very long jobs, roughly matching the runtimes of jobs from the trace.

LANL Mustang

Microsoft Fiddle

LANL Mustang, Microsoft Fiddle, and Alibaba traces. The figures display the average number of job arrivals throughout a day. The interquartile range is shown as the shaded region.

LANL Mustang

Microsoft Fiddle

Optimal dynamic and static offline schedules for the last day of the LANL Mustang, Microsoft Fiddle, and the second to last day of the Alibaba trace. The left y axis shows the number of servers of the static and dynamic offline optima at a given time (black). The right y axis shows the number of jobs (i.e., the load) at a given time (red).

We have seen traces from very different real-world use cases. The Microsoft Fiddle trace is based on a heterogeneous server architecture, and the Alibaba trace receives heterogeneous loads. The PMR and valley lengths of the days used in our analysis are shown in <ref>. Interestingly, as shown in <ref>, their TPMR, peak distances, and valley lengths are mostly similar.

§.§ Assumptions

We impose a couple of assumptions to simplify our analysis. First, and already mentioned, our analysis is inherently limited by the traces we used as a basis for our experiments. While we examine a wide variety of traces, the high variability in traces indicates they are a fundamental limitation to any estimation of real-world performance.

Another common limitation of models is that the interarrival times of jobs are on the order of seconds or smaller <cit.>. However, this is not a limitation of our analysis as we are using a general Poisson process with an appropriate mean arrival rate in our delay model.

In the context of high-performance computing, jobs typically have a gang schedulinggang scheduling requirement, i.e., a requirement that related jobs are processed simultaneously even though they are run on different hardware <cit.>. For simplification, we assume this requirement always to be satisfied. However, this is not a substantial limitation as the scheduling of jobs within a time slot is not determined by the discussed algorithms and instead left to the server operator. Nevertheless, in principle, the gang scheduling requirement may render some schedules infeasible if the processing time on servers exceeds the length of a time slot when gang scheduling constraints are considered.

There are also some limitations resulting from the design of our model. As was mentioned previously, we assume that the jobs arrive at the beginning of a new time slot rather than at random times throughout the time slot. Moreover, we assumed that for every job, a server type exists that can process this job within one time slot. In other words, there exists no job running longer than $\delta$. We have seen in <ref> that this assumption is violated in most practical scenarios. In <ref>, we described how this assumption can be removed. The same approach can also be used to remove the assumption that jobs must arrive at the beginning of a time slot.

§.§ Alternatives to Right-Sizing Data Centers

To determine the benefit of dynamically right-sizing data centers, we must first describe the alternative strategies to managing a data center. We will then use these approaches as a point of reference in our analysis.

Most data centers are statically provisioned; that is, the configuration of active servers is only changed rarely (often manually) and remains constant during most periods <cit.>. To support the highest loads, the data centers are peak-provisioned, i.e., the number of servers is chosen such that they suffice to process all jobs even during times where most jobs arrive. Moreover, as a safety measure, data centers are typically provisioned to handle much higher loads than the loads encountered in practice <cit.>.

characteristic     LANL Mustang     Microsoft Fiddle     Alibaba
duration     5 years     30 days     8 days

number of jobs     20 million     120 thousand     14 million

median interarrival time     0 seconds     8 seconds     0 seconds

PMR     621.94     89.43     3.93

TPMR     2.5     1.68     1.77

mean peak distance     100 minutes     105 minutes     89 minutes

mean valley length     120 minutes     115 minutes     74 minutes

diurnal pattern     yes     -     yes

Characteristics of the LANL Mustang, Microsoft Fiddle, and Alibaba traces.

Naturally, traces with a high PMR or long valleys are more likely to benefit from alternatives to static provisioning. Therefore another widely used alternative is valley fillingvalley filling, which aims to schedule lower priority jobs (i.e., some batch jobs) during valleys. In an ideal scenario, this approach can achieve $\text{PMR} \approx 1$, which would allow for efficient static provisioning. Crucially, this approach requires a large number of low-priority jobs which may be processed with a significant delay (requiring a considerable minimum perceptible delay $\delta_i$ for a large number of jobs of type $i$), and thus in most cases, valleys cannot be eliminated entirely.  <cit.> showed that dynamic-right sizing can be combined with valley filling to achieve a significant cost reduction. The optimal balancing of dynamic right-sizing and valley filling is mainly determined by the change to the PMR.  <cit.> showed that cost savings of 20% are possible with a PMR of 2 and a PMR of approximately 1.3 can still achieve cost savings of more than 5%. Generally, the cost reduction vanishes once the PMR approaches $1$, which may happen between 30% to 70% mean background load <cit.>. The results when dynamic right-sizing is used together with valley filling can be estimated from previous results.

§.§ Performance Metrics

Let $OPT$ denote the dynamic offline optimum and $OPT_s$ denote the static offline optimum. In our analysis, the normalized costnormalized cost of an online algorithm is the ratio of the obtained cost and the dynamic optimal offline cost, i.e. $NC(ALG) = c(ALG) / c(OPT)$. Further, we base our estimated cost reductioncost reduction on an optimal offline static provisioning:
\begin{align*}
CR(ALG) = \frac{c(OPT_s) - c(ALG)}{c(OPT_s)}.
\end{align*}
Note that this definition is similar to the definition of regret, but expressed relative to the overall cost. We refer to $SDR = c(OPT_s) / c(OPT)$ as the static/dynamic ratiostatic/dynamic ratio, which is closely related to the potential cost reductionpotential cost reduction $PCR = CR(OPT)$.

§.§ Previous Results

<cit.> showed that the cost reduction is directly proportional to the PMR and inversely proportional to the normalized switching cost. Additionally,  <cit.> showed that, as one would expect, the possible cost reduction decreases as the delay cost assumes a more significant fraction of the overall hitting costs. In practice, this can be understood as the effect of making the model more conservative.

trace     PMR     mean valley length (hours)
Facebook 2009-0     2.115     2.565

Facebook 2009-1     1.913     1.522

Facebook 2010     1.549     1.435

LANL Mustang     6.575     1.167

Microsoft Fiddle     3.822     2.125

Alibaba     1.339     2.792

PMR and mean valley length of the traces used in our analysis. Note that the valley lengths are typically shorter than the normalized switching cost of our model.

§.§ Model Parameters

We now describe how we parametrized our model in our case studies. In our models, we strive to choose conservative estimates to under-estimate the cost savings from dynamically right-sizing data centers. This approach is similar to the study by  <cit.>. <Ref> gives an overview of the used parameters producing the results of subsequent sections.

Energy We use the linear energy consumption model from <ref> in our experiments. In their analysis,  <cit.> choose energy cost and energy consumption such that the fixed energy cost (i.e., the energy cost of a server when idling) is $1$ and the dynamic energy cost is $0$ as, on most servers, the fixed costs dominate the dynamic costs <cit.>. We investigate this model and an alternative model. In the alternative model, we estimate the power consumption of a server with 1 kW during peak loads and with 500 W when idling to yield a conservative estimate (as cooling costs are included). According to the U.S. Energy Information Administration (EIA), the average cost of energy in the industrial sector in the United States during April 2021 was 6.77 cents per kilowatt-hour <cit.>. We use this as a conservative estimate as data centers typically use a more expensive portfolio of energy sources. If the actual carbon cost of the used energy were to be considered, which is the case in some data centers as discussed in <ref>, energy costs are likely to be substantially higher.

Revenue Loss According to measurements, a 500 ms increase in delay results in a revenue loss of 20% or 0.04%/ms <cit.>. Thus, scaling the delay measured in ms by 0.1 can be used as a slight over-approximation of revenue loss.  <cit.> choose the minimal perceptible delay as 1.5 times the time to run a job, which is a very conservative estimate if valley filling is assumed a viable alternative. In our model, we choose the minimal perceptible delay as 2.5 times the time to run a job which is equivalent as we also added the processing time of a job to the delay. In the case of valley filling, jobs are typically processed with a much more significant delay. Similar to  <cit.>, we also estimate a constant network delay of 10 ms.

Switching Cost We mentioned in <ref> that in practice, the switching cost should be on the order of operating a server between an hour to several hours. To obtain a conservative estimate, we choose $\beta$ such that the normalized switching cost times the length of a time slot equals 4 hours.

Time Slot Length We choose a time slot length of 1 hour. We further assume that the average processing time of jobs is $\delta / 2$ unless noted otherwise.

parameter     model 1     model 2
time slot length     1 hour     1 hour

energy cost     $c=1$     $c=0.0677$

energy consumption     $\Phi_{\text{min}}=1, \Phi_{\text{max}}=1$     $\Phi_{\text{min}}=0.5, \Phi_{\text{max}}=1$

revenue loss     $\gamma = 0.1, \delta_i = 2.5 \eta_i$     $\gamma = 0.1, \delta_i = 2.5 \eta_i$

normalized switching cost     4 hours     4 hours

Models used in our case studies. $\eta_i$ is the processing time of jobs of type $i$.

§ UNI-DIMENSIONAL ALGORITHMS

The results of this section are based on the final day of the LANL Mustang, Facebook, and the second to last day of the Alibaba trace. We begin by discussing the general features of the traces. Then, we compare the uni-dimensional online algorithms with respect to their achieved normalized cost, cost reduction, and runtime.

Fractional vs. Integral Cost For all traces, the ratio of the fractional and the integral costs is 1 for a precision of at least $10^{-3}$. This is not surprising due to the large number of servers used in each model.

Ratio of static and dynamic optima

Correlation of the PMR and the static/dynamic ratio

Ratio of static and dynamic offline optima for each trace. The LANL Mustang and Microsoft Fiddle traces have a significantly higher PMR than the remaining traces. Generally, we observe a strong correlation of PMR and the static/dynamic ratio.

Normalized cost

Cost reduction

Effect of the ratio of static and dynamic optima on the cost reduction and normalized cost achieved by the memoryless algorithm.

Dynamic vs. Static Cost The dynamic and static costs differ significantly depending on the trace. The ratio of dynamic and static optimal costs for each trace is shown in <ref>.

<Ref> shows a strong positive correlation between the average cost reduction achieved by the memoryless algorithm and the ratio of the static and dynamic optima. As $OPT_s / OPT$ is directly linked to the PMR, this also indicates a strong correlation between cost reduction and the PMR. Even under our very conservative estimates of parameters, we achieve a significant cost reduction when the ratio of the static and dynamic offline optimum exceeds 1.5. Similar to  <cit.>, we observe that cost savings increase rapidly as the PMR increases.

We also observe in <ref> that as the static/dynamic ratio increases, the normalized costs achieved by the memoryless algorithm increases too but not as much as the potential energy savings, resulting in the observed significant cost reduction.

Facebook 2009-0

Facebook 2009-1

Facebook 2010

LANL Mustang

Normalized costs of uni-dimensional online algorithms. For the Facebook 2010 trace, the second model results in an optimal schedule constantly using all servers, explaining the disparate performance compared to the first model. Further, Probabilistic and Rand-Probabilistic perform very poorly in this setting and are therefore not shown for model 2. Generally, Probabilistic and Rand-Probabilistic achieve similar results. Interestingly, we observe that LCP outperforms Int-LCP when the potential cost reduction is small. In contrast, Int-LCP and the probabilistic algorithms outperform Memoryless and LCP significantly when the potential cost reduction is large. <Ref> compares the schedules obtained by LCP and Int-LCP in greater detail.

Facebook 2009-0

Facebook 2009-1

Facebook 2010

LANL Mustang

Cost reduction of uni-dimensional online algorithms. For the Facebook 2010 trace, the second model results in an optimal schedule constantly using all servers, explaining the disparate performance compared to the first model. Further, Probabilistic and Rand-Probabilistic perform very poorly in this setting and are therefore not shown for model 2. Results are mainly determined by the normalized cost and the potential cost reduction (or static/dynamic ratio). We observe that the achieved cost reduction is dominated by the potential cost reduction (see <ref>).

Weak positive correlation of achieved normalized cost and cost reduction. Intuitively, one would expect a strong negative correlation, i.e., the achieved cost reduction increases as the normalized cost approaches 1. Here, we find a positive correlation as the achieved cost reduction is dominated by the potential cost reduction (see <ref>).

Normalized Cost In the application of right-sizing data centers, we are interested in the cost associated with integral schedules. <Ref> shows the normalized cost of each algorithm. For fractional algorithms, we consider the cost of the associated integral schedule obtained by ceiling each configuration. Notably, LCP and Int-LCP perform differently depending on the trace and used model. We explore this behavior in <ref>.

Cost Reduction <Ref> shows the achieved cost reduction. In general, we observe in <ref>, <ref>, and <ref> that the achieved cost reduction is dominated by the potential cost reduction (which is mainly influenced by the PMR, see <ref>). When the potential cost reduction is small, algorithms with a smaller normalized cost in a particular setting, achieve a significantly higher cost reduction.

Facebook 2009-0

LANL Mustang

Comparison of the schedules obtained by LCP and Int-LCP under our first model. We observe that LCP is much “stickier” than Int-LCP, which is beneficial when the potential cost reduction is small (i.e., for the Facebook 2009-0 trace) but detrimental when the potential cost reduction is large (i.e., for the LANL Mustang trace). In our experiments, we observe that the memoryless algorithm tends to behave similarly to LCP (i.e., is more “sticky”), whereas the probabilistic algorithms tend to behave similarly to Int-LCP (i.e., are less “sticky”).

Runtimes of uni-dimensional online algorithms.

Runtime <ref> shows the distribution of runtimes (per iteration) of the online algorithms using the Facebook 2009-1 trace. The memoryless algorithm, LCP, and Int-LCP are very fast, even as the number of time slots increases. The runtime of Probabilistic and Rand-Probabilistic is slightly dependent on the used trace and model but generally good. However, when resulting schedules are thight around the upper bound of the decision space, as is the case for the Facebook 2010 trace under our second model, the probabilistic algorithms perform take multiple minutes per iteration. Rand-Probabilistic is significantly slower than Probabilistic as due to the relaxation the integrals need to be computed in a piecewise fashion. The runtime of RBG grows linearly with time and is shown in <ref> for the first four time slots.

Cost profiles of uni-dimensional online algorithms. Note that integral algorithms seem to prefer revenue loss and switching costs over energy costs, whereas fractional algorithms prefer energy costs. However, this is likely because integral algorithms balance energy costs and revenue loss more accurately than the ceiled schedules of fractional algorithms.

Cost Makeup An interesting aspect of the (integral) schedules obtained by the online algorithms is the makeup of their associated costs to understand whether an algorithm systematically prefers some cost over another. We measure this preference of an algorithm as the normalized deviation, i.e., the cost of the algorithm minus the mean cost among all algorithms divided by the standard deviation of costs among all algorithms. We then average the results between all traces. <Ref> shows the cost profiles for each algorithm. We observe that fractional algorithms prefer energy cost over revenue loss and switching cost, while integral algorithms prefer revenue loss and switching cost over energy cost. This is likely because fractional algorithms cannot balance energy cost and revenue loss optimally. When fractional schedules are ceiled, this results in an additional energy cost while reducing revenue loss. In absolute terms, the deviations make up less than 1% of the overall costs of each type.

We have seen that even in our conservative model, significant cost savings with respect to the optimal static provisioning in hindsight can be achieved in practical settings when the PMR is large enough or the normalized switching cost is less than the typical valley length. Due to the conservative estimates of our model, it is likely that in practice, much more drastic cost savings are possible. For example, when energy costs are higher, or the switching costs are on the order of operating a server in an idle state for one hour rather than four hours.

§ MULTI-DIMENSIONAL ALGORITHMS

Now, we turn to the discussed multi-dimensional algorithms. We begin by analyzing the simplified settings, SLO and SBLO, from <ref>. Then, we analyze the gradient-based methods from <ref>.

§.§ Smoothed Load Optimization

Recall that for SLO, during a single time slot a server can process at most one job. Hence, we cannot use dynamic job durations to model the different runtimes of jobs on servers with two GPUs and servers with eight GPUs. Instead, we use a simplified model based on our second model, which is described in <ref>. Note that we disregard revenue loss and that we assume, servers operate at full utilization (if they are active). Overall, we obtain the operating costs $c = (243.720, 219.348)$ and switching costs $\beta = (487.440, 663.672)$.

cost     simplified model
operating cost     servers with eight GPUs have $0.9$ times the energy consumption (per processed job) as servers with two GPUs due to improved cooling efficiency

switching cost     servers with eight GPUs have $1.3$ times the switching cost as servers with two GPUs due to an increased associated risk

Simplified model used in our case studies of SLO and SBLO. The model parameters are based on our second model described in <ref>.

The achieved normalized cost and cost reduction of lazy budgeting are shown in <ref> and <ref>, respectively. The dynamic offline optimal schedule primarily uses 8-GPU-servers and only uses 2-GPU-servers for short periods. The static offline optimal schedule uses 122 8-GPU-servers and no 2-GPU-servers as they have a larger operating cost, which would have to be paid throughout the entire day. The lazy budgeting algorithms primarily use 2-GPU-servers due to their lower switching cost and stick with 8-GPU-servers once they were powered up. <Ref> and <ref> show the schedules obtained by the online algorithms in comparison to the offline optimal schedule. The runtime of the deterministic and randomized variants is shown in <ref>.

Normalized cost

Cost reduction

Performance of lazy budgeting (SLO) for the Microsoft Fiddle trace when compared against the offline optimum. The results of the randomized algorithm are based on five individual runs.

SLO (deterministic)

SLO (randomized)

SBLO ($\epsilon = 1/4$)

Comparison of the schedules obtained by lazy budgeting for the Microsoft Fiddle trace and the offline optimum. For SLO, the deterministic algorithm is shown in blue and one result of the randomized algorithm is shown in red. The lazy budgeting algorithm for SBLO is shown in green. Note that the lazy budgeting algorithms prefer the 2-GPU-servers initially due to their low switching costs. For SLO, the randomized algorithm appears to be less “sticky” than the deterministic algorithm, resulting in a better normalized cost.

Runtime of lazy budgeting algorithms.

§.§ Smoothed Balanced-Load Optimization