利用机器人末端执行器对曲面工件进行跟踪，能够得到曲面工件的轮廓，精确的轮廓轨迹可以为机加工提供初始轨迹. 机器人在运动过程中由于自身刚度不足^[1-2]和外界环境的变化^[3]导致接触力不稳定同时波动较大，利用机器人力控制算法能够有效地保持接触力的稳定.

诸多学者对跟踪场景中的力控制算法进行研究，并在传统的力控制算法中加入了智能算法，如Tung等^[4]在控制器中加入了自适应模糊控制算法，在线优化了控制规则；Abu-Mallouh等^[5]使用神经网络，补偿了跟踪时系统滞后造成的误差和库仑摩擦力；Li等^[6]利用神经网络，对未知的跟踪环境进行分类；Ye等^[7]采用模糊控制方法，在线调整机器人阻抗模型中的阻尼刚度参数. 智能算法^[8-10]应用在跟踪实验时，往往需要先获取先验知识和模型，先验知识和模型的准确程度决定了控制算法的好坏；强化学习不需要人为地选取先验模型，在与环境反复交互的过程中，能够自主地发现最优行为，减轻人工作业的难度，达到较好的控制效果. 强化学习在机器人实际应用中已经取得了较好的效果^[11-15].

强化学习与机器人结合时，使用的算法通常可以分为无模型强化学习和有模型强化学习^[16]，无模型强化学习算法虽然具有一定的通用性，但是需要同外界环境进行的交互次数较多，在实际应用中效率较低. 有模型强化学习算法构建了模型，提高了数据利用率，同时在交互过程中模型会不断地修正，因此在实际应用中比无模型强化学习交互的次数少，成为研究的热点. 有模型强化学习中，基于概率动力学模型的强化学习算法（PILCO）的收敛速度快，在实际中有较好的应用^[17-19]，机器人恒力跟踪场景可以重复实验，有利于满足该算法的要求，因此可以将强化学习应用到机器人恒力跟踪场景中.

本文对机器人曲面恒力跟踪场景中机器人末端执行器的受力进行分析，对接触过程进行建模. 利用概率动力学模型，对机器人的接触状态进行拟合和预测，强化学习在概率动力学模型中搜索策略，推理机器人的输出控制参数. 在实验中，将强化学习的输入状态改为一段时间内的状态平均值，以减少接触状态下信号的干扰，最终得到稳定的跟踪力.

机器人恒力跟踪实验平台如图1所示，六维力传感器安装在机器人的末端，传感器坐标系{T}、工作台坐标系{S}和机器人基坐标系{B}的姿态相同. 机器人始终保持v_s的速度沿着V_x方向运动，如图2所示，当机器人末端执行器和工件接触后，机器人末端执行器会根据受力的大小在V_y方向微小地移动，从而调整机器人和工件的接触状态，其中机器人运动坐标系{V}相对于传感器坐标系{T}的姿态始终保持不变. 当末端执行器与曲面工件接触时，如图3所示，末端执行器受到曲面的法向力F_n和切向力F_τ，恒定的F_n能够反映出曲线的轮廓. 为了得到F_n，需要将曲面接触力坐标系{C}中的力映射到已知的传感器坐标系{T}中，曲面接触力坐标系{C}中心与传感器坐标系{T}中心重合，曲面接触力坐标系的X轴方向与曲面的切向方向相同，Y轴方向始终垂直于曲面轮廓，由图3的受力分析结果可得，F_n和F_τ为

(1) $\left. \begin{gathered} {F_{\text{τ}} } = { - ^T}{F_x}\cos \theta { - ^T}{F_y}\sin \theta {\text{，}} \\ {F_{\rm n}} = { - ^T}{F_x}\sin \theta { + ^T}{F_y}\cos \theta {\text{。}} \\ \end{gathered} \right\} $

式中：^TF_x和^TF_y分别为传感器坐标系{T}坐标系下X和Y方向的力由六维力传感器测得； $\theta $为X_T与X_C之间的夹角.

式(1)中夹角 $\theta $未知，因此需要对 $\theta $进行估计， $\theta $等于曲面的切线倾斜角. 当机器人末端与曲面接触不分离时，在一个采样周期内，通过得到V_x方向和V_y方向的位移可以求得该周期内角度的正切值，反正切运算后可得一个采样周期内曲面的切线倾斜角：

