基于改进动态系统稳定估计器的机器人技能学习方法

近年来, 机器人产业高速发展, 整体规模持续增长, 在制造业和服务业等众多领域都有广泛应用. 随着工业产品迭代速度日益增长, 个性化需要与日俱增, 传统依靠手工编程完成特定任务的方法难以适应新的需求. 因此, 迫切需要开发简单实用, 且可以灵活适用于多种任务的机器人技能学习方法.

机器人示教学习(Learning from demonstration, LfD)的灵感最初来源于人类的模仿学习, 近几年获得学术界和工业界的广泛关注[1–4]. 机器人通过观察用户演示来学习新技能, 同时将新技能泛化到不同场景下的相似任务中去, 一般包含演示、学习、复现三个阶段. 演示阶段需要解决的问题是如何向机器人进行技能演示, 常用的方法有视觉示教、动觉示教、遥操作示教和双臂示教. 学习阶段是对技能进行建模, 抽象的技能符号由示教数据具体表示, 然后利用示教数据训练模型参数. 复现阶段的性能主要体现在稳定性、复现精度、泛化能力和抗干扰性能4个方面.

动态系统(Dynamical system, DS)是对机器人技能进行建模的常用方法[5], 该方法将规划和执行集成到一起, 并将所有可能的解决方案嵌入到模型中以实现目标[6]. 在非线性DS基础上发展起来的动态运动原语只要进行一次演示就可以对运动进行建模[7], 动态运动原语描述的运动模型由非线性DS和线性DS组成, 其中非线性部分保证了轨迹复现的相似度, 线性部分则确保了模型全局稳定性, 两者的切换依靠相位变量平稳进行. 尽管动态运动原语提供了一种有效而精确的方法来对复杂的动态进行编码, 但是单变量编码丢失了各自由度之间的相关信息, 而且该方法本质上仍依赖于时间, 在面对时间扰动时需要用启发式方法重置相位变量[8].

为弥补动态运动原语的缺陷, 文献[9]提出了动态系统的稳定估计器(Stable estimator of dynamical systems, SEDS). 它首先利用高斯混合模型(Gaussian mixture models, GMM)和高斯混合回归(Gaussian mixture regression, GMR)的概率学习方法对轨迹进行初步拟合. 概率学习方法是轨迹编码中的常用方法, 它可以保留演示的固有可变性[10], 但是无法确保训练得到的动态系统具有全局稳定性. 因此SEDS在后续优化中加入了稳定性约束, 确保机器人在不受扰动的情况下能够到达目标点. 然而, 过于严格的稳定性约束可能会在学习过程中限制模型的精度. 针对SEDS方法中稳定性和精度难以平衡的问题, 文献[11]利用微分同胚变换改进了SEDS, 称作τ$\tau$-SEDS. 该方法在保证系统稳定性的同时, 很大程度上克服了模型的精度问题. 但是增加了模型复杂度, 导致学习更加耗时, 仅限于离线学习. DS方法的快速稳定学习[12]同时考虑了估计DS的稳定性、准确性和学习速度三个因素. DS方法快速稳定学习的快速学习能力在很大程度上方便了它的实际应用, 但在稳定性和准确性方面还不够优越. 文献[13]提出了一种基于流形浸没和淹没的学习方法来解决精度与稳定性的矛盾, 该方法保证了有效提取动力学特征和稳定形式的高精度, 而且能处理交叉运动的情形, 但模型复杂度较高.

