Finite-horizon stochastic dynamic programming obtains its boundary condition from the terminal reward. What value equation remains when the process has no terminal date? Discounting makes the infinite reward stream bounded and turns the Bellman equations into fixed-point problems.
The undiscounted expected total reward of policy is
One drawback of this model is that we could easily encounter values that are or , even in a setting as simple as a single-state MDP which loops back into itself and where the accrued reward is nonzero.
Therefore, it is often more convenient to work with an alternative formulation which guarantees the existence of a limit: the expected total discounted reward of policy is defined to be:
for and when , in which case, .
Finally, another possibility for the infinite-horizon setting is the so-called average reward or gain of policy defined as:
We won’t be working with this formulation in this course due to its inherent practical and theoretical complexities.
Extending the previous notion of optimality from finite-horizon models, a policy is said to be discount optimal for a given if:
Furthermore, the value of a discounted MDP, , is defined by:
More often, we refer to by simply calling it the optimal value function.
As for the finite-horizon setting, the infinite horizon discounted model does not require history-dependent policies, since for any there exists a with identical total discounted reward:
Random Horizon Interpretation of Discounting¶
Can geometric termination give the discount factor a probabilistic meaning rather than treating it only as an algebraic device? The use of discounting can be motivated both from a modeling perspective and as a means to ensure that the total reward remains bounded. From the modeling perspective, we can view discounting as a way to put more or less weight on the immediate rewards vs. the long-term consequences. There is also another interpretation which stems from that of a finite horizon model but with an uncertain end time. More precisely:
Let denote the expected total reward obtained by using policy when the horizon length is random. We define it by:
Under this geometric horizon the expected length is , often called the effective horizon of the discounted problem: corresponds to planning about 100 steps ahead on average. The construction can be simulated directly. After collecting each reward, flip a coin that comes up heads with probability : continue on heads and end the episode on tails. The undiscounted return of such an episode is an unbiased estimate of . In practice we rarely simulate this way. Truncating episodes at random throws away the rewards that come after the coin lands on tails, which increases the variance of the estimate, and in learning it discards signal that could have been used. The equivalence is still useful to keep in mind, because it lets an algorithm trade a discount factor for a stopping probability and back.
Drawing an MDP as a graph, and the cemetery state¶
Small MDPs are often drawn as directed graphs, as Puterman does for his toy examples. Each state is a node. Each action available in a state contributes arrows out of that node, one for every possible next state, labeled with the probability . The action sets may differ across states: in Figure 1(a), state has two actions and state only one, as allowed by the state-dependent sets . Without the actions, such a picture is the transition graph of a Markov chain. Fixing a decision rule chooses (or randomizes over) the arrows out of each node and leaves exactly that Markov chain, whose transition matrix is the of the next section. A controlled system is a family of Markov chains, one per decision rule, and folds the choice of action into the transition probabilities.
Figure 1:(a) A two-state MDP drawn as a graph. In , action keeps the system in and action moves it to ; in only is available. (b) The same MDP with a cemetery state . Every original transition probability is multiplied by , and from every state the system moves to with probability , whatever the action. is absorbing and pays no reward.
The graph gives a third reading of discounting. Add a new state , sometimes called a cemetery or sink state. From every state and every action , move to with probability , and otherwise follow the original dynamics, so that the new transition probabilities are for and . Once in the system stays there forever and collects reward zero. In this auxiliary MDP, the undiscounted expected total reward of any policy equals its discounted value in the original MDP, because the probability of still being alive at step is . This is the construction of Theorem 1 written as a state rather than as a random variable, and it is how discounted problems are reduced to total-reward problems with absorbing states, sometimes called stochastic shortest-path problems.
Not every process has an absorbing state. A system can run forever, and for its long-run behavior to be well defined we need some structure. The usual assumption is ergodicity: the chain settles into a stationary distribution that spreads its mass over many states, and long-run averages are averages under that distribution. This is the setting of the average-reward criterion mentioned at the start of the chapter, which we do not pursue.
Vector Representation in Markov Decision Processes¶
How do transition kernels, rewards, and policies become matrices and vectors that expose the Bellman equations as operator equations?
Let V be the set of bounded real-valued functions on a discrete state space S. This means any function satisfies the condition:
where the notation represents the sup-norm (or -norm) of the function .
When working with discrete state spaces, we can interpret elements of V as vectors and linear operators on V as matrices, allowing us to leverage tools from linear algebra. The sup-norm ( norm) of matrix is defined as:
where represents the -th component of the matrix .
For a Markovian decision rule , we define:
For a randomized decision rule , these definitions extend to:
In both cases, denotes a reward vector in , with each component representing the reward associated with state . Similarly, is a transition probability matrix in , capturing the transition probabilities under decision rule .
For a nonstationary Markovian policy , the expected total discounted reward is given by:
Using vector notation, and writing for the -step transition matrix of the nonstationary policy (with ), this can be expressed as:
Each term has a direct reading: is the probability of being in state at time when starting from , so is the expected discounted reward collected steps ahead. The matrix products average over every path the system can take under the policy, which is what a Monte Carlo estimate would approximate by simulation.
Factoring out of every term after the first gives a recursive relationship:
where is the policy shifted by one step. Unrolling this recursion recovers the series above, so the series is the solution of the recursion.
For a stationary policy with constant decision rule , the total expected reward simplifies to:
This last expression is called a Neumann series expansion, and it’s guaranteed to exist under the assumptions of bounded reward and discount factor strictly less than one.
Note that for any induced matrix norm (i.e., a norm satisfying for all vectors ) and any matrix , the spectral radius is bounded by:
This inequality provides a practical way to verify the convergence condition by checking the simpler condition rather than trying to compute the eigenvalues directly.
We can now verify that is invertible and the Neumann series converges.
Norm of the transition matrix: Since is a stochastic matrix (each row sums to 1 and all entries are non-negative), its -norm is:
Norm of the scaled matrix: Using the homogeneity property of norms, we have:
Bounding the spectral radius: Since the spectral radius is bounded by the matrix norm:
Verifying convergence: Since by assumption, we have:
This strict inequality guarantees that is invertible and the Neumann series converges.
Therefore, the Neumann series expansion converges and yields:
The middle sum is indexed from to match the statement of Theorem 2, while the last one keeps the convention used for the value function. They are the same series, shifted by one index.
Consequently, for a stationary policy, can be determined as the solution to the linear equation:
which can be rearranged to:
We can also characterize as the solution to an operator equation. More specifically, define the transformation by
for any . Intuitively, takes a value function as input and returns a new value function that combines immediate rewards () with discounted future values ().
Therefore, we view as an operator mapping elements of to : i.e., . The fact that the value function of a policy is the solution to a fixed-point equation can then be expressed with the statement:
This is a fixed-point equation: the value function is a fixed point of the operator .
Solving Operator Equations¶
Which iterative and Newton-like methods solve a fixed point when direct matrix inversion is unavailable or inappropriate?
The operator equation we encountered in MDPs, , is a specific instance of a more general class of problems known as operator equations. These equations appear in various fields of mathematics and applied sciences, ranging from differential equations to functional analysis.
Operator equations can take several forms, each with its own characteristics and solution methods:
Fixed Point Form: , where . This form is common in fixed-point problems, such as our MDP equation: we seek a fixed point such that .
General Operator Equation: , where . Here, and can be different spaces. We seek an that satisfies the equation for a given .
Nonlinear Equation: , where . A special case of the general operator equation where we seek roots or zeros of the operator.
Variational Inequality: Find such that for all . Here, is a closed convex subset of , and (the dual space of ). These problems often arise in optimization, game theory, and partial differential equations.
Figure 2:The same iteration on four unknowns: a number, a pendulum trajectory, the value of a one-state MDP and the water value function of the reservoir. Each method draws a line through with some slope and jumps to where that line meets . Slope zero gives successive approximation and value iteration. The slope of gives Newton’s method and policy iteration.
Successive Approximation Method¶
For equations in fixed-point form, the basic tool is successive approximation, also called fixed-point iteration: start from any guess and apply the operator repeatedly. Whether this works depends on the operator, and the standard sufficient condition is that it be a contraction.
In other words, a contraction is a Lipschitz map whose Lipschitz constant is strictly less than one. The property is attached to a norm: an operator can be a contraction in one norm and fail to be one in another. For MDPs the norm that works is the sup-norm , and from here on without a subscript means the sup-norm. When we later say that or is a contraction, the precise statement is that it is a sup-norm contraction. Results proved in the sup-norm do not carry over to other norms for free; that a Bellman operator composed with a projection need not be a contraction in a weighted norm is one of the technical roots of the instability that Sutton calls the deadly triad.
The theorem is an existence and uniqueness statement, and its proof is constructive: the sequence used to show that a fixed point exists is the algorithm that computes it. Uniqueness matters, because nothing about an equation alone rules out many solutions. The theorem gives a sufficient condition only. An operator that is not a contraction may still have a fixed point, and other fixed-point theorems cover some of those cases.
Figure 3:What a contraction shrinks. Each dot is a value function of a two-state MDP, and each step applies the Bellman optimality operator (defined later in this chapter) to all of them at once. The dashed box, whose side is the cloud’s diameter in the sup-norm, shrinks by the factor at every step, and the whole cloud collapses onto the unique fixed point .
The stopping test only measures the size of the last step, a quantity computed from the iterates themselves. On its own it says nothing about the distance to the fixed point. Turning it into such a guarantee takes extra work, and this is what the a posteriori bound in Theorem 3 provides for a contraction: stopping when guarantees . When is close to one, a small step can still leave the iterate far from . For MDPs we will choose the tolerance with this in mind, so that value iteration returns an -optimal policy (Proposition 2). Since the error shrinks by at least a factor per step, the convergence rate is linear.
However, the contraction mapping condition is not the only one that can lead to convergence. For instance, if is nonexpansive (i.e., Lipschitz continuous with Lipschitz constant 1) and is a Banach space with certain geometrical properties (e.g., uniformly convex), then under additional conditions (e.g., has at least one fixed point), the successive approximation method can still converge, albeit potentially more slowly than in the contraction case.
In practice, when dealing with specific problems like MDPs or differential equations, the properties of the operator often naturally align with one of these convergence conditions. For example, in discounted MDPs, the Bellman operator is a contraction in the supremum norm, which guarantees the convergence of value iteration.
For policy evaluation, successive approximation replaces the linear solve with the loop . A direct solve of costs on the order of operations, while each iteration costs one matrix-vector product, and the loop never needs to form or factor a matrix. Here is a contraction with factor , by the same norm argument we used for the Neumann series.
The second general tool for operator equations is the Newton-Kantorovich method, which extends Newton’s method to function spaces. We postpone it to Newton-Kantorovich and Policy Iteration: applied to the Bellman optimality equation it yields policy iteration, just as successive approximation applied to the same equation yields value iteration.
Optimality Equations for Infinite-Horizon MDPs¶
How does actionwise maximization turn fixed-policy evaluation into the Bellman optimality equation?
Recall that in the finite-horizon setting, the optimality equations are:
where is the value function at time step for state , is the set of actions available in state , is the reward function, is the discount factor, and is the transition probability from state to state given action .
Intuitively, we would expect that by taking the limit of to infinity, we might get the nonlinear equations:
which are called the optimality equations or Bellman equations for infinite-horizon MDPs.
We can adopt an operator-theoretic perspective by defining operators on the space of bounded real-valued functions on the state space . For a deterministic Markov rule , define the policy-evaluation operator:
Terminology varies here. Bertsekas and Sutton and Barto call the Bellman operator for , and part of the RL theory literature calls it the Bellman consistency operator. In this book we keep the name Bellman operator for the operator with the maximum, defined next, which is closer to Bellman’s own usage, and call the policy-evaluation operator because evaluating a policy is what it is used for.
The Bellman optimality operator is then:
where is the set of Markov deterministic decision rules, is the reward vector under decision rule , and is the transition probability matrix under decision rule .
Note that while we write , we do not implement the above operator by enumerating all decision rules. Rather, the fact that we compare policies based on their value functions in a componentwise fashion means that maximizing over the space of Markovian deterministic rules reduces to the following update in component form:
The two forms of , a maximum over decision rules and a maximum over actions in each state, are equal. The following proposition makes this precise.
For convenience, we define the greedy selector that extracts a decision rule from a value function:
Equivalently, is any decision rule with . Puterman calls such a rule -improving. A -improving rule, one that is greedy with respect to the optimal value function, is called conserving. In the RL literature both are simply called greedy policies: given , the policy can act myopically because already summarizes the future. This selector will be useful for expressing algorithms succinctly:
Value iteration: iterate , then extract at termination
Policy iteration: with solving
Optimality of Deterministic Stationary Policies¶
The definition of optimality at the start of the chapter compares with every policy in : history-dependent, nonstationary and randomized. The comparison is componentwise, so must be at least as good as every other policy in every starting state. Searching that class directly is hopeless. A history-dependent rule with depends on a history that grows with , and a nonstationary policy needs one rule per stage, which for an infinite horizon means infinitely many. Even restricted to deterministic Markov policies, enumerating every policy and evaluating each one is out of reach. The next result shows that nothing is lost by searching a much smaller set.
See Puterman (1994), Section 6.2, for the proof. The contraction property is shown in the proof of Proposition 2 below. Two consequences are worth stating. First, the supremum over is attained, by a deterministic stationary policy, so we may restrict the search to the set of deterministic Markov decision rules applied at every stage. Second, the optimal value function is unique but the optimal policy need not be: when several actions attain the maximum in some state, each choice gives an optimal policy. Note also that the theorem is specific to the discounted infinite-horizon case. In a finite-horizon problem the optimal policy depends on the stage in general.
Q-Functions¶
The optimality equations can also be written in terms of a function of a state and an action. The optimal Q-function (called Q-factors in the operations research literature) is the unique solution of
As with , the right-hand side defines an operator on functions of , and later chapters write it with the same symbol, . It is a sup-norm contraction with factor when , so exists and is unique. The two optimal functions determine each other:
and a conserving rule can be read off directly as . Probabilistically, conditions on both the state and the first action, while averages the action out. Optimality is still defined through , because that is how policies are ordered.
When and are known, as in the examples of this chapter, there is no reason to prefer : computing is cheaper and the greedy rule follows from it. The Q-function matters in model-free RL, where we observe sampled transitions and rewards but not the transition probabilities. Extracting a greedy action from requires the expectation , which needs the model. Extracting it from only requires a maximum over the action argument, which works for a learned , a neural network for instance, without knowing . This is why the Q-function form is the starting point of the simulation-based chapters. The convenience has a price. A Q-function has one value per state-action pair instead of one per state, and it keeps the detail that averages away, so learning typically takes more capacity, more samples and better coverage of the actions than learning would.
Algorithms for Solving the Optimality Equations¶
What computational trade-off separates value iteration, policy evaluation, and policy improvement?
The optimality equations are operator equations. Therefore, we can apply general numerical methods to solve them. Applying the successive approximation method to the Bellman optimality equation yields a method known as “value iteration” in dynamic programming. Following Puterman and the older operations research literature, we reserve this name for the iteration with the max. The same iteration applied to the affine operator is iterative policy evaluation, although RL texts often call any iteration of this form value iteration. A direct application of the blueprint for successive approximation yields the following algorithm:
The termination criterion in this algorithm is based on a specific bound that provides guarantees on the quality of the solution. This is in contrast to supervised learning, where we often use arbitrary termination criteria based on computational budget or early stopping when the learning curve flattens. This is because establishing implementable generalization bounds in supervised learning is challenging.
However, in the dynamic programming context, we can derive various bounds that can be implemented in practice. These bounds help us terminate our procedure with a guarantee on the precision of our value function and, correspondingly, on the optimality of the resulting policy.
This detail is easy to miss. The value iteration in Sutton & Barto (2018) stops when the change falls below a small threshold , without relating to the quality of the returned policy. The threshold above is chosen so that the answer comes with a certificate: you fix the accuracy in advance, the loop runs until the test passes, and the greedy policy it returns is then guaranteed to be -optimal (Proposition 2). Two caveats apply. The guarantee concerns the exact arithmetic; floating-point rounding is not part of the analysis, and with close to one a very small may be below what float32 or even float64 can resolve. And whether a certificate matters depends on the application. Nobody needs one to solve a grid world, but an economist who reports a policy’s monetary value, or an engineer who must document a safety margin, does.
The value function is only an instrument. What we want in the end is a policy, and given (or any approximation of it) we obtain one by the greedy step : for each state, compute the one-step lookahead for every action and keep the best. With arrays of shape this is a single argmax along the action axis. Nothing is simulated forward, because already summarizes the future.
A Sharper Stopping Rule: the Span Seminorm¶
A good stopping rule stops as early as possible while still delivering the requested accuracy. The sup-norm test is not the best one available. A sharper test measures the change between iterates with the span seminorm
It is a seminorm rather than a norm because for every constant : the span ignores a uniform shift of all values. Let be the last change made by value iteration, with smallest and largest components and . Two properties of give two-sided bounds. It is monotone ( componentwise implies ), and it shifts constants by : . From we get , and repeating the argument,
The same argument with in place of , for , shows that . Subtracting,
So stopping when also returns an -optimal greedy policy (Puterman (1994), Section 6.6). Since , this test is never later than the sup-norm test, and it is often much earlier: value iteration tends to spend many iterations shifting all values by nearly the same amount, which the span ignores. On a random MDP with 30 states, 4 actions, and rewards around 5, the sup-norm test with needs 242 iterations and the span test needs 11, and both return an optimal policy. The bounds also give a better value estimate than itself: the midpoint is within of .
Bellman contraction laboratory¶
Use the controls below to change the discount factor, transition persistence, reward asymmetry, and starting value. Before moving , predict how it will change the slope of the error envelope. The middle panel compares the observed error and Bellman residual with the contraction bound; the text below reports the final greedy policy.
Figure 4:Static preview of value-iteration traces and a geometric contraction bound. The online book provides controls for , transitions, rewards, and the initial value.
Policy Iteration¶
Value iteration searches the space of value functions and reads a policy off at the end, so it reaches the policy indirectly. Policy iteration searches the space of policies directly. Bellman called this idea approximation in policy space, and Howard (1960) turned it into the algorithm below. Each iteration has two steps. The first evaluates the current policy, which is the linear problem we already know how to solve, so any policy-evaluation routine can be reused as a subroutine. The second improves the policy by acting greedily with respect to the value just computed.
Recall that and are built by averaging over the actions the rule selects: and . The tie-breaking instruction in step 4 (“keep if it is among the maximizers”) matters for the stopping test. Without it, the algorithm could cycle among equally good greedy rules and never see the same rule twice.
Two remarks follow from the proof. First, the stopping rule is exact. When the same rule appears twice, solves the optimality equation and the policy is optimal, not merely -optimal. This argument counts policies, so it is specific to finite state and action spaces; for continuous spaces, convergence needs the Newton-Kantorovich machinery of Newton-Kantorovich and Policy Iteration. Second, the proof never used that is greedy in every state. Any with , strictly in at least one state, improves the policy. Improving a single state per iteration is enough for convergence. This matters later, because many RL algorithms contain a policy-iteration loop whose improvement step is only partial: a small change at each step, perhaps only in the states that were visited.
In practice policy iteration converges in very few iterations, often a handful, even for large problems. The reason is that it is Newton’s method in disguise, a connection made by Puterman & Brumelle (1979) and developed in Newton-Kantorovich and Policy Iteration. Each iteration costs more than a value-iteration step, though. Both algorithms pay about operations for a greedy step, but policy iteration adds a direct policy evaluation, a linear solve of order .
Modified Policy Iteration¶
Policy evaluation can also be done by successive approximation, and nothing forces us to run that inner loop to convergence. Stopping it after a few sweeps gives modified policy iteration (MPI) Puterman & Shin (1978). One way to see what this does is through the Neumann series. Write the exact evaluation of the greedy rule as a correction to the current guess:
The step is the Bellman residual , propagated through the inverse. MPI truncates the series after terms:
where the second form, sweeps of policy evaluation started from , is how it is implemented. The integer is a hyperparameter that interpolates between the two algorithms we have seen. With the update is , which is value iteration. With it is exact evaluation, which is policy iteration. In between, each iteration costs matrix-vector products and no linear solve.
The stopping test is the one from value iteration, and it gives the same guarantee, provided the returned policy is greedy with respect to (one more greedy step), exactly as value iteration returns the policy greedy with respect to its last iterate. Started from a with , the iterates increase monotonically to and the returned policy is -optimal (Puterman (1994), Section 6.5). The number of sweeps can also change from one iteration to the next.
Why should a truncated evaluation work at all? In Newton’s method, the inverse Hessian (or inverse Jacobian) rescales the gradient step, a step-size schedule adapted to the local curvature. Here plays that role for the Bellman residual, and value iteration replaces it by the identity. Truncating the Neumann series replaces it by an approximation, in the same way that quasi-Newton and inexact Newton methods use approximate Hessians or approximate linear solves. Newton-Kantorovich and Policy Iteration makes this precise.
Exact Scheduling MDP for Inference Serving¶
How can a large request-level system be reduced to a finite MDP, and which predictive distinctions disappear in that reduction?
The inference examples have so far treated the scheduling rule as fixed and the GPU clock as the action. In a different decision interface, the clock is fixed and the action is the phase that receives the next unit of service. Prefill processes a waiting prompt and moves its request into the decode batch. Decode advances every request in that batch by one step of token generation. Serving either phase delays the other.
The trade-off comes from batching. A decode step costs about the same time and power whether it advances one request or six, because generating a token is limited by reading the model weights from memory, not by arithmetic. On the measured NVIDIA L4 run used below, one decode stream produces a token every 57.0 ms; with four concurrent streams each stream still produces a token every 59.5 ms, and with eight every 60.8 ms. Six streams therefore finish requests almost six times as fast as one. A scheduler that decodes as soon as one request is ready keeps the batch small and wastes that capacity. A scheduler that keeps prefilling to grow the batch leaves the requests already in it waiting. The MDP decides where to switch.
An exact request-level Markov state would contain every prompt length, generated token count, cache allocation, and waiting time. For computation, these quantities are aggregated into the finite state
where counts waiting prefill jobs and counts active decode jobs. The actions are
An action is masked when its phase is empty. Prefill is also masked at , which represents the cache limit in this abstraction. The whole policy is therefore one table, which the figure below draws in full.
Time advances in decision periods of 0.1 seconds. A Bernoulli arrival with probability per period is calibrated from the load-normalized version of the same five-minute Azure trace used in the modeling chapter. The phase rates and powers come from one measured NVIDIA L4 run of Qwen/Qwen2.5-7B-Instruct served by vLLM 0.28.0, at the 1,125 MHz requested clock level; its batch-balanced median realized graphics clock was 939.375 MHz.
Prefill. The measured prefill rate and the trace’s mean prompt length give the probability that one prefill action completes a prompt. A completed prompt moves from to .
Decode. The single-stream time per output token and the trace’s mean output length give the probability that a lone decode stream finishes its response in one period. With active streams, each stream finishes independently with probability , where is the measured slowdown of one stream in a batch of , interpolated between the measured batch sizes 1, 4, and 8. The number of completions is binomial, so its mean grows almost linearly in .
Arrivals and drops. After service, an arrival joins the prefill queue. An arrival that finds is dropped.
These choices define a complete transition matrix on 49 states.
The one-step cost counts the requests in the system, dropped arrivals, and energy:
By Little’s law, the long-run average of divided by the arrival rate is the mean time a request spends in the system, so the first term is a latency cost. A dropped request costs as much as holding one request for 100 periods, or 10 seconds. That is about twice the mean time a request spends in the system under the optimal policy, so dropping a request is never a cheap way to shorten the queue. For prefill, decode, and idle, respectively, the measured profile gives joules per decision period. Each value is the measured phase-power summary at the requested 1,125 MHz level multiplied by 0.1 seconds, rather than a direct request-level energy measurement. Only the ratio enters the stage cost.
With , cost-minimizing value iteration applies
Iteration stops when . The final Bellman residual is checked independently and must be below 10-8.
The figure compares the value-iteration policy with two transparent rules. Prefill priority prefills whenever a prompt waits and the batch has room. Decode priority decodes whenever a decode job is active. Each replay applies one policy for 60 seconds from an empty system. All three replays use the same random numbers for arrivals and service, so they see the same requests, and every difference between them comes from the decisions.
Figure 5:Exact value iteration on the measured-L4-calibrated inference-scheduling MDP, compared with prefill priority and decode priority. Each policy is drawn over all 49 states. The panel beside it reports the policy’s exact discounted cost from the empty system and its smallest gap over all states. The replay samples queue transitions from the reduced model at its fixed measured-profile clock; these trajectories are simulated, not direct vLLM observations.
Figure 6:The exact policy, prefill priority, and decode priority over all 49 states, with each policy’s exact discounted cost from the empty system. The interactive version adds controls for playing, stepping through, and scrubbing an episode sampled from the reduced transition model under each policy. The queue trajectory is simulated rather than observed directly from vLLM.
The optimal policy is a switching curve. When one or two prompts wait, it prefills until four requests are decoding and then decodes. When three or more wait, it fills the batch to five before decoding. It decodes whenever no prompt waits, and idles only in the empty state. Decode priority is the opposite extreme. It prefills only when the batch is empty, so it never runs more than one stream, and a lone stream finishes fewer requests per period (0.059) than arrive (0.071). Its prompt queue therefore fills, and arrivals are dropped.
How can a reader check that value iteration returned the optimal policy rather than a plausible one? For a fixed policy , the policy evaluation equation is a linear system with 49 unknowns, and solving it gives the exact discounted cost of in every state. An optimal policy satisfies in every state for every stationary policy , so each heuristic must lose everywhere, not only on average.
The exact rows come from the linear solves. The simulated rows average 2,000 seeded episodes of 300 periods from the empty system. Value iteration required 2,435 sweeps on this high-discount problem. Its independently recomputed Bellman residual is , and evaluating its policy by a linear solve reproduces its value function to within 10-8 in every state. Prefill priority is worse than the optimal policy by at least 3.02 in every state and disagrees with it in 8 states, all on the switching curve, where prefill priority keeps filling the batch to six. Decode priority is worse by at least 131.82 in every state. In the simulated episodes it drops 4.35 arrivals per episode and roughly doubles the time a request spends in the system. Prefill priority drops slightly fewer arrivals than the optimal policy, 0.14 against 0.20 per episode. The optimal policy minimizes the total discounted cost, not each term of it, and accepts a few more drops for less time spent waiting in the batch.
The reduced scheduling MDP has the following transition parameters and value-iteration certificate.
Download the dynamic-programming certificate (CSV)
The certificate concerns the supplied 49-state transition matrix. The result is exact for the stated finite MDP, not for vLLM. Aggregating request ages and lengths removes distinctions that can affect head-of-line waiting and cache release; all prompts and responses are treated as having the trace’s mean length. The calibrated transition kernel is stationary, the clock is fixed, temperature is absent from the state, the batch is capped at six, and the action set excludes mixed prefill-decode batches. Measurements from one L4 deployment calibrate the phase rates, batch slowdown, and powers; they do not establish scheduling performance across L4 systems. The request-level replay from the earlier chapters remains a model-audit tool, not part of the optimality proof.
Inspect the scheduling MDP, value iteration, and the policy check
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132def make_scheduling_mdp( *, arrival_probability: float = 0.07, prefill_completion_probability: float = 0.14, decode_completion_probability: float = 0.06, decode_batch_efficiency: Sequence[float] | None = None, gamma: float = 0.99, decision_period_s: float = 0.1, action_energy: Sequence[float] = (0.85, 1.0, 0.25), drop_penalty: float = 100.0, energy_weight: float = 0.1, ) -> SchedulingMDP: """Construct the 49-state prefill/decode scheduling abstraction. ``decode_completion_probability`` is the probability that a decode stream running alone finishes its response during one decision period. With ``d`` active streams each one finishes independently with probability ``1 - (1 - q) ** efficiency[d]``, so a decode step completes about ``d * efficiency[d]`` times as many jobs as a lone stream. """ cap = SCHEDULING_QUEUE_CAP side = cap + 1 if not 0.0 <= arrival_probability <= 1.0: raise ValueError("arrival_probability must lie in [0, 1]") if not 0.0 <= prefill_completion_probability <= 1.0: raise ValueError("prefill_completion_probability must lie in [0, 1]") if not 0.0 <= decode_completion_probability <= 1.0: raise ValueError("decode_completion_probability must lie in [0, 1]") if not 0.0 < gamma < 1.0: raise ValueError("gamma must lie strictly between zero and one") if decision_period_s <= 0.0: raise ValueError("decision_period_s must be positive") if drop_penalty < 0.0 or energy_weight < 0.0: raise ValueError("drop_penalty and energy_weight must be nonnegative") energy = np.asarray(action_energy, dtype=float) if energy.shape != (3,) or np.any(energy < 0.0): raise ValueError("action_energy must contain three nonnegative values") efficiency = ( np.ones(side, dtype=float) if decode_batch_efficiency is None else np.asarray(decode_batch_efficiency, dtype=float) ) if efficiency.shape != (side,) or np.any(efficiency <= 0.0): raise ValueError("decode_batch_efficiency must hold one positive value per batch size") job_probability = 1.0 - (1.0 - decode_completion_probability) ** efficiency job_probability[0] = 0.0 states = np.array([(p, d) for p in range(side) for d in range(side)], dtype=int) number_states = states.shape[0] transitions = np.zeros((number_states, 3, number_states), dtype=float) drop_transitions = np.zeros_like(transitions) costs = np.full((number_states, 3), np.inf, dtype=float) valid = np.zeros((number_states, 3), dtype=bool) drops = np.zeros((number_states, 3), dtype=float) maximum_energy = max(float(np.max(energy)), 1e-12) for state_index, (p_value, d_value) in enumerate(states): p_initial = int(p_value) d_initial = int(d_value) for action in range(3): is_valid = ( action == 2 or (action == 0 and p_initial > 0 and d_initial < cap) or (action == 1 and d_initial > 0) ) if not is_valid: continue valid[state_index, action] = True if action == 0: service_outcomes = ( (p_initial - 1, d_initial + 1, prefill_completion_probability), (p_initial, d_initial, 1.0 - prefill_completion_probability), ) elif action == 1: service_outcomes = tuple( (p_initial, d_initial - completed, probability) for completed, probability in enumerate( _binomial_probabilities( d_initial, float(job_probability[d_initial]) ) ) ) else: service_outcomes = ((p_initial, d_initial, 1.0),) for p_base, d_after, service_probability in service_outcomes: for arrival, arrival_prob in ( (0, 1.0 - arrival_probability), (1, arrival_probability), ): dropped = arrival == 1 and p_base >= cap p_after = min(cap, p_base + arrival) next_index = _state_index(p_after, d_after) probability = service_probability * arrival_prob transitions[state_index, action, next_index] += probability if dropped: drop_transitions[state_index, action, next_index] += probability drops[state_index, action] += probability costs[state_index, action] = ( p_initial + d_initial + drop_penalty * drops[state_index, action] + energy_weight * energy[action] / maximum_energy ) mdp = SchedulingMDP( states=states, transitions=transitions, stage_cost=costs, valid_actions=valid, drop_probability=drops, arrival_probability=float(arrival_probability), prefill_completion_probability=float(prefill_completion_probability), decode_completion_probability=float(decode_completion_probability), action_energy=energy, gamma=float(gamma), decision_period_s=float(decision_period_s), description=( "Exactly solved reduced abstraction with capped prefill and decode counts; " "measured single-stream decode speed and measured batch scaling, so a " "larger decode batch finishes more jobs per step; it is not an exact " "reproduction of vLLM." ), decode_batch_efficiency=efficiency, decode_job_completion_probability=job_probability, drop_transitions=drop_transitions, drop_penalty=float(drop_penalty), energy_weight=float(energy_weight), ) mdp.validate() return mdp
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75def solve_scheduling_mdp( mdp: SchedulingMDP | None = None, *, tolerance: float = 1e-10, maximum_iterations: int = 100_000, ) -> DPSolution: """Solve the reduced discounted scheduling model by value iteration.""" model = make_scheduling_mdp() if mdp is None else mdp model.validate() if tolerance <= 0.0 or maximum_iterations <= 0: raise ValueError("tolerance and maximum_iterations must be positive") value = np.zeros(model.states.shape[0], dtype=float) iterations = 0 for iterations in range(1, maximum_iterations + 1): continuation = np.einsum("sak,k->sa", model.transitions, value) q_values = model.stage_cost + model.gamma * continuation q_values = np.where(model.valid_actions, q_values, np.inf) updated = np.min(q_values, axis=1) if np.max(np.abs(updated - value)) <= tolerance: value = updated break value = updated continuation = np.einsum("sak,k->sa", model.transitions, value) q_values = np.where( model.valid_actions, model.stage_cost + model.gamma * continuation, np.inf, ) bellman_update = np.min(q_values, axis=1) residual = float(np.max(np.abs(bellman_update - value))) policy = np.argmin(q_values, axis=1).astype(int) labels = tuple(SCHEDULING_ACTIONS[int(action)] for action in policy) return DPSolution( value=value, policy=policy, policy_labels=labels, iterations=iterations, bellman_residual=residual, mdp=model, ) def decode_priority_policy(mdp: SchedulingMDP) -> np.ndarray: """Decode whenever a decode job is active, else prefill, else idle.""" p_value, d_value = mdp.states[:, 0], mdp.states[:, 1] return np.where(d_value > 0, 1, np.where(p_value > 0, 0, 2)).astype(int) def prefill_priority_policy(mdp: SchedulingMDP) -> np.ndarray: """Prefill whenever the batch has room and a prompt waits, else decode.""" p_value, d_value = mdp.states[:, 0], mdp.states[:, 1] return np.where( (p_value > 0) & (d_value < SCHEDULING_QUEUE_CAP), 0, np.where(d_value > 0, 1, 2), ).astype(int) def policy_value(mdp: SchedulingMDP, policy: Sequence[int]) -> np.ndarray: """Exact discounted cost-to-go of a stationary policy by a linear solve.""" policy_array = np.asarray(policy, dtype=int) number_states = mdp.states.shape[0] if policy_array.shape != (number_states,): raise ValueError("policy has an incompatible shape") rows = np.arange(number_states) if np.any(~mdp.valid_actions[rows, policy_array]): raise ValueError("policy selects an invalid action") transition = mdp.transitions[rows, policy_array] cost = mdp.stage_cost[rows, policy_array] return np.linalg.solve(np.eye(number_states) - mdp.gamma * transition, cost)
Newton-Kantorovich and Policy Iteration¶
Value iteration is successive approximation applied to the Bellman optimality equation. This section develops the other general tool for operator equations, the Newton-Kantorovich method, and shows that applying it to the same equation gives the policy iteration algorithm of Policy Iteration. This explains why policy iteration needs so few iterations, and why modified policy iteration, which solves the Newton system only approximately, still works.
Newton-Kantorovich Method¶
The Newton-Kantorovich method is a generalization of Newton’s method from finite dimensional vector spaces to infinite dimensional function spaces: rather than iterating in the space of vectors, we are iterating in the space of functions.
Newton’s method is often written as the familiar update:
which makes it look as though the essence of the method is “take a derivative and invert it.” But the real workhorse behind Newton’s method (both in finite and infinite dimensions) is linearization.
At each step, the idea is to replace the nonlinear operator by a local surrogate model of the form
where is a linear map capturing how small perturbations in the input propagate to changes in the output. This is a Taylor-like expansion in Banach spaces: the role of the derivative is precisely to provide the correct notion of such a linear operator.
To find a root of , we impose the condition that the surrogate vanishes at the next iterate:
Solving this linear equation gives the increment . In finite dimensions, is the Jacobian matrix; in Banach spaces, it must be the Fréchet derivative.
But what exactly is a Fréchet derivative in infinite dimensions? To understand this, we need to generalize the concept of derivative from finite-dimensional calculus. In infinite-dimensional spaces, there are several notions of differentiability, each with different strengths and requirements:
1. Gâteaux (Directional) Derivative
We say that the Gâteaux derivative of at in a specific direction is:
This quantity measures how the function changes along the ray . While this limit may exist for each direction separately, it doesn’t guarantee that the derivative is linear in . This is a key limitation: the Gâteaux derivative can exist in all directions but still fail to provide a good linear approximation.
2. Hadamard Directional Derivative
Rather than considering a single direction of perturbation, we now consider a bundle of perturbations around . We ask how the function changes as we approach the target direction from nearby directions. We say that has a Hadamard directional derivative if:
This is a stronger condition than Gâteaux differentiability because it requires the limit to be uniform over nearby directions. However, it still doesn’t guarantee linearity in .
3. Fréchet Derivative
The strongest and most natural notion: is Fréchet differentiable at if there exists a bounded linear operator such that:
This definition directly addresses the inadequacy of the previous notions. Unlike Gâteaux and Hadamard derivatives, the Fréchet derivative explicitly requires the existence of a linear operator that provides a good approximation. Key properties:
must be linear in (unlike the directional derivatives above)
The approximation error is , uniform in all directions
This is the “true” derivative: it generalizes the Jacobian matrix to infinite dimensions
Notation: or
Relationship:
In the context of the Newton-Kantorovich method, we work with an operator where both and are Banach spaces. The Fréchet derivative is the best linear approximation of near , and it’s exactly this linear operator that we use in our linearization .
Now apart from those mathematical technicalities, Newton-Kantorovich has in essence the same structure as that of the original Newton’s method. That is, it applies the following sequence of steps:
Linearize the Operator: Given an approximation , we consider the Fréchet derivative of , denoted by . This derivative is a linear operator that provides a local approximation of near .
Set Up the Newton Step: The method then solves the linearized equation for a correction :
This equation represents a linear system where is chosen so that the linearized operator equals zero.
Update the Solution: The new approximation is then given by:
This correction step refines , bringing it closer to the true solution.
Repeat Until Convergence: We repeat the linearization and update steps until the solution converges to the desired tolerance, which can be verified by checking that is sufficiently small, or by monitoring the norm .
The convergence of Newton-Kantorovich does not hinge on being a contraction over the entire domain (as it could be the case for successive approximation). The convergence properties of the Newton-Kantorovich method are as follows:
Local Convergence: Under mild conditions (e.g., is Fréchet differentiable and is invertible near the solution), the method converges locally. This means that if the initial guess is sufficiently close to the true solution, the method will converge.
Global Convergence: Global convergence is not guaranteed in general. However, under stronger conditions (e.g., is analytic and satisfies certain bounds), the method can converge globally.
Rate of Convergence: When the method converges, it typically exhibits quadratic convergence. This means that the error at each step is proportional to the square of the error at the previous step:
where is the true solution and is some constant. This quadratic convergence is significantly faster than the linear convergence typically seen in methods like successive approximation.
Newton-Kantorovich Applied to Bellman Optimality¶
We now apply the Newton-Kantorovich framework to the Bellman optimality equation. Let
The problem is to find such that , or equivalently . The operator is piecewise affine, hence not globally differentiable, but it is directionally differentiable everywhere in the Hadamard sense and Fréchet differentiable at points where the maximizer is unique.
We consider three complementary perspectives for understanding and computing its derivative.
Perspective 1: Max of Affine Maps¶
In tabular form, for finite state and action spaces, the Bellman operator can be written as a pointwise maximum of affine maps:
where is the transition matrix associated with action . Each is affine in . The operator therefore computes the upper envelope of a finite set of affine functions at each state.
At any , let the active set at state be
Then the Hadamard directional derivative exists and is given by
If the active set is a singleton, this expression becomes linear in , and is Fréchet differentiable at , with
where is the greedy policy at .
Perspective 2: Envelope Theorem¶
Consider now a value function approximated as a linear combination of basis functions:
At a node , define the parametric maximization
Define
so that . Since is linear in , we can apply the envelope theorem (Danskin’s theorem): if the optimizer is unique or selected measurably, then
We do not need to differentiate the optimizer itself. The result extends to the subdifferential case when ties occur, where the Jacobian becomes set-valued.
This result is useful when solving the collocation equation . Newton’s method requires the Jacobian , and this expression allows us to compute it without involving any derivatives of the optimal action.
Perspective 3: The Implicit Function Theorem¶
The third perspective applies the implicit function theorem to understand when the Bellman operator is differentiable despite containing a max operator. The maximization problem defines an implicit relationship between the value function and the optimal action, and the implicit function theorem tells us when this relationship is smooth enough to differentiate through.
The Bellman operator is defined as
The difficulty is that the max operator encodes a discrete selection: which action achieves the maximum. To apply the implicit function theorem, we reformulate this as follows. For each action , define the action-value function:
The optimal action at satisfies the optimality condition:
Now suppose that at a particular , action is a strict local maximizer in the sense that there exists such that
This strict inequality is the regularity condition needed for the implicit function theorem. It ensures that the optimal action is unique at and remains so in a neighborhood of .
To see why, consider any perturbation with small. Since is linear in , we have:
The perturbation term is bounded: . Therefore, for , the strict gap ensures that
Thus remains the unique maximizer throughout the neighborhood .
The implicit function theorem now applies: in this neighborhood, the mapping is constant (and hence smooth), taking the value . This allows us to write
as an explicit formula that holds throughout the neighborhood. Since is an affine (hence smooth) function of , we can differentiate it:
More precisely, for any perturbation :
This is the Fréchet derivative:
where is the greedy policy.
The role of the implicit function theorem: It guarantees that when the maximizer is unique with a strict gap (the regularity condition), the argmax function is locally constant, which removes the non-differentiability of the max operator. Without this regularity condition (specifically, at points where multiple actions tie for optimality), the implicit function theorem does not apply, and the operator is not Fréchet differentiable. The active set perspective (Perspective 1) and the envelope theorem (Perspective 2) provide the tools to handle these non-smooth points.
Connection to Policy Iteration¶
We return to the Newton-Kantorovich step:
Suppose for the greedy policy . Then
which is exactly policy evaluation for . Recomputing the greedy policy from yields the next iterate.
Thus, policy iteration is Newton-Kantorovich applied to the Bellman optimality equation. At points of nondifferentiability (when ties occur), the operator is still semismooth, and policy iteration corresponds to a semismooth Newton method. The envelope theorem is what justifies the simplification of the Jacobian to , bypassing the need to differentiate through the optimizer. This completes the equivalence.
The Semismooth Newton Perspective¶
The three perspectives we developed above (the active set view, the envelope theorem, and the implicit function theorem) all point toward a deeper framework for understanding Newton-type methods on non-smooth operators. This framework, known as semismooth Newton methods, was developed precisely to handle operators like the Bellman operator that are piecewise smooth but not globally differentiable. The connection between policy iteration and semismooth Newton methods has been rigorously developed in recent work Gargiani et al. (2022).
The classical Newton-Kantorovich method assumes the operator is Fréchet differentiable everywhere. The derivative exists, is unique, and varies continuously with the base point. But the Bellman operator violates this assumption at any value function where multiple actions tie for optimality at some state. At such points, the implicit function theorem fails, and there is no unique Fréchet derivative.
Semismooth Newton methods address this by replacing the notion of a single Jacobian with a generalized derivative that captures the behavior of the operator near non-smooth points. The most commonly used generalized derivative is the Clarke subdifferential, which we can think of as the convex hull of all possible “candidate Jacobians” that arise from limits approaching the non-smooth point from different directions.
For the Bellman residual , the Clarke subdifferential at a point can be characterized explicitly using our first perspective. Recall that at each state , we defined the active set . When this set contains multiple actions, the operator is not Fréchet differentiable. However, it remains directionally differentiable in all directions, and the Clarke subdifferential consists of all matrices of the form
In words, the generalized Jacobian is the set of all matrices where is any policy that selects an action from the active set at each state. When the maximizer is unique everywhere, this set reduces to a singleton, and we recover the classical Fréchet derivative. When ties occur, the set has multiple elements: precisely the convex combinations mentioned in Perspective 1.
The semismooth Newton method for solving proceeds by selecting an element at each iteration and solving
What this tells us is that any choice from the Clarke subdifferential yields a valid Newton-like update. In the context of the Bellman equation, choosing where is any greedy policy corresponds exactly to the policy evaluation step in policy iteration. The freedom in selecting which action to choose when ties occur translates to the freedom in selecting which element of the subdifferential to use.
Under appropriate regularity conditions (specifically, when the residual function is BD-regular or CD-regular), the semismooth Newton method converges locally at a quadratic rate Gargiani et al. (2022). This means that near the solution, the error decreases quadratically:
This theoretical result explains an empirical observation that has long been noted in practice: policy iteration typically converges in very few iterations, often just a handful, even when the state and action spaces are enormous and the space of possible policies is exponentially large.
The semismooth Newton framework also suggests a spectrum of methods interpolating between value iteration and policy iteration. Value iteration can be interpreted as a Newton-like method where we choose at every iteration, ignoring the dependence of on entirely. This choice guarantees global convergence through the contraction property but sacrifices the quadratic local convergence rate. Policy iteration, at the other extreme, uses the full generalized Jacobian , achieving quadratic convergence but at the cost of solving a linear system at each iteration.
Between these extremes lie methods that use approximate Jacobians or solve the Newton system only approximately. Modified policy iteration is one of them: it keeps the Jacobian but replaces its inverse by the first terms of the Neumann series. Another natural variant is to choose for some scalar . This leads to the update
This is relaxed value iteration with step size . Taking damps the step: the update remains a global contraction, but with the larger modulus , so it is slower than standard value iteration. Taking lengthens the step, which is over-relaxation. The sup-norm argument then gives a contraction only for , with modulus , which exceeds . Near the solution, however, the rate is set by the eigenvalues of . The constant mode now contracts at , so over-relaxation is faster whenever has no eigenvalues close to -1. For every , the update requires only pointwise operations rather than solving a linear system. The Newton perspective thus unifies existing algorithms and generates new ones by systematically exploring different approximations to the generalized Jacobian.
The connection to semismooth Newton methods places policy iteration within a broader mathematical framework that extends far beyond dynamic programming. Semismooth Newton methods are used in optimization (for complementarity problems and variational inequalities), in PDE-constrained optimization (for problems with control constraints), and in economics (for equilibrium problems). The Bellman equation, viewed through this lens, is simply one instance of a piecewise smooth equation, and the tools developed for such equations apply directly.
Policy Iteration as a Newton Method¶
Comparing the Newton step above with Algorithm 3, the linearization happens in the policy evaluation step 3, where we solve , and the improvement step recomputes the point at which the operator is linearized. The linear system can be solved by direct methods (Gaussian elimination), by successive approximation, or by Krylov methods such as GMRES. Solving it only approximately, with a few sweeps of successive approximation, is Algorithm 4: an inexact Newton method whose direction comes from a truncated Neumann series.
Example: Sailing to a Windward Mark in a Shifting Wind¶
How should a sailor beat to a mark upwind when the wind direction keeps changing, and do the rules of thumb that racing sailors learn on the water follow from the Bellman equation?
A sailboat cannot sail straight into the wind. To reach a mark that lies upwind, it sails at an angle to the wind on one tack, with the wind coming over one side of the boat, and then turns through the wind onto the other tack. The boat is on starboard tack when the wind comes over its right side and on port tack when it comes over its left side. Each change of tack costs time. When the wind shifts, one tack points closer to the mark (it is lifted) and the other points further away (it is headed). Racing sailors are taught to tack when headed and to avoid sailing out to the laylines, the lines from which the mark can be reached on a single tack, too early. Vanderbei (1996) posed this problem as a stochastic dynamic program on a grid. It later became the “sailing” benchmark for online planning in Péret & Garcia (2004), Kocsis & Szepesvári (2006) and Bonet & Geffner (2012), where value iteration supplies the exact optimal values against which sampled planners are scored.
States. The lake is a grid of waypoints. The state is : the boat’s waypoint , with pointing east and north; the direction the wind blows from; and the current tack . The tack must be part of the state because the cost of the next leg depends on it. There are states. The start is the waypoint and the mark is , 16 rows to the north. The mark is an absorbing state with zero cost.
Actions. Each leg goes to one of the 8 neighbouring waypoints, so the action is a heading . Let be the angle in degrees between the wind direction and the heading. The heading points into the wind and is not allowed, and neither is a heading that leaves the lake.
Costs. The cost of a leg is its duration in minutes. Following Vanderbei (1996), a leg between two neighbouring waypoints in a straight line takes minutes close-hauled, minutes on a beam reach and minutes on a broad reach. Dead downwind takes minute, the default of the benchmark code of Bonet & Geffner (2012). A diagonal leg is times longer, and we let it take times as long, which reproduces the range of leg costs, 1 to 8.6 minutes, reported by Kocsis & Szepesvári (2006). Let be the tack after the leg: starboard if the wind comes from the right of the heading, port if it comes from the left, and unchanged when sailing dead downwind. A change of tack wastes minutes Vanderbei, 1996. The cost is
where for N, E, S, W and for the diagonals.
Transitions. The boat moves deterministically to the neighbouring waypoint in direction . At the end of the leg the wind keeps its direction or shifts by according to a Markov chain , so
We compare three chains. In the steady chain the wind never changes. The oscillating chain uses teaching values for a breeze that swings around north: from north it shifts left or right with probability 0.1 each; from any other direction it moves back towards north with probability 0.2, further away with probability 0.1, and otherwise persists. The shifty chain is the default matrix of the benchmark code of Bonet & Geffner (2012), rotated so that the wind it favours blows from the north. It keeps the wind’s direction with probability 0.2 to 0.4 only.
Objective. The goal is the least expected time to the mark,
where is the first leg that ends at the mark. This is a stochastic shortest-path problem: there is no discount, and the horizon is random because the race ends when the mark is reached. A discount factor would make a minute late in the race cheaper than a minute at the start, which is not what a racer wants, so we keep . The Bellman operator
is then not a contraction in the sup norm, and the argument of Proposition 2 does not apply. A different argument does. Some policy reaches the mark with probability one (a proper policy), and every leg costs at least one minute, so a policy that fails to reach the mark has infinite expected cost. Under these two conditions Bertsekas & Tsitsiklis (1991) show that is the unique fixed point of and that value iteration converges to it from any starting function. We run Algorithm 2 with in place of from and stop when minutes. We also run Algorithm 3 from a baseline policy described below.
Figure 7:The optimal policy for the oscillating breeze. (a) Wind from the north, boat on starboard tack. Each arrow is the optimal heading from that waypoint; red arrows change tack. Contours give the expected minutes to the mark (triangle) and the dashed lines are the laylines of the mean wind. (b) The wind has backed to the north-west, heading every starboard boat; almost everywhere the policy tacks onto port, which now points straight north. (c) One race, in which the optimal boat (blue) and the layline sailor (orange) see the same wind on each leg; circles mark tacks. This race is a large win for the optimal boat, 90.9 against 106.9 minutes, while the expected gap is 2.5 minutes. (d) Sup-norm error of value iteration from , its residual, and the error of policy iteration started from the layline sailor. In the online book the replay races the two boats leg by leg with the wind arrow and the arrow field for the current wind, shows the distribution of finishing times over 4000 races, and animates the value-iteration sweeps.
The baseline is a layline sailor who ignores the shifts when deciding where to tack. This sailor stays on starboard tack until reaching the port layline of the mean (north) wind, then sails port tack towards the mark, and does the mirror image on the right side of the course. On each leg the sailor takes the allowed heading on that tack that ends closest to the mark. A second baseline, the certainty-equivalent planner, replans every leg as if the wind it sees now will never change. The table gives the expected time from the start with the wind from the north. All values are exact policy evaluations. For the optimal policy and the layline sailor, the mean of 4000 simulated races per chain agrees with them within four standard errors.
| Wind | Optimal (value iteration) | Certainty-equivalent | Layline sailor | Value-iteration sweeps | Policy-iteration improvements |
|---|---|---|---|---|---|
| Steady wind from the north | 93.5 | 93.5 | 93.5 | 39 | 10 |
| Oscillating breeze (teaching values) | 84.1 | 85.6 | 86.6 | 42 | 6 |
| Shifty breeze (mdp-engine matrix) | 72.4 | 74.3 | 75.2 | 45 | 5 |
Expected time in minutes from the start to the windward mark, 16 rows upwind on a 21 by 21 lake, with the wind from the north at the start. Certainty-equivalent plans each leg as if the wind seen now will never change; the layline sailor tacks only on the laylines of the mean wind. Value iteration starts from zero and stops when the sup-norm residual is below 1e-09 minutes; policy iteration starts from the layline sailor.
In a steady wind the three policies coincide. The best course is one tack, 16 diagonal legs of minutes plus one 3-minute tack, or 93.5 minutes, and the build script checks that value iteration returns exactly this number. Shifts make the race shorter on average for every policy, because a lifted boat can sail straight north at 4 minutes per row instead of . The optimal policy uses them best. In the oscillating breeze it needs 84.1 minutes against 86.6 for the layline sailor; in the shifty breeze, 72.4 against 75.2. The gain is about 3% of the race, and it is a gain in expectation. Over 4000 paired races in the oscillating breeze, the optimal boat finishes first in 62% of them and ties in 8%, so the layline sailor still wins about three races in ten. The certainty-equivalent planner sits in between, at 85.6 and 74.3 minutes, so reacting to the wind seen now captures part of the gain and planning for the wind still to come captures the rest.
Two pieces of sailing lore can be read off the policy. The first is tack on the headers. Take the open water below the mark and at least three columns inside both laylines, and look at every state in which the wind has shifted from north, on either tack. In the oscillating breeze the optimal policy ends the next leg on the lifted tack in all 900 such states; in the shifty breeze it does so in 93% of them. A header leaves the headed boat sailing across the course with no progress towards the mark, so paying 3 minutes to tack is worth it even though the shift may soon reverse. Panel (b) of Figure 7 shows this as a field of red arrows. The second is stay away from the laylines. In the mean wind a starboard boat eight rows below the mark tacks as soon as it is four columns left of the mark, halfway to the layline (panel (a)). The reason is the one sailors give. On the layline itself, any shift hurts: a header forces an extra tack, and a lift means the boat has sailed further than needed. Staying near the middle of the course keeps both options open. For a planner that assumes the wind seen now will persist, like the certainty-equivalent planner, every tack point inside the laylines is equally good, so nothing pushes it away from them. Neither rule is written into the model; both come out of the minimization in .
The two algorithms reach in different ways. Value iteration from computes, at sweep , the least expected cost of legs with no penalty for where the boat ends. For a boat far from the mark, the cheapest legs are cheap reaching and running legs that do not lead to the mark, so stays well below there. The values settle near the mark first, and the region within one minute of spreads outwards from the mark until it covers the lake at sweep 37. The value-iteration mode of the replay shows this. The sup-norm error stays above one minute for the first 36 sweeps and then collapses. In the oscillating breeze the residual reaches exactly zero after 42 sweeps. Under the optimal policy every course reaches the mark within 33 legs whatever the wind does, so value iteration ends like backward induction with a horizon of a few dozen legs rather than converging geometrically. Policy iteration started from the layline sailor needs 6 improvement steps. Each step evaluates a whole policy, the linear system of Algorithm 3 (solved here by iterating that policy’s own Bellman operator), and the error falls from 14.4 minutes to zero.
The model is a teaching version of the real problem. Real wind shifts are usually much smaller than , and they happen in time, not at waypoints; the boat’s speed also depends on wind strength, which is fixed here. The leg times and tack delay are the published values of the benchmark, not measurements of any boat, and the oscillating chain, the lake and the two baselines are our choices. No measured sailing data are used. What carries over is the structure: a cost-to-go that depends on the wind as well as the position, a value function that is the unique fixed point of the Bellman operator without any discount, and sailing rules of thumb recovered as properties of the optimal policy.
Example: Maintaining a Bridge Deck¶
When should a highway agency repair a concrete bridge deck that has not failed yet, and how much does a condition-based rule save over a fixed schedule?
Every highway bridge in the United States is inspected, usually every two years, and the inspector gives the deck a condition rating from 9 (excellent) to 0 (failed). The ratings are reported each year to the National Bridge Inventory (NBI), where the deck rating is item 58. Decks get worse slowly and at random, through cracking, corrosion of the reinforcing steel and spalling. An agency can do nothing, seal and patch the deck, overlay or rehabilitate it, or replace it. Work costs money now; a deck in poor condition costs money later, through the risk of a failure that closes the bridge. The trade-off repeats every year without a natural end, so we model it as an infinite-horizon discounted MDP and solve it with the algorithms of this chapter.
States. The state is the deck rating, observed once a year. Ratings 9 to 5 are kept apart and ratings 4 and below are merged into one poor state, so . Merging the lowest ratings follows the deck model of Saifullah et al. (2024), whose costs we use below.
Actions. In every state : a minor repair (sealing and patching), a major repair (rehabilitation or an overlay), or a new deck.
Transitions. The do-nothing kernel is estimated from NBI records by counting. We use the deck histories assembled by Liu et al. (2023), which contain the item-58 rating of 88,596 decks for every year from 1993 to 2019 and are mirrored in their public code repository. No history in that sample ever goes up, so it holds only decks that received no rating-raising work. With the number of years in which a deck went from rating to rating , the maximum-likelihood estimate is
computed from 2,303,496 year-to-year pairs. The diagonal is large: a deck stays at rating 9 with probability 0.85 per year, at 8 with 0.92, at 7 with 0.97, at 6 with 0.97 and at 5 with 0.98. These correspond to mean stays of 6.5, 12.9, 30.9, 36.4 and 55.7 years. No deck improves without work, so the poor state is absorbing under “nothing”.
The effects of work are teaching assumptions, in the spirit of the component model of Andriotis & Papakonstantinou (2020). Let be the rating right after the work: a minor repair raises the rating by one but cannot fix a poor deck, a major repair raises it by two up to at most 8, and replacement gives a new deck, . The deck then deteriorates for the rest of the year as if nothing had been done:
Costs. Costs are in US dollars per square meter of deck. The work costs are , , and , from Table 7 of Saifullah et al. (2024). Their Table 6 gives an annual failure probability of 0.001 for ratings 9 and 8, 0.005 for 7, 6 and 5, and 0.01 for the poor state, and they charge a failure ten times the rebuilding cost. We take and leave out the further cost they accrue while a failed bridge is out of service. The one-period cost is the work cost plus the expected failure cost of the repaired deck,
A deck left alone thus costs 26.5 per year at ratings 9 and 8, 132.5 at ratings 7 to 5, and 265 when poor.
Discount. We use , a 3% annual rate, as Saifullah et al. (2024) do.
The optimal value is the smallest expected discounted cost from rating , and it solves the Bellman equation with a minimum in place of the maximum,
We solve it twice. Value iteration (Algorithm 2, with a minimum) starts from and stops when with dollars per square meter. With the error shrinks by about 3% per sweep, and the rule is met after 728 sweeps. Policy iteration (Algorithm 3) starts from “do nothing everywhere”, solves the linear system at each step, and stops after three evaluations. Its second policy already repairs every deck rated 7 or worse; the third changes the action at two ratings and is optimal. This is the Newton behavior described above: a handful of linear solves instead of hundreds of contractions. The two methods return the same policy, and the tests check it against all stationary deterministic policies.
We compare the optimal policy with two rules that need no model. The reactive rule does nothing until the deck is poor and then replaces it. The schedule applies a minor repair every years whatever the rating, and replaces the deck when it is poor; we report the best between 1 and 40. Because the schedule depends on the year, we evaluate it exactly on the Markov chain over the pair (year modulo , rating).
Figure 8:(a) Do-nothing transition probabilities between deck ratings, counted from NBI records of 88,596 decks, 1993 to 2019; entries below 0.005 are left blank. (b) Optimal action at each rating for failure costs from 2 to 40 times the replacement cost; the framed column is the base case. (c) Expected discounted cost from each rating in the base case for the optimal policy, the best schedule (a minor repair every 9 years) and the reactive rule. In the online book the figure is a replay: two groups of 12 decks start from the same ratings and use the same random numbers for 50 years, one under the optimal policy and one under a baseline, with crews marking each repair and counters of the discounted cost so far.
In the base case the optimal policy does nothing at ratings 9 and 8, applies a minor repair at ratings 7, 6 and 5, and a major repair when the deck is poor. It never replaces a deck. The policy is a control limit: work starts at the first rating where the failure risk jumps, from 26.5 to 132.5 dollars per square meter per year, and the work gets heavier as the deck gets worse. A minor repair at rating 7 costs 400 and moves the deck back to 8, which saves about 106 per year for as long as it stays there, a mean of 12.9 years. For a new deck the expected discounted cost is 1,864 dollars per square meter, against 2,951 for the best schedule and 3,200 for the reactive rule.
| Failure cost (× replacement) | Optimal action at ratings 9, 8, 7, 6, 5, ≤4 | Optimal | Best schedule | Reactive |
|---|---|---|---|---|
| 2 | nothing, nothing, nothing, nothing, nothing, nothing | 892 | 1,141 (every 40 yr) | 1,090 |
| 5 | nothing, nothing, minor, minor, nothing, major | 1,841 | 2,248 (every 40 yr) | 2,223 |
| 10 (base) | nothing, nothing, minor, minor, minor, major | 2,311 | 3,623 (every 9 yr) | 4,112 |
| 20 | nothing, nothing, minor, minor, minor, major | 3,242 | 5,615 (every 6 yr) | 7,889 |
| 40 | nothing, nothing, minor, major, minor, replace | 5,101 | 8,821 (every 4 yr) | 15,444 |
Expected discounted cost in USD per m² of deck (γ = 0.97), averaged over the first-year NBI ratings of the decks in the sample. Schedule is a minor repair every k years, with the best k, plus a replacement whenever the deck is poor; reactive does nothing until the deck is poor and then replaces it. Do-nothing transitions: GitHub mirror of NBI records.
The table shows that the answer is driven by the price of failure. When a failure costs only twice the replacement cost, the best policy is never to touch the deck. The poor state then costs 53 per year, and no repair pays for itself at a 3% discount rate. The control limit is not guaranteed either. With a failure cost of five times replacement, the optimal policy repairs decks rated 7 and 6 but leaves decks rated 5 alone. A minor repair moves a deck from 5 to 6, where the failure probability is the same, so it only delays the arrival in the poor state. With a failure cost of 40 times replacement, a deck rated 6 gets a major repair and a deck rated 5 a minor one, for the same reason: a major repair lifts 6 to 8, out of the higher-risk band, while from 5 it would only reach 7. Monotone optimal policies follow from conditions on the costs and transitions Puterman, 1994, and the step-shaped failure table breaks them. In every row the optimal policy beats both rules, and the gap grows with the failure cost, to 3,720 dollars per square meter against the best schedule and 10,344 against the reactive rule at 40 times replacement.
The replay shows the same comparison on sample paths. In the base case, after 50 years the 12 decks under the optimal policy have cost 1,863 dollars per square meter on average, in discounted terms, against 3,196 under the reactive rule and 2,880 under the best schedule. The schedule wastes money on decks that are still rated 9 or 8, while the reactive rule lets decks sit at ratings 7 to 5, where the yearly risk is five times higher, for decades.
Several parts of this model are real and several are not. The do-nothing kernel comes from real inspection records, but it is pooled over all decks regardless of age, traffic, climate and material, and it is estimated from a sample that excludes every deck whose rating went up. Decks that deteriorated fast are the ones most likely to have been repaired, so the mean stays above are probably too long. The work costs and failure probabilities are published values that their authors call indicative, the effects of the repairs are assumptions, and the cost has no term for traffic delays or posting a weight limit. Inspections are treated as exact, while real ratings are noisy, which turns the problem into the partially observed one studied by Andriotis & Papakonstantinou (2020). Finally, an agency manages thousands of decks under one budget, which couples the per-deck problems that we solved separately here.
Example: Controlling an Invasive Plant in a River Network¶
Which stretch of a river should be treated each year to keep an invasive plant in check, and how large can the river get before an exact answer is out of reach?
Tamarisk, or saltcedar, is a shrub that has invaded the banks of many rivers in the western United States. It crowds out native trees such as willows and cottonwoods, and its seeds travel with the water, so an infestation upstream reseeds everything below it. Land managers can kill the plants in a stretch of river, or kill them and plant natives in their place, but each treatment costs money and only part of the network can be treated each year. Hall et al. (2012) formulated this problem as an MDP on a river network, and Alkaee Taleghan et al. (2015) used it as a benchmark for planning with a simulator. We build a small version of their model, solve it exactly with value iteration and policy iteration, and then watch the exact methods run out of room as the river grows.
States. The river is a binary tree of reaches, the stretches of river between two confluences. Water flows from the leaves toward reach 1, the outlet. Each reach has one site for a tree, as in the smallest instances of Alkaee Taleghan et al. (2015), and the site is empty (), holds a native tree () or holds tamarisk (). The state is , so .
Actions. For each reach the manager chooses : do nothing, eradicate, or restore. Eradication kills a tamarisk plant with probability . Restoration first eradicates and then plants a native tree in the empty site, which takes with probability . Without a budget the joint action ranges over choices. With the budget of Alkaee Taleghan et al. (2015), at most one reach is treated per year,
Transitions. A year has three steps. First the treatment acts on each reach independently, with a kernel . Then every plant dies with probability 0.2, with kernel . Last, the survivors reseed the empty sites. Each surviving plant releases 100 seeds, and a seed from reach lands in reach with probability
where and count the downstream and upstream steps on the path from to , and the weight that remains flows out of the outlet. The values 0.5 and 0.1 are the downstream and upstream dispersal parameters of Alkaee Taleghan et al. (2015); the geometric form is our simplification of the model they cite. Each empty reach also receives on average 1 tamarisk seed and 4 native seeds from outside the network. We treat the seeds that reach as Poisson with means
An empty site stays empty if no seed arrives. Otherwise one seed takes it, and tamarisk wins with probability , where 1.5 is one of the values of tamarisk’s competitive advantage examined by Hall et al. (2012). An occupied site keeps its plant. Call the resulting distribution of reach next year. Given the survivors, the reaches fill independently, so
with and Kronecker products of matrices, and a dense matrix. The coupling between reaches enters only through .
Costs and discount. Following Alkaee Taleghan et al. (2015), each invaded reach costs 1.0 for the damage plus 0.1 for the tamarisk site it holds, and the treatments cost for nothing, eradicate and restore:
The costs are in abstract units. The goal is to minimize the expected discounted cost .
Solution. Value iteration (Algorithm 2, with a minimum in place of the maximum) computes
and stops when with . The factored form of means that no transition matrix is ever stored. Each sweep computes once and then applies the small factors of to for every action. Policy iteration (Algorithm 3) starts from “do nothing”, solves exactly at each step, and improves greedily. On the 7-reach river, with states and 15 actions, value iteration stops after 148 sweeps and policy iteration after 4 improvements. The two value functions agree to within .
Figure 9:One simulated episode on a 7-reach river that starts fully invaded. Each segment is a reach, colored by what grows there; water flows down to the outlet. A ring marks the reach treated that year, R for restore and E for eradicate. Top: the optimal policy computed by value iteration. Bottom: the heuristic that restores the invaded reach farthest upstream. Both rows use the same random numbers, and the titles give the discounted cost accumulated before that year. In the online book the replay runs year by year, and the reader can change the starting river, the comparison policy and the random seed.
Every policy in the table below is evaluated exactly, by solving its policy-evaluation equations, so the numbers are expectations rather than averages of simulated episodes.
| Policy | Every reach invaded | Two headwater reaches | Outlet only | Uniform start |
|---|---|---|---|---|
| Optimal (value iteration) | 36.89 | 8.67 | 4.29 | 13.30 |
| Restore the most upstream invaded reach | 38.33 | 9.24 | 5.03 | 14.53 |
| Restore the most downstream invaded reach | 40.33 | 9.23 | 5.04 | 14.93 |
| Do nothing | 61.71 | 25.03 | 15.30 | 34.18 |
| Restore every invaded reach (no budget) | 19.83 | 7.34 | 4.81 | 9.05 |
| Optimal with no budget | 18.98 | 6.45 | 4.02 | 8.18 |
Expected discounted cost (γ = 0.9) on the 7-reach river, computed exactly by solving the policy-evaluation equations. The first four policies respect the budget of one treated reach per year; the last two do not. Uniform start averages over all 2187 starting states.
The optimal policy treats some invaded reach in every one of the 2,059 states that contain tamarisk, yet it agrees with the “most upstream first” rule in only 678 of them and with the “most downstream first” rule in 672. From a fully invaded river it restores reach 2, which sits between two headwater reaches and the outlet, rather than a headwater reach or the outlet itself. When a single tamarisk plant stands among natives, it eradicates instead of restoring, because the native seed rain will usually refill the site and the cheaper treatment suffices. These choices lower the expected cost from a fully invaded river from 38.33 under the best heuristic to 36.89, and from 14.53 to 13.30 averaged over all starting states. The budget matters more than the choice of reach. Restoring every invaded reach each year, which the budget forbids, costs 19.83 from a fully invaded river, and the optimal policy without a budget costs 18.98.
The same solvers run on rivers of 1 to 8 reaches (Figure 10). The number of value-iteration sweeps barely changes, from 140 to 150, because it depends on and rather than on the size of the model. Policy iteration needs 3 or 4 improvements at every size. What grows is the work per sweep. The states grow as , and a dense table of would hold numbers with the budget and without it. For 7 reaches without a budget that is about 1010 numbers, or 84 GB. The factored backup avoids this table, but the coupling matrix alone has entries: 43 million for 8 reaches, and about , or 28 GB, for 10. On the machine that built the figures, the 8-reach model took about 10 s to build and a few seconds to solve, and the 7-reach model without a budget took several times as long to solve as with it. Timings vary from run to run; the growth does not. The same blow-up hits the actions. The budget keeps the joint action set at , but without it the minimization in each backup ranges over joint actions.
Figure 10:Left: the number of states, the number of joint actions with the one-treatment budget, and the number of entries in a dense table of without the budget, against the number of reaches . Right: seconds to build the model and solve it on one machine. Times below 0.1 s are dominated by overhead.
The model is a teaching version of the benchmark. The parameter values for mortality, seed output, dispersal, treatment success, costs and the discount factor come from Alkaee Taleghan et al. (2015), and the competitive advantage from Hall et al. (2012), but the river, the one-site reaches, the geometric dispersal kernel, the Poisson seed rain and the order of events within a year are our choices, and no field data enter the model. Real rivers have hundreds of reaches and many sites per reach. Alkaee Taleghan et al. (2015) describe the full dynamics as a dynamic Bayesian network too complex for exact inference, and plan from a simulator instead. The example still shows the point that the rest of the book builds on. Value iteration and policy iteration converge in a number of steps that hardly depends on the size of the river, but each step touches every state, and the number of states is exponential in the number of reaches. Managing a real river needs value functions that are represented compactly and backups that are estimated from samples, which is the subject of approximate Bellman equations and fitted Q-iteration.
Summary and Outlook¶
Discounting converts an infinite stream of rewards into a bounded value and makes the Bellman maps contractions in the sup norm. Value iteration applies successive approximation; policy iteration alternates a fixed-policy linear solve with greedy improvement, and modified policy iteration interpolates between the two by truncating that solve. Value iteration’s stopping rule comes with a certificate of -optimality, and the span seminorm gives a sharper one. The inference-scheduling example used that linear solve as a check: no heuristic policy beat the value-iteration policy in any of the 49 states. It also showed the price of state reduction: a smaller MDP averages over distinctions that can still affect future completions.
The hard maximum in the Bellman operator selects one action and is nondifferentiable at ties. Can the decision rule remain stochastic while the operator becomes smooth? Entropy-regularized dynamic programming answers by placing a convex regularizer on the action distribution.
Self-checks¶
Solution to Exercise 1
At most . The discounted Bellman operator is a -contraction in the sup norm.
Solution to Exercise 2
.
Solution to Exercise 3
Policy evaluation solves . Policy improvement computes action values from that and chooses a greedy action in each state.
Solution to Exercise 4
Two states can have the same numbers of waiting and decoding requests while their active requests have different remaining output lengths. The reduced state discards those lengths, so its transition kernel averages over them.
- Puterman, M. L. (1994). Markov Decision Processes: Discrete Stochastic Dynamic Programming. John Wiley & Sons.
- Sutton, R. S., & Barto, A. G. (2018). Reinforcement Learning: An Introduction (2nd ed.). MIT Press. http://incompleteideas.net/book/the-book-2nd.html
- Howard, R. A. (1960). Dynamic Programming and Markov Processes. MIT Press.
- Puterman, M. L., & Brumelle, S. L. (1979). On the Convergence of Policy Iteration in Stationary Dynamic Programming. Mathematics of Operations Research, 4(1), 60–69.
- Puterman, M. L., & Shin, M. C. (1978). Modified Policy Iteration Algorithms for Discounted Markov Decision Problems. Management Science, 24(11), 1127–1137.
- Gargiani, M., Zanelli, A., Liao-McPherson, D., Summers, T. H., & Lygeros, J. (2022). Dynamic Programming Through the Lens of Semismooth Newton-Type Methods. IEEE Control Systems Letters, 6, 2996–3001. 10.1109/LCSYS.2022.3181213
- Vanderbei, R. J. (1996). Optimal Sailing Strategies. Statistics and Operations Research Program, Princeton University. https://vanderbei.princeton.edu/sail/sail.html
- Péret, L., & Garcia, F. (2004). On-Line Search for Solving Markov Decision Processes via Heuristic Sampling. In R. L. de Mántaras & L. Saitta (Eds.), Proceedings of the 16th European Conference on Artificial Intelligence (ECAI) (pp. 530–534). IOS Press.
- Kocsis, L., & Szepesvári, C. (2006). Bandit Based Monte-Carlo Planning. Machine Learning: ECML 2006, 4212, 282–293. 10.1007/11871842_29
- Bonet, B., & Geffner, H. (2012). Action Selection for MDPs: Anytime AO* Versus UCT. Proceedings of the AAAI Conference on Artificial Intelligence, 26(1), 1749–1755. 10.1609/aaai.v26i1.8369
- Bertsekas, D. P., & Tsitsiklis, J. N. (1991). An Analysis of Stochastic Shortest Path Problems. Mathematics of Operations Research, 16(3), 580–595. 10.1287/moor.16.3.580
- Saifullah, M., Papakonstantinou, K. G., Bhattacharya, A., Stoffels, S. M., & Andriotis, C. P. (2024). Multi-agent deep reinforcement learning with centralized training and decentralized execution for transportation infrastructure management. arXiv preprint arXiv:2401.12455.
- Liu, H., Nehme, J., & Lu, P. (2023). An application of convolutional neural network for deterioration modeling of highway bridge components in the United States. Structure and Infrastructure Engineering, 19(6), 731–744. 10.1080/15732479.2021.1979597
- Andriotis, C. P., & Papakonstantinou, K. G. (2020). Deep reinforcement learning driven inspection and maintenance planning under incomplete information and constraints. arXiv preprint arXiv:2007.01380.
- Hall, K., Alkaee Taleghan, M., Albers, H. J., Crowley, M., & Dietterich, T. G. (2012). Managing Invasive Species in a River Network. Third International Conference on Computational Sustainability (CompSust’12). https://www.cs.ubc.ca/~crowley/papers/compsust2012.pdf