(2) $\left. \begin{aligned} & \tilde \theta = {\rm arctan}\left({{{v_{\rm{s}}}\Delta t}}/{{\Delta y}}\right) {\text{，}} \\ & \Delta y = kv\Delta t {\text{。}} \end{aligned} \right\}$

式中： $\widetilde \theta $为一个采样周期内的曲线倾斜角；v_s为机器人沿着V_x方向的速度； $\Delta t$为机器人沿着V_x方向每个周期的时间； $\Delta y$为机器人沿着V_y方向的偏移量，其大小与上位机在每个周期内给机器人的偏置电压成正比；v为偏置电压； $k$为常数. 由于工件曲面没有突变的轮廓，则 $\widetilde \theta $的曲线较平滑，然而实验时 $\widetilde \theta $包含较多的噪声和干扰，需要对 $\widetilde \theta $进行滤波以便满足实际的需要，滤波选用平均值滤波：

(3) ${\theta _{{t}}} =\left\{ \begin{aligned} &\sum\limits_{t=1}^i {{{\tilde \theta }_{{t}}}} /t, \quad t \leqslant i \leqslant w; \\ & \sum\limits_{t=i - w + 1}^i {{{\tilde \theta }_{{t}}}} /w, \quad w \leqslant t \leqslant i. \\ \end{aligned} \right.$

式中： ${\theta _{{t}}}$为t时刻滤波后的曲面切线倾斜角， ${\tilde \theta _{{i}}}$为i时刻的曲面切线倾斜角，w为滤波窗口常数. 滤波后得到 $\theta $的平均误差小于5°，满足实际应用的需要.

机器人末端执行器与环境是刚性接触，接触过程中振动较大，造成F_n波动很大，需要利用机器人力控制算法调整机器人与环境的接触状态，从而达到期望的接触效果. 机器人力控制算法有阻抗控制、导纳控制、力位混合控制等^[20]，在机器人曲面跟踪场景中期望的接触力是恒定的，显式力控制比较适用于该场景^[21]. 显式力控制是将末端执行器与接触面之间的实际接触力和期望的接触力进行比较，根据力误差调整机器人状态，以达到期望的力. 显式控制分为基于力的显式力控制和基于位置的显式力控制^[22]，基于力的显式力控制根据期望与实际的力之差调节机器人的关节力矩以实现恒力控制；基于位置的显式力控制将期望与实际的力之差转化为机器人的偏移位移，从而改变机器人当前的位置，间接调节接触力. 由于基于位置的显式力控制比较容易实现，选择该算法.

基于位置的显式力控制的控制流程如图4所示.将F_n与期望力F_d进行比较后得到力误差，力控制器将力误差转化为机器人的位移调整量，控制器根据位移调整量的大小发送一个相对位置命令后，机器人高精度的位置控制器及时调整机器人的位移，从而间接控制机器人与环境接触力. 在系统是非线性的情况下，基于位置的显式力控制器也容易于实现，具有较好的鲁棒性.

力控制器通常选择PD、PI或PID控制器^[23]，由于PD算法不会产生积分失控，容易获得较好的调节品质. 选择PD控制器，仅在一维方向上进行分析，控制器设计为

(4) $\left. \begin{gathered} \Delta {u_{{t}}} = {k_{\rm{p}}}(\Delta {f_{{t}}} - \Delta {f_{{{t - 1}}}}) + {k_{\rm{d}}}\Delta {f_{{t}}} {\text{，}} \\ \Delta {f_{{t}}} = {F_{\rm{n}}} - {F_{\rm{d}}} {\text{，}} \\ {u_t} = {u_{t - 1}} + \Delta {u_{{t}}} {\text{，}} \\ \Delta {x_t} = g(\Delta {u_t}). \\ \end{gathered} \right\}$

式中： $\Delta {u_{{t}}}$为t时刻增加的输出量，实验中输出量为偏置电压； ${u_t}$为t时刻的输出量；k_p和k_d为比例参数； $\Delta {f_{{t}}}$为t时刻法向力与期望力之差； $\Delta {x_t}$为t时刻机器人位移改变量，与输出量有关，用g表示，当机器人系统参数不变时，两者的关系是确定不变的.

在机器人与环境接触过程中，当k_p和k_d不协调时，机器人与环境会产生剧烈的碰撞或者过度的挤压，造成接触力突然过大或过小，有可能损坏机器人和配套设施. 需要限制比例参数的范围，将控制器设计为

(5) $\Delta {u_{{t}}} = ({k_{\rm{p}}} + {k_{{\rm{\Delta p}}}})(\Delta {f_{{t}}} - \Delta {f_{{{t - 1}}}}) + {k_{\rm{d}}}\Delta {f_{{t}}}.$

式中：k_p和k_d设置为定值，保证机器人和环境始终处于接触状态，同时接触力不会过大或过小；k_Δp为偏移参数，用于调节机器人的位置.

k_Δp的优化过程是比较繁琐的，强化学习可以很好地解决这种序贯决策问题，它可以根据机器人与环境接触的状态调整参数以达到最优奖励值，从而得到最优策略. 在有模型强化学习中，基于概率动力学模型的PILCO算法在实际应用中迭代次数较小，数据利用效率较高^[24]，比较适用于机器人恒力跟踪场景. 在该场景中，可以将PILCO算法分为概率动力学模型学习、强化学习最优策略的学习和恒力跟踪实验的策略应用3个部分.

在经典的控制中，一般需要使用函数逼近方法来获得模型. 在选择高次多项式或非参数模型的情况下，当参数过多或模型较复杂时会导致过度拟合，模型会被覆盖，开始产生拟合噪声^[25]. PILCO的模型是概率动力学模型，在模型拟合时将不确定性纳入考虑中，即在模型拟合时将模型的不确定性考虑成与模型不相关的噪声，由于考虑了模型的不确性，与传统的函数逼近方法相比，得到的模型偏差较小.

假设动态系统t时刻的状态和t−1时刻的状态、输出相关，模型为

(6) ${{ x}_{{t}}} = f({{ x}_{{{t - 1}}}},{{ u}_{{{t - 1}}}}).$

在标准的高斯回归模型中，假设数据模型输入值为数组，其中包括状态值和参数输入值，即 ${{X}} = [ \left( {{x_1},{u_1}} \right),\left( {{x_2},{u_{\rm{2}}}} \right), \cdots ,({x_{{n}}},{u_{{n}}})\}]$，模型训练目标值为 ${{Y}} \!=[{y_1},{y_2}, \cdots ,{y_{{n}}}]$，其中 ${y_{{t}}}{\rm{ }} \!=\! {x_{{t}}} - {x_{{{t - 1}}}} + {\xi _{{t}}},$ ${x_{{t}}} \in {{X}},$ ${y_{{t}}} \in {{Y}},$ $\xi $为独立的高斯噪声，即 $\xi$ 满足： ${\rm N}(0,{\sigma_\xi^2}) $分布. 为了得到训练目标的值，根据模型在t−1时刻的值对t时刻进行一步的预测，概率满足高斯分布：

(7) $\left. \begin{split} &p({x_{{t}}}|{x_{{{t - 1}}}},{u_{{{t - 1}}}}) = N\left({x_{{t}}}|{\mu _{{t}}},{\sigma ^2_{\xi{t}}}\right),\\ &{\mu _t} = {x_{{{t - 1}}}} + {E}[{y_{{t}}}],\\ &{\sigma ^2_{\xi{t}}} = {{\mathop{\rm var}} }[{y_{{t}}}]. \end{split} \right\}$

令 $\tilde x = ({x,u})$，内核函数选择平方指数协方差：

(8) $k({{ x}_{\rm{p}}},{{ x}_{\rm{q}}}) = {\alpha ^2}\exp \;\left[ - \frac{1}{2}{({{ x}_{\rm{p}}} - {{ x}_{\rm{q}}})^{\rm T}}{{{\varLambda }}^{ - 1}}({{ x}_{\rm{p}}} - {{ x}_{\rm{q}}})\right]{\text{.}} $

式中：x_p、x_q为属于 $\tilde x$的任意元素； $\alpha $为函数f信号方差； ${\varLambda }$为特征长度l_i的平方对应的对角矩阵，即 ${{\varLambda }} = {\rm diag} \;([{l_1}^2, \cdots ,{l_{{n}}}^2])$. 假设超参数为θ_p，包括特征长度l₁,···,l_n，信号方差α²和噪声方差 ${\sigma _{\xi{}}}$，在θ_p和输入X下Y概率的对数的似然边际为

(9) $ \begin{split} \log p({{Y}}|{{X}},{{\theta }_{\rm p}}) =& - \frac{1}{2}{{{Y}}^{\rm T}}{({{{K}}_\theta } + {\delta _\varepsilon }^2{{I}})^{ - 1}}{{Y}} - \\&\frac{1}{2}\log |{{{K}}_\theta } + {\delta _\varepsilon }^2{{I}}{|^{ - 1}} - \frac{D}{2}\log\; (2{\text π} ). \end{split} $

式中：D为输入空间的维数； ${{{K}}_{\theta {\rm p}}}$为关于θ_p的核函数；I为单位矩阵. 超参数值可以通过Ⅱ型超参数最大似然估计(ML-Ⅱ)获得^[25]：

(10) $ {\hat {\bf\theta }} \in \arg \mathop {\max }\limits_{ \theta_{\rm p}} \log p({{Y}}|{{X}},{{\theta }}_{\rm p}). $

在得到超参数后，对输入参数进行高斯过程回归预测：当输入测试参数为x*时，根据先验概率得到后验概率分布：

(11) $P\left( {{{ y}^*}|{{ x}^*}} \right) = N({{ k}_*}^{\rm T}{\left( {{{K}} + {\sigma _\varepsilon }^2{{I}}} \right)^{ - 1}}{ Y},{{ k}_{**}} - {{ k}_*}^{\rm T}{({{K}} + {\sigma _\varepsilon }^2{{I}})^{ - 1}}{{ k}_*}).$

式中：y*为预测输出值；k_*=k(X, x*), k_**=k(x*, x*).

强化学习是以马尔科夫决策过程为理论基础的. 马尔科夫决策过程由一个元组构成M=(S，A，P，r，ρ₀，γ)，其中S为系统状态集合；A为动作集合；P为状态的转移概率，概率动力学模型已经求出；r为回报函数；ρ₀为初始状态；γ为折扣因子，一般取值为[0, 1]^[26]. 强化学习主体(agent)在与环境交互过程中，通过选择不同策略，即根据不同的状态选择合适的动作得到不同的奖励，最终得到最优解. 强化学习算法的种类有很多，大体可以分为基于策略的强化学习和基于价值的强化学习. 基于策略的强化学习比较适合动作为连续集的场景，适合机器人跟踪场景. 基于策略的强化学习结构如图5所示，主体的策略π更新过程是通过改变策略参数φ实现的，φ的更新的步长与目标函数（该场景中设置为强化学习期望回报）的梯度有关. 强化学习期望回报与奖励函数和状态概率分布有关，奖励函数用于衡量机器人的实际状态到期望状态的距离. 在机器人跟踪领域中，主体目标是寻找到最优的策略，最优策略是对长期的策略进行评估，可以利用3.1节中概率动力模型中学习到的模型，对长期的回报进行近似的推断，在得到强化学习中的目标函数后，根据目标函数梯度改变输入的策略，直至更新的策略满足最优解.

强化学习的期望回报为

(12) $\left. \begin{aligned} {V^\pi } =& \sum\limits_{t = 1}^T {E[c({{ x}_{{t}}})|{\text π}]} ,P({{ x}_0}) = N({\mu _0},{\sigma _0^2}); \\ {V^\pi } = &\sum\limits_{t = 1}^T {\int {c({{ x}_{{t}}})p({{ x}_{{t}}})} } {\rm d}{{ x}_{{t}}} . \\ \end{aligned} \right\}$

式中： ${V^\pi }$为在策略 $\pi $下的期望回报；初始状态 ${{ x}_0}$满足高斯分布； $p({{ x}_{{t}}})$为在状态x_t下的状态的概率分布； $c({{ x}_t})$为奖励函数，

(13) $c({ x}) = 1 - \exp\; \left.\left( - \frac{1}{2}\left\| {{ x} - {{ x}_{{\rm{target}}}}} \right\|_W^2\right/{\delta _{\rm{c}}}^2\right) {\text{，}}$

其中，||x−x_target||_W为x和x_target之间的马氏距离，W为加权因子，x_target为目标状态， ${\delta _{\rm{c}}}$为缩放因子. 在设定目标状态后，奖励函数比较容易得到，而式(12)中x_t时间的状态分布 $p({{ x}_{{t}}})$是未知的，需要对 $p({{ x}_t})$进行推理：在t−1时刻，首先由 $p({{ x}_{{{t}} - 1}})$的状态分布推理出动作分布 $p({u_{{{t - 1}}}})$，其次将 $p({{ x}_{{{t}} - 1}})$与 $p({u_{{{{t} - 1}}}})$结合得到联合分布 $p({{ x}_{{{{t} - 1}}}},{u_{{{{t} - 1}}}})$， $p({{ x}_{{{{t} - 1}}}},{u_{{{{t} - 1}}}})$可以得到状态改变量 $p({ \varDelta} _{t})$，其中 ${ { \varDelta} _{{t}}} = {{ x}_{{t}}} - {{ x}_{{{{t} - 1}}}}$，最后得到 $p({{ x}_t})$. $p({{ x}_{{t}}})$的期望和方差^[25]为

(14) $\left. \begin{aligned} {\mu _{{t}}} =& {\mu _{{{{t} - 1}}}} + {\mu _{\bf \Delta} } , \\ {\sigma _{{t}}} =& {\sigma _{{{t - 1}}}} + {\sigma _\Delta } + \operatorname{cov}\; [{{ x}_{{{t - 1}}}},{{ \varDelta} _{{t}}}] + \operatorname{cov}\; [{{ \varDelta} _{{t}}},{{ x}_{{{t - 1}}}}]. \end{aligned} \right\}$

在策略和状态关系通过策略参数线性化或者非线性化后，最优策略由最优策略参数确定，当策略满足以下条件时即为最优策略：

(15) ${{\text π} ^*} \in \arg \min {V^{\text π}}({{ x}_0}).$

为了得到最优策略参数，需要求解期望回报相对于策略参数的梯度，令： ${{ \varepsilon} _t}{\bf =} {E}[c({{ x}_t})]$，则

(16) $\frac{{{\rm d}{\varepsilon _{{t}}}}}{{{\rm d}\varphi }} = \frac{{{\rm d}{\varepsilon _{{t}}}}}{{{\rm d}p({x_{{t}}})}}\frac{{{\rm d}p({x_{{t}}})}}{{{\rm d}\varphi }} = \frac{{\partial {\varepsilon _{{t}}}}}{{\partial {\mu _{{t}}}}}\frac{{{\rm d}{\mu _{{t}}}}}{{{\rm d}\varphi }} + \frac{{\partial {\varepsilon _{{t}}}}}{{\partial {\sigma_{\xi{t}}}}}\frac{{{\rm d}{\sigma_{\xi{t}}}}}{{{\rm d}\varphi }}.$

式中： $\displaystyle{{\partial {\varepsilon _{{t}}}}}/{{\partial {\mu _{\rm{t}}}}}$， $\displaystyle{{\partial {\varepsilon _t}}}/{{\partial{\sigma_{\xi{t}}}}}$在给定 $c({x_{{t}}})$后可以解析得到； $\displaystyle{{{\rm d}{\mu _t}}}/{{{\rm d}\varphi }}$、 $\displaystyle{{{\rm d}{\sigma_{\xi{t}}}}}/{{{\rm d}\varphi }}$可以通过t−1时刻 $\displaystyle{{{\rm d}{\mu _{t - 1}}}}/{{{\rm d}\varphi }}$、 $\displaystyle{{{\rm d}{\sigma_{\xi(t-1)}}}}/{{{\rm d}\varphi }}$的值递推出来：

(17) $\left. \begin{gathered} \frac{{{\rm d}{\mu _{{t}}}}}{{{\rm d}\varphi }} = \frac{{\partial {\mu _{{t}}}}}{{\partial {\mu _{{{t - 1}}}}}}\frac{{{\rm d}{\mu _{{{t - 1}}}}}}{{{\rm d}\varphi }} + \frac{{\partial {\mu _t}}}{{\partial {\sigma_{\xi(t-1)}}}}\frac{{{\rm d}{\sigma_{\xi(t-1)}}}}{{{\rm d}\varphi }} + \frac{{\partial {\mu _{{t}}}}}{{\partial \varphi }}, \\ \frac{{{\rm d}{\sigma_{\xi{t}}}}}{{{\rm d}\varphi }} = \frac{{\partial {\sigma_{\xi{t}}}}}{{\partial {\mu _{{{t - 1}}}}}}\frac{{{\rm d}{\mu _{{{t - 1}}}}}}{{{\rm d}\varphi }} + \frac{{\partial {\sigma_{\xi{t}}}}}{{\partial{\sigma_{\xi(t-1)}}}}\frac{{{\rm d}{\sigma_{\xi(t-1)}}}}{{{\rm d}\varphi }} + \frac{{\partial {\sigma_{\xi{t}}}}}{{\partial \varphi }} . \end{gathered} \right\}$

在得到所有梯度后，利用共轭梯度法可得策略参数更新步长，详细推导见文献[25].通过不停地迭代，直至策略参数满足式(15)最优策略条件.

由于机器人与环境接触过程具有不确定性，在相同的偏移参数k_Δp下，相同时刻得到的状态x会不同，直接利用某时刻的状态进行计算会导致采集数据误差较大. 初始状态下，将接触过程平均分成30个时间段，假设在不同时间段内策略对各段接触的影响是相互独立的，提取每个时间段内力信号，求出F_n与F_d之差Δf，Δf的平均力 $f_{{\rm{ave}}}^{}$、Δf的方差 $f_{\operatorname{var} }^{}$和每段时间周期内截止时刻对应的角度 ${\theta _{}}$，构成机器人的状态值. 在更新策略 $\pi $后，得到30组新的偏移参数(k_Δp=[k¹_Δp，k²_Δp，k³_{Δp，···，}k³⁰_Δp])和角度，下一次实验时以得到的角度为时间节点，对偏移参数进行更新.

机器人的状态值和偏移参数之间的关系需要非线性化处理，由于RBF网路有较好的输入和输出的映射功能，可以任意精度地逼近任意的非线性函数，同时具有唯一的最佳逼近特征，学习和收敛的速度较快. 在实验中，选择高斯内核的RBF网络，将偏移参数和状态x关系非线性化为

(18) $\left. \begin{gathered} \pi ({ x},{ \varphi} ) = \sum\nolimits_{i = 1}^n {{w_{{i}}}} {\phi _{{i}}}({ x}), \\ {\phi _{{i}}}({ x}) = \exp\; \left( - \frac{1}{2}{({ x} - {{ \mu} _{{{\rm r}i}}})^{\rm T}}{{{\varGamma }}^{ - 1}}({ x} - {{ \mu} _{{{\rm r}i}}})\right). \end{gathered} \right\}$

式中：n为节点； $w$为权重； ${{ \varGamma }}$为权重矩阵； ${ \mu }$为高斯基函数的中心位置； $\varphi $为策略参数，其中包括 $\left\{ {{{ w}_{{i}}},{{\varGamma }},{{ \mu} _{\rm{r}}}} \right\}$. 算法结构如图6所示，由概率动力学模型模块和强化学习模块构成. 概率动力学模型模块输入由第j次实验的状态（包括平均力、方差和角度）和输入偏移参数构成，输出为j+1次实验的状态，强化学习模块接收第j+1次实验状态和角度后由RBF网络输出为j+1次实验的偏移参数.

为了方便计算，需要对RBF的输入值作归一化处理，使输入值为[−1, 1]. 实验中，期望目标值x_target由平均力和力的方差构成，即x_target=[f^ave_target,f ^var_target]，其中平均力和力的方差的目标都为0.

实验装置如图7所示，由机器人、控制器、倍福模块、六维力传感器和上位机构成. 机器人采用安川MH24；六维力传感器选用德国的ME-FKD40；力信号通过倍福模块采集并传送给上位机，利用上位机运行力控制器(式(5))中的算法，将输出信号传送给倍福模块，并由倍福模块转化为电压信号；控制器内置软件接收倍福−10~10 V的模拟信号产生位移偏移，位移偏移的方向与模拟信号电压符号相关，位移偏移与电压绝对值成正比；控制器可以将机器人的关节角度传送给上位机. 实验时，上位机采集六维力传感器中F_x和F_y的力信号，离线运行算法后，得到偏移参数后进行下一次实验.

实验使用的曲面工件如图8所示. 曲面工件的轮廓通过三维软件中样条曲线命令拉伸而成，机器人的初始轨迹如图8的虚线所示，初始轨迹的运行程序通过手动示教设定，其中A点到B点的初始轨迹为直线. 当机器人到达A点时，开始采集力信号，内置软件开始运行；机器人始终以1 mm/s的速度沿V_x方向移动，在V_y方向开始偏移，位移偏移的大小是根据显式力控制算法来实现的，其中F_d为30 N；k_p=0.004，k_d=0.000 02；k_Δp为[0，0.02]. 当机器人达到B点时，内置软件停止运行，即偏移停止，离开B点后当力信号为0时，力信号采集结束.

