Uczenie maszynowe i sztuczne sieci neuronowe/Wykład 12

Uogólnienie koncepcji MDP na przypadki ciągłe

Klasycznym przykładem problemu, który można sformułować w postaci ciągłego procesu decyzyjnego Markowa jest sterowanie odwróconego wahadła. Możemy sobie wyobrazić, że mamy wózek o masie [math]M[/math] który może poruszać się po 1-wymiarowym torze i na tym wózku na zawiasie o jednym stopniu swobody umieszczony jest nieważki pręt o długości [math]l[/math] na którego końcu znajduje się punktowa masa [math]m[/math]. Zadaniem naszym jest utrzymanie wahadła w pionie poprzez przykładanie do wózka odpowiedniej siły [math]F[/math]. Przestrzeń stanu dla tego problemu jest 4-wymiarowa, bo mamy [math](x,\dot{x}, \theta , \dot{\theta })[/math]

Najprostszym pomysłem na zastosowanie do tego przykładu (lub podobnych) metod uczenia z wymuszaniem (RL) jest dyskretyzacja przestrzeni stanów i przestrzeni akcji i zastosowanie algorytmów typu iteracja funkcji wartościującej lub iteracja strategii. Pomysł ten napotyka w praktyce na bardzo poważny problem ze skalowaniem znany jako "klątwa wymiarów". Polega ona na tym, że jeśli każdy z [math]n[/math] wymiarów przestrzeni próbkujemy na [math]k[/math] poziomów to otrzymujemy siatkę o [math]k^{n}[/math] węzłów. Widać zatem, że złożoność obliczeniowa będzie rosła wykładniczo z wymiarami przestrzeni.

Dyskretyzacja akcji

Zazwyczaj w interesujących problemach przestrzeń akcji jest mniej wymiarowa niż przestrzeń stanów (pomyślmy o tym, że sterowanie obiektami w wielu grach da się zrealizować za pomocą klawiszy strzałek).

Zanim przejdziemy do konkretnego algorytmu, wprowadzimy najpierw koncepcję symulatora/modelu dla MDP. Model jest sposobem na opis [math]P_{sa}[/math]. Zakładamy, że model to "czarna skrzynka" (kawałek kodu/program), która otrzymuje na wejściu stany i akcje a na wyjściu zwraca nowe stany.

[math](s,a) \rightarrow [model] \rightarrow s'[/math]

W naszym przykładzie z wahadłem modelem jest układ równań opisujących dynamikę wózka i wahadła:

[math]\left( M + m \right) \ddot{x} - m \ell \ddot{\theta }\cos \theta + m \ell \dot{\theta }^2 \sin \theta = F[/math]

[math]\ell \ddot{\theta }- g \sin \theta = \ddot{x} \cos \theta [/math]

Możemy go przepisać jako układ równań pierwszego rzędu w którym stan jest 4-wymiarowym wektorem:

[math] s = \left[ \begin{array}{ccc} x \\ \dot{x} \\ \theta \\ \dot{\theta }\end{array} \right] [/math]

Nasz symulator mógłby działać np. rozwiązując ten układ równań metodą Eulera pierwszego rzędu:

[math] s_{t+1} = s_{t} + \Delta t \dot{s} [/math]

Jeśli nie umiemy wypisać równań dynamiki układu ale mamy fizyczny model i np. operatora tego modelu to można spróbować rozwiązać go przez obserwację wielu realizacji zachowania układu (model, operator), wg poniższego algorytmu:

• Obserwujemy [math]m[/math] realizacji stanów i akcji, które doprowadziły do przejść między stanami:

  [math] s_{0}^{(1)} \left.\begin{array}{c}a_{0}^{(1)} \\ \rightarrow \\ \end{array}\right. s_{1}^{(1)} \left.\begin{array}{c}a_{1}^{(1)} \\ \rightarrow \\ \end{array}\right. s_{2}^{(1)} \cdots \left.\begin{array}{c}a_{T-1}^{(1)} \\ \rightarrow \\ \end{array}\right. s_{T}^{(1)} [/math]

  [math] \vdots [/math]

  [math] s_{0}^{(m)} \left.\begin{array}{c}a_{0}^{(m)} \\ \rightarrow \\ \end{array}\right. s_{1}^{(m)} \left.\begin{array}{c}a_{1}^{(m)} \\ \rightarrow \\ \end{array}\right. s_{2}^{(m)} \cdots \left.\begin{array}{c}a_{T-1}^{(m)} \\ \rightarrow \\ \end{array}\right. s_{T}^{(m)} [/math]

  Te realizacje dostarczają nam zbioru uczącego [math]\left\lbrace (s_{t},a_{t}), s_{t+1}\right\rbrace [/math]
• Za pomocą jednego z algorytmów uczenia z nauczycielem (np. regresja liniowa, regresja liniowa z jądrem, sieci warstwowe z algorytmem wstecznej propagacji błędów itd.) wytwarzamy przybliżenie odwzorowania

  [math]s_{t_+1} = f(s_{t},a_{t})[/math]

  Dla ustalenia uwagi załóżmy, że w naszym przykładzie z wahadłem będzie to odwzorowanie liniowe postaci:

  [math]s_{t+1} = A s_{t} + B a_{t}[/math]

  gdzie [math]A[/math] jest macierzą ([math]4 \times 4[/math]), zaś [math]B[/math] jest wektorem (u nas 4 wymiarowym). Zadaniem algorytmu uczącego byłoby wyznaczenie [math]A[/math] i [math]B[/math] takich, że:

  [math]\arg \min _{A,B} \sum _{i=1}^{m}\sum _{t=0}^{T-1}||s_{t+1}^{(i)} - (A s_{t}^{(i)} +B a_{t}^{(i)})||^{2} [/math]

• Teraz mamy dopasowany model, w wersji deterministycznej:

  [math]s_{t+1} = As_{t} + Ba_{t}[/math]

  Możemy też rozważać wersję stochastyczną postaci:

  [math]s_{t+1} = As_{t} + Ba_{t} + \epsilon _{t}[/math]

  gdzie [math]\epsilon \sim N(0,\sigma ^{2})[/math]

Estymacja funkcji wartościującej

W przypadku stanów dyskretnych podawaliśmy funkcję wartościującą [math]V(s)[/math] w postaci tablicy wartości: dla każdego stanu mieliśmy przypisaną jakąś liczbę. Funkcję [math]V[/math] można było znaleźć rozwiązując równania Bellmana (dla [math]n[/math] stanów mieliśmy układ [math]n[/math] równań liniowych z [math]n[/math] niewiadomymi). W przypadku ciągłym nie da się zastosować tego podejścia bezpośrednio. Musimy dopasować jakąś funkcję ciągłą, która przybliżałaby [math]V[/math]. Można zastosować podejście analogiczne jak przy regresji. Najpierw trzeba wybrać jakieś cechy [math]\phi (s)[/math], którymi będziemy charakteryzować stan [math]s[/math]. W najprostszym przypadku może to być sam stan: [math]\phi (s) = s[/math] (w przypadku wahadła [math][x,\dot{x}, \theta , \dot{\theta }]^{T}[/math]). Ale może okazać się korzystne wzbogacenie wektora cech o jakieś dodatkowe składowe np.: [math]\phi (s) = [1, x, \dot{x}, \dot{x} ^{2}, \theta x, \dot{\theta }^{2}, \dots ]^{T}[/math] Wówczas funkcję wartościującą można przybliżyć przez:

[math]V(s) = \beta ^{T} \phi (s) [/math]

gdzie [math]\beta [/math] to wektor parametrów regresji liniowej (Jest to równanie analogiczne do tego, które stosowaliśmy przy regresji liniowej [math]\rightarrow [/math] tam mieliśmy [math]y = h_{\theta }(x) = \theta ^{T} \phi (x)[/math], tu żeby nie przeciążać notacji wektor parametrów regresji nazwaliśmy [math]\beta [/math].)

Przypomnijmy, że najważniejszym punktem algorytmu iteracji funkcji wartościującej było uaktualnianie funkcji wartości zgodnie z równaniem:

[math]V(s) = R(s) + \gamma \max _{a} \sum _{s^{\prime }} P_{sa}V(s^{\prime })[/math]

co bardziej ogólnie można wyrazić:

[math]V(s) = R(s) + \gamma \max _{a} E_{s^{\prime } \sim P_{sa}}[V(s^{\prime })][/math]

Widać, że przejście do przypadku ciągłych stanów wymaga znalezienia sposobu wyznaczania wartości oczekiwanej: [math]E_{s^{\prime } \sim P_{sa}}[V(s^{\prime })][/math].

Po tym wstępie przyjrzyjmy się algorytmowi dopasowywanej funkcji wartościującej:

• Losowo próbkuj [math]m[/math] stanów [math] s^{1}, s^{(2)}, \dots , s^{(m)} \in S[/math]
• Inicjalizuj [math]\beta = 0[/math]
• Powtarzaj {

  Dla [math]i = 1, \dots , m \lbrace [/math]

  Dla każdej akcji [math]a \in A[/math] {

  Próbkuj [math]s_{1}^{\prime }, s_{2}^{\prime }, \dots , s_{k}^{\prime } \sim P_{s^{(i)}a}[/math] (używając modelu/symulatora MDP)

  przypisz [math]q(a) = \frac{1}{k} \sum _{j=1}^{k} R(s^{(i)}) + \gamma V(s_{j}^{\prime })[/math]

  // wielkość [math]q(a)[/math] jest estymatorem [math]R(s) + \gamma E_{s^{\prime } \sim P_{sa}}[V(s^{\prime })][/math]

  }

  Przypisz [math]y^{(i)} = \max _{a} q(a)[/math]

  // tak obliczone [math]y^{(i)}[/math] jest estymatorem [math]R(s) + \gamma \max _{a} E_{s^{\prime } \sim P_{sa}}[V(s^{\prime })][/math]

  }

  // w wersji dyskretnej alg. iteracji f. wartościującej uaktualnialiśmy

  //[math]V(s^{(i)}) = y^{(i)}[/math]

  //w tym algorytmie, ponieważ [math]V(s)[/math] nie jest tablicą ale funkcją ciągłą

  //dopasowujemy parametry f. wartościującej (jeśli jest liniowa to np. za pomocą regresji liniowej,

  //można w tym kroku zastosować jakiś inny algorytm uczenia z nauczycielem)

  Przypisz [math]\beta = \arg \min _{\beta } \frac{1}{2} \sum _{i=1}^{m}(\beta ^{T} \phi (s^{(i)}) - y(i))^{2}[/math]

  }

Zauważmy, że w przypadku gdy model/symulator MDP jest deterministyczny to można uprościć powyższy algorytm kładąc [math]k=1[/math]. Jest tak ponieważ wartość oczekiwana w przypadku deterministycznym jest dokładnie równa wartości obliczonej jako [math]V(s^{\prime })[/math] ([math]s^{\prime }[/math] jest jednoznaczne mając [math]s[/math] i [math]a[/math]). W przypadku stochastycznym trzeba empirycznie stwierdzić do jakich stanów [math]s^{\prime }[/math] przejdziemy i jaka jest średnia [math]V(s^{\prime })[/math].

Jeśli chodzi o funkcję nagrody [math]R(s)[/math] zwykle ją zadajemy zgodnie z intuicją wynikającą z problemu. W przypadku wahadła można przyjąć np.:

[math] R = \left\lbrace \begin{array}{ r l} -1 & \text{gdy wahadło się przewróci } \\ 0 & \text{w przeciwnym razie} \end{array} \right.[/math]

Jak wyznaczyć strategię?

Algorytm estymowanej funkcji wartościującej oblicza przybliżoną postać funkcji [math]V^{*}[/math]. Akcje podejmujemy wtedy w następujący sposób:

[math]a = \arg \max _{a} E_{s^{\prime } \sim P_{sa} } [V(s^{\prime })][/math]

W ogólnym przypadku stochastycznym dla policzenia wartości oczekiwanej w powyższym wzorze konieczne jest pobranie próbki stanów [math]s^{\prime } \sim P_{sa}[/math] i obliczenie wartości średniej [math]V(s^{\prime })[/math]. Istnieją jednak dwa ważne przypadki kiedy można uprościć obliczenia wartości oczekiwanej:

• dla modelu/symulator deterministycznego do obliczenia wartości oczekiwanej wystarczy wyliczyć [math]V(s^{\prime })[/math] i wybrać tą akcję, która prowadzi do stanu [math]s^{\prime }[/math] o największej funkcji wartościującej. Zatem mamy:

  [math]a = \arg \max _{a} V(f(s,a))[/math]

• dla stochastycznego symulatora/modelu postaci:

  [math]s_{t+1} = f(s_{t},a_{t}) + \epsilon _{t}[/math]

gdzie [math]\epsilon _{t}[/math] jest próbką szumu gusowskiego [math]\sim N(0,\sigma ^{2})[/math]. W tym przypadku można zastosować przybliżenie:

[math]E_{s^{\prime }\sim P_{sa}}[V^{*}(s^{\prime })] \approx V^{*}(E[s^{\prime }]) = V^{*}(f(s,a))[/math]

I w tym przybliżeniu akcję wybieramy jako:

[math]a = \arg \max _{a} V(f(s,a))[/math]

Np. w naszym przykładzie

[math]V^{*}(s) = \beta ^{T} \phi (s) [/math]

Przykładowa implementacja kontrolera wahadła

Poniższy przykład zaczerpnięto z: http://mechanistician.blogspot.com/2009/05/lec17-fitted-value-iteration_31.html:

# -*- coding: utf-8 -*-
import matplotlib
matplotlib.use('TkAgg')
import sys
from pylab import *
from numpy import *

#To add more artificial features just append them to feats.
def get_feat_vec(x,x_dot,theta,theta_dot):
    feats=[1.0,
           x,
           x_dot,
           theta,
           theta_dot,
           x_dot**2,
           theta_dot**2,
           theta*x]
    return array(feats).reshape(1,len(feats))

#Displays the cart pole graph.
class show_cart_params: pass
show_cart_args=show_cart_params()
show_cart_args.first=True
show_cart_args.trial=0
def show_cart(x,x_dot,theta,theta_dot):
    figure(1)
    length=3
    x1=x
    y1=0
    x2=x+length*sin(theta)
    y2=length*cos(theta)
    cart_linex1,cart_liney1=[x-.4,x+.4],[-.25,-.25]
    cart_linex2,cart_liney2=[x+.4,x+.4],[-.25,0]
    cart_linex3,cart_liney3=[x+.4,x-.4],[0,0]
    cart_linex4,cart_liney4=[x-.4,x-.4],[0,-.25]
    base_linex,base_liney=[x,x],[-.5,-.25]
    axis([-3,3,-0.5,3.5])
    global show_cart_args
    
    if show_cart_args.first:
        show_cart_args.first=False 
        show_cart_args.pend_line,=plot([x1,x2],[y1,y2],c='black')
        show_cart_args.cart_line1,=plot(cart_linex1,cart_liney1,c='cyan')
        show_cart_args.cart_line2,=plot(cart_linex2,cart_liney2,c='cyan')
        show_cart_args.cart_line3,=plot(cart_linex3,cart_liney3,c='cyan')
        show_cart_args.cart_line4,=plot(cart_linex4,cart_liney4,c='cyan')
        show_cart_args.base_line,=plot(base_linex,base_liney,c='cyan')

    show_cart_args.pend_line.set_data([x1,x2],[y1,y2])
    show_cart_args.cart_line1.set_data(cart_linex1,cart_liney1)
    show_cart_args.cart_line2.set_data(cart_linex2,cart_liney2)
    show_cart_args.cart_line3.set_data(cart_linex3,cart_liney3)
    show_cart_args.cart_line4.set_data(cart_linex4,cart_liney4)
    show_cart_args.base_line.set_data(base_linex,base_liney)
    if show_cart_args.trial!=num_failures:
        title('trial #%d'%num_failures)
        show_cart_args.trial=num_failures
    draw()

#Returns whether or not the system is out of bounds.
def is_terminal(x,theta):
    return -2.4>x or x>2.4 or -twelve_degrees>theta or theta>twelve_degrees

#Simulates cart pole dynamics.
def cart_pole(action,x,x_dot,theta,theta_dot):   
    gravity=9.8
    masscart=1.0
    masspole=.3
    total_mass=masspole+masscart
    length=.7
    polemass_length=masspole*length
    force_mag=10.0
    tau=.02
    fourthirds=1.3333333333333

    action_flip_prob=0.0
    force_noise_factor=0.0
    no_control_prob=0.0

    if action_flip_prob>random.random():
        action=1-action

    if action>0:
        force=force_mag
    else:
        force=-force_mag

    force=force* \
           (1-force_noise_factor+random.random()*2*force_noise_factor)

    if no_control_prob>random.random():
        force=0

    costheta=cos(theta)
    sintheta=sin(theta)

    temp=(force+polemass_length*theta_dot*theta_dot*sintheta)/total_mass
    thetaacc=(gravity*sintheta-costheta*temp)/(length*(fourthirds \
            -masspole*costheta*costheta/total_mass))
    xacc=temp-polemass_length*thetaacc*costheta/total_mass

    new_x=x+tau*x_dot
    new_x_dot=x_dot+tau*xacc
    new_theta=theta+tau*theta_dot
    new_theta_dot=theta_dot+tau*thetaacc

    return get_feat_vec(new_x,new_x_dot,new_theta,new_theta_dot)

