Jarekz: /* Przykładowa implementacja kontrolera wahadła */

2015-05-21T18:15:19Z

Przykładowa implementacja kontrolera wahadła

Jarekz: Utworzono nową stronę "==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..."

2015-05-21T18:14:00Z

Utworzono nową stronę "==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..."

Nowa strona

==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>
[[Plik:Cart-pendulum.png]]

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:

<ul>

<li>
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>

<li>
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>

<li>
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>

</ul>

==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:

<ul>

<li>
Losowo próbkuj <math>m</math> stanów <math> s^{1}, s^{(2)}, \dots , s^{(m)} \in S</math>
<li>
Inicjalizuj <math>\beta = 0</math>

<li>
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>
}

</ul>

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:

<ul>

<li>
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>

<li>
dla stochastycznego symulatora/modelu postaci:

::<math>s_{t+1} = f(s_{t},a_{t}) + \epsilon _{t}</math>

</ul>

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:
<source lang = py>
# -*- 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()

</source>

@@ Linia 198: / Linia 198: @@
 ===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:
-<source lang = py>
+<source lang = python>
 # -*- coding: utf-8 -*-
 import matplotlib

Uczenie maszynowe i sztuczne sieci neuronowe/Wykład 12 - Historia wersji

Jarekz: /* Przykładowa implementacja kontrolera wahadła */

Jarekz: Utworzono nową stronę "==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..."