α-min: A Compact Approximate Solver For Finite-Horizon POMDPs

1 Introduction

Most Partially Observable Markov Decision Process (POMDP) algorithms focus on providing near optimal solutions to infinite-horizon POMDP problems. Near-optimal performance comes at the cost of providing solutions requiring many alpha vectors even when solving small size POMDP problems [Poupart et al., 2011]. In applied fields such as conservation biology, POMDP solutions are often too complex to be analyzed and communicated to managers to be implemented [Tulloch et al., 2015]. Indeed, deriving simple management rules has proven difficult when the POMDP solution has a large number of alpha vectors [Nicol and Chad`es, 2012]. As a result, the solutions must be explored by simulations in order to derive general rules of thumb [Chad`es et al., 2011]. In the best case, these rules are tested via simulations and the loss of performance is reported. This practice is time consuming and does not offer performance guarantees. There is a need to provide an alternative for POMDP users who require simple approximate solutions. To tackle this issue, we develop α-min, a finite-horizon POMDP solver that calculates a good policy given a limit on the number of α-vectors.

Section 2 provides an overview of POMDPs. Section 3 introduces the principles of our approach that relies on calculating α-vectors that minimize the gap between a tight upper bound and a current lower bound approximating the optimal value function. In section 4, we provide a MILP formulation of our approach and present the properties of a first algorithm ϵ-min that finds approximate solutions given a fixed maximum gap. From ϵ-min, we then derive α-min that finds good approximate solutions given a fixed number of α-vectors. Section 5 assesses the performance of α-min on four small benchmark problems, on a larger randomly generated problem, and on a novel computational sustainability problem that seeks to best allocate resources to protect a meta-population of threatened Sumatran tigers. Finally, we discuss the results and future works in Section 6.

Table 3 contains the main notations of this paper. Proofs and proof sketches are available as supplementary material 1 .

2 POMDP Overview

POMDPs are a convenient model for solving sequential decision-making optimization problems when the decisionmaker does not have complete information about the current state of the system [Sigaud and Buffet, 2013]. Formally, a discrete finite-horizon POMDP is specified as a tuple {S, A, O, τa, Ωa, R, H}, where

• H = {0, ..., T -1}, T ∈ N, is the time horizon. Elements of H are called time-steps and T is the number of time-steps.
• S, ∀t ∈ H, st ∈ S is the state of the system at t.
• A, ∀t ∈ H, at ∈ A is the taken action at t.
• O, ∀t ∈ H, zt ∈ O is the observation at t.

1 https://sites.google.com/site/ijcaialphamin/home

