Action-Constrained Markov Decision Processes With Kullback-Leibler Cost
Ana Bušić Sean Meyn

























=1  

  
  



 
  



#1#1
           [t]#1
       #1

#1read notes: #1  

#1SPM:  #1   

#1#1
#1#1











#1<#1>
 


𝒰
𝒰



c_KL

w
W
W̌


w


@wd
@hta
@htb
@box
@measure#1@box=#1@wd=@box @hta=@box

#1@measure#1
  @box=height @hta width@ @wd
  @htb=@box @htb@hta @htb@
  height @box width@ #1@htb1[-1]@box




W̌



DD-r
x^∘

_B_ζ 


diag 
T
τ_s

H^e
G


#1Λ_#1



H^∘
H^∘  

 ↰   
 ▿ 

IPD
SPD
CPD

 







Π

⊃

#1#1
 







⊓⊔
501em=0pt=0






#1[.9][c].85
=0pt
=12pt plus .1pt      
=6pt plus 1.5pt minus 1pt 
 #1




#1

 §.§ 

#1


 §.§ 








T


-

w⟶
v⟶
#1#1⟹
d⟶

#1#1⟶
#1#1


tr  

rank 
deg
sign
supp 



[1.6mm]|||



det
lim suplim sup
lim inflim inf
arg  min

arg  max









A
K
S
U
V
W

X
X_u
X_n
()

Y

Z











 




R_+
R
 


 
I

I


T
Z


Z

Z_+


Q

C




 
π̌ 
P̌
Ž
p̌
μ̌


Ř





A
B
C
D
E
F
G
H
I
J
K
L
M
N
O
P
Q
R
S
T
U
V
W
X
Y
Z

a
b
c
d
e
f
g
h
i
j
k
l
m
n
o
p
q
r
s
t
u
v
w
x 
z






#1



 







a
b
d
e
m
q
r
v 
x
y
u
w
 

  

z  

A
B  

C 
D 
E   
F     
G     
H     
I      
   
     
L
M
N  
Q  
R  
R̅  
 


 
 

S  
T  
U  
 

N
Q
 
 
W


X
X

Q

q̃

y

ỹ




Y   
Y


 
to 0pt1.25pt
 

V  
W  
  
Z  


Z



w
y
z
Q
S
U
W






θ
Φ


η

π_0 

 

ϕ 
ψ


ρ
ζ

Υ 



  





Ψ 

Σ

P






ĉ
ẑ

f̂
ĝ
m̂
ĥ
p̂
q̂



Q

v̂
ŵ
ŷ


F
Î
Ĵ

N̂

Q
T
W

η̂

θ̂
ψ̂
π̂_0
ν̂
μ̂
α̂
β̂
γ̂

λ̂
Λ̂



A
B
C
C
E
H
I
M
P
Q

d
g
k
u
x

P
E








==
Φ
N
P
Q
T
W






μ̃
π̃_0
X̃
m̃
θ̃
ẽ
σ̃
τ̃






A
B
C
D
E
F
G
H
J
I
K
L
M
N
O
P
Q
R
S
T
U
V
W
X
Y
Z






 


Var 
 
 










#1#20cm2#1#2 

#1#20pt#1#2

   
 #1#2#3#1#2#3
 


#1#20#1#2
1#1#2
3#1#2
3#1#2

1ddt
1ddt
3ddt
3ddt


1ddr
1ddr
3ddr
3ddr


1ddy
1ddy
3ddy
3ddy


1d^to 2pt+dt
1d^to 2pt+dt
3d^to 2pt+dt
3d^to 2pt+dt


1ddα
1ddα
3ddα
3ddα

1ddζ
1ddζ
3ddζ
3ddζ

1d^+dζ
1d^+dζ
3d^+dζ
3d^+dζ



112
112
312
312

113
113
313
313

114
114
314
314
 
 134
134
334
334

116
116
316
316

118
118
318
318



dist=

:=

μ^Leb
M



P
P 
E

#1_#1P

i.o.


a.s.
vec  



1/ N

#1,#2,1#2 ∑_#1^#2

   enters  