SEDS中另一个值得注意的问题是混合高斯分量个数的选取, 但是对于该问题的相关研究较少. 过多或者过少的分量个数选取都会导致模型无法有效提取演示的动力学特征, 因此该问题具有一定的研究价值. 通常用于确定有限混合模型的最佳分量的方法是贝叶斯信息准则[14], 然而这种模型选择方法存在一些明显的缺陷[15], 常常过高估计模型分量的个数, 导致过拟合. 贝叶斯非参数模型是一种定义在无限维参数空间上的贝叶斯模型, 其利用在适当数量的模型分量密度上产生后验分布来调整模型大小, 因此可以根据数据自适应聚类个数, 其中狄利克雷过程混合模型是最常用的贝叶斯非参数模型之一[16-18]. Figureoa等[19]提出了一种物理一致的贝叶斯非参数混合模型, 该方法可以自动估计最佳的混合分量个数, 并且将相似性测度融入先验信息, 提高了复现和泛化的精度. 但是该方法使用吉布斯采样计算模型的后验概率, 计算复杂度较高.

鉴于SEDS存在的上述缺陷, 本文提出了改进的SEDS (Improved SEDS, i-SEDS). 该方法有效地解决了SEDS中稳定性和精度无法兼顾的问题, 并且可以自动确定合适的分量个数. 仿真以及Franka-panda协作机器人的实验结果验证了本文方法的有效性和优越性. 本文的主要贡献有: 1)使用狄利克雷过程高斯混合模型(Dirichlet process GMM, DPGMM)代替GMM拟合演示, 并利用变分推断(Variational inference, VI)训练模型, 该模型可以根据演示数据自动确定合适的混合分量个数. 仿真分析超参数对DPGMM的影响, 为超参数的选择提供了指导意义; 2)采用参数化的李雅普诺夫函数修改了原SEDS中的稳定性约束条件, 从而提高了学习轨迹的精度, 解决了稳定性和精度难以兼顾的问题.

1.   问题描述

将机器人离散运动公式化为由自治动态系统驱动的控制律. 考虑一个状态变量ξ∈Rd$\xi \in {\mathbb{R}}^d$, ξ$\xi$通常表示笛卡尔空间中末端执行器的坐标. 有时ξ$\xi$也可用于定义机器人系统的广义坐标(例如, 机械手的关节角度或轮式机器人的姿势)

式中, f(ξ)$f(\xi )$是一个非线性连续可微函数, 用于编码特定行为, 例如拾取搬运物品, 伸手拿杯子等. 从初始位置ξ0${\xi }_0$开始, 机器人运动ξt,t∈[0,+∞)${\xi }_t,t\in [0,+\mathrm{\infty })$由式(1)积分得到, 运动最后收敛至单一吸引子ξ∗${\xi }^{\ast }$, f(ξ∗)=0$f({\xi }^{\ast })=0$. 本文需要估计该吸引子上全局渐近稳定的动态系统f(ξ)$f(\xi )$.

通常, 将基于一组M条演示轨迹的数据{ξt,m,ξ˙t,m}Tm,Mt=0,m=1${\left\{{\xi }_{t,m},{\dot{\xi }}_{t,m}\right\}}_{t=0,m=1}^{T^m,M}$对f(ξ)$f(\xi )$进行估计, Tm$T^m$表示第m$m$条演示轨迹的序列长度, 估计得到的DS为f^(ξ)$\hat{f}(\xi )$. 一般采用回归算法进行估计, 如高斯过程回归[20]、局部加权映射回归算法[21]和混合高斯回归等. 但这些标准的回归技术在训练过程中没有考虑DS的渐近稳定性, f^(ξ)$\hat{f}(\xi )$很可能是不稳定的或收敛至伪吸引子. 因此SEDS方法将不稳定DS的建模转换为非线性优化问题, 构造了全局稳定的运动估计, 具体流程如图1所示.

图 1  基于SEDS的机器人示教学习流程图

Fig. 1  Flow chart of LfD based on SEDS