#################################################
#             Symulacja
#################################################

gamma=0.995
feat_vec=get_feat_vec(0,0,0,0)
n=shape(feat_vec)[1]
twelve_degrees=0.2094384
m=100000
w_vec=zeros((n,1))
feat_mat=ones((m,n))
diff=sys.maxint
new_diff=diff-1
epochs=0

#generate a bunch of state samples to fit the value function
x=random.uniform(-1.1,1.1);x_dot=0.0;theta=0.0;theta_dot=0.0
for i in range(m):
    if random.random()>.5:
        action=1
    else:
        action=0    
        
    feat_vec=cart_pole(action,x,x_dot,theta,theta_dot)
    feat_mat[i,:]=feat_vec

    if is_terminal(feat_vec[0,1],feat_vec[0,3]):
        x=random.uniform(-1.1,1.1);x_dot=0.0
        theta=random.uniform(-twelve_degrees,twelve_degrees);theta_dot=0.0
    else:
        x=feat_vec[0,1];x_dot=feat_vec[0,2]
        theta=feat_vec[0,3];theta_dot=feat_vec[0,4]

#Iteratively score each of the generated states and use the
#resulting feature,target tuple to estimate the parameters
#of a linear function.
while 10>epochs and diff>new_diff:
    diff=new_diff
    epochs+=1

    y_vec=zeros((m,1))   
    for i in range(m):
        feat_vec=feat_mat[i,:].reshape(1,n)
        x=feat_vec[0,1];x_dot=feat_vec[0,2]
        theta=feat_vec[0,3];theta_dot=feat_vec[0,4]

        # POBRANIE NAGRODY DLA AKTUALNEGO STANU
        if is_terminal(x,theta):
            R=-1.0
        else:
            R=0.0         
        
        # OBLICZENIE FUNKCJI WARTOŚCIUJĄCEJ DLA AKCJI = 0 
        feat_vec0=cart_pole(0,x,x_dot,theta,theta_dot)
        score1=R+gamma*dot(feat_vec0,w_vec)
        # OBLICZENIE FUNKCJI WARTOŚCIUJĄCEJ DLA AKCJI = 1 
        feat_vec1=cart_pole(1,x,x_dot,theta,theta_dot)
        score2=R+gamma*dot(feat_vec1,w_vec)
        
        if score1>score2:
            action=0
            score=score1
        elif score2>score1:
            action=1
            score=score2
        else:
            if random.random()>.5:
                action=1
                score=score2
            else:
                action=0
                score=score1

        y_vec[i,0]=score
    
    result=linalg.lstsq(feat_mat,y_vec) # tu fitujemy model liniowy
    w_new=result[0]
    new_diff=result[1]
    w_vec=w_new+0#copy

###### Ilustracja działania nauczonego kontrolera ##########

ion()
time=0
time_steps_to_failure={}
num_failures=0
time_at_start_of_current_trial=0
x=random.uniform(-1.1,1.1);x_dot=0.0;theta=0.0;theta_dot=0.0
display_started=False
max_num_failures_to_run=250

#Run the learned controller a few times to plot and compare
#against the discretized version.
while max_num_failures_to_run>num_failures:
    
    feat_vec0=cart_pole(0,x,x_dot,theta,theta_dot)
    score1=dot(feat_vec0,w_vec)
    
    feat_vec1=cart_pole(1,x,x_dot,theta,theta_dot)
    score2=dot(feat_vec1,w_vec)
    
    if score1>score2:
        action=0
    elif score2>score1:
        action=1
    else:
        if random.random()>.5:
            action=1
        else:
            action=0

    feat_vec=cart_pole(action,x,x_dot,theta,theta_dot)
    x=feat_vec[0,1];x_dot=feat_vec[0,2]
    theta=feat_vec[0,3];theta_dot=feat_vec[0,4]
    time+=1

    if display_started:
        show_cart(x,x_dot,theta,theta_dot)

    if is_terminal(x,theta):
        num_failures+=1
        if num_failures>max_num_failures_to_run-2:
            display_started=True

        time_steps_to_failure[num_failures]= \
            time-time_at_start_of_current_trial
        time_at_start_of_current_trial=time

        x=random.uniform(-1.1,1.1);x_dot=0.0;theta=0.0;theta_dot=0.0

#Plot learning curve:
figure(2)
plot(time_steps_to_failure.values(),'k')
title('Time steps per trial')
ioff()
show()