(#1)(#1)  
ε

 · 



 
theoremTheorem[section]
corollary[theorem]Corollary
proposition[theorem]Proposition
lemma[theorem]Lemma




#1Lemma <ref>
#1Prop. <ref>
#1Theorem <ref>
#1Corollary <ref>


#1Section <ref>
 
#1Chapter <ref>

#1Appendix <ref>
 
  
  




κ


a

c

f
g
h
l
m
n

p

y
A
B
C
E
F
M
P
Q
T


n
ρ


ρ
μ
ν
α
β

α
η
 



Φ
σ

Σ


τ





σ̂
τ̂
P̂
Φ̂

F̂



















     
rmnumx
romannumx(rmnumx)rmnumx






anum
alphanum(anum)anum







 

 

 




   


Δ

γ
Ψ


 

γ

 
 









        


figures/

#1#2



    < g r a p h i c s >
 





S 
 
 
V


H
V
g_p 

#1Fig. <ref>


P^tot

I

R_+
R














 







  

 
 
 
 
 
 

 
  
  
  

 
 
 
 
 


 















Action-Constrained Markov Decision Processes 
  

  With Kullback–Leibler CostFunding from the ANR under grant ANR-16-CE05-0008, and NSF under awards EPCN 1609131, CPS 1646229 is gratefully acknowledged.

    Ana BušićInria and DI ENS, École normale supérieure, CNRS, PSL Research University, Paris, France
    (ana.busic@inria.fr, <http://www.di.ens.fr/ busic/>).
  Sean MeynDepartment of Electrical and Computer
Engineering at the University of Florida, Gainesville (meyn@ece.ufl.edu, <http://www.meyn.ece.ufl.edu/>). 

    December 30, 2023
============================================================================================================================================================================================================================================================================================================================











This paper concerns computation of optimal policies in which the one-step reward function contains a cost term that models Kullback-Leibler divergence with respect to nominal dynamics. This technique was introduced by Todorov in 2007, where it was shown under general conditions that the solution to the average-reward optimality equations reduce to a simple eigenvector problem. Since then many authors have sought to apply this technique to control problems and models of bounded rationality in economics. 

A crucial assumption is that the input process is essentially unconstrained. For example, if the nominal dynamics include randomness from nature (e.g., the impact of wind on a moving vehicle), then the optimal control solution does not respect the exogenous nature of this disturbance. 

This paper introduces a technique to solve a more general class of action-constrained MDPs. The main idea is to solve an entire parameterized family of MDPs, in which the parameter is a scalar weighting the one-step  reward function. The approach is new and practical even in the original unconstrained formulation.





















  
*Keywords:  
 Markov decision processes,  Computational methods.



§ INTRODUCTION


 


Consider a Markov Decision Process (MDP) with finite state space , general action space ,  and  one-step reward  function ×→.   Two standard optimal control criteria are finite-horizon:

    ^*_T(x) = max∑_t=0^T[  (X(t),U(t))  | X(0)=x ]

where T≥ 0 is fixed,   and average reward:

    η^*(x) =max{lim inf_T→∞1/T∑_t=0^T-1[  (X(t),U(t))    | X(0)=x ]  } .

where ={X(t): t≥ 0}, ={U(t): t≥ 0} denote the state and input sequences.
 

In either case, the maximum is over all  admissible input sequences; it is obtained as deterministic state feedback under general conditions. In the average-reward framework the optimal policy is typically stationary:  U(t) = ϕ^*(X(t)) for a mapping ϕ^*→,  and η^*(x) does not depend upon the initial condition x   (see <cit.>).   


 
A special class of MDP models was introduced by <cit.>, for which either   optimal control problem has an attractive solution.  The reward   function is assumed to be the sum of two terms: 

    (x,μ)     = 
    (x) 
    		-  D(μ P_0(x,)).
 




The first term is a function   → that is completely unstructured.   The second is a “control cost”,  defined using   Kullback–Leibler  (K-L) divergence (also known as relative entropy).   The control cost is based on deviation from nominal (control-free) behavior;  modeled by a nominal transition matrix P_0:

    D(μ P_0(x,))∑_x'μ(x') log(μ(x') /P_0(x,x')).
 


It is shown that the solution with respect to the average reward criterion is obtained as the solution to the following eigenvector problem:  let (λ,v) denote the Perron-Frobenius eigenvalue-eigenvector pair for the positive matrix with entries  (x,x') = exp((x)) P_0(x,x'),  x,x'∈.
The eigenvector property
v = λ v implies that the “twisted” matrix 

    (x,x') 
    =
    1/λv(x')/v(x)(x,x')   ,    x,x'∈  .

is a transition matrix on .  This transition matrix defines the dynamics of the model under optimal control.  A similar model was introduced in the earlier work of  <cit.>, but without the complete solution reviewed here. 



























l.28


.95TodCar.pdf

Optimal hill climb
 


 

 

Since the publication of <cit.> there has been significant theoretical advancement, with proposed applications to economics <cit.>, distributed control <cit.>, and neuroscience <cit.>.


It is appealing to imagine that rational economic agents are solving an eigenvector problem to maximize their utility.  However, a careful look at the controlled dynamics (<ref>) suggests a limitation of this MDP formulation:  how can this transformation respect exogenous disturbances from nature?  
An essential assumption in this prior work is that for each x,  and any pmf μ,  it is possible to choose the action so that P(x,x') = μ(x').  This is equivalent to the assumption that  the action space 
 consists of all probability mass functions on , and the controlled transition matrix is entirely determined by the input as follows:

    {X(t+1)=x'| X(t) =x, U(t) = μ} = μ(x') ,     x,x'∈, μ∈  .

This modeling assumption presents a significant limitation, as pointed out in <cit.>:  “It prevents us from modeling systems subject to disturbances outside the actuation space”.


 
 
f:TodCar is based on an example of <cit.>.  Reaching the parking spot at the top of the hill in minimum time (or  minimal fuel) is formulated as a total cost problem, similar to (<ref>).   
The figure has been modified to indicate that wind and rain influence the behavior of the car on the track.    The optimal solution cannot take the form (<ref>) when this additional randomness is included in the model, since this would mean our control action would modify the weather.  







  
Contributions

In this paper the K-L cost framework is broadened to include constraints on the pmf μ appearing in (<ref>).  
The new approach to computation   is based on the solution of an entire family of MDP problems,  parameterized by a scalar ζ appearing as a weighting factor in the one-step reward function.  Letting X(t) denote the state, and R(t) denote the randomized policy at time t, this one-step reward is of the form

    (X(t),R(t)) = ζ(X(t)) -  (X(t),R(t) )

in which  denotes relative entropy with respect to nominal dynamics (see (<ref>)).   

The main results of the paper are contained in Theorems <ref>  and <ref>, 
with parallel results for the total- and average-reward control problems.   In each case, it is shown that   the solution to an entire family of MDPs can be obtained through the solution of a single ordinary differential equation (ODE).   

The ODE solution is most elegant in the average-reward setting.  For each ζ, the solution to the average-reward optimization problem is based on a relative value function h^*_ζ→.  For the MDP with d states, each function is viewed as a vector in ^d with entries {h_ζ^*(x^i) : 1≤ i≤ d}.  A vector field ^d→^d is constructed so that these functions solve the ODE 

    h^*_ζ  =  (h^*_ζ) ,
        
    with boundary condition h^*_0≡ 0.


One step in the construction of  is differentiating each side of the dynamic programming equations; a starting point of the 50 year old sensitivity theory of  <cit.>, and  more recent <cit.>.

  More closely related is the sensitivity theory surrounding Perron-Frobenius eigenvectors that appears in the theory of large deviations <cit.>.
 The goals of this prior work are different, and we are not aware of comparable algorithms that simultaneously solve the family of control problems.

  The optimal control formulation is far more general than in the aforementioned work <cit.>, as it allows for inclusion of exogenous randomness in the MDP model.   The dynamic programming equations become significantly  more complex in this generality, so that in particular, the Perron-Frobenious computational approach used in  prior work is no longer applicable.   
 









 

In addition to its value as a computational tool,  there is a significant benefit to solve the entire collection of optimal control problems for a range of the parameter ζ.  For example, this provides a means to understand the tradeoff between state cost and control effort.     Simultaneous computation of the optimal policies is also an essential ingredient of the distributed control architecture introduced in  <cit.>.   

The ODE algorithm
is easily implemented for  problems of moderate size.   In this paper an example is provided in which the the size of the state space d is greater than 1,000;   the action space is an open subset ^d-1 since actions correspond to randomized decision rules.   The optimal solutions for the desired range of ζ were obtained in less than one hour using a standard laptop running Matlab.


 



 



  
 
The remainder of the paper is organized as follows.   
s:design describes the new Kullback–Leibler cost criterion and numerical techniques for the MDP solutions.    This is applied to a path-finding problem in s:drone.
Conclusions and topics for future research are contained in s:conc.






§ MDPS WITH KULLBACK–LEIBLER COST
 


 
 


 §.§ MDP model
   

The dynamics of the MDP are assumed of the form (<ref>),  where the action space consists of a convex subset of probability mass functions (pmf) on .  An explanation of the one-step reward  (<ref>) will be provided after a few preliminaries.


A transition matrix P_0 is given that describes  nominal (control-free) behavior.   It is assumed to be irreducible and aperiodic.   It follows that P_0 admits a unique invariant pmf, denoted π_0.  For any other transition matrix, with unique invariant pmf π, the
 Donsker-Varadhan rate function is denoted,

    K(P P_0) = ∑_x,x'π(x) P(x,x')   log(P(x,x') /P_0(x,x'))
  
under the usual convention that “0log(0) = 0”.  
It is called a “rate function” because it defines the relative entropy rate between two stationary Markov chains, see
<cit.>.












    



As in <cit.>, the rate function is used here 
to model the cost of deviation from the nominal transition matrix P_0. 
The two control objectives surveyed in the introduction will be specialized as follows, based on the utility function    → and a scaling parameter ζ≥ 0.
For  the finite-horizon optimal control problem,  


    ^*_T(x,ζ) = max∑_t=0^T_x[ ζ(X(t))   - (X(t),P(t)) ]  ,

where   
the expectation is conditional on X(0)=x, and 

    (x, P) = D(P(x,) P_0(x,))∑_x'  P(x,x') log(P(x,x') /P_0(x,x'))

for any x∈  and transition matrix P.

The average reward optimization problem is analogous:

    η^*(ζ) =max( lim inf_T→∞1/T∑_t=0^T-1_x  [ ζ(X(t))   - (X(t),P(t))  ]   ) .

In each case, the  maximum is over all transition matrices {P(t)}. 
The average reward optimization problem can be cast as
the solution to the convex optimization problem,

    η^*_ζ =max_π,P{ζπ() -  K(P P_0)  :  π   P  = π}

where  the maximum is over all transition matrices.    

In this context, the one-step reward   appearing in (<ref>, <ref>) is a function of pairs (x,P):

    (x,P)  ζ(x) -   (x, P)

for any x∈  and transition matrix P.

There is practical value to considering a parameterized family of reward functions.  For one, it is useful to understand the sensitivity of the control solution to the relative weight given to utility    and the penalty on control action.   This is well understood in classical linear control theory – consider for example the celebrated symmetric root locus in linear optimal control <cit.>.















  
*Nature & nurture   

Exogenous randomness from nature imposes additional constraints in the optimal control problem  (<ref>) or (<ref>).   

It is assumed that the state space is the cartesian product of two finite sets:  = ×, and the state is similarly expressed X(t)=(X_u(t), X_n(t)).   At a given time t it is assumed that X_n(t+1) is conditionally independent of the input at time t,  given the value of X(t).   This is formalized by the  following conditional-independence assumption: 

    P(x,x') = R(x, x_u') Q_0(x,x_n') ,   x=(x_u,x_n) ∈,  x_u'∈, x_n'∈
 
The matrix R defines the randomized decision rule for X_u(t+1) given X(t).  
The matrix Q_0 is fixed and models the distribution of   X_n(t+1) given X(t)=x,
 and each are subject to the pmf constraint:
∑_x_u' R(x, x_u')=∑_x_n'  Q_0(x,x_n') =1
for each x. 

Subject to the constraint (<ref>), the two optimal control problems (<ref>, <ref>) are transformed to  the final forms considered in this paper:

    ^*_T(x,ζ)    =   max∑_t=0^T  _x[  (X(t),R(t))   ]
    
    η^*(ζ)    =   max{lim inf_T→∞1/T∑_t=0^T-1_x[  (X(t),R(t)) ]  }

where in each case the maximum is over sequences of randomized decision rules {R(0),…,R(T)},     

    (x,R)     ζ(x) -  (x,R) 
     
     and  (x,R)   ∑_x'  P(x,x') log(P(x,x') /P_0(x,x'))  
     =   ∑_x'_u  R(x,x'_u) log(R(x,x'_u) /R_0(x,x_u'))





 §.§ Notation




For any transition matrix P, an invariant pmf  is interpreted as a row vector, so that invariance can be expressed π P=π.   Any function f→ is interpreted as a d-dimensional column vector,  and we use the standard notation Pf  (x) =   ∑_x'P(x,x')f(x'),  x∈.
The fundamental matrix is   the inverse, 

    Z = [I - P + 1⊗π]^-1

where 1⊗π is a matrix in which each row is identical, and equal to π.  
If P is irreducible and aperiodic, then it can be expressed as the power series  
Z = ∑_n=0^∞  [P - 1⊗π]^n,




with [P - 1⊗π]^0  I  (the d× d identity matrix),
and [P - 1⊗π]^n  = P^n - 1⊗π  for n≥ 1.  

Any function g×→ is regarded as an unnormalized log-likelihood ratio:  Denote for x,x'∈,

    P_g(x,x')  P_0(x,x')exp(  g( x' | x )  -  g(x)    ),

in which g( x' | x ) is the value of g at (x,x')∈×, and g(x) is the normalization constant,

    g(x)     
    log( ∑_x' P_0(x,x')exp(  g(x'|  x)     ) )

The rate function can be expressed in terms of  its  invariant pmf π_g, the bivariate pmf  Π_g(x,x') = π_g(x) P_g(x,x'), and the log moment generating function (<ref>):

    K(P_g P_0)    = ∑_x,x'Π_g(x,x')    [  g( x' | x )  -  g(x) ]
    
       = ∑_x,x'Π_g(x,x')      g( x' | x )  -     ∑_x π_g(x) g(x)
  

The unusual notation is introduced because  g( x' | x ) will take the form of a conditional expectation in all of the results that follow:  given any function h→ we denote

    h(x'_u|  x)  =  ∑_x_n' Q_0(x,x_n') h(x_u',x_n')   .

In this case the transformation only transforms the dynamics of _u:

    P_h(x,x')  = R_h(x,x_u') Q_0(x,x_n')  ,  
    R_h(x,x'_u)  R_0(x,x'_u)exp( h( x'_u | x )  -  g(x)    )  .






  


 §.§ ODE for finite time horizon



Here an ODE is constructed to compute the   value functions {^*_τ (x,ζ)  : 1≤τ≤ T , ζ≥ 0}.    To aide exposition it is helpful to first look at the general problem:   Assume that the state space  is finite, the action space  is general,  and let {P_u(x,x')} denote the controlled transition matrix.   The  one-step reward on state-action pairs is of the form (x,u) = ζ(x) - c(x,u),  where c×→_+.   Assume that  c(x,u)≡ 0  for a unique value u=u_0.







For each 1≤τ≤ T denote, as in (<ref>),

    ^*_τ(x,ζ) = max∑_t=0^τ_x[ w(X(t), U(t)) ]
  
where the maximum is over all admissible  
inputs {U(t)=ϕ_t(X(0),…,X(t))}.
Each value function can be regarded as the maximum over functions   {ϕ_t}  (subject to measurability conditions and 
hard constraints on the input).  It is assumed that the maximum (<ref>) is finite for each (x ,ζ). 
 
The dynamic programming equation (principle of optimality) holds:  for τ≥ 1,

    ^*_τ(x,ζ) = max_u {ζ(x)  -  c(x,u)     + ∑_x' P_u(x,x') ^*_τ-1 (x')  }
  
Assume that a maximizer  ϕ^*_τ-1,ζ(x) exits for each τ,ζ, and x.

A crucial observation is that for each x, the value function appearing in (<ref>) is the maximum
of functions that are affine in ζ.  It follows that ^*_τ(x,ζ) is convex as a function of ζ, and hence absolutely continuous.  Consequently, the right derivative   H^*_τ(x,ζ) ^*_τ(x,ζ) exists everywhere.  A recursive equation follows from 
(<ref>):

    H^*_τ(x,ζ) = (x)       + ∑_x'_τ-1,ζ(x,x') H^*_τ-1 (x',ζ)
  
where _τ-1,ζ(x,x') = P_u^*(x,x') with u^* =  ϕ^*_τ-1,ζ(x).   

In matrix notation this becomes H^*_τ =  _τ-1,ζ, where
 _0,ζ =I, and
 for any 1≤τ≤ T,

    _τ-1,ζ = I +  _τ-1,ζ+_τ-1,ζ_τ-2,ζ
    			+ _τ-1,ζ_τ-2,ζ⋯_0,ζ

 This is similar to a truncation of the  power series representation of the fundamental matrix (<ref>).

Denote ^*_ζ(x) ={^*_k(x,ζ)  : 0≤ k≤ T}, regarded as a vector in ^||×(T+1), parameterized by the non-negative constant ζ.  The following result follows from the preceding arguments: 

  
The family of functions {^*_ζ} solves the ODE
^*_ζ =( ^*_ζ ),  ζ≥ 0,  
with boundary condition ^*_0 = 0.   The vector field can be described in block-form as follows, with T+1 blocks:

    ^*_k(,ζ)  =_k( ^*_ζ ) ,   0≤ k≤ T  .

The identity _0() = holds for any .
For k≥ 1, the right hand side depends on  its argument only through the associated policy:  for any sequence of functions =(_0,…,_T),

    _k()     =  Z_k-1
    where    
      Z_k-1   = 
     I +  P_k-1+P_k-1P_k-2
    			+ P_k-1P_k-2⋯ P_0 	 
    			
    
    			P_i(x,x')    = 	P_ϕ_i(x) (x,x') ,   all i, x, x',
    ϕ_i(x)     =  _u   {  -  c(x,u)     + ∑_x' P_u(x,x') _i (x')  } ,      1≤ i, k≤ T.


The theorem provides valuable computational tools for models of moderate 
cardinality and moderate time-horizon.  
Two questions remain:


  * What is ϕ_i for the problem under study in this paper?

  * Can a tractable ODE be constructed in infinite-horizon optimal control problems?

The answer to the second question is the focus of Section <ref>.  

The answer to (i) is contained in the following.   For any function →, denote

    R_(x,) = _R  {(x,R) +∑_x' P(x,x') (x') } ,   x∈  ,

subject to the constraint that P depends on R via (<ref>), and with  defined in (<ref>).
 


  
For any function  the maximizer R_ is unique and can be expressed

    R_(x,x_u') = R_0(x,x_u')  exp(  (x'_u|  x)  -   Λ(x)    )

where  (x'_u|  x)  =  ∑_x_n' Q_0(x,x_n') (x_u',x_n')  for each x∈, x_u'∈,
and Λ(x) is a normalizing constant, defined so that R_(x,) is a pmf for each x.






















 Given the form of the reward  and the constraint on P, 
the optimization problem of interest here can be written, for each x, as
 
    R_ (x,) = _μ{μ(  ) - D(μμ_0)  }

where the variable μ() represents R(x,),   μ_0=R_0(x,), and

    μ(  ) =  ∑_x' = (x'_u,x'_n) R(x,x'_u)Q_0(x,x'_n) (x'_u,x'_n) = ∑_x'_uμ(x'_u) (x'_u|  x)

The proposition is a consequence of this combined with   Theorem 3.1.2  of  <cit.>  (i.e., convex duality between relative entropy and the log moment generating function).

 

It follows from the proposition that the vector field is smooth in a neighborhood of the optimal solution {^*_ζ : ζ≥ 0}.   These results are central to the average-reward case considered next.  

 


 §.§ Average reward formulation



We consider now the case of average reward (<ref>), subject to the structural constraint (<ref>).    The associated average reward optimization equation (AROE) is expressed as follows: 

    max_R{(x,R)
    
    +∑_x' P(x,x') h^*_ζ(x') } = h^*_ζ(x) + η^*(ζ)

 In which   η^*(ζ) is the optimal average reward, and h^*_ζ is the relative value function.  The maximizer defines a transition matrix:

    _ζ =_P {ζπ() -  K(P P_0)  : π   P  = π}

Recall that the relative value function is not unique, since  a new solution is obtained by adding a non-zero constant;  the normalization  h^*_ζ()=0 is imposed, where ∈ is a fixed state.



   The proof of t:IPD (i) is a consequence of t:phiR.  The second result is obtained on combining   Lemmas B.2–B.4 of   <cit.>.


There exist optimizers {_ζ, _ζ : ζ∈},  and solutions to
the AROE    {h^*_ζ,η^*(ζ): ζ∈} with the following properties: 




  * The optimizer _ζ can be obtained from the relative value function h^*_ζ as follows:  

    _ζ(x,x')
     P_0(x,x')exp(  h_ζ(x'_u|  x)  -  h_ζ(x)    )

where  for x∈, x_u'∈,

    h_ζ(x'_u|  x)  =  ∑_x_n' Q_0(x,x_n') h^*_ζ(x_u',x_n'),

and h_ζ(x) is the normalizing constant (<ref>) with   h=h_ζ. 
 


  * {_ζ, _ζ, h^*_ζ,η^*(ζ): ζ∈} are continuously differentiable in the parameter ζ.


 

Representations for the derivatives in t:IPD (ii),  in particular the derivative of h_ζ^* with respect to ζ,  lead to a representation for the ODE   used to compute the    transition matrices {_ζ}. 

 
 

 
It is convenient to generalize the problem slightly here:  let {h_ζ^∘ : ζ∈} denote a family of functions on , continuously differentiable in the parameter ζ.  They are not necessarily relative value functions, but we maintain the structure established in  t:IPD for the  family of transition matrices.  Denote,

    h_ζ (x'_u|  x)   =  ∑_x_n' Q_0(x,x_n') h^∘_ζ(x_u',x_n'),  
       
    x∈, x_u'∈

and then define as in (<ref>),

    P_ζ(x,x')
     P_0(x,x')exp(  h_ζ(x'_u|  x )  -  h_ζ(x)    )
 
The function   h_ζ→ is a normalizing constant, exactly as in (<ref>):

    h_ζ^∘(x)     
    log( ∑_x' P_0(x,x')exp(  h_ζ(x'_u|  x)     ) )


 




We begin with a general method to construct  a family of functions {h_ζ^∘ : ζ∈} based on an ODE.     The ODE is expressed,  

    h_ζ^∘  =  (h_ζ^∘) ,    ζ∈,

with boundary condition h_0^∘≡ 0.  A particular instance of the method will result in h_ζ^∘ = h_ζ^* for each ζ.
 Assumed given is a mapping  from transition matrices to functions on .   Following this, the vector field  is obtained through the following two steps:  For a function h→,  


  * Define a new transition matrix via (<ref>),

    P_h(x,x')  P_0(x,x')exp(  h(x_u' | x)  -  h(x)    ),   x,x'∈,

in which h(x_u' | x) = ∑_x_n' Q_0(x,x_n') h(x_u',x_n'),  and h(x) is a normalizing constant.   



  * Compute = (P_h), and define (h) =.
It is assumed   that the functional  is constructed so that  ()=0 for any h.



   




 
We now specify the functional  , whose domain consists of  transition matrices that are irreducible and aperiodic.    For  any transition matrix  P in this domain, the fundamental matrix Z is obtained using (<ref>), and then =(P) is defined as

    (x) = ∑_x' [ Z(x,x')-Z(,x') ]  (x'),     x∈

The function    is a solution to Poisson's equation,

    P  = - + ,     where π()  ∑_xπ(x) (x).






  


Consider the ODE 
(<ref>)
with boundary condition h_0^∘≡ 0,
and with  =(P) defined using (<ref>) for each transition matrix P
that is irreducible and aperiodic.
The solution to this ODE exists, and the resulting functions { h^∘_ζ : ζ∈}    
coincide with the relative value functions {h^*_ζ: ζ∈}.  Consequently,   _ζ = P_h_ζ for each ζ.



 

The proof requires validation of the representation H^*_ζ = (_ζ) for each ζ, where h^*_ζ is the relative value function, 
_ζ is defined in (<ref>), and

    H^*_ζ = h^*_ζ
 
Substituting the maximizer _ζ in the form (<ref>) into the AROE gives the fixed point equation ζ   + h^*_ζ = h^*_ζ + η^*(ζ).





Differentiating each side   then gives,

    +  _ζ H^*_ζ= H^*_ζ + η^*(ζ) .

This is Poisson's equation, and it follows that _ζ() =  η^*(ζ).  Moreover,  since h^*_ζ()=0 for every ζ,  we must have H^*_ζ()=0 as well.  Since the solution to Poisson's equation with this normalization is unique, we conclude that  (<ref>) holds, and hence H^*_ζ = (_ζ) as claimed.  



l^∙

X^L
ω






§ EXAMPLE
 

 
We consider a variant of the example of  <cit.> in which a UAV 
(unmanned aerial vehicle) 
 needs to reach a target subject to energy costs, and subject to disturbances from wind.   The location of the UAV at time t is denoted L_t, and evolves according to the controlled linear dynamics:

    L_t+1 = L_t + W_t + U_t

where ={U_t} is the control sequence, and  ={W_t}  models the impact of the wind. There are d_L locations across a two-dimensional grid. 

Wind is location-dependent: It is assumed that the  wind profile over the region is determined by a stochastic process ={ N_t} and a function  such that for each t,

    W_t=(L_t, N_t ).

The process   is assumed to be Markovian with finite state space {1,…,d_N}, and state transition matrix denoted Q_0.  This is the
 nature component of the MDP model,  with state process X_t=(L_t,N_t),  t≥ 0.

A nominal model is described by a randomized policy in which U_t=0 with high probability.   The specific form used in the experiments 
was constructed as follows.   On denoting  L^+_t = L_t + W_t, a transition matrix R_0^L is constructed with the interpretation

    R_0^L(l^+, l')  = { L_t+1 =  l' | L^+_t  = l^+ } ,   t≥ 0  .

The nominal randomized strategy is the d_L× d_L matrix,

    R_0(x, u)  = { U_t = u| X_t  = x } = R_0^L(l+(l,n) , l+(l,n)+u ) ,        x=(l,n)  .

The overall transition matrix is the product:

    P_0(x,x') = R_0^L(l+w(n,l), l') Q_0(n,n'),     x=(n,l), x'=(n',l')  .



  





The goal of the control problem is to reach a target location  and remain there.  To ensure that the set { (,n) :  1≤ n≤ d_N} is absorbing, a separate rule is imposed on R_0 for these states:   W_t + U_t  =0 if L_t =.

The reward function  is taken to be a scaled negative cost:  = -  c,  where c→_+,
with c()=0 and c(l)>0 for l≠.    The optimal steady-state mean is zero in this model, and the relative value function is the negative of the cost to go:

    -h^*(x) = 
    J^*(x) min_x[∑_t=0^τ_∙{ζ c(L_t) +(X_t,R(t) ) }]

where τ_∙  (unknown a-priori) is the first hitting time to .   An example is illustrated in f:CostToGo15x15n2n4 — the details are provided in the following.


 

  

 


  
Details of the numerical experiment

The  set of locations  is taken to be a rectangular grid of the form = { (i,j) : 1≤ i≤ d_a,  1≤ j ≤ d_o}, in which   d_a, d_o≥ 2 and d_L=d_a× d_o  (the subscripts are  meant to suggest latitude and longitude).  The function c appearing in (<ref>) was taken to be the indicator function,  c(l)={l≠},  with =(d_a,d_o).


The values  d_a=d_o=15, and d_N=5 are fixed throughout.
The size of the state space is thus d_a× d_o× d_N= 1,125,  and the action space is a subset of the simplex in ^1,125.  




The transition matrix for nominal control was taken   of the following form:

    R_0^L(l, l')   = κ(l)   exp{-1/2σ^2_u l'-l^2 }  ,   l,l'∈ ,

where κ(l) >0 is chosen so that R_0^L(l, ) is a pmf on  for each l∈.   The value σ^2_u = 1/2 was used in the numerical results that follow.











The Markov chain  was taken to be a skip-free symmetric random walk on the integers { 1,…,d_N}.   For a fixed δ_n∈(0,1) the probability of transition is Q_0(n,n+ 1) =Q_0(n,n-1) =δ_n,
where addition is modulo d_N,  
and Q_0(n,n) =1-δ_n for any n.
Recall that this means

    {N_t+1 = n+1 | N_t = n}
    =
    {N_t+1 = n-1 | N_t = n} =  δ_n.

The value δ_n=0.05 was chosen in these experiments.



Recall that →^2 is used to defined the wind process .  
For each value of n,  the function (,n) can be interpreted as a vector field on .   For each n,  a slowly varying continuous function was constructed on the two-dimensional rectangle [1,d_a]×[1,d_o]. The function (,n) was  taken to be its quantization to the lattice .
 The values were restricted to the set of pairs {(i,j): |i|≤ 1, |j|≤ 1}.  

 

The family of optimal policies was obtained using the ODE method,   and the solution for three values of ζ is illustrated in  f:Quiver15by15n2n4.pdf.    Each of the arrows shown is proportional to  the conditional expectation:

    v(l,n) [L_t+1 - L_t | L_t = l  N_t=n]

in which l∈ is the position on the grid.  The figure shows only the values n=2 and n=4 (the most interesting to view because of obvious spatial variability). 











  











 
If the position l=(l_a,l_o)  is far from the boundary of , say,  min(l_a,l_o)≥ 4 and  min(d_a-l_a,d_o-l_o)≥ 4, then

    [U_t | L_t = l  N_t=n]≈ 0 and 
    v(l,n) ≈(l,n),    ζ=0






For the case ζ=1 the vector field is transformed so that vectors near the target state point in this direction; for ζ=2 this behavior is  more apparent.   For states far from the target the control effort seems to be lower – most likely the optimal policy waits for more favorable weather that will push the UAV in the North-East direction.













The eigenvalues of _ζ are shown in f:EigenvalueComparison15x15
for ζ=0,1,2.    Most of the eigenvalues are driven near zero for ζ=2. 
Those three that are independent of ζ are  the three eigenvalues of Q_0,
{0.9095, 
    0.9655, 1   }.




While the vector field and eigenvalues change significantly when ζ is doubled from 1 to 2, the cost to go  J^* defined in (<ref>) grows relatively slowly with ζ.  Shown on the right hand side of f:CostToGo15x15n2n4 are comparisons for these two values of ζ.  One plot with n=2 and the other n=4.
The plot on the far right shows J_ζ^*(l,ζ) for 0≤ζ≤ 1 and l=(1,1)  (the location farthest from ).   

These plots are easily obtained because of the nature of the algorithm:   the optimal policy and value function are generated for any range of ζ of interest.



   















§ CONCLUSIONS
 



The ODE approach for solving MDPs has  simple structure for the class of models considered in this paper.  We are currently looking at approaches to approximate dynamic programming as has been successful in the unconstrained model  <cit.>.


It is likely that the ODE has special structure for other classes of MDPs, such as the  “rational inattention” framework of <cit.>.  The computational efficiency of this approach  will depend in part on numerical properties of the ODE, such as its sensitivity for complex models.   Applications to distributed control were the original motivation for this work, with particular attention to “demand dispatch” <cit.>.  It is believed that this paper will offer  new computational tools in this ongoing research.

 


abbrv




 