SEDS使用GMM-GMR回归方法从演示中估计初始DS. 首先利用具有K$K$个分量的GMM估计联合概率分布p(ξ,ξ˙|Θ)=∑Kk=1πkN(ξ,ξ˙|μk,Σk)$p(\xi ,\dot{\xi }|\mathrm{\Theta })={\sum }_{k=1}^K{\pi }_k\mathcal{N}\left(\xi ,\dot{\xi }|{\mu }_k,{\mathrm{\Sigma }}_k\right)$, 其中πk${\pi }_k$、μk${\mu }_k$和Σk${\mathrm{\Sigma }}_k$分别是各高斯分量的权重、均值和协方差, 即Θk={πk,μk,Σk}${\mathrm{\Theta }}_k=\{{\pi }_k,{\mu }_k,{\mathrm{\Sigma }}_k\}$, Θ={Θ1,⋯,ΘK}$\mathrm{\Theta }=\{{\mathrm{\Theta }}_1,\cdots ,{\mathrm{\Theta }}_K\}$, 然后GMR部分根据联合概率分布推导出后验概率p(ξ˙|ξ,Θ)$p(\dot{\xi }|\xi ,\mathrm{\Theta })$, 对其取均值得到:

式中, Ak=Σξ˙ξk(Σξk)−1,$A_k={\mathrm{\Sigma }}_k^{\dot{\xi }\xi }{\left({\mathrm{\Sigma }}_k^{\xi }\right)}^{-1},$ bk=μξ˙k−Akμξk,$b_k={\mu }_k^{\dot{\xi }}-A_k{\mu }_k^{\xi },$hk(ξ)=πkp(ξ|k)/∑Ki=1πip(ξ|i)$h_k\left(\xi \right)={\pi }_{k}p\left(\xi |k\right)/{\sum }_{i=1}^K{\pi }_{i}p\left(\xi |i\right)$. Akξ+bk$A_k\xi +b_k$表示第k$k$个线性DS, hk(ξ)$h_k\left(\xi \right)$是状态相关的非线性系数, ∑Kk=1hk(ξ)=1${\sum }_{k=1}^K{h}_k(\xi )=1$, 且hk(ξ)>0$h_k(\xi )>0$.

上述方法有足够的灵活性来建模各种运动, 但是无法确保DS的渐近稳定性, 为此, SEDS给出了在保证系统稳定性下的参数学习方法, 即通过最小化带非线性项约束的目标函数来估计参数Θ$\mathrm{\Theta }$.

式中, ≺0$\prec 0$表示矩阵的负定, T=∑Mm=1Tm$\mathcal{T}={\sum }_{m=1}^M{T}^m$是全部训练数据的总数量. 目标函数J(Θ)$J(\mathrm{\Theta })$的选取有均方误差(Mean square error, MSE)和对数似然(Log likelihood, LL)两种.

式中, 对数ln也可换成以其他常数为底的对数, 上式的约束是f^$\hat{f}$全局稳定收敛的充分条件, 这些条件由Lyapunov第二稳定性方法推导得到, 见定理1.

定理1. 如果存在径向无界且满足式(5)的Lyapunov函数V(ξ):Rd↦R$V(\xi ):{\mathbb{R}}^d\mapsto \mathbb{R}$, 则DS在ξ∗${\xi }^{\ast }$点处全局渐近稳定.

SEDS方法可以有效地从演示中估计出具有全局渐近稳定性的DS, 是LfD领域中的常用方法, 但是该方法存在两个问题: 1)用于估计初始DS的GMM-GMR方法需要人工指定混合分量个数K$K$, 当状态ξ$\xi$维度升高或演示轨迹变得复杂时, 通过经验确定超参数K$K$将变得很困难; 2)在SEDS中, V(ξ)$V(\xi )$是二次李雅普诺夫函数(Quadratic Lyapunov function, QLF), 即V(ξ)=(ξ−ξ∗)T(ξ−ξ∗)$V(\xi )=(\xi -{\xi }^{\ast }{)}^{\mathrm{T}}(\xi -{\xi }^{\ast })$. 在几何学上, QLF只允许L2$L_2$范数(即∥ξ−ξ∗∥2$\Vert \xi -{\xi }^{\ast }{\Vert }_2$)距离是单调递减的轨迹[11]. 这导致了SEDS在包含高曲率或非单调(暂时远离吸引子)的高度非线性运动中表现不佳, 不能同时保证稳定性和精度性.