在k_Δp的范围内随机选取30个参数作为初始随机策略，同时参数分布符合高斯分布，实验流程如图9所示. 在执行策略后，将得到的实验数据归一化处理，方便后续处理，利用采集的数据求出模型的超参数，得到概率动力模型. 在该模型下，不断更新策略的参数改变策略，求出最优策略，继续实验直至接触力与期望力的误差小于3 N.

在初始参数下，得到的力信号如图10(a)所示，期望值30 N内层的虚线范围为±2 N，外层的虚线范围为±3 N. 在随机策略下，接触过程中前期波动较大，后期波动小. 机器人末端执行器在与曲面工件轮廓前期接触时，x方向挤压力波动较大，因此F_n波动较大；机器人末端执行器在与曲面工件轮廓后期接触时，x方向挤压力相对较小同时波动不大，因此F_n的变化相对前期较小，所以相对平稳.

将第1次的实验结果进行离线处理，得到30组状态和偏移参数的值，输入PILCO算法，得到最优策略；重复实验8次，得到的法向力信号如图10、11所示. 在7次迭代时，前期的法向力稳定且收敛. 接触过程的后期由于相对稳定，收敛速度相对较快，在迭代第8次后，法向力与期望值的误差绝对值绝大部分小于2 N，到达期望的效果.

迭代时力误差绝对值的最大值|Δf |_max、力误差绝对值的平均值 $\left|\Delta \bar f \right| $和力误差的标准差Δf_s在迭代过程中的变化如表1所示. 力误差绝对值的最大值在第1次迭代后基本没有什么变化，由于在机器人和曲面接触的过程中冲击阶段时间较短，在冲击阶段的调节的效果比较有限，力误差绝对值的最大值变化不多，而跟踪阶段力误差绝对值的平均值和力误差的标准差的总体趋势是下降的.

在迭代过程中，PILCO算法对每一个时间段内的奖励函数的值(见式(13))进行预测，奖励函数的值能够反映机器人当前时间段的状态与期望目标状态的偏差. 如图12所示，预测奖励函数的值为实线，上下的范围为误差范围，实际的奖励函数的值为带星号的实线，当奖励函数的值接近0时，接触状态接近期望的状态. 第1次迭代过程时，由于采集的数据较少，PILCO算法对机器人与环境接触的状态进行推理时，预测的状态和实际的状态有一定差距；当采集的数据增加后，在第8次迭代中，实际的奖励函数的值在预测的奖励函数的值的范围之内，预测的状态真实地反映了实际中的状态，证明了策略的准确性.

将PILCO算法和机器人恒力跟踪的模糊迭代算法进行比较，其中恒力跟踪的模糊迭代算法应用的原理是根据机器人与环境的接触力的大小离线迭代机器人的轨迹，同时在迭代过程中利用模糊算法补偿机器人的轨迹. 2种算法的效果对比如图13所示. PILCO算法的迭代次数为8次，模糊迭代算法迭代了18次；PILCO算法中法向力与期望力的误差大部分为±2 N，模糊迭代算法中为±3 N；PILCO算法在冲击阶段的效果较差，模糊迭代算法可以对每个过程都进行调节，因此在冲击阶段的控制效果较好. 当机器人与环境的接触力第1次达到期望力至机器人离开环境这段时间内，PILCO算法和模糊迭代算法的力误差绝对值的最大值、力误差绝对值的平均值和力误差的标准差的对比如表2所示. PILCO算法和模糊迭代算法的力误差标准差相差不大，PILCO算法的力误差绝对值平均值比模糊迭代算法减小了29%.

机器人末端执行器对曲面工件轮廓进行恒力跟踪实验时，对末端执行器的受力进行分析，使用机器人显式力控制对机器人的位移进行调整. 为了得到恒定的接触力，通过有模型强化学习PILCO算法，对接触过程中的状态和输入参数的关系进行学习，自主地得到输出参数. 在实验中，将强化学习的输入状态改为一段时间内的状态平均值，以减少接触状态下信号的干扰. 实验结果显示，接触力的误差大部分小于2 N. 将PILCO算法和模糊迭代算法的效果进行对比可知，PILCO算法的迭代次数较少，跟踪力的稳定性效果更好，力误差绝对值的平均值比模糊迭代算法减少了29%.