Iadine Chad`es CSIRO iadine.chades@csiro.au

• τa is the transition matrix for action a. Elements are τa (st, st+1).
• Ωa is the observation matrix for action a. Elements are Ωa(st+1, zt+1).
• R is the reward matrix. Elements are R (at, st).

For sake of clarity, we define the following notation:

• For action a ∈ A and for observation z ∈ O, let Ma,z be the matrix of dimension S × S such that Ma,z(st+1, st) = Ωa(z, st+1)τa(st+1, st).
• For every a ∈ A, the vector ra = R(a, .) corresponds to the row of the matrix R corresponding to the action a.

The optimal decision at time t may depend on the complete history of past actions and observations. Because it is neither practical nor tractable to use the history of the actionobservation trajectory to compute an optimal solution, belief states (also called beliefs), i.e. probability distributions over states, are used to summarize and overcome the difficulties of imperfect detection [ ˚Astr¨om, 1965]. A POMDP can be cast into a fully observable Markov decision process defined over the continuous belief state space, B.

Solving exactly a finite-horizon POMDP means finding an optimal policyΠ0 = ∪t∈Hπt, where πt : B →Amaps belief states at time t to actions. Π0 maximizes the expected sum of rewards E[∑t∈H rat · bt] over the time horizon H (· denotes the scalar product). For each time-step t, for a given belief state bt and a given policy Πt = ∪t′∈{t,...,T-1}πt′ the expected sum E[∑t′∈{t,...,T-1} rat · bt] is also referred to as the value function Vt,Πt (bt). A value function allows us to rank strategies by assigning a real value to each belief bt. An optimal policy Πt is a policy such that, ∀ bt ∈ B, ∀ Π′t , Vt,Πt (bt) ≥ Vt,Π′t (bt). Several strategies can be optimal and share the same optimal value function Vt, which can be computed using Bellman's principle of optimality [Bellman, 1957]: ∀bt ∈ B,

where the belief bt+1 can be computed as follows : ∀st+1 ∈ S,

Equation 1 can be rewritten Vt = BL(Vt+1) where BL is the Bellman operator [Shani et al., 2012], sometimes also called backup operator [Pineau et al., 2006]. While various algorithms from the operations research and artificial intelligence literature have been developed over the past years, exact resolution of POMDPs is intractable: finite-horizon POMDPs are PSPACE-complete [Papadimitriou and Tsitsiklis, 1987] and infinite-horizon POMDPs are undecidable [Madani et al., 2003].

3 α-min principle

Equation 1 can be solved by directly manipulating α-vectors [Smallwood and Sondik, 1973]. For every t ∈ H, there exists a finite set Γt of vectors of dimension |S| (the so-called αvectors) which define entirely Vt such as |Γt| is minimal and:

This formulation is equivalent to equation 1. Given that ΓT-1 = {ra| a ∈ A, ra is not dominated}, one can in theory build a set of α-vectors defining the value function Vt from Γt+1 at any time-step t ∈ H′ = {0, ..., T - 2}, because every αt ∈ Γt can be written:

where at ∈ A and αat,zt+1 t+1 , zt+1 ∈ O are elements of Γt+1.

The set of α-vectors P(Γt+1) given by Equation 4 is such that Γt ⊆ P(Γt+1) (α-vectors of P(Γt+1) potentially belong to Γt). In the case of exact resolution, for a given Γt+1, a natural way to compute Γt is first to compute P(Γt+1) entirely, and then prune the dominated vectors that are not useful for representing the value function. In practice, this approach is computationally expensive [Sigaud and Buffet, 2013, sections 7.3 and 7.4].

Point-based algorithms are recent approximate approaches to solve POMDPs [Shani et al., 2012]. Value functions are updated according to a subset of beliefs ˜Bt ⊆ B sampled to be as relevant as possible to get good approximations of Γt at each time-step. Every alpha vector αbt t of the current approximation ˜Γt ⊆ Γt, corresponding to the belief bt ∈ ˜Bt, is generated as follows:

where ˜Γt+1 is an approximation of Γt+1.

Definition 1. We call α the vector-function such that for every bt ∈ B, α(bt) = argmaxαt∈P(˜Γt+1) αt · bt. Note that for every bt ∈ B, we have α(bt) · bt = BL( ˜ Vt+1)(bt).

Using the point-based approach on a finite-horizon, one can build a lower approximation of the optimal value function at each time-step. We start by setting ˜ΓT-1 = ΓT-1 and for every t ∈ H′ , we build ˜Γt = {α(bt)|bt ∈ ˜Bt}. We call ˜ Vt the corresponding value function (according to Equation 3 where Γt, Vt and B are respectively replaced by ˜Γt, ˜ Vt and ˜Bt). We have ˜Γt ⊆ Γt, t ∈ H, thus every ˜ Vt is potentially only a lower bound of Vt. If ˜Γt+1 = Γt+1 (i.e. ˜ Vt+1 = Vt+1), then the maximal error on Vt is by definition gapt = maxbt∈B(BL( ˜ Vt+1)(bt) - ˜ Vt(bt)). By induction from ˜ VT-1 = VT-1, which does not have any error, and according to Lemma 1, the maximal error on any ˜ Vt is ∑t′∈{t,...,T-1} gapt′ . Thus the maximal error on V0 at t = 0 is bounded by gap = ∑t∈H gapt for any b0. This measure of error is usual in approximation approaches [Hansen, 1998; Hauskrecht, 2000].

Lemma 1. Let t ∈ H′ and ˜ Vt+1 be an approximate lower representation of Vt+1 such that: ∀ bt+1 ∈ B, ˜ Vt+1(bt+1) ≤

Vt+1(bt+1) ≤ ˜ Vt+1(bt+1) + ϵ. Then, ∀ bt ∈ B, ˜ Vt(bt) ≤ Vt(bt) ≤ ˜ Vt(bt) + ϵ, where ˜ Vt = BL(˜ Vt+1) and Vt = BL(Vt+1).

Typical point-based methods sample belief states by simulating interactions with the environment and then updating the value function over a selection of those sampled belief states [Pineau et al., 2006; Shani et al., 2006]. Our approach consists, for every given time-step t, in searching iteratively for belief states b∗t of Equation 6, in order to improve ˜ Vt .

The belief point set expansion consisting of adding b∗t to ˜Bt and the improvement of ˜ Vt consisting of adding α(b∗t ) to ˜Γt aim to reduce the current gap as much as possible iteratively (Figure 1).

Several algorithms have been proposed to approximately solve POMDPs. For example, Roy et al. [2005] aimed to reduce the dimensionality of the belief space using principal components analysis, rather than minimize the size of the solution itself. Poupart and Boutilier [2004] sought to reduce the dimensionality of the belief space and to apply the "bounded policy iteration" method [Poupart and Boutilier, 2003] based on bounded stochastic finite state controllers. In both cases [Poupart and Boutilier, 2004; Roy et al., 2005], the methods do not provide an intrinsic performance guarantee, i.e. no gap is provided between the approximate and the optimal solutions. One must use an external algorithm to compute this gap.

Unlike point-based sampling approaches [Shani et al., 2012] or other approximation approaches [Hauskrecht, 2000; Poupart and Boutilier, 2004; Roy et al., 2005], our method generates iteratively new non-dominated α-vectors by solving Equation 6 directly.

4 AMILP approach

Given a time-step t ∈ H′ , solving Equation 6 exactly is an NP-hard optimization problem. This is a direct consequence of the NP-hardness of the exact backup operation [Littman et al., 1995]. Indeed, otherwise we could use Equation 6 to construct an exact representation of the optimal value function Vt given the exact value function Vt+1, in a time which is a polynomial in the instance and the number of α-vectors needed to describe Vt.

Equation 6 can be formulated as the following quadratic program QP :

where the variables are αt and bt. The first term αt · bt of the first constraints corresponds to the α-vector formulation of BL( ˜ Vt+1)(bt) and ˜αt · bt corresponds to the α-vector formulation of ˜ Vt(bt). Note that ˜αt are not variables but coefficients of the |˜Γt| first constraints.the first type.

The problem QP is difficult to solve, not only because it is a non-concave quadratic program, but mainly because P(˜Γt+1) is not known. In the remainder of the paper, we will demonstrate how we can reformulate QP into a MILP that can be solved efficiently.

Because the last constraint αt ∈ P(˜Γt+1) is not expressed as a set of linear inequalities, we have still to find a set of such linear inequalities describing Convt, the convex hull of P(˜Γt+1). Since we are in a maximization case, we can focus on describing Ct, the polyhedron composed of hyperplanes of Convt with positive normal vectors and hyperplanes αt,s = 0, s ∈ S, where αt,s are the components of αt. Note that extreme points of Ct with strictly positive coordinates correspond to non-dominated α-vectors of P(˜Γt+1), while interior points of P(˜Γt+1) correspond to dominated αvectors.

Let us describe Ct using an infinite number of constraints: Ct = {αt ∈ R|S| | αt · bt ≤ αt(bt) · bt, bt ∈ B,αt ≥ 0}. We can approximate Ct with a finite number of constraints by considering the convex polyhedron ˜ C n t = {αt ∈ R|S| | αt ·˚bit ≤ αt(˚bit ) ·˚bit , i = 0, ..., n-1, αt ≥ 0} where ˚b0t , ...,˚bn-1 t are n ≥ |S| beliefs of B. We have Ct ⊆ ˜ C n t .

For a given C n t , instead of solving QP, we can now solve the easier quadratic program QPn(˚b1t , ...,˚bn-1 t ), which can provide an approximation of QP with guaranteed performance (Proposition 2):

Proposition 2. Let ˆαnt , ˆ˜αnt , ˆbnt , and ˆ gn t be an optimal solution of QPn.

This formulation allows us to decide whether to expand the belief set by adding ˆbnt to ˜Bt in the case where δn t is small enough, or to construct, from the current solution of QPn, a new relevant belief˚bnt defining a new constraint of ˜ C n t , in order to get a better approximation ˜ C n+1 t of Ct.

Our aim is then to find a good description of Ct, i.e. using few hyperplanes. It is well known in linear programming, and particularly in polyhedral approaches, that the best possible hyperplanes describing a convex polyhedron correspond to its so-called facets [Mahjoub, 2014]: the set of these facets is indeed the minimal set of hyperplanes needed to describe a convex polyhedron and is always finite. The problem of finding relevant facets is known as separation problem [Mahjoub, 2014]. Gr¨otschel et al. [1981] showed that the cost of optimization on a given polyhedron does not depend on the number of constraints of the system describing the polyhedron, but rather on the separation problem associated with this system. In our case, optimizing over Ct means finding a new relevant non-dominated α-vector and its associated belief. Unfortunately, since QPn is quadratic, one cannot directly apply the results of Gr¨otschel et al. [1981] to solve QPn.

An important step in proving the convergence of our proposed algorithms is to solve the following separation problem, formally: compute a new facet F n t of Ct from any current solution ˆαnt of QPn which is not already in Ct. This can be done by using the algorithmGenerateFacet(ˆαnt , ˜Γt+1) (Proposition 3).

Proposition 3. Given a set ˜Γt+1 of α-vectors at time-step t representing the value function ˜ Vt+1, and a vector αt ∈ R|S| , one can decide if αt belongs to Ct or not, and if not, generate a facet of Ct. We call the corresponding algorithm GenerateFacet(αt, ˜Γt+1).

Each time we solve QPn and δn t is not satisfying, we can generate a new facet F n t from ˆαnt using GenerateFacet(ˆαnt , ˜Γt+1). We then can set ˚bnt = bF n t t to add the corresponding constraint to ˜ C n t (which becomes ˜ C n+1 t ), where bF n t t is the belief corresponding to the facet F n t . In doing this, ˜ C n t converges to Ct (Proposition 4).

Proposition 4. ˜ C n t converges to Ct since at each iteration we compute a new facet of Ct, and the convex polyhedron Ct has a finite number of facets. Thus, there exists n∗ ∈ N such that the solution ˆαn∗ t of QPn∗ belongs to Ct, i.e. δn∗ t = 0. Therefore, for a given precision ϵp > 0 there exists nϵp ∈ N such that δnϵp t ≤ ϵp.

QPn is however still not easy to solve, because the objective function is clearly not concave, due to the term αt · bt. But, since ˆαnt necessarily belongs to one of the extreme points of ˜ C n t , it satisfies the System 7:

Let β1, ..., β|S| ∈ [0, 1] such that ˆbnt = ∑j∈1,...,|S| βj˚bt ij and ∑j∈1,...,|S| βj = 1. Such βj and˚bij t always exist if we assume that there are at least |S| beliefs˚bit affinely independent (for example, by setting the |S| first˚bit to the extreme points of the simplex B). αt · bt can be easily re-written ∑j=1,...,|S| βij (αt(˚bij t ) ·˚bij t ), which is a linear expression. This technique of linearization is similar to considering interpolations of BL( ˜ Vt+1) as an upper bound [Poupart et al., 2011]. Our approach differs from [Poupart et al., 2011] because we calculate a tight upper bound in order to guarantee a near best possible improvement of the lower bound instead of computing upper bound and lower bound with two independent heuristics.

Thus, we can finally reformulate QPn into a MILP (Proposition 5).

Proposition 5. Let b0t , ..., bn-1 t be n ≥ |S| beliefs. Problem QPn(b0t , ..., bn-1 t ) can be reformulated as the following MILP, called MILPn(b0t , ..., bn-1 t ), containing only continuous variables except n that are 0-1 variables.

where M ≥ T ×Rmax with Rmax the maximum reward over the matrix R. M is then an upper bound of the term αt ·˚bit .

Any optimal solution of MILPn is also an optimal solution of QPn. We keep the same notation for the solutions: (ˆαnt , ˆ˜αnt , ˆbnt , ˆ gn t ).

We now have all the tools to solve Equation 6 with bounded error (Proposition 6). The procedure to find the corresponding near optimal belief (FBB) is given in Algorithm 1.

Proposition 6. Given a time-step t, ˜ Vt, ˜ Vt+1 and a representation ˜Bt of the belief space B, one can solve approximately the optimization problem corresponding to Equation 6 to within specified precision ϵp, using Algorithm 1, called FBB (Find Best Belief).

Algorithm 1 Find the best belief for expanding ˜Bt according to Equation 6 to within specified precision ϵp

1: procedure FBB(t, ˜Bt, ˜Γt, ˜Γt+1, ϵp) 2: n← 0 3: for s ∈ S do 4: ˚bnt ← es ▷ es, s ∈ S are extreme points of B 5: n ← n + 1 6: while ∆> ϵp do 7: (ˆαnt ,ˆbnt , ˆ gn t ) ← Solve MILPn(˚b0t , ...,˚bn-1 t ) 8: ∆ ← δn t = ˆαnt · ˆbnt -α(ˆbnt ) · ˆbnt 9: F n t ←GenerateFacet(ˆαnt , ˜Γt+1) 10: ˚bnt ← bF n t t 11: n ← n + 1 return ˆ gn t ,ˆbt

Proposition 7. In the worst case, Algorithm 1 requires O(PN ×2N+|S|+1 ) operations, where N designates the maximum number of facets needed to describe Ct and PN is a polynomial in N and the size of the POMDP.

Theorem 8. In finite time, one can solve approximately any finite-horizon POMDP with a specified arbitrary maximum gap ϵ, using Algorithm 2 by setting ϵp ≤ ϵ.

Algorithm 2, which we call ϵ-min, aims to provide a compact solution under the constraint of respecting a gap less than or equal to a fixed parameter ϵ uniformly spread across each time-step. ϵ-min is attractive because one can specify a required maximum gap, but ϵ-min may lead naturally to the calculation of a large number of α-vectors per time-step, since the algorithm adds new α-vectors until the required maximum gap is reached. This behavior is not suitable when looking for simple solutions.

However, we can easily adapt ϵ-min to constrain the maximum number of α-vectors to use per time-step. We call this algorithm α-min (Algorithm 3).

Theorem 9. In finite time, one can solve approximately any finite-horizon POMDP with a specified arbitrary maximum number N of α-vectors per time-step using Algorithm 3 by setting ϵp small enough. Algorithm 3 provides a maximum gap between the solution it computes and an optimal solution.

1: procedure ϵ-MIN(H, A, S, O, ϵ, ϵp) 2: ˜ΓT-1 ← {ra| a ∈ A, ra is not dominated} 3: gapT-1 = 0 4: for t ∈ H′ do 5: ˜Bt = binit t ▷ binit t is an arbitrary belief 6: gapt ← ∞ 7: for t = T - 2, T - 3, ..., 0 do 8: while gapt > ϵ T-1 do 9: (gapt, b∗t ) ←FBB(t, ˜Bt, ˜Γt, ˜Γt+1, ϵp T -1 ) 10: ˜Bt ← ˜Bt ∪ {b∗t } 11: ˜Γt ← ˜Γt ∪ {α(b∗t )} 12: gap ← ∑ t∈H′ gapt 13: return {˜Γt, t ∈ H}, gap

Algorithm 3 α-min: Solve POMDP with a maximum number N of α-vectors, to within precision ϵp

1: procedure α-MIN(H, A, S, O, N, ϵp) 2: Lines 2 to 7 of Algorithm 2 3: while |˜Γt| < N and gapt > ϵp T-1 do 4: Lines 9 to 13 of Algorithm 2

5 Experiments

We first assess the performance of α-min on four small finitehorizon POMDP problems from the literature 2 and a larger randomly generated problem (random30). We compare its performance to a leading infinite-horizon POMDP solver Sarsop [Kurniawati et al., 2008]. Sarsop results were obtained by solving the POMDPs over an infinite-horizon with γ = 0.999 and a maximum computational time of 1000s. Sarsop lower bounds (LB) were calculated as the expected sum of rewards cumulated over T time-steps by simulation of the infinite policy using a γ=1. α-min results were obtained using a fixed number of α-vectors set arbitrarily with a maximum computational time of 1000s per time-step on a 94.4 GB, 3.47GHz, 19 cores computer and CPLEX 12.5. Overall, the performance of α-min is encouraging with surprisingly good gaps obtained considering the small quantity of α-vectors (Table 1). The lower bounds of Sarsop and α-min are close. Table 2 shows the behavior of the lower bounds and gaps generated by α-min as the allowed number of α-vectors is increased for the problemmilos-aaai97.

Finally, we illustrate the benefits of using α-min to best manage or monitor four sub-populations of Sumatran tigers with declining connectivity over a 10-year time horizon [Linkie et al., 2006]. We model this problem as a non-stationary finitehorizon POMDP with 16 states representing the status extinct or extant of each sub-population; 13 actions representing the decisions of doing nothing, managing and/or monitoring each sub-population; 16 observations representing the difficulty of detecting Sumatran tigers in each sub-population (absent or present); and 10 transition matrices representing the declining connectivity between subpopulations over time. The probabilities of going extinct with and without management were derived based on tiger census estimates [Chad`es et al., 2008; Linkie et al., 2006]. Detection probabilities were derived based on [McDonald-Madden et al., 2011], and projected fragmentation scenarios followed [Linkie et al., 2006]. Interested readers can refer to the supplementary material 3 for the POMDP files corresponding to this problem and to [Chad`es et al., 2008; McDonald-Madden et al., 2011;

2 http://pomdp.org/examples/index.shtml

3 https://sites.google.com/site/ijcaialphamin/home

Regan et al., 2011] for limitations and advantages of using POMDPs in conservation problems. For sake of illustration, Figure 2 presents the policy graph obtained for T=10 and |α|=7 by α-min. In this case, the solution is guaranteed to be at most 10% from the optimal strategy (gap < LB 10 ). The CPU time required to solve the problem was 378 seconds. Note that because this is a non-stationary finite-horizon POMDP, it is not possible to provide a Sarsop solution.

6 Discussion

We proposed two new algorithms, ϵ-min and α-min, to approximately solve finite-horizon POMDPs. Both algorithms rely on finding at each time-step α-vectors to minimize the gap between the optimal value function and its current approximation, as formally expressed by Equation 6. ϵ-min calculates an approximate solution given a maximum gap, which can lead to the calculation of very many α-vectors. α-min is more directly applicable to computational sustainability as it provides good POMDP solutions given a fixed number of α-vectors. α-min can easily be adapted to find approximate solutions given a CPU time limit, or both a time limit and a maximum number of α-vectors.

α-min greedily expands the set of beliefs and the set of αvectors by adding iteratively new beliefs and α-vector which aim to reduce the current gap as much as possible, until N α-vectors have been added. One interesting possible future direction would be to improve α-min in order to generate a set of α-vectors of a given cardinality N with the guarantee that no strictly-better set of the same cardinality exists (to within a given precision).

Unlike most point-based approaches, our algorithms assume that we do not know the initial belief b0. Taking advantage of b0 could be explored in future work.

Our method can also be adapted to the infinite-horizon case, since it is a point-based approach. In our case, it was particularly interesting to propose a solver for finite-horizon POMDPs as it allows us to generate non-stationary policies, given that, in the context of computational sustainability, the transition matrices and rewards might change over time.

Finally, the complexity bound of Proposition 7 could probably be improved. However finding an "efficient" complexity bound, e.g. polynomial in the instance, is unlikely given the non-approximability results for POMDPs in general and for finite-horizon POMDPs in particular [Lusena et al., 2001].