针对问题1), 本文利用贝叶斯非参数模型能根据数据自动确定最佳聚类个数的特性, 采用了DPGMM拟合演示数据, 并推导GMR估计初始DS; 针对问题2), 本文选择参数化QLF (Parametrized QLF, P-QLF)来确保估计DS的稳定性, 即V(ξ)=(ξ−ξ∗)TP(ξ−ξ∗)$V(\xi )=(\xi -{\xi }^{\ast }){}^{\mathrm{T}}P(\xi -{\xi }^{\ast })$. 其中P$P$的作用是将简单的QLF重塑为 “椭圆” 形式, 放松了约束条件, 从而能够精确复现不收缩(高曲率且朝向目标具有非单调性)的演示轨迹.

2.   DPGMM-GMR算法

DPGMM-GMR算法由DPGMM的建模过程、利用变分推断进行模型求解和推导GMR估计初始DS三部分组成.

2.1   狄利克雷过程高斯混合模型

狄利克雷过程(Dirichlet process, DP)一般用DP(G0,α)$DP(G_0,\alpha )$表示, 其中G0$G_0$是连续的基分布, α$\alpha$是集中参数, 用来控制从DP中抽取的分布的离散化程度[22]. 假设从DP中随机抽取一个样本分布G$G$, 然后再从G$G$中独立地抽取N个随机变量{θ∗n}Nn=1$\{{\theta }_n^{\ast }{\}}_{n=1}^N$:

折棍过程(Stick breaking process, SBP)通常用于建模和推断DP[23]. 给定两个无限随机变量v={vk}∞k=1$v=\{v_k{\}}_{k=1}^{\mathrm{\infty }}$和{θ∗k}∞k=1$\{{\theta }_k^{\ast }{\}}_{k=1}^{\mathrm{\infty }}$, 其中vk$v_k$服从Beta分布, 即vk∼Beta(1,α)$v_k\sim \mathrm{B}\mathrm{e}\mathrm{t}\mathrm{a}(1,\alpha )$, 则SBP可以表示为:

式中, πk(v)=vk∏k−1i=1(1−vi)${\pi }_k\left(\mathit{v}\right)=v_k{\prod }_{i=1}^{k-1}(1-v_i)$, ∑∞k=1πk(v)=1${\sum }_{k=1}^{\mathrm{\infty }}{\pi }_k\left(v\right)=1$.

通过SBP, 可以将独立于G0$G_0$绘制的原子θk${\theta }_k$视为混合模型分量的分布参数, 该模型包含无限数量的分量, 并且每个分量的权重为πk(v)${\pi }_k(v)$, 从而建立了狄利克雷过程混合模型.

为了便于表述, 将全部演示数据表示为X={ξt,ξ˙t}Tt=1$X=\{{\xi }_t,{\dot{\xi }}_t{\}}_{t=1}^{\mathcal{T}}$, 其中xt={ξt,ξ˙t}$x_t=\{{\xi }_t,{\dot{\xi }}_t\}$. 假设每个观测值xt$x_t$都对应一个参数θ∗t${\theta }_t^{\ast }$, 服从概率密度p(xt|θ∗t)$p(x_t|{\theta }_t^{\ast })$. 整个参数集{θ∗t}Tt=1$\{{\theta }_t^{\ast }{\}}_{t=1}^{\mathcal{T}}$服从一个共同的DP先验, 由于分布G是离散的, 因此{θ∗t}Tt=1$\{{\theta }_t^{\ast }{\}}_{t=1}^{\mathcal{T}}$可以划分到相同簇内, 以πk(v)${\pi }_k(v)$的概率共享相同的θk${\theta }_k$. 然后引入指示变量Z={zt}Tt=1$Z=\{z_t{\}}_{t=1}^{\mathcal{T}}$, 也称作隐变量. 对于每个观测xt$x_t$, 都对应一个指示变量zt$z_t$, 用来指示数据xt$x_t$所属分量. 因此狄利克雷过程混合模型满足下式概率分布:

式中, Cat(⋅)$\mathrm{C}\mathrm{a}\mathrm{t}(\cdot )$表示分类分布.

在狄利克雷过程混合模型的基础上明确定义基分布G0$G_0$为Normal-Wishart分布, 即θk∼NIW(θ0)${\theta }_k\sim \mathcal{N}\mathcal{I}\mathcal{W}({\theta }_0)$, 就得到了DPGMM, 这意味着每个高斯分量的均值和精度都服从一个独立的Normal-Wishart分布, 根据式(8), DPGMM的公式化表示如下:

式中, μk和Λk${\mu }_k和{\mathrm{\Lambda }}_k$分别表示第k$k$个高斯分量的均值和精度, m0和β0$m_0和{\beta }_0$分别是μk${\mu }_k$的高斯先验的均值和方差, W0和ν0$W_0和{\nu }_0$分别是Λk${\mathrm{\Lambda }}_k$的Wishart先验的精度矩阵和自由度.

式(9)是给定的共轭先验分布, 式(10)是模型的似然函数. DPGMM的概率结构如图2所示. 然而在无限维的设定下, 实际中模型不易处理. 因此本文采用一种常见的策略, 即基于DP的截断SBP[24]. 预定义最大的混合分量数K^$\hat{K}$, 并且令q(vK^=1)=1$q(v_{\hat{K}}=1)=1$. 即对于所有k>K^$k>\hat{K}$, πk(v)${\pi }_k(v)$等于零.

图 2  DPGMM概率结构图

Fig. 2  DPGMM probability structure diagram

2.2   变分推断

第2.1节已为所有随机变量明确定义了共轭先验和模型的似然函数, 本节采用变分贝叶斯方法得出参数可靠的后验分布. 将隐变量和未知参数表示为Θ={v,Z,μ,Λ}$\mathrm{\Theta }=\{v,Z,\mu ,\mathrm{\Lambda }\}$, 超参数表示为Ξ={α,m0,β0,W0,ν0}$\mathrm{\Xi }=\{\alpha ,m_0,{\beta }_0,W_0,{\nu }_0\}$. 由于实际后验p(Θ|X,Ξ)$p(\mathrm{\Theta }|X,\mathrm{\Xi })$难以计算, 因此变分推断方法引入任意分布q(Θ)$q(\mathrm{\Theta })$以逼近p(Θ|Ξ,X)$p(\mathrm{\Theta }|\mathrm{\Xi },X)$. 在此假设下, 导出模型的对数似然值lnp(X)$\ln p(X)$, 即:

其中

KL(q||p)$KL(q||p)$表示变分后验 q(Θ)$q(\mathrm{\Theta })$和实际后验p(Θ|X,Ξ)$p(\mathrm{\Theta }|X,\mathrm{\Xi })$之间的KL(Kullback-Leibler)散度, 刻画了两个分布之间的相似程度. 为降低求解难度, 一般采用平均场理论限制概率分布q(Θ)$q(\mathrm{\Theta })$的范围, 然后寻找这个类别中使得KL散度最小化的概率分布[25]. 因此需要将q(Θ)$q(\mathrm{\Theta })$分解为:

同时, 考虑DPGMM的共轭先验配置, 根据共轭的定义, 期望后验q(Θ)$q(\mathrm{\Theta })$与先验p(Θ)$p(\mathrm{\Theta })$具有相同的函数形式[26], 可得:

指示变量zt$z_t$满足[27]:

式中, ⟨⋅⟩$⟨\cdot ⟩$表示数学期望.

采用变分贝叶斯期望最大化(Expectation maximization, EM)算法. 对该模型参数进行迭代求解. 该算法类似于EM算法, 在E步中, 计算指示变量的期望; 在M步中, 该期望值用于重新计算其他参数的变分分布.

期望通过变分推断得到DPGMM模型参数的后验分布后, 可以得到如下的模型联合概率密度: