Matlab Procedures for Markov Decision Processes

Summary: We first summarize certain essentials in the analysis of homogeneous Markov chains with costs and rewards associated with states, or with transitions between states. Then we consider three cases: a. Gain associated with a state; b. One-step transition gains; and c. Gains associated with a demand under certain reasonable conditions. Matlab implementations of the results of analysis provide machine solutions to a variety of problems.

In order to provide the background for Matlab procedures for Markov decision processes, we first summarize certain essentials in the analysis of homogeneous Markov chains with costs and rewards associated with states, or with transitions between states. References are to Pfeiffer: PROBABILITY FOR APPLICATIONS.  

1. Some cost and reward patterns 
   Consider a finite, ergodic (i.e., recurrent, positive, aperiodic, irreducible) homogeneous Markov chain X_N, with state space E={1,2,⋯,M{. To facilitate use of matrix computation programs, we number states from one, except in certain examples in which state zero has a natural meaning. Associated with this chain is a cost or reward structure belonging to one of the three general classes described below. The one-step transition probability matrix is P=[p(i,j)], where p(i,j)=P(X_n+1=j|X_n=i). The distribution for X_n is represented by the row matrix π(n)=[p₁(n),p₂(n),⋯,p_M(n)], where p_i(n)=P(X_n=i). The long-run distribution is represented by the row matrix π=[π₁,π₂,⋯,π_M].
   1. Type 1. Gain associated with a state. A reward or gain in the amount g(i) is realized in the next period if the current state is i. The gain sequence{G_n:1≤n{ of random variables G_n+1=g(X_n) is the sequence of gains realized as the chain X_N evolves. We represent the gain function g by the column matrix g=[g(1)g(2)⋯g(M)]^T.
   2. Type 2. One-step transition gains. A reward or gain in the amount g(i,j)=g_ij is realized in period n+1 if the system goes from state i in period n to state j in period n+1. The corresponding one-step transition probability is p(i,j)=p_ij. The gain matrix is g=[g(i,j)]. The gain sequence{G_n:1≤n{ of random variablesG_n+1=g(X_n,X_n+1) is the sequence of gains experienced as the chain X_N evolves.
   3. Type 3. Gains associated with a demand. In this case, the gain random variables are of the form
      G_n+1=g(X_n,D_n+1)

      (1)
      If n represents the present, the random vector U_n=(X₀,X₁,⋯,X_n) represents the “past” at n of the chain X_N. We suppose {D_n:1≤n{ is iid and each {D_n+1,U_n{ is independent. Again, the gain sequence{G_n:1≤n{ represents the gains realized as the process evolves. Standard results on Markov chains show that in each case the sequence G_N={G_n:1≤n{ is Markov.
   A recurrence relation. We are interested in the expected gains in future stages, given the present state of the process. For any n>0 and any m>1, the gain G_n^(m) in the mperiods following period n is given by
   G_n^(m)=G_n+1+G_n+2+⋯+G_n+m

   (2)
   If the process is in state i, the expected gain in the next period q_i is
   q_i=v_i⁽¹⁾=E[G_n+1|X_n=i]

   (3)
   and the expected gain in the next m periods is
   v_i^(m)=E[G_n^(m)|X_n=i]

   (4)
   We utilize a recursion relation that allows us to begin with the transition matrix P and the next-period expected gain matrixq=[q₁q₂⋯q_m]^T and calculate the v_i^(m) for any m>1. Although the analysis is somewhat different for the various gain structures, the result exhibits a common pattern. In each case
   v_i^(m)=E[G_n^(m)|X_n=i]=E[G_n+1|X_n=i]+

   m−1

   ∑

   k=1

    E[G_n+k+1|X_n=i]

   (5)
   =q_i+

   m−1

   ∑

   k=1

    ∑E[G_n+k+1|X_n+1=j]p(i,j)

   (6)
   =q_i+∑E[

   m−1

   ∑

   k=1

    G_n+k|X_n=j]p(i,j)

   (7)
   We thus have the fundamental recursion relation
   (*)v_i^(m)=q_i+∑v_j^(m−1)p(i,j)

   (8)
   The recursion relation (*) shows that the transition matrix P=[p(i,j)] and the vector of next-period expected gains, which we represent by the column matrixq=[q₁,q₂,⋯,q_M]^T, determine the v_i^(m). The difference between the cases is the manner in which the q_i are obtained.
   • Type 1: . q_i=E[g(X_n)|X_n=i]=g(i)
   • Type 2: . q_i=E[g(X_n,X_n+1)|X_n=i]=E[g(i,X_n+1)|X_n=i]=∑_j∈Eg(i,j)p(i,j)
   • Type 3: . q_i=E[g(X_n,D_n+1)|X_n=i]=E[g(i,D_n+1)]=∑_kg(i,k)P(D=k)
   • Matrix form: . For computational purposes, we utilize these relations in matrix form. If
     v(n)=[v₁⁽ⁿ⁾v₂⁽ⁿ⁾⋯v_M⁽ⁿ⁾]^Tandq=[q₁q₂⋯q_M]^T

     (9)
     then
     (*)v(m+1)=q+Pv(m)forallm>0,withv(0)=0

     (10)
   Examination of the expressions for q_i, above, shows that the next-period expected gain matrix q takes the following forms. In each case, g is the gain matrix. In case c, p_D is the column matrix whose elements are P(D=k).
   • Type 1: q=g
   • Type 2: q=thediagonalofPg^T
   • Type 3: q=gp_D
2. Some long-run averages 
   Consider the average expected gain for m periods
   E[

   1

   m

    G_n^(m)]=

   1

   m

    

   m

   ∑

   k=1

    E[G_n+k]=

   1

   m

    

   m

   ∑

   k=1

    E{E[G_n+k|X_n−1]{

   (11)
   Use of the Markov property and the fact that for an ergodic chain

   1

   m

    

   m

   ∑

   k=1

    p^k(i,j)→π_jasm→∞

   (12)
   yields the result that,
   limE[

   1

   m

    G_n^(m)]=∑P(X_n−1=i)∑q_jπ_j=∑q_jπ_j=g

   (13)
   and for any state i,
   limE[

   1

   m

    G_n^(m)|X_n=i]=lim

   1

   m

    v_i^(m)=g

   (14)
   The expression for g may be put into matrix form. If the long-run probability distribution is represented by the row matrix π=[π₁π₂⋯π_M] and q=[q₁q₂⋯;q_M]^T, then
   g=πq

   (15)
3. Discounting and potentials 
   Suppose in a given time interval the value of money increases by a fraction a. This may represent the potential earning if the money were invested. One dollar now is worth 1+a dollars at the end of one period. It is worth (1+a)ⁿ dollars after n periods. Set α=1/(1+a), so that 0<α≤1. It takes α dollars to be worth one dollar after one period. It takes αⁿdollars at present to be worth one dollar n periods later. Thus the present worth or discounted value of one dollar n periods in the future is αⁿ dollars. This is the amount one would need to invest presently, at interest rate a per unit of time, to have one dollar n periods later. If f is any function defined on the state space E, we designate by f the column matrix[f(1)f(2)⋯f(M)]^T. We make an exception to this convention in the case of the distributions of the state probabilities π(n)=[p₁(n)p₂(n)⋯p_M(n)], their generating function, and the long-run probabilities π=[π₁π₂⋯π_M], which we represent as row matrices. It should be clear that much of the following development extends immediately to infinite state space E=N={0,1,2,⋯{. We assume one of the gain structures introduced in Sec 1 and the corresponding gain sequence {G_n:1≤n{. The value of the random variable G_n+1 is the gain or reward realized during the period n+1. We let each transition time be at the end of one period of time (say a month or a quarter). If n corresponds to the present, then G_n+k is the gain in period n+k. If we do not discount for the period immediately after the present, then α^k−1G_n+k is the present value of that gain. Hence

   ∞

   ∑

   k=1

    α^k−1G_n+kisthetotaldiscountedfuturegain

   (16)
   Definition. The α-potential of the gain sequence {G_n:1≤n{ is the function φ defined by
   φ(i)=E[

   ∞

   ∑

   n=0

    αⁿG_n+1|X₀=i]∀i∈E

   (17)
   Remark. φ(i) is the expected total discounted gain, given the system starts in state i. We next define a concept which is related to the α-potential in a useful way. Definition. Theα-potential matrixR^α for the process X_N is the matrix
   R^α=

   ∞

   ∑

   n=0

    αⁿPⁿwithP⁰=I

   (18)

   THEOREM 1: 3.1

   Let X_N be an ergodic, homogeneous Markov chain with state space E and gain sequence {G_n:1≤n{. Let q=[q₁q₁⋯q_M]^T where q_i=E[G_n+1|X_n=i]fori∈E. Forα∈(0,1), let φ be the α-potential for the gain sequence. That is,

   φ(i)=E[

   ∞

   ∑

   n=0

    ,αⁿ,G_n+1,|,X₀,=,i]∀i∈E

   (19)

   Then, if R^α is the α-potential matrix for X_N, we have

   φ=R^αq

   (20)

   — □

   THEOREM 2: 3.2

   Consider an ergodic, homogeneous Markov chain X_N with gain sequence {G_n:1≤n{ and next-period expected gain matrix q. For α∈(0,1), the α-potential φ and the α-potential matrix R^α satisfy

   [I−αP]R^α=Iand[I−αP]φ=q

   (21)

   If the state space E is finite, then

   R^α=[I−αP]⁻¹andφ=[I−αP]⁻¹q=R^αq

   (22)

   — □

   EXAMPLE 1: A numerical example

   Suppose the transition matrix P and the next-period expected gain matrix q are

   P=[

   0.2	0.3	0.5
   0.4	0.1	0.5
   0.6	0.3	0.1

   ]andq=[

   2

   5

   3

   ]

   (23)

   For α=0.9, we have

   R^0.9=(I−0.9P)⁻¹=[

   4.4030	2.2881	3.3088
   3.5556	3.1356	3.3088
   3.6677	2.2881	4.0441

   ]andφ=R^0.9q=[

   30.17

   32.72

   30.91

   ]

   (24)

   — □
   The next example is of class c, but the demand random variable is Poisson, which does not have finite range. The simple pattern for calculating the q_i must be modified.

   EXAMPLE 2: Costs associated with inventory under an (m,M) order policy.

   If k units are ordered, the cost in dollars is

   c(k)={

   0	k=0
   10+25k	0<k≤M

   (25)

   For each unit of unsatisfied demand, there is a penalty of 50 dollars. We use the (m,M) inventory policy described in Ex 23.1.3. X₀=M, and X_n is the stock at the end of period n, before restocking. Demand in period n is D_n. The cost for restocking at the end of period n+1 is

   g(X_n,D_n+1)={

   10+25(M−Xn)+50max{Dn+1−M,0{	0≤Xn<m
   50max{Dn+1−Xn,0{	m≤Xn≤M

   (26)

   For m=1,M=3, we have

   g(0,D_n+1)=85+50I_[3,∞)(D_n+1)(D_n+1−3)

   (27)

   g(i,D_n+1)=50I_[3,∞)(D_n+1)(D_n+1−i)1≤i≤3

   (28)

   We take m=1,M=3 and assume D is Poisson (1). From Ex 23.1.3 in PA, we obtain

   P=[

   0.0803	0.1839	0.3679	0.3679
   0.6321	0.3679	0	0
   0.2642	0.3679	0.3679	0
   0.0803	0.1839	0.3679	0.3679

   ]

   (29)

   The largest eigenvalue |λ|≈0.26, so n≥10 should be sufficient for convergence. We use n=16.

   Taking any row of P¹⁶, we approximate the long-run distribution by

   π=[0.28580.28470.26320.1663]

   (30)

   We thus have

   q₀=85+50E[I_[3,∞)(D_n+1)(D_n+1−3)]=85+50

   ∞

   ∑

   k=4

    (k−3)p_k(thetermfork=3iszero)

   (31)

   For the Poisson distribution

   ∞

   ∑

   k=n

    kp_k=λ

   ∞

   ∑

   k=n−1

    p_k. 

   Hence q₀=85+50[

   ∞

   ∑

   k=3

    p_k−3

   ∞

   ∑

   k=4

    p_k]≈85+50[0.0803−3×0.0190]=86.1668 q₁=50

   ∞

   ∑

   k=3

    (k−1)p_k=50[

   ∞

   ∑

   k=2

    p_k−

   ∞

   ∑

   k=3

   p_k]=50p₂=9.1970 q₂=50

   ∞

   ∑

   k=3

    (k−2)p_k=50[0.2642−2×0.0803]=5.1809 q₃=q₀−85=1.1668 so that

   q=[

   86.1668

   9.1970

   5.1809

   1.1668

   ]

   (32)

   Then, for α=0.8 we have

   R^0.8=(I−0.8P)⁻¹=[,

   1.9608	1.0626	1.1589	0.8178
   1.4051	2.1785	0.8304	0.5860
   1.1733	1.2269	2.1105	0.4893
   0.9608	1.0626	1.1589	1.8178

   ]andφ=R^0.8q=[

   185.68

   146.09

   123.89

   100.68

   ]

   (33)

   Recall that we are dealing with expected discounted costs. Since we usually consider starting with a full stock, the amount φ(M)=φ(3)=100.68 is the quantity of interest. Note that this is smaller than other values of φ(j), which correspond to smaller beginning stock levels. We expect greater restocking costs in such cases.
4. Evolution of chains with costs and rewards
   1. No discounting A generating function analysis of
      (*)v(m+1)=q+Pv(m)forallm>0,withv(0)=0

      (34)
      shows
      v(n)=ng1+v+transients,wherev=B₀q

      (35)
      Here g=πq, is a column matrix of ones, P₀=P^∞, and B₀ is a matrix which analysis also shows we may approximate by
      B₀=B(1)≈I+P+P²+⋯+Pⁿ−(n+1)P₀

      (36)
      The largest |λ_i|<1 gives an indication of how many terms are needed.

      EXAMPLE 3: The inventory problem (continued)

      We consider again the inventory problem. We have

      P=[

      0.0803	0.1839	0.3679	0.3679
      0.6321	0.3679	0	0
      0.2642	0.3679	0.3679	0
      0.0803	0.1839	0.3679	0.3679

      ]andq=[,

      86.1668

      9.1970

      5.1819

      1.1668

      ,]

      (37)

      The eigenvalues are 1,0.0920+i0.2434,0.0920−0.2434 and 0. Thus, the decay of the transients is controlled by|λ^*|=(0.0920²+0.2434²)^1/2=0.2602. Since |λ^*|⁴≈0.0046, the transients die out rather rapidly. We obtain the following approximations

      P₀≈P¹⁶≈[

      0.2852	0.2847	0.2632	0.1663
      0.2852	0.2847	0.2632	0.1663
      0.2852	0.2847	0.2632	0.1663
      0.2852	0.2847	0.2632	0.1663

      ]sothatπ≈[0.28520.28470.26320.1663]

      (38)

      The approximate value of B₀ is found to be

      B₀≈I+P+P²+P³+P⁴−5P₀≈[

      0.4834	−0.3766	−0.1242	0.0174
      0.0299	0.7537	−0.5404	−0.2432
      −0.2307	−0.1684	0.7980	−0.3989
      −0.5166	−0.3766	−0.1242	1.0174

      ,]

      (39)

      The value g=πq≈28.80 is found in the earlier treatment. From the values above, we have

      v=B₀q≈[

      37.6

      6.4

      −17.8

      −47.4

      ]sothatv(n)≈[

      28.80n+37.6

      28.80n+6.4

      28.80n−17.8

      28.80n−47.4

      ]+transients

      (40)

      The average gain per period is clearly g≈28.80. This soon dominates the constant terms represented by the entries in v.
   2. With discounting Let the discount factor be α∈(0,1). If n represents the present period, then G_n+1= the reward in the first succeeding period G_n+k= the reward in thekth succeding period. If we do not discount the first period, then
      G_n^(m)=G_n+1+αG_n+2+α²G_n+3+⋯+α^m−1G_n+m=G_n+1+αG_n+1^(m−1)

      (41)
      Thus
      v_i^(m)=E[G_n^(m)|X_n=i]=q_i+α

      M

      ∑

      j=1

       p(i,j)v_j^(m−1)

      (42)
      In matrix form, this is
      v(n)=q+αPv(n−1)

      (43)
      A generating function analysis shows that
      v_i⁽ⁿ⁾=v_i+transients1≤i≤M

      (44)
      Hence, the steady state part of
      v_i^(m)=q_i+α

      M

      ∑

      j=1

       p(i,j)v_j^(m−1)isv_i=q_i+α

      M

      ∑

      j=1

       p(i,j)v_j1≤i≤M

      (45)
      In matrix form,
      v=q+αPvwhichyieldsv=[I−αP]⁻¹q

      (46)
      Since the q_i=E[G_n+1|X_n=i] are known, we can solve for the v_i. Also, since v_i^(m) is the present value of expected gain m steps into the future, the v_i represent the present value for an unlimited future, given that the process starts in state i.
5. Stationary decision policies 
   We suppose throughout this section that X_N is ergodic, with finite state space E={1,2,⋯,M{. For each state i∈E, there is a class A_i of possible actions which may be taken when the process is in state i. A decision policy is a sequence of decision functions d₁,d₂⋯ such that
   Theactionatstagenisd_n(X₀,X₁,...,X_n).

   (47)
   The action selected is in A_i whenever X_n=i. We consider a special class of decision policies. Stationary decision policy
   d_n(X₀,X₁,⋯,X_n)=d(X_n)withdinvariantwithn

   (48)
   The possible actions depend upon the current state. That is d(X_n)∈A_i whenever X_n=i. Analysis shows the Markov character of the process is preserved under stationary policies. Also, if E is finite and every stationary policy produces an irreducible chain, then an optimal stationary policy is optimal. We suppose each policy yields an ergodic chain. Since the transition probabilities from any state depend in part on the action taken when in that state, the long-run probabilities will depend upon the policy. In the no-discounting case, we want to determine a stationary policy that maximizes the gain g=πq for the corresponding long-run probabilities. We say “a stationary policy,” since we do not rule out the possibility there may be more than one policy which maximizes the gain. In the discounting case, the goal is to maximize the steady state part v_i in the expressions
   v_i⁽ⁿ⁾=v_i+transients1≤i≤M

   (49)
   In simple cases, with small state space and a small number of possible actions in each state, it may be feasible to obtain the gain or present values for each permissible policy and then select the optimum by direct comparison. However, for most cases, this is not an efficient approach. In the next two sections, we consider iterative procedures for step-by-step optimization which usually greatly reduces the computational burden.
6. Policy iteration method– no discounting 
   We develop an iterative procedure for finding an optimal stationary policy. The goal is to determine a stationary policy that maximizes the long run expected gain per period g. In the next section, we extend this procedure to the case where discounting is in effect. We assume that each policy yields an ergodic chain. Suppose we have established that under a given policy (1) v_i⁽ⁿ⁾=q_i+∑p(i,j)v_j⁽ⁿ⁻¹⁾, where q_i=E[G_n+1|X_n=i] In this case, the analysis in Sec 4 shows that for large n, (2) v_i⁽ⁿ⁾≈ng+v_i, where g=∑π_iq_i We note that v_i and g depend upon the entire policy, while q_i and p(i,j) depend on the individual actions a_i. Using (1) and (2), we obtain
   ng+v_i=q_i+∑p(i,j)[(n−1)g+v_j]=q_i+(n−1)g+∑p(i,j)v_j

   (50)
   From this we conclude (3) g+v_i=q_i+∑p(i,j)v_j for all i∈E Suppose a policy d has been used. That is, action d(i) is taken whenever the process is in state i. To simplify writing, we drop the indication of the action and simply write p(i,j) for p_ij(d(i)), etc. Associated with this policy, there is a gain g. We should like to determine whether or not this is the maximum possible gain, and if it is not, to find a policy which does give the maximum gain. The following two-phase procedure has been found to be effective. It is essentially the procedure developed originally by Ronald Howard, who pioneered the treatment of these problems. The new feature is the method of carrying out the value-determination step, utilizing the approximation method noted above.
   1. Value-determination procedure We calculate g=∑π_iq_i=πq. As in Sec 4, we calculate
      v=B₀q≈[I+P+P²+⋯+P^s−(s+1)P₀]qwhereP₀=limPⁿ

      (51)
   2. Policy-improvement procedure We suppose policy d has been used through period n. Then, by (3), above,
      g+v_i=q_i+∑p(i,j)v_j

      (52)
      We seek to improve policy d by selecting policy d^*, with d^*(i)=a_ik^*, to satisfy
      q_i^*+∑p_ij^*v_j=max{q_i(a_ik)+∑p_ij(a_ik)v_j:a_ik∈A_i{,1≤i≤M

      (53)
   Remarks
   • In the procedure for selecting d^*, we use the “old” v_j.
   • It has been established that g^*≥g, with equality iff g is optimal. Since there is a finite number of policies, the procedure must converge to an optimum stationary policy in a finite number of steps.
   We implement this two step procedure with Matlab. To demonstrate the procedure, we consider the following

   EXAMPLE 4: A numerical example

   A Markov decision process has three states: State space E={1,2,3{.

   Actions: State 1: A₁={1,2,3{ State 2: A₂={1,2{ State 3: A₃={1,2{

   Transition probabilities and rewards are:

   TABLE 1

   p1j(1): [1/3 1/3 1/3]	g1j(1): [1 3 4]
   p1j(2): [1/4 3/8 3/8]	g2j(2): [2 2 3]
   p1j(3): [1/3 1/3 1/3]	g3j(3): [2 2 3]
   p2j(1): [1/8 3/8 1/2]	g2j(1): [2 1 2]
   p2j(2): [1/2 1/4 1/4]	g2j(2): [1 4 4]
   p3j(1): [3/8 1/4 3/8]	g3j(1): [2 3 3]
   p3j(2): [1/8 1/4 5/8]	g3j(2): [3 2 2]

   Use the policy iteration method to determine the policy which gives the maximum gain g.

   A computational procedure utilizing Matlab We first put the data in an m-file. Since we have several cases, it is expedient to include the case number. This example belongs to type 2. Data in file md61.m

   type = 2
   PA = [1/3 1/3 1/3; 1/4 3/8 3/8; 1/3 1/3 1/3; 0 0 0;
        1/8 3/8 1/2; 1/2 1/4 1/4; 0 0 0;
    3 4; 2 2 3; 2 2 3; 0 0 0;    
        3/8 1/4 3/8; 1/8 1/4 5/8]
   GA = [
   1    % Zero rows are separators
   2 1 2; 1 4 4; 0 0 0;
   
   2 3 3; 3 2 2]
   1 2]
   A = [1 2 3 0 1 2 0
    
   

   md61
   type = 2
   PA = 0.3333    0.3333    0.3333
       0.2500    0.3750    0.3750
            0         0         0
       0.3333    0.3333    0.3333
   0.1250    0.3750    0.5000
       0.3750    0.2500    0.3750
   0.5000    0.2500    0.2500
            0         0         0
       0.1250    0.2500    0.6250
    
   GA =  1     3     4
        2     2     3
        2     2     3
        0     0     0
        2     1 
      2
        1     4     4
        0     0     0
        2     3     3
        3     2     2
    
   
   A =   1     2     3     0     1     2     0     1     2
   

   Once the data are entered into Matlab by the call for file “md61,” we make preparation for the policy-improvement step. The procedure is in the file newpolprep.m

   % file newpolprep.m
   % version of 3/23/92
   disp('Data needed:')
   transition probabilities for states and actions')
   disp('  Matrix GA of
   disp('  Matrix PA of
     gains for states and actions')
   disp('  Type number to identify GA matrix types')
   sition gains')
   disp('    Type 3.  Gains associated
   disp('    Type 1.  Gains associated with a state')
   disp('    Type 2.  One-step tra
   n with a demand')
   disp('  Row vector of actions')
   disp('  Value of alpha (= 1 for no discounting)')
   counting) ');
   PA = input('Enter matrix PA of transi
   c  = input('Enter type number to show gain type ');
   a  = input('Enter value of alpha (= 1 for no di
   stion probabilities ');
   GA = input('Enter matrix GA of gains ');
   if c == 3
   t('Enter row vector of demand probabilities ');
   end
   if c 
   PD = inp
   u== 1
   
   QA = GA';
   elseif c ==2
    QA = diag(PA*GA');     % (qx1)
   else
    QA = GA*PD';
   
   end
   

   A  = input('Enter row vector A of possible actions ');  % (1xq)
   m  = length(PA(1,:));
   q  = length(A);
   o approximate P0 ');
   s  = input('Enter s, the power of P  
   n  = input('Enter n, the power of P  
   tin the approximation of V ');
    col is index; second is A; third is QA
   DIS = ['    Index      Action 
   QD = [(1:q)' A' QA];    % Firs
   t   Value'];
   disp(DIS)
   disp(QD)
   if a == 1
    call = 'Call for newpol';
   else
    call = 'Call for newpold';
   
   end
   (call)
   dis
   p
   

   newpolprep   % Call for preparatory program in file npolprep.m
   Data needed:
    transition probabilities for states and actions
    Matrix GA o
    Matrix PA o
   ff gains for states and actions
    Type number to identify GA matrix types
    gains
      Type 3.  Gains associated with 
      Type 1.  Gains associated with a state
      Type 2.  One-step transitio
   na demand
    Row vector of actions
    Value of alpha (= 1 for no discounting)
   Enter type number to show gain type 2
   gains GA
   Enter row vector A of possible action
   Enter value of alpha (=1 for no discounting) 1
   Enter matrix PA of transition probabilities  PA
   Enter matrix GA of
    s  A
   Enter n, the power of P  to approximate P0  16
   Enter s, the power of P  in the approximation of V  6
      Index      Action    Value
   00    1.6250
      6.0000    2.0
   1.0000    1.0000    2.6667
      2.0000    2.0000    2.3750
      3.0000    3.0000    2.3333
      4.0000         0         0
      5.0000    1.0
   0000    2.5000
      7.0000         0         0
      8.0000    1.0000    2.6250
      9.0000    2.0000    2.1250
    
   
   Call for newpol
   

   The procedure has prompted for the procedure newpol (in file newpol.m)

   % file:  newpol.m  (without discounting)
   % version of 3/23/92
   icy as row matrix of indices ');
   D  = A(d);          
   d  = input('Enter po
   l     % policy in terms of actions
   probabilities for policy
   Q  = QA(d',:);           % selects
   P  = PA(d',:);           % selects
     next-period gains for policy
   P0 = P^n;                % Display to check convergence
   pected gain per period
   C = P + eye(P);
   for j=2:s
   PI = P0(1,:);            % Long-run distribution
   G = PI*Q                 % Long-run e
   x
    C = C + P^j;           % C = I + P + P^2 + ... + P^s
   end
   V = (C - (s+1)*P0 )*Q;  % B = B0*Q
   disp(' ')
   sp('Policy in terms of actions')
   disp(D)
   disp('Val
   disp('Approximation to P0; rows are long-run dbn')
   disp(P0)
   disp('Policy in terms of indices')
   disp(d)
   d
   iues for the policy selected')
   disp(V)
   disp('Long-run expected gain per period G')
   disp(G)
    Action   Test Value'];
   disp(DIS)
   disp(T)
   disp('Check current policy a
   T = [(1:q)' A' (QA + PA*V)];  % Test values for determining new policy
   DIS =['    Index   
    gainst new test values.')
   disp('--If new policy needed, call for newpol')
   
   disp('--If not, use policy, values V, and gain G, above')
   

   newpol
   Enter policy as row matrix of indices [2 6 9]  % A deliberately poor choice
    
   Approximation to P0; rows are long-run dbn
   0.2642    0.2830    0.4528
   ms of indices
       2     6    
   0.2642    0.2830    0.4528
      0.2642    0.2830    0.4528
    
   Policy in te
   r 9
    
   Policy in terms of actions
   
   2     2     2
   

   Long-run expected gain per period G
      2.2972
    Action   Test Value
      1.0000  
   Index   
      1.0000    2.7171
      2.4168
      3.0000    3.0000
   2.0000    2.0000
        2.3838
      4.0000         0         0
    2.0000    2.5677
      7.0000  
   5.0000    1.0000    1.6220
      6.0000  
           0         0
      8.0000    1.0000    2.6479
   cy against new test values.
   -
   9.0000    2.0000    2.0583
    
   Check current pol
   i-If new policy needed, call for newpol
       % New policy is needed
   newpol
    
   Enter policy as row matrix of 
   --If not, use policy and gain G, above 
    indices [1 6 8]
    
   Approximation to P0; rows are long-run dbn
      0.3939    0.2828    0.3232
    1     6     8
    
   Policy in te
   0.3939    0.2828    0.3232
      0.3939    0.2828    0.3232
    
   Policy in terms of indices
     
    rms of actions
       1     2     1
    
   Values for selected policy
      0.0526
     -0.0989
      0.0223
    
   .0000    2.6587
      2.0000    2.0000
   Long-run expected gain per period G
      2.6061
    
      Index     Action   Test Value
      1.0000    
   1    2.3595
      3.0000    3.0000    2.3254
      4.0000         0         0
      5.0000    1.0000    1.6057
    2.0000    2.1208
    
   Check cur
   6.0000    2.0000    2.5072
      7.0000         0         0
      8.0000    1.0000    2.6284
      9.0000  
    rent policy against new test values.
   --If new policy needed, call for newpol    
   
   --If not, use policy, values V, and gain G, above
   

   Since the policy selected on this iteration is the same as the previous one, the procedure has converged to an optimal policy. The first of the first three rows is maximum; the second of the next two rows is maximum; and the first of the final two rows is maximum. Thus, we have selected rows 1, 5, 6, corresponding to the optimal policy d^*∼(121). The expected long-run gain per period g=2.6061.

   The value matrix is

   v=[,

   v1

   v2

   v3

   ,]=[

   0.0527

   −0.0989

   0.0223

   ,]

   (54)

   We made a somewhat arbitrary choice of the powers of P used in the approximation of P₀ and B₀. As we note in the development of the approximation procedures in Sec 4, the convergence of Pⁿ is governed by the magnitude of the largest eigenvalue other than one. We can always get a check on the reliability of the calculations by checking the eigenvalues for P corresponding to the presumed optimal policy. For the choice above, we find the eigenvalues to be 1, -0.1250, 0.0833. Since0.125⁴≈0.0002, the choices of exponents should be quite satisfactory. In fact, we could probably use P⁸ as a satisfactory approximation to P₀, if that were desirable. The margin allows for the possibility that for some policies the eigenvalues may not be so small. — □
7. Policy iteration– with discounting 
   It turns out that the policy iteration procedure is simpler in the case of discounting, as suggested by the evolution analysis in Sec 4. We have the following two-phase procedure, based on that analysis.
   1. Value-determination procedure. Given the q_i and transition probabilities p(i,j) for the current policy, solve v=q+αPv to get for the corresponding v_i
      v=[I−αP]⁻¹q

      (55)
   2. Policy-improvement procedure Given {v_i:1≤i≤M{ for the current policy, determine a new policy d^*, with each d^*(i)=a_i* determined as that action for which
      q_i^*+α

      M

      ∑

      j=1

       p^*(i,j)v_j=max{q_i(a_ik)+

      M

      ∑

      j=1

       p_ij(a_ik)v_ja_ik∈A_i{

      (56)
   We illustrate the Matlab procedure with the same numerical example as in the previous case, using discount factor a=0.8. The data file is the same, so that we call for it, as before.

   EXAMPLE 5

   Assume data in file md61.m is in Matlab; if not, call for md61. We use the procedure newpolprep to prepare for the iterative procedure by making some initial choices.

   newpolprep
   Data needed:
    Matrix PA of transition probabilities for states and actions
    Matrix GA of gains for states and actions
    Type 1.  Gains associated with a state
    Type number to identify GA matrix types
    
       Type 2.  One-step transition gains
      Type 3.  Gains associated with a demand
   
    Row vector of actions
   r no discounting)
    Value of alpha (= 1 f
   o

   Enter type number to show gain type 2
   Enter value of alpha (= 1 for no discounting) 0.8
   Enter matrix PA of transition probabilities PA
   of possible actions A
   Enter
   Enter matrix GA of gains GA
   Enter row vector A
     n, the power of P  to approximate P0 16
   ion of V 6
    
      Index      Action  Test Value
      1.0
   Enter s, the power of P  in the approxima
   t000    1.0000    2.6667
      2.0000    2.0000    2.3750
   3.0000    3.0000    2.3333
   5000
      7.0000         0     
   4.0000         0         0
      5.0000    1.0000    1.6250
      6.0000    2.0000    2
   .    0
      8.0000    1.0000    2.6250
      9.0000    2.0000    2.1250
    
   
   Call for newpold
   

   The procedure for policy iteration is in the file newpold.m.

   % file newpold.m  (with discounting)
   % version of 3/23/92
   cy as row matrix of indices ');
   D = A(d);
   P = PA(d',
   d = input('Enter pol
   i:);        % transition matrix for policy selected
   policy selected
    
   V = (eye(P) - a*P)\Q;            % Present values
   Q = QA(d',:);        % average next-period gain for
     for unlimited future
   T = [(1:q)' A' (QA + a*PA*V)];   % Test values for determining new policy
   disp(' ')
   disp('    Policy in terms of actions')
   disp(D)
   dis
   disp('Approximation to P0; rows are long-run dbn')
   disp(P0)
   disp('    Policy in terms of indices')
   disp(d)
   p('  Values for policy')
   disp(V)
   eeded, call for newpold')
   disp('--If not, us
   DIS =['    Index     Action    Test Value'];
   disp(DIS)
   disp(T)
   disp('Check current policy against new test values.')
   disp('--If new policy 
   ne policy and values above')

   newpold
   Enter policy as row matrix of indices [3 5 9]  % Deliberately poor choice
    
   Approximation to P0; rows are long-run dbn
   2828    0.3232
      0.3939    0
      0.3939    0.2828    0.3232
      0.3939    0
   ..2828    0.3232
    Policy in terms of indices
    
   3     5     9
    
     10.3778
      9.6167
     10.1722
   Policy in terms of actions
       3     1     2
    
    Values for policy
    
       2.0000   10.3872
      3.0000 
   Index     Action  Test Value
      1.0000    1.0000   10.7111
      2.000
   0   3.0000   10.3778
      4.0000         0         0
      5.0000    1.0000    9.6167
    10.7133
      9.0000    2.0000 
   6.0000    2.0000   10.6089
      7.0000         0         0
      8.0000    1.0000 
      10.1722
    
   Check current policy against new test values.
   --If new policy needed, call for newpold
   oximation to P0; rows are long-run db
   --If not, use policy and values above
    
   newpold
   Enter policy as row matrix of indices [1 6 8]
    
   App
   rn
      0.3939    0.2828    0.3232
      0.3939    0.2828    0.3232
      0.3939    0.2828    0.3232
    
      Policy in terms of indices
   
   1     6     8
    
      Policy in terms of actions
       1     2     1
    
    Values for policy
     13.0844
     12.9302
   
   13.0519
   

      Index     Action  Test Value
      1.0000    1.0000   13.0844
      3.0000    3.0000   12.7511
      2.0000    2.0000   12.7865
   4.0000         0         0
      7.0000         0         0
   5.0000    1.0000   12.0333
      6.0000    2.0000   12.9302
      8.0000    1.0000   13.0519
   icy needed, call for newpold
   9.0000    2.0000   12.5455
    
   Check current policy against new test values.
   --If new po
   l--If not, use policy and values above

   Since the policy indicated is the same as the previous policy, we know this is an optimal policy. Rows 1, 6, 8, indicate the optimal policy to be d^*∼(1,2,1). It turns out in this example that the optimal policies are the same for the discounted and undiscounted cases. That is not always true. The v matrix gives the present values for unlimited futures, for each of the three possible starting states.