[PYTHON] Simulation d'estimation de trajectoire à l'aide de SLAM à base de graphes

1.Tout d'abord

Selon "Probability Robotics" [^ 1], SLAM est à peu près divisé en __ SLAM__ en ligne </ font> et __ SLAM hors ligne__ </ font> 2 Il est classé en deux méthodes. __ SLAM en ligne__ </ font> signifie que la probabilité postérieure de la carte et de la posture est calculée à chaque fois.

\begin{align}
    p(\ x_t,\ m\ |\ z_{1:t},\ u_{1:t}\ ) \tag{1.1}\\
\end{align}

Comme vous pouvez le voir à partir de la formule $ (1.1) $, SLAM en ligne </ font> a la dernière attitude vis-à-vis des observations passées et ne contrôle que $ x_t $ __ __ Est en cours de calcul. SLAM en ligne </ font> est un __ algorithme en temps réel __ pour calculer séquentiellement les dernières attitudes. L'algorithme d'estimation de l'auto-position par le filtre de Kalman étendu [^ 2] et le filtre de particules [^ 3] que j'ai créé précédemment correspond à ce SLAM en ligne </ font>.

D'autre part, dans __ SLAM__ </ font> hors ligne (parfois appelé complete SLAM), non seulement la dernière posture $ x_t $, mais toute la trajectoire du robot $ x_ {1 La probabilité postérieure est calculée pour: t} $.

\begin{align}
    p(\ x_{1:t},\ m\ |\ z_{1:t},\ u_{1:t}\ ) \tag{1.2}\\
\end{align}

SLAM hors ligne </ font> est un algorithme batch car le calcul est exécuté à un moment précis.

2 Graph-Based SLAM Le SLAM courant de nos jours est le SLAM hors ligne utilisant des graphiques, et le SLAM hors ligne utilisant des graphiques est généralement appelé Graph-Based SLAM [^ 4]. Cette fois, j'aimerais implémenter ce SLAM basé sur des graphes en Python et vérifier son fonctionnement. Les articles auxquels il est fait référence lors de l'implémentation sont "Un didacticiel sur le SLAM basé sur un graphique [^ 5]" et "Une explication sur le SLAM basé sur un graphique [^ 6]".

3 Posture du robot et modèle de mouvement

Le modèle de posture et de mouvement du robot s'applique au "Chapitre 5 Robot Motion" dans le livre "Probabilistic Robotics" [^ 1], veuillez donc vous référer au livre pour plus de détails. Par conséquent, une brève explication sera donnée ici. Tout d'abord, la posture du robot $ \ boldsymbol {x} $ à un instant discret $ t $ est définie comme suit.

\begin{align}
    \boldsymbol{x}_{t} =
       \begin{bmatrix}
        x_{t} \\
        y_{t} \\
        \theta_{t} \\
       \end{bmatrix} \tag{3.1}\\
\end{align}

Au temps $ t $, le robot reçoit __ vitesses de traduction $ v_ {t} $ __ et __ vitesses de rotation $ \ omega_ {t} $ _, et au temps $ t '$ posture $ \ boldsymbol x { On suppose que t '} $ est obtenu. Cependant, dans l'environnement réel, _translation speed $ v {t} $ __ et _rotation speed $ \ omega {t} $ __ sont bruyants, donc ils peuvent être exprimés par les formules suivantes.

\begin{align}
    \begin{bmatrix}
        x_{t'} \\
        y_{t'} \\
        \theta_{t'} \\
    \end{bmatrix}
    = 
    \begin{bmatrix}
        x_{t} \\
        y_{t} \\
        \theta_{t} \\
    \end{bmatrix}
    +
    \begin{bmatrix}
        -\frac{v_{t}}{\omega_{t}}sin\theta_{t}
           + \frac{v_{t}}{\omega_{t}}sin(\theta_{t}+\omega_{t}\Delta t) \\
        \frac{v_{t}}{\omega_{t}}cos\theta_{t}
           - \frac{v_{t}}{\omega_{t}}cos(\theta_{t}+\omega_{t}\Delta t) \\        
        \omega_{t}\Delta t + \gamma_{t}\Delta t \\
    \end{bmatrix}

    \tag{3.2}\\
\end{align}

cependant,

\begin{align}
    \begin{bmatrix}
        v_{t} \\
        \omega_{t} \\
        \gamma_{t} \\
    \end{bmatrix}
    = 
    \begin{bmatrix}
        v^\ast_{t} \\
        \omega^\ast_{t} \\
        0 \\
    \end{bmatrix}
    +
    \begin{bmatrix}
        \varepsilon_{\alpha_1 {v_{t}^\ast}^2 + {\alpha_2 \omega_{t}^\ast}^2} \\
        \varepsilon_{\alpha_3 {v_{t}^\ast}^2 + {\alpha_4 \omega_{t}^\ast}^2} \\
        \varepsilon_{\alpha_5 {v_{t}^\ast}^2 + {\alpha_6 \omega_{t}^\ast}^2} \\
    \end{bmatrix}

    \tag{3.3}\\
\end{align}

Est. Ici, $ v_ {t} ^ \ ast $ et $ \ omega_ {t} ^ \ ast $ représentent respectivement les vraies valeurs de la vitesse de translation et de la vitesse de rotation. De plus, $ \ varepsilon $ représente le bruit, et $ \ alpha_1 $ à $ \ alpha_6 $ sont des paramètres de bruit spécifiques au robot. Ici, $ \ varepsilon_b $ signifie une variable d'erreur avec une moyenne de zéro et un écart type de $ b $.

3.1 Définition du système de coordonnées

Lors de la mise en œuvre de SLAM basé sur des graphiques, les trois systèmes de coordonnées suivants sont définis.

__ Système de coordonnées mondial __

Un système de coordonnées basé sur un certain point du monde. Les points de repère fixes ne se déplacent pas dans le système de coordonnées mondial. world_system.png

__Robot système de coordonnées __

Un système de coordonnées basé sur un certain point du robot. Le point de référence se déplace avec le robot. Par conséquent, le repère stationnaire se déplacera à mesure que le robot se déplace. En outre, les points de repère sont essentiellement observés dans ce système de coordonnées de robot. robot_system.png

4 Modèle d'observation de référence

Pour le modèle d'observation de point de repère, "2.2 Observation" dans "Explication du SLAM basé sur un graphique [^ 6]" est appliqué, veuillez donc vous y référer. Par conséquent, une brève description sera donnée ici. En observant le repère $ L_c $ dans le système de coordonnées du robot à un instant discret $ t $, on suppose que des informations sur les trois repères suivants peuvent être obtenues.

--Distance au repère $ L_c $: $ d_ {c, t} $ --Direction où le repère $ L_c $ a été observé: $ \ varphi_ {c, t} $ --Landmark $ L_c $ Direction à laquelle vous êtes face: $ \ psi_ {c, t} $

observation.png

La formule est la suivante.

\begin{align}
    \begin{bmatrix}
        d_{c,t} \\
        \varphi_{c,t} \\
        \psi_{c,t} \\
    \end{bmatrix}
    = 
    \begin{bmatrix}
        d^\ast_{c,t} \\
        \varphi^\ast_{c,t} \\
        \psi^\ast_{c,t} \\
    \end{bmatrix}
    +
    \begin{bmatrix}
        \varepsilon_{\beta_1 {d^\ast_{c,t}}^2} \\
        \varepsilon_{\beta_2} \\
        \varepsilon_{\beta_3} \\
    \end{bmatrix}

    \tag{4.1}\\
\end{align}

$ d ^ \ ast_ {c, t} $, $ \ varphi ^ \ ast_ {c, t} $, $ \ psi ^ \ ast_ {c, t} $ représentent les vraies valeurs, respectivement, de $ \ beta_1 $ $ \ Beta_3 $ est un paramètre de bruit spécifique au robot.

5 Problème SLAM basé sur un graphique

5.1 Calcul de la posture relative du robot

Pour résoudre le problème, considérons un graphe avec des nœuds (sommets) et des arêtes (arêtes). Les nœuds correspondent à la posture du robot à chaque instant, et le bord est placé entre deux nœuds qui observent le même repère.

rel_pose.png

A ce moment, on peut voir qu'il existe les deux méthodes suivantes pour calculer la posture relative du robot à chaque instant. ① Différence de posture estimée du robot entre les nœuds à chaque instant: $ \ Delta \ boldsymbol x_ {t, t '} $

Cela peut être fait en calculant simplement la différence d'attitude.

\begin{align}
    \Delta \boldsymbol{x}_{t,t'} = \boldsymbol{x}_{t'} - \boldsymbol{x}_t
    \tag{5.1}\\
\end{align}
② Vecteur d'observation à chaque instant où le même repère a été observé: $ \ Delta \ boldsymbol x ^ {L} _ {c, t, t '} $

La figure ci-dessous montre le résultat de la conversion d'un repère en un système de coordonnées basé sur le repère basé sur le résultat de l'observation.

robot_system.png

A partir de cette figure, la posture relative du robot par rapport au repère $ \ boldsymbol x ^ L_ {c} $ est exprimée par la formule suivante.

\begin{align}
    \boldsymbol{x}^L_{c}
    = 
    \begin{bmatrix}
        x^L_{c} \\
        y^L_{c} \\
        \theta^L_{c} \\
    \end{bmatrix}
    = 
    \begin{bmatrix}
        d_{c} \, cos \left( \pi + \varphi_{c} - \psi_{c} \right) \\
        d_{c} \, sin \left( \pi + \varphi_{c} - \psi_{c} \right) \\
        \frac{\pi}{2} - \psi_{c} \\
    \end{bmatrix}

    \tag{5.2}\\
\end{align}

Ceci est calculé pour les temps $ t $ et $ t '$, et la différence est calculée.

\begin{align}
    \Delta \boldsymbol{x}^{L}_{c,t,t'} = \boldsymbol{x}^L_{c,t'} - \boldsymbol{x}^L_{c,t}
    \tag{5.3}\\
\end{align}

Considérons maintenant la différence entre $ \ Delta \ boldsymbol x_ {t, t '} $ et $ \ Delta \ boldsymbol x ^ L_ {c, t, t'} $.

\begin{align}
    \boldsymbol{e}_{c,t,t'} = \Delta \boldsymbol{x}_{t,t'} - \Delta \boldsymbol{x}^L_{c,t,t'}   \tag{5.4}\\
\end{align}

$ \ Delta \ boldsymbol x_ {t, t '} $ lorsque cette erreur $ \ boldsymbol e_ {c, t, t'} $ vaut zéro, c'est-à-dire $ \ boldsymbol x_ {t} $ et $ \ boldsymbol x_ { On peut dire que t '} $ est la posture estimée optimale. Cependant, puisque chaque arête contient une erreur et que chaque nœud est également associé à d'autres arêtes, il est peu probable que toutes les erreurs $ \ boldsymbol e_ {c, t, t '} $ soient nulles. C'est impossible. Par conséquent, il faut considérer une méthode de calcul de la posture estimée dans laquelle toutes les erreurs $ \ boldsymbol e_ {c, t, t '} $ peuvent être dites optimales.

5.2 Calcul de la matrice d'information

Pour le modèle de calcul de la matrice d'information, reportez-vous à "3.1.2 $ \ Sigma_ {c, t, t '}, \ Omega_ {c, t, t'} $ calcul" dans "Explication du SLAM basé sur un graphe [^ 6]". Il sera appliqué, veuillez donc vous y référer. Par conséquent, une brève description sera donnée ici.

\begin{align}
    \boldsymbol{\Sigma}_{c,t}
    = 
    \begin{bmatrix}
        \left(\beta_1 d^\ast_{c,t}\right)^2 & 0 & 0 \\
        0 & \left(d^\ast_{c,t} \ sin \beta_2\right)^2 & 0 \\
        0 & 0 & {\beta_2}^2 + {\beta_3}^2 \\
    \end{bmatrix}

    \tag{5.5}\\
\end{align}

Ensuite, convertissez $ \ boldsymbol \ Sigma_ {c, t} $ et $ \ boldsymbol \ Sigma_ {c, t '} $ calculés dans le système de coordonnées de mesure en système de coordonnées universel. La matrice de rotation suivante est utilisée pour la conversion.

\begin{align}
    \boldsymbol{R}_{c,t}
    = 
    \begin{bmatrix}
        c & -s & 0 \\
        s & c & 0 \\
        0 & 0 & 1 \\
    \end{bmatrix}

    \tag{5.6}\\
\end{align}

cependant,

\begin{align}
    c = cos(\theta_t + \varphi_t)       \tag{5.7}\\
    s = sin(\theta_t + \varphi_t)       \tag{5.8}\\
\end{align}

Est.

  • Le signe de la formule $ (10) $ et "s" dans "Explication of graph-based SLAM [^ 6]" dans la référence est différent. J'ai une direction positive dans le sens antihoraire, mais les références ont-elles une direction positive dans le sens horaire?

Enfin, la matrice de covariance $ \ boldsymbol \ Sigma_ {c, t, t '} $ et la matrice d'information $ \ boldsymbol \ Omega_ {c, t, t'} $ dans le système de coordonnées mondial sont les suivantes.

\begin{align}
    \boldsymbol{\Sigma}_{c,t,t'} &= \boldsymbol{R}_{c,t}\,\boldsymbol{\Sigma}_{c,t}\,\boldsymbol{R}_{c,t}^T
                     + \boldsymbol{R}_{c,t'}\,\boldsymbol{\Sigma}_{c,t'}\,\boldsymbol{R}_{c,t'}^T   \tag{5.9}\\
    \boldsymbol{\Omega}_{c,t,t'} &= \boldsymbol{\Sigma}_{c,t,t'}^{-1}    \tag{5.10}\\
\end{align}

5.3 Méthode d'optimisation

C'est finalement le problème d'optimisation qui est au cœur du sujet principal. Tout d'abord, le robot crée toutes les paires de deux arêtes des résultats d'observation obtenus à chaque instant. Soit l'ensemble de toutes les paires $ \ mathcal {C} $ et définissons une fonction d'évaluation pour trouver la somme des carrés des erreurs.

\begin{align}
    \boldsymbol{F}(\boldsymbol{x}) = \displaystyle \sum_{ \boldsymbol{e}_{c,t,t'} 
       \in \mathcal{C} }^{  } \ \underbrace{\boldsymbol{e}_{c,t,t'}^T \, \boldsymbol{\Omega}_{c,t,t'} \, \boldsymbol{e}_{c,t,t'}}_{\boldsymbol{F}_{t,t'}}   \tag{5.11}\\
\end{align}

Le but ultime est de trouver $ \ hat {\ boldsymbol {x}} $ qui minimise $ \ boldsymbol {F} (\ boldsymbol {x}) $ dans cette équation ci-dessus.

La méthode Gauss-Newton et la méthode Levenberg-Marquardt sont généralement utilisées comme méthode d'approche pour trouver la valeur minimale $ \ hat {\ boldsymbol {x}} $. Pour la méthode Gauss-Newton et la méthode Levenberg-Marquardt, les sites suivants peuvent être utiles, veuillez donc jeter un œil si vous êtes intéressé.

Cette fois, nous allons le résoudre par la méthode de Gauss-Newton de la même manière que "Un tutoriel sur le SLAM basé sur des graphes [^ 5]".

Tout d'abord, $ \ boldsymbol F_ {c, t, t '} $ et $ \ boldsymbol x_ {t'} $, respectivement, $ \ Delta \ boldsymbol x_ {t} $, $ \ Delta \ boldsymbol x_ {t ' } $ Prenons le cas du changement. Exprimé comme une expression

\begin{align}
    \boldsymbol{e}_{c}(\boldsymbol{x}_{t} + \Delta\boldsymbol{x}_{t},\,\boldsymbol{x}_{t'} + \Delta\boldsymbol{x}_{t'}) 
      \ &=\ \boldsymbol{e}_{c}(\boldsymbol{x}_{t,t'}+\Delta\boldsymbol{x}) \tag{5.12}\\
      \ &\simeq\ \boldsymbol{e}_{c,t,t'} + 
         \left.\frac{\partial \boldsymbol{e}_{c,t,t'}}
                    {\partial \boldsymbol{x}}
         \right|_{\boldsymbol{x}=\boldsymbol{x}_{t,t'}}\Delta\boldsymbol{x}
   \tag{5.13}
\end{align}

Il est exprimé comme. Cependant, la matrice Jacobi utilise $ \ boldsymbol {J} _c $

\boldsymbol{J}_{c,t,t'} = \left.\frac{\partial \boldsymbol{e}_{c,t,t'}}{\partial\boldsymbol{x}}\right|_{\boldsymbol{x}=\boldsymbol{x}_{t,t'}}   \tag{5.14}

Il est exprimé comme. Donc,

\begin{align}
    \boldsymbol{F}(\boldsymbol{x}+\Delta\boldsymbol{x}) &= \sum_{ \boldsymbol{e}_{c,t,t'} \in \mathcal{C} } \boldsymbol{F}(\boldsymbol{x}_{t,t'}+\Delta\boldsymbol{x})   \tag{5.15}\\
   &\simeq \sum_{ \boldsymbol{e}_{c,t,t'} \in \mathcal{C} } \ 
      (\boldsymbol{e}_{c,t,t'} + \boldsymbol{J}_{c,t,t'}\Delta{\boldsymbol{x}} )^T
      \, \boldsymbol{\Omega}_{c,t,t'} \,
      (\boldsymbol{e}_{c,t,t'} + \boldsymbol{J}_{c,t,t'}\Delta{\boldsymbol{x}} )   \tag{5.16}\\
   &= \sum_{ \boldsymbol{e}_{c,t,t'} \in \mathcal{C} } \ 
      \underbrace{\boldsymbol{e}_{c,t,t'}^T \boldsymbol{\Omega}_{c,t,t'} \boldsymbol{e}_{c,t,t'}}_{c_{c,t,t'}}
     + 2 \, \underbrace{\boldsymbol{e}_{c,t,t'}^T \boldsymbol{\Omega}_{c,t,t'} \boldsymbol{J}_{c,t,t'}}_{\boldsymbol{b}_{c,t,t'}} \Delta{\boldsymbol{x}} + \Delta{\boldsymbol{x}}^T \underbrace{\boldsymbol{J}_{c,t,t'}^T \boldsymbol{\Omega}_{c,t,t'} \boldsymbol{J}_{c,t,t'}}_{\boldsymbol{H}_{c,t,t'}}\Delta{\boldsymbol{x}}   \tag{5.17}\\

   &= \sum_{ \boldsymbol{e}_{c,t,t'} \in \mathcal{C} } \ 
      c_{c,t,t'} + 2 \, \boldsymbol{b}_{c,t,t'} \Delta{\boldsymbol{x}}
     + \Delta{\boldsymbol{x}}^T \boldsymbol{H}_{c,t,t'}\Delta{\boldsymbol{x}}   \tag{5.18}\\

   &= \underbrace{\sum_{ \boldsymbol{e}_{c,t,t'} \in \mathcal{C} } c_{c,t,t'}}_{\boldsymbol{c}}
     + 2 \underbrace{\left(\sum_{ \boldsymbol{e}_{c,t,t'} \in \mathcal{C} } \, \boldsymbol{b}_{c,t,t'} \right)}_{\boldsymbol{b}} \Delta{\boldsymbol{x}} 
     + \Delta{\boldsymbol{x}}^T \underbrace{\left( \sum_{ \boldsymbol{e}_{c,t,t'} \in \mathcal{C} }  \boldsymbol{H}_{c,t,t'} \right)}_{\boldsymbol{H}} \Delta{\boldsymbol{x}}   \tag{5.19}\\

   &= \boldsymbol{c}
    + 2 \boldsymbol{b} \Delta{\boldsymbol{x}}
    + \Delta{\boldsymbol{x}}^T \boldsymbol{H} \Delta{\boldsymbol{x}}   \tag{5.20}\\
\end{align}

$ \ Boldsymbol {c} $ dans l'expression $ (5.20) $ est une constante. De plus, lorsque le point de variation de la formule $ (5.20) $ est calculé,

\begin{align}
   2 \boldsymbol{b} + 2 \boldsymbol{H} \Delta{\boldsymbol{x}} = 0   \tag{5.21}\\
\end{align}

Par conséquent, $ \ Delta \ boldsymbol {x} $ est

\begin{align}
   \Delta{\boldsymbol{x}} = -\boldsymbol{H}^{-1} \boldsymbol{b}   \tag{5.22}\\
\end{align}

Sera. Mettez à jour $ \ boldsymbol {x} $ avec le $ \ Delta \ boldsymbol {x} $ calculé.

\begin{align}
   \boldsymbol{x} \leftarrow \boldsymbol{x}+\Delta \boldsymbol{x}   \tag{5.23}\\
\end{align}

Calculez $ \ boldsymbol {F} (\ boldsymbol {x}) $ de la formule $ (5.11) $ en utilisant la mise à jour $ \ boldsymbol {x} $, et la formule $ (5.12) $ jusqu'à ce que le résultat du calcul soit suffisamment petit. Répétez le calcul à partir de. Soit $ \ boldsymbol {x} $ lorsque le résultat du calcul est suffisamment petit, la solution optimale $ \ hat {\ boldsymbol {x}} $.

6 Mise en œuvre

J'ai implémenté SLAM basé sur des graphes en Python. Le script est publié sur GitHub. ・ GitHub --graph_based_slam.py

Le résultat de l'exécution du script est le suivant. (Cliquez sur l'image pour accéder à YouTube.) extended_kalman_filter

  • Une explication détaillée sera ajoutée plus tard.

7 Enfin

J'ai essayé d'implémenter le SLAM basé sur des graphiques, qui est populaire récemment, mais honnêtement, il était difficile de mettre en œuvre et d'écrire cet article. Parfois, cela ne fonctionnait pas, mais je suis heureux d'avoir réussi à le faire fonctionner comme prévu. C'est encore un article mince, mais je le publierai ici. Je pense qu'il y a de nombreux endroits où l'explication est insuffisante, mais si vous pouvez le souligner, je vais l'ajouter. De plus, le commentaire de M. Ueda a été très utile pour sa mise en œuvre.

  • "Explication du SLAM basé sur un graphique [^ 6]"

Ensuite, je voudrais exécuter SLAM sur la machine réelle ♪

Les références

[^ 1]: "Probability Robotics (Premium Books Version)", My Navi Publishing, 2005. [^ 2]: Visualisation de l'opération d'estimation de l'auto-position par filtre de Kalman étendu [^ 3]: Visualisation de l'opération d'estimation d'auto-position par filtre à particules [^ 4]: Masahiro Tomono: Reconnaissance environnementale des robots mobiles - Construction de cartes et estimation de l'auto-position