双足机器人具有人类的外观特征以及步行能力，因而能适应人类的生活与工作环境，具有广泛的应用前景^[1]. 在实际应用中，由于外部环境的不确定性和时变性，双足机器人不可避免地受到各种外部扰动的影响^[2-4]. 从已公开的一些仿真以及实验结果来看^[5]，目前双足机器人的抗扰动能力依然较弱，鲁棒性仍然较差. 传统的基于零力矩点（zero moment point，ZMP）的准静态步行控制虽然有较好的实用性，但控制方式过于保守，要求机器人的ZMP点时刻落在支撑多边形内，与人类的步行方式不符^[6].

近年来，国内外学者开始关注双足机器人动态步行的鲁棒控制问题^[7-9]. 相比准静态步行，动态步行^[10]不关注机器人的局部稳定性，通常基于庞卡莱映射方法，从动力学模型及其解的角度，分析步态的全局稳定性. 针对外部作用力下的动态步行双足机器人，Veer等^[11]假设作用力信息已知，提出混杂零动力学控制方法，并给出稳定动态步行的解析条件. Wang等^[12]考虑机器人在双足支撑期受到外部作用力的情形，将机器人的抗扰动控制问题转化为动态平衡优化问题，并提出基于递归神经网络的运动/力混合控制方法. 为了实现非平整路面下的动态稳定步行，Hamed等^[13]提出连续控制器设计方法，并基于庞卡莱分析，将控制器的参数设计问题转化为双线性矩阵不等式（bilinear matrix inequality，BMI）优化问题. Sun等^[14]将非平整路面步行控制问题统一转化为虚拟斜坡步行控制问题，并基于倒立摆模型，提出神经网络控制方法. 总体来看，国内外学者关于动态步行鲁棒控制的研究主要考虑非平整路面^[14-15]、限时不变作用力^[16-17]以及周期作用力^[11,18]等外部扰动情形. 在实际中，机器人的外部扰动，如风力以及人为牵引力等，可能具有时变性和随机性. 在此情形下，对机器人外部扰动的幅值的估计难免会存在误差. 目前，国内外学者较少关注不确定性扰动下双足机器人动态步行的鲁棒控制问题. 研究^[19]表明，当存在不确定性扰动时，动力学系统的状态轨迹会密集地分布于相空间，而不会收敛到某一轨道，这使得其长期动力学行为的描述没有意义，从而难以利用常规的庞卡莱映射方法对其进行稳定性分析.

针对以上问题，本研究将随机振动系统的特定庞卡莱映射方法拓展到不确定性扰动下双足机器人的稳定性分析. 为了实现不确定性扰动下双足机器人的鲁棒控制，提出自适滑模控制器. 考虑到双足机器人在实际应用中常会遭遇非平整路面，进一步将该控制器拓展到非平整路面的鲁棒控制.

如图1所示，三维 (3-D)双足机器人由躯干、2条大腿以及2条小腿构成. 其中，小腿长 ${L_1} = 0.275$ m，质量 ${m_1} = 0.875$ kg；大腿长 ${L_2} = 0.275$ m，质量 ${m_2} = 0.875$ kg；躯干宽 $W = 0.150$ m，高度 ${L_3} = 0.050$ m，质量 ${m_3} = 5.500$ kg. 该机器人包含9个自由度，每个自由度都可被独立驱动. 其中， ${q_0}$、 ${q_1}$、 ${q_2}$为支撑腿踝关节自由度， ${q_3}$、 ${q_8}$为膝关节自由度， ${q_4}$、 ${q_5}$为支撑腿髋关节自由度， ${q_6}$、 ${q_7}$为摆动腿髋关节自由度. 根据机器人与地面的接触形式，其运动状态可分为单足支撑相以及双足支撑相. 为了简化问题，假设偏航自由度 ${q_0}$受到摩擦力的影响，在整个单足支撑相保持静止，并假设机器人的双足支撑相为瞬时刚性碰撞，在碰撞过程中位姿不突变，也不发生反弹现象.

在单足支撑相阶段，假设机器人受到不确定性扰动的影响，如图2所示. 为了便于分析，假设不确定性扰动水平作用于质心. 根据拉格朗日建模法，机器人在单足支撑相的动力学模型的表达式为

(1) ${{{D}}_v}({{q}}){\ddot{{q}}} + {{H}_v}({{q}},{\dot{{q}}}) = {{u}} + {{{E}}^{\rm T}}{{{F}}_{\rm ex}}.$

式中： ${{q}} = {[{q_1},\cdots,{q_8}]^{\rm T}}$为坐标向量； $v \in \{ 1,2\} $， $v = 1$表示左脚支撑， $v = 2$表示右脚支撑； ${{{D}}_v}({{q}})$为惯性矩阵； ${{{H}}_v}({{q}},{\dot{ q}})$为科氏力和重力项； ${{u}}$为输入力矩； ${{E}} = {{\partial {{{P}}_{\rm cm}}({{q}})} / \partial }{{q}}$为雅可比矩阵， ${{{P}}_{\rm cm}}({{q}})$为质心位置； ${{{F}}_{\rm ex}} = {{f}}\xi (t)$为不确定性作用力，此处假设为有界随机激励形式， ${{f}}$为常数向量， $\xi (t)$为时间段 $[0,T]$内的有界随机过程， $T$为步行周期.

根据文献[20]、[21]， $\xi (t)$的每个样本可估计为

(2) $\eta (t) = \sum\limits_{i=1}^N {{A_i}\cos\; ({w_i}t + {\phi _i})}. $

式中： ${w_i}$为非负的独立分布于区间 $[{w_{\rm l}},{w_{\rm r}}]$的随机变量； ${\phi _i}$为均匀分布于 $[0,2\text{π} )$的随机相位； $N$为充分大的正整数； ${A_i} = \sqrt {2{S_i}\Delta {w_i}} $， $\Delta {w_i} = ({w_{\rm r}} - {w_{\rm l}})/N$为频率增量， ${S_i}$为功率谱密度函数，表达式为

(3) ${S_i} = \frac{1}{{2\pi }}\left(\frac{{{\sigma ^2}}}{{4{{({w_i} - \varOmega )}^2} + {\sigma ^4}}} + \frac{{{\sigma ^2}}}{{4{{({w_i} + \varOmega )}^2} + {\sigma ^4}}}\right).$

其中， $\varOmega $为中心频率， $\sigma $为谱宽参数. 由于随机过程 $\xi (t)$的每个样本均为 $N$项谐波函数的和，一旦生成即可视作确定性函数.

在单足支撑相阶段，定义状态向量 ${{{x}}_v} = {[{{q}},{\dot{ q}}]^{\rm T}}$，则式（1）可改写为如下状态空间形式：

(4) ${{\dot{{x}}}_v} = {{{f}}_v}({{{x}}_v}) + {{{g}}_v}({{{x}}_v}){{u}} + {{{f}}_{\rm ex}}.$

式中： ${{{f}}_v}({{{x}}_v}) = \left[ {\begin{array}{*{20}{c}} {{\dot{ q}}} \\ { - {{D}}_v^{ - 1}{{{H}}_v}} \end{array}} \right]$， ${{{g}}_v}({{{x}}_v}) = \left[ {\begin{array}{*{20}{c}} {{{\bf 0}_{8 \times 8}}} \\ { - {{D}}_v^{ - 1}} \end{array}} \right]$， ${{{f}}_{\rm ex}} = $ $\left[ {\begin{array}{*{20}{c}} {{{\bf 0}_{8 \times 1}}} \\ {{{D}}_v^{ - 1}{{{E}}^{\rm T}}{{{F}}_{\rm ex}}} \end{array}} \right]$.

在双足支撑相阶段，选取广义坐标向量 ${{{q}}_{\rm{e}}} = $ ${[{x_{\rm st}},{y_{\rm st}},{z_{\rm st}},{{{q}}_0},{{{q}}_1},\cdots,{{{q}}_8}]^{\rm T}}$. 根据前述双足支撑相的假设，可以求得机器人在双足支撑相的碰撞方程为

(5) $\left[ {\begin{array}{*{20}{c}} {{\dot{ q}}_{\rm e}^ + } \\ {{{{F}}_{v,{\rm{sw}}}}} \end{array}} \right] = \mathop {\left[ {\begin{array}{*{20}{c}} {{{{D}}_{v,{\rm e}}}}&{ - {{E}}_{v,{\rm{sw}}}^{\rm T}} \\ {{{{E}}_{v,{\rm{sw}}}}}&{\bf 0} \end{array}} \right]}\nolimits^{ - 1} \left[ {\begin{array}{*{20}{c}} {{{{D}}_{v,{\rm e}}}{\dot{ q}}_{\rm e}^ - } \\ {\bf 0} \end{array}} \right].$

式中： ${\dot{ q}}_{\rm e}^ - $、 ${\dot{ q}}_{\rm e}^ + $分别为碰撞前后的广义坐标速度， ${{{D}}_{v,{\rm{e}}}}$为双足支撑相的惯性矩阵， ${{{E}}_{v,{\rm{sw}}}}$为摆动脚位置 ${{{P}}_{\rm{sw}}} = {[{x_{\rm{sw}}},{y_{\rm{sw}}},{z_{\rm{sw}}}]^{\rm T}}$对广义坐标向量的雅可比矩阵， ${{{F}}_{v,{\rm{sw}}}}$为摆动腿所受的脉冲碰撞力.

在碰撞结束后，支撑腿发生切换，关节坐标也相应地进行交换，交换顺序为 $[{{{q}}_{1,{\rm{sw}}}},{{{q}}_{2,{\rm{sw}}}},{{{q}}_8},\cdots,{{{q}}_3}]$ $ \to $ $[{{{q}}_1},\cdots,{{{q}}_8}]$，如图1所示. 结合碰撞方程（式（5））以及坐标交换过程，机器人双足支撑相的整体动力学模型表达式^[17]为

(6) ${{x}}_{{\rm{Nex}}(v)}^ + = {{{\varDelta}} _v}({{x}}_v^ - ).$

式中： ${{x}}_v^ - = {[{{q}}_v^ - ,{\dot{ q}}_v^ - ]^{\rm T}}$为上一步的末状态， ${{x}}_{{\rm{Nex}}(v)}^ + = $ $ {[{{q}}_{{\rm{Nex}}(v)}^ + ,{\dot{ q}}_{{\rm{Nex}}(v)}^ + ]^{\rm T}}$为下一步的初始状态， $ {{{\varDelta}} _v}(\cdot) $为机器人双足支撑相的离散映射，标签函数 ${\rm{Nex}}(v)$的表达式为

(7) ${\rm{Nex}}(v) = \left\{ {\begin{array}{*{20}{l}} {1,\;\;\;\;\;\;v = 2}; \\ {2,\;\;\;\;\;\;v = 1} . \end{array}} \right.$

结合式（4）、（6），机器人在不确定性扰动下的动力学模型可表示为

(8) $\sum = \left\{ {\begin{array}{*{20}{l}} {{{{\dot{ x}}}_1} = {{{f}}_1}({{{x}}_1}) + {{{g}}_1}({{{x}}_1}){{u}} + {{{f}}_{\rm ex}},\;\;\;\;\;{{{x}}_1} \notin {{{S}}_1}}; \\ {{{x}}_2^ + = {{{\varDelta}} _1}({{x}}_1^ - ),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{{x}}_1} \in {{{S}}_1}}; \\ {{{{\dot{ x}}}_2} = {{{f}}_2}({{{x}}_2}) + {{{g}}_2}({{{x}}_2}){{u}} + {{{f}}_{\rm ex}},\;\;\;\;{{{x}}_2} \notin {{{S}}_2}}; \\ {{{x}}_1^ + = {{{\varDelta}} _2}({{x}}_2^ - ),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{{x}}_2} \in {{{S}}_2}}. \end{array}} \right.$

式中： ${{{S}}_v} = \{ {{{x}}_v}|{z_{\rm{sw}}}({{{x}}_v}) = 0,\;{\dot z_{\rm{sw}}}({{{x}}_v}) < 0\} $ 为切换曲面，此时机器人与地面发生碰撞，单足支撑相结束.

在不确定性扰动作用下，动力学系统的状态轨迹会密集地分布于相空间，而不会收敛到某一轨道，这使得其长期动力学行为的描述没有意义^[19]，从而难以利用常规的庞卡莱映射方法对此类系统进行稳定性分析.

可将特定庞卡莱映射方法^[19,21]用于分析随机振动系统在有限时间段 $[0,{T_0}]$的稳定性：若在 $t = 0$时从相空间某有界集合 $X$出发的任意状态轨迹能在 $t = {T_0}$时返回集合 $X$，则集合 $X$为时间段 $[0,{T_0}]$内在Lyapunov意义下稳定的轨道集合. 基于特定庞卡莱映射方法，随机振动系统的稳定性分析被转化为确定性周期系统的稳定性分析. 本研究将特定庞卡莱映射方法拓展到不确定性扰动下双足机器人在有限步数 $n$步内的稳定性分析. 为了便于分析，假设 $n$为偶数，且 $n \geqslant 2$.

记随机过程 $\xi (t)$在第 $k$步的样本为 ${\eta _k}$，则机器人在第 $k$步的动力学方程可以表示为

(9) ${{{D}}_v}({{q}}){\ddot{ q}} + {{{H}}_v}({{q}},{\dot{ q}}) = {{u}} + {{{E}}^{\rm T}}{{f}}{\eta _k}.$

为了应用特定庞卡莱映射方法，对机器人的不确定性扰动作周期化处理. 对于任意正整数 $l$，令

(10) $ {\eta _k}(t) = {\eta _{(k + ln)}}(t);\;\;t \in [0,T], \;1 \leqslant k \leqslant n. $

这种周期化处理是合理的，因为研究只关注机器人在 $n$步内的稳定性. 基于式（9）、（10），机器人不确定性扰动系统被转化为如下周期扰动系统：

(11) $\sum = \left\{ {\begin{aligned} & {{{{\dot{ x}}}_1} = {{{f}}_1}({{{x}}_1}) + {{{g}}_1}({{{x}}_1}){{u}} + {{f}}{\eta _1},\;\;\;\;{{{x}}_1} \notin {{{S}}_{\eta ,1}}} ;\\ & {{{x}}_2^ + = {{{\varDelta}} _1}({{x}}_1^ - ),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{{x}}_1} \in {{{S}}_{\eta ,1}}}; \\ & {{{{\dot{ x}}}_2} = {{{f}}_2}({{{x}}_2}) + {{{g}}_2}({{{x}}_2}){{u}} + {{f}}{\eta _2},\;\;\;{{{x}}_2} \notin {{{S}}_{\eta ,2}}} ;\\ & {{{x}}_1^ + = {{{\varDelta}} _2}({{x}}_2^ - ),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{{x}}_2} \in {{{S}}_{\eta ,2}}}; \\ & \;\;\;\;\;\;\;\;\;\;\;\qquad\qquad\qquad \vdots \\ & {{{{\dot{ x}}}_2} = {{{f}}_2}({{{x}}_2}) + {{{g}}_2}({{{x}}_2}){{u}} + {{f}}{\eta _n},\;\;\;{{{x}}_2} \notin {{{S}}_{\eta ,n}}} ;\\ & {{{x}}_1^ + = {{{\varDelta}} _2}({{x}}_2^ - ),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{{x}}_2} \in {{{S}}_{\eta ,n}}} . \end{aligned}} \right.$

式中： ${{{S}}_{\eta ,k}} = \{ {{{x}}_v}|{z_{\rm{sw}}}({{{x}}_v},{\eta _k}) = 0,\;{\dot z_{\rm{sw}}}({{{x}}_v},{\eta _k}) < 0\} $为扰动 ${{f}}{\eta _k}$所对应的切换曲面， $k = 1, \cdots ,n$.

系统（式（11））可以视为具有 $n$个连续运动相的周期混杂系统，因此可采用确定性方法对其进行稳定性分析. 定义系统（式（11））第 $k$个连续运动相对应的切换时间映射 ${{{T}}_k}:{{{x}}_v} \to {\bf{R}}$为

(12) ${{{T}}_k}({{{x}}_0}) = \inf\; \{ t \geqslant 0|{{{\varPhi}} _k}(t,{{{x}}_0}) \in {{{S}}_{\eta ,k}}\} .$

式中： ${{{\varPhi }}_k}(t,{{{x}}_0})$为式（9）的解，并满足 ${{{\varPhi }}_k}(0,{{{x}}_0}) = {{{x}}_0}$.

定义其第 $k$个连续运动相对应的状态映射 ${{{P}}_k}:{{{S}}_{\eta ,(k - 1)}} \to {{{S}}_{\eta ,k}}$为

(13) ${{{P}}_k}({{x}}) = {{{\varPhi}} _k}({{{T}}_k}({{{\varDelta}} _{{\rm{Nex}}(v)}}({{x}})\;),{{{\varDelta}} _{{\rm{Nex}}(v)}}({{x}}))\;.$

从而，特定庞卡莱映射可被定义为 ${{P}}:{{{S}}_{\eta ,n}}$ $ \to {{{S}}_{\eta ,n}}$，表达式为

(14) ${{P}}({{x}}) = {{{P}}_1} \circ \cdots \circ {{{P}}_n}({{x}}).$

由特定庞卡莱映射 ${{P}}$的构建可以看出，特定庞卡莱映射实为转化后的周期扰动系统（式（11））的常规庞卡莱映射，因而可基于常规庞卡莱映射方法进行处理^[19]. 记特定庞卡莱映射P的第 $l$次返回映射为 ${({{P}})^l}$. 当 $l \to \infty $时，若 ${({{P}})^l}$收敛，则机器人在有限步数 $n$步内是在Lyapunov意义下稳定的.

由上述分析可看出，特定庞卡莱映射方法是由常规庞卡莱映射方法发展而来的，但特定庞卡莱映射方法只关注系统在有限时间内的稳定性，并须先对外部扰动作周期化处理，将不确定性系统转化为确定性周期系统，然后便可基于常规庞卡莱映射方法进行处理. 若系统本身为确定性周期系统，则可直接利用常规庞卡莱映射方法进行稳定性分析.

针对不确定性扰动作用下的3-D双足机器人提出自适应滑模控制器，基于该控制器的机器人控制系统框图如图3所示. 相比以往的滑模控制器^[22-23]，该控制器无须准确估计外部扰动的幅值信息，对幅值估计误差具有一定的适应性. 自适应滑模控制器的具体设计如下：

1）符号说明. $\left\| \cdot \right\|$表示对矩阵或向量的每个元素取绝对值； ${({{a}})_i}$表示向量 ${{a}}$的第 $i$个元素；对于等维度的非负向量 ${{a}}$、 ${{b}}$，若对于任意 $i$，有 ${({{a}})_i} \leqslant {({{b}})_i}$或 ${({{a}})_i} < {({{b}})_i}$，记为 ${{a}} \leqslant {{b}}$或 ${{a}} < {{b}}$；对于等维度向量 ${{c}}$、 ${{d}}$， ${{c}}/{{d}}$表示将 ${{c}}$的每个元素对应除以 ${{b}}$的元素.

2）控制目标. 机器人各关节的参考轨迹记为 ${{{q}}_{\rm d}}$，关节控制误差 ${\tilde{ q}} = {{q}} - {{{q}}_{\rm d}}$. 定义向量

(15) ${ s} = \dot {\tilde{ q}} + \lambda {\tilde{ q}}.$

记其元 ${({{s}})_i}$为 ${s_i}$， $\lambda > 0$. 若向量 ${{s}} \equiv {\bf 0}$，则机器人的状态收敛到理想关节轨迹. 但在实际中，选取 ${{s}} \equiv {\bf 0}$作为控制目标易导致颤振问题^[22]. 为了避免该问题，本研究设计控制器，使机器人的状态轨迹收敛到某个关于常数 $\phi $的薄边界层B = $ \{ ({{ q}},{\dot{ q}})|$ $\left\| { s} \right\|$ $ < \phi { I}_{8 \times 1}\} $，即对于任意 $i$，当 $\left| {{s_i}} \right| > \phi $时，机器人状态满足如下滑模条件：

(16) $\frac{1}{2}\frac{\rm d}{{{\rm d}t}}s_i^2 = {\dot s_i}{s_i} < 0.$

此时，若定义Lyapunov函数 $V = {{{s}}^{\rm T}}{{s}}$，可知 $\dot V < 0$.

3）控制器设计. 设机器人第 $k$步的外部扰动 ${{f}}{\eta _k}$的幅值估计为 ${\hat{ F}}$， $\left\| {{ f}{\eta _k}} \right\| < {\hat{ F}}$恒成立，并定义向量 ${{\varPsi}} = $ $\left\| {{{D}}_v^{ - 1}{{{E}}^{\rm T}}} \right\|{\hat{ F}}/\phi $，则机器人第 $k$步的控制器可设计为

(17) ${{u}} = {{{D}}_v}{\ddot {{q}}_{\rm d}} + {{{H}}_v}({{q}},{\dot{ q}}) - {{{D}}_v}\lambda \dot {\tilde{ q}} - {{{D}}_v}{{Ks}}.$

式中： ${{K}} = {\rm diag} \;({k_1}, \cdots ,{k_8})$， ${k_i} = {({{\varPsi}} )_i}$.

控制器（式（17））的收敛性的简单证明如下. 首先，对式（15）求导，并结合式（9）、（17），可以得到

(18) ${\dot{ s}} = {{D}}_v^{ - 1}{{{E}}^{\rm T}}{{f}}{\eta _k} - {{Ks}}.$

令 ${{M}} = {{D}}_v^{ - 1}{{{E}}^{\rm T}}{{f}}{\eta _k}$，取 ${m_i} = {({{M}})_i}$，可得 ${{\dot s}_i} = {m_i} - {k_i}{s_i}$，从而有

(19) $\frac{1}{2}\frac{\rm d}{{{\rm d}t}}s_i^2 = {{\dot s}_i}{s_i} = {m_i}{s_i} - {k_i}s_i^2.$

此外，由于 $\left\| {{M}} \right\| \leqslant \left\| {{{D}}_v^{ - 1}{{{E}}^{\rm T}}} \right\|\left\| {{{f}}{\eta _k}} \right\| < \left\| {{{D}}_v^{ - 1}{{{E}}^{\rm T}}} \right\|{\hat{ F}}$ $ = \phi {{\varPsi}} $，有

(20) $\left| {{m_i}} \right| < \phi \left| {{k_i}} \right|.$

所以，当 $\left| {{s_i}} \right| > \phi $时，有 $\left| {{m_i}} \right| < \left| {{k_i}{s_i}} \right|$，从而有

(21) ${m_i}{s_i} < \left| {{k_i}{s_i}} \right|\left| {{s_i}} \right| = {k_i}s_i^2.$

再将其代入式（19），有 ${{\dot s}_i}{s_i} < 0$，满足滑模条件（式（16）），证明了控制系统是收敛的.

在实际中，对机器人外部扰动的幅值估计难免会存在误差. 因此，控制器须对幅值估计误差有一定的适应性，特别是当外部扰动的估计幅值小于实际幅值时，控制系统不至于发散失稳. 设机器人第 $k$步外部扰动 ${{f}}{\eta _k}$的实际幅值为 ${{{F}}_{\rm m}}$， ${{{F}}_{\rm m}}$的最大元素为 ${{{f}}_{\rm m}}$，幅值估计 ${\hat{ F}}$的最小元素为 ${{\hat{f}}_{\rm m}}$，且令 ${{{f}}_{\rm r}} = {{{f}}_{\rm m}}/{{\hat{f}}_{\rm m}}$，则可以类似证明，基于控制器（式（17）），机器人的状态轨迹会收敛到边界层：

(22) ${B_1} = \{ ({{ q}},{\dot{ q}})|\left\| {{s}} \right\| < \phi {{{f}}_{\rm r}}{{ I}}_{8 \times 1}\}. $

因此，当外部扰动的幅值估计在一定范围内偏大或偏小时，控制系统仍然收敛. 实际上，当外部扰动的幅值估计偏小时，根据式（19）， ${{s}}$的模会增大，状态轨迹开始远离参考轨迹，而此时控制器（式（17））的反馈输入 ${{Ks}}$也会相应增大，从而控制系统依然能够保持收敛性；当外部扰动的幅值估计偏大时， ${{s}}$的模会减小，反馈输入 ${{Ks}}$也相应减小，从而系统不至于收敛到 ${{s}} = {\bf 0}$，导致发生颤振.

考虑到双足机器人在实际应用中常会遭遇非平整路面，进一步将自适应滑模控制器拓展到非平整路面的鲁棒控制. 在非平整路面上行走时，机器人容易失稳摔倒的重要原因是其与地面碰撞后会产生较大的速度突变，导致系统状态远离预定的参考轨迹^[7]. 为了实现非平整路面下的鲁棒控制，本研究提出碰撞速度不变性条件，基于落地速度控制进行在线轨迹规划以避免碰撞带来过大的速度突变，再基于自适应滑模控制器对机器人进行反馈控制.

1）碰撞速度不变性条件. 对于碰撞方程（式（5）），若机器人的落地速度 ${\dot {{P}}_{\rm{sw}}}({{q}}_{\rm{e}}^ - ,{\dot{ q}}_{\rm{e}}^ - ) = {\bf 0}$，则碰撞后机器人速度不发生突变，且碰撞作用力为零. 实际上，机器人的碰撞方程（式（5））可以改写为

(23) $ \left. {\begin{aligned} & {{\dot{ q}}_{\rm e}^ + - {{ D}}_{v,{\rm e}}^{ - 1}{{E}}_{v,{\rm{sw}}}^{\rm T}{{{F}}_{v,{\rm{sw}}}} = {\dot{ q}}_{\rm e}^ - } , \\ & \qquad \quad {{{{E}}_{v,{\rm{sw}}}}{\dot{ q}}_{\rm e}^ + = {\bf 0}} . \end{aligned}} \right\}$

从而有

(24) $ - {{{E}}_{v,{\rm{sw}}}}{{ D}}_{v,{\rm e}}^{ - 1}{{E}}_{v,{\rm{sw}}}^{\rm T}{{{F}}_{v,{\rm{sw}}}} = {{{E}}_{v,{\rm{sw}}}}{\dot{ q}}_{\rm e}^ - .$

将 ${\dot {{P}}_{\rm{sw}}}({{q}}_{\rm e}^ - ,{\dot{ q}}_{\rm e}^ - ) = {{{E}}_{v,{\rm{sw}}}}{\dot{ q}}_{\rm e}^ - = {\bf 0}$代入上式，可以得到

(25) ${{{E}}_{v,{\rm{sw}}}}{{ D}}_{v,{\rm e}}^{ - 1}{{E}}_{v,{\rm{sw}}}^{\rm T}{{{F}}_{v,{\rm{sw}}}} = {\bf 0}.$

${{{E}}_{v,{\rm{sw}}}}$、 ${{ D}}_{v,{\rm e}}^{ - 1}$都是满秩矩阵，所以矩阵 ${{{E}}_{v,{\rm{sw}}}}{{ D}}_{v,{\rm e}}^{ - 1}{{E}}_{v,{\rm{sw}}}^{\rm T}$可逆，从而有

(26) ${{{F}}_{v,{\rm{sw}}}} = {\bf 0}.$

将其代入式（23），可得 ${\dot{ q}}_{\rm e}^ - = {\dot{ q}}_{\rm e}^ + $，即碰撞后机器人速度不发生突变.

2）基于落地速度控制的在线轨迹规划. 为了控制落地速度，根据机器人预估落地位置进行在线轨迹规划，并以支撑腿关节 ${{{q}}_{\rm st}} = {[{{{q}}_1}, \cdots ,{{{q}}_5}]^{\rm T}}$以及摆动脚位置 ${{{P}}_{\rm{sw}}}$作为规划目标. 将支撑腿关节轨迹设计为如下贝塞尔多项式函数：

(27) ${{{q}}_{\rm st}}(t) = \sum\limits_{k = 0}^3 {{{{\alpha}} _k}\frac{{3!}}{{k!(3 - k)!}}{\tau ^k}{{(1 - \tau )}^{3 - k}}} .$

式中： $\tau = t/T$； ${{{\alpha}} _k}$为贝塞尔多项式系数，可通过机器人的初始状态以及末状态来确定.

将摆动脚的轨迹设计为

(28) $ {{{P}}_{\rm{sw}}}(t) = \left\{ {\begin{aligned} & {\sum\limits_{k = 0}^5 {{{{\beta}} _k}\frac{{5!}}{{k!(5 - k)!}}{\tau ^k}{{(1 - \tau )}^{5 - k}}} ,\;\;t \leqslant T} ;\\ & {{{{P}}_{\rm{sw,f}}} + \frac{1}{2}{{{a}}_{\rm s}}{{(t - T)}^2},\;\;\;T < t \leqslant T + {t_0}}; \\ & {{{{P}}_{\rm c}}{\rm{ + }}{{{v}}_{\rm s}}(t - T - {t_0}),\;\;\;\;\;\;\;\;\;t > T + {t_0}} . \end{aligned}} \right.$

式中： ${{{ \beta}} _k}$为贝塞尔多项式系数， ${{{P}}_{\rm{sw,f}}}$为预估落地位置， ${{{a}}_{\rm s}}$为常向量， ${t_0}$为常数， ${{{v}}_{\rm s}} = {{{a}}_{\rm s}}{t_0}$， ${{{P}}_{\rm c}} = {{{P}}_{\rm{sw,f}}} + $ $0.5{{{a}}_{\rm s}}t_0^2 $. 当 $t \leqslant T$时， ${{{P}}_{\rm{sw}}}(t)$为贝塞尔多项式函数，系数 ${{{\beta}} _k}$可通过摆动脚的初始位置 ${{{P}}_{{\rm{sw}},0}}$及速度 ${\dot {{P}}_{{\rm{sw}},0}}$、预设的中间位置 ${{{P}}_{{\rm{sw}},{\rm m}}}$、预估落地位置 ${{{P}}_{\rm{sw,f}}}$及速度 ${\dot {{P}}_{\rm{sw,f}}}$来确定. 为了控制落地速度，设定 ${\dot {{P}}_{\rm{sw,f}}} = {\bf 0}$. 若由于路面高度估计误差以及外部扰动的影响，机器人在预期时间段 $[0,T)$内未与地面发生碰撞，则摆动脚先以恒定的加速度 ${{{a}}_{\rm s}}$试探性下落，直至达到限定速度 ${{{v}}_{\rm s}}$，此后便以恒定速度 ${{{v}}_{\rm s}}$下落，直至发生碰撞.

3）关节角度计算与控制. 为了应用自适应滑模控制器（式（17）），须计算摆动腿关节 ${{{q}}_{\rm{sw}}} = {[{{{q}}_6},{{{q}}_7},{{{q}}_8}]^{\rm T}}$参考轨迹的角度信息、角速度信息以及角加速度信息. 以左脚支撑为例，摆动腿关节的运动学模型如图4所示. 图中， ${{{P}}_{\rm l,h}}$为左髋关节位置， ${{{P}}_{\rm r,h}}$为右髋关节位置， ${{{P}}_{\rm r0}}$、 ${{{P}}_{\rm r1}}$、 ${{{P}}_{\rm r2}}$为3个单位正交向量， ${{{P}}_{\rm r3}}$为右髋关节到摆动脚的向量，长度为 $L$， ${\gamma _1}$为向量 ${{{P}}_{\rm r1}}$与 ${{{P}}_{\rm r3}}$的夹角. 设 ${{{P}}_{\rm r3}}$在 ${{{P}}_{\rm r0}}$上的投影向量为 ${\lambda _1}{{{P}}_{\rm r0}}$， ${{{P}}_{\rm r3}}$在 ${{{P}}_{\rm r2}}$上的投影向量为 ${\lambda _2}{{{P}}_{\rm r2}}$，则根据几何关系有

(29) $\left. {\begin{aligned} & \quad \quad{{{{q}}_6} = - {\rm{atan}}({{{\lambda _1}}}/{{{\lambda _2}}})}, \\ & \quad{{{{q}}_7} = - ({\pi }/{2} - {\gamma _1} + {\gamma _2})} ,\\ & {{{{q}}_8} = \pi - {\rm{acos}}\;\left(\frac{{L_1^2 + L_2^2 - {L^2}}}{{2{L_1}{L_2}}}\right)} . \end{aligned}} \right\}$

直接利用式（29）对时间求导来获得 ${{\dot{ q}}_{\rm{sw}}}$、 ${\ddot {{q}}_{\rm{sw}}}$会产生较大的计算量. 为了提高计算效率，采用解析法求解 ${{\dot{ q}}_{\rm{sw}}}$、 ${\ddot {{q}}_{\rm{sw}}}$. 首先，利用式（27）、（28）对时间求导可得 ${{\dot{ q}}_{\rm st}}$、 ${\ddot {{q}}_{\rm st}}$、 ${\dot {{P}}_{\rm{sw}}}$、 ${\ddot {{P}}_{\rm{sw}}}$. 然后，根据运动学关系对 ${{{P}}_{\rm{sw}}}$求导，有

(30) ${\dot {{P}}_{\rm{sw}}} = {{{E}}_{1,{\rm{sw}}}}{\dot{ q}} = \left[ {\begin{array}{*{20}{c}} {{{{J}}_1}},&\!\!\!\!\!\!{{{{J}}_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{{\dot{ q}}}_{\rm st}}} \\ {{{{\dot{ q}}}_{\rm{sw}}}} \end{array}} \right],$

从而

(31) ${{\dot{ q}}_{\rm{sw}}} = {{J}}_2^ - ({\dot {{P}}_{\rm{sw}}} - {{{J}}_1}{{\dot{ q}}_{\rm st}}).$

再对式（30）求导，可得

(32) ${\ddot {{P}}_{\rm{sw}}} = {\dot {{E}}_{1,{\rm{sw}}}}{\dot{ q}} + \left[ {\begin{array}{*{20}{c}} {{{{J}}_1}},&\!\!\!\!\!\!{{{{J}}_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{{\ddot{ q}}}_{\rm st}}} \\ {{{{\ddot{ q}}}_{\rm{sw}}}} \end{array}} \right].$

从而有

(33) ${\ddot {{q}}_{\rm{sw}}} = {{J}}_2^ - ({\ddot {{P}}_{\rm{sw}}} - {\dot {{E}}_{1,{\rm{sw}}}}{\dot{ q}} - {{{J}}_1}{\ddot {{q}}_{\rm st}}).$

至此，已获得各个关节参考轨迹的角度信息、角速度信息以及角加速度信息，从而可利用控制器（式（17））对机器人进行鲁棒控制.

以下仿真实验假设机器人在第 $k$步的外部扰动 ${{f}}{\eta _k}$水平作用于质心，其中， ${{f}} = {[{{{f}}_x},0,0]^{\rm T}}$， ${{{f}}_x}$为扰动强度. 仿真实验的其他初始参数设置如下.

1）初始步态参数设置. 在无扰动情形下，机器人的行走步态左右对称，以左脚支撑相为例，机器人各个关节的初始状态和末状态分别为 ${{x}}_1^ + = {[{{q}}_1^ + ,{\dot{ q}}_1^ + ]^{\rm T}}$、 ${{x}}_1^ - = {[{{q}}_1^ - ,{\dot{ q}}_1^ - ]^{\rm T}}$，摆动脚的中间位置为 ${{{P}}_{{\rm{sw}},{\rm m}}}$，且

$\begin{aligned} {{q}}_1^ + = & [ - {\rm{0}}{\rm{.315\;7}},\;{\rm{0}}{\rm{.010\;3}},\; {\rm{0}}{\rm{.776\;9}},\;{\rm{0}}{\rm{.170\;0}},\; - {\rm{0}}{\rm{.144\;4}},\; \\ & - 0.145\;0,\;0.521\;7,\;{\rm{0}}{\rm{.762\;5}}{]^{\rm T}} ,\; \end{aligned} $

$\begin{aligned} {\dot{ q}}_1^ + = & [ - 1.936\;6,\;{\rm{0}}{\rm{.551\;6}},\; {\rm{1}}{\rm{.766\;6}},\; - {\rm{0}}{\rm{.453\;5}},\; - 1.185\;8,\; \\ & - 1.124\;9,\;0.634\;1,\;0.256\;9{]^{\rm T}},\; \end{aligned} $

$\begin{aligned} {{q}}_1^ - = &[ - 0.652\;9,\; - 0.013\;5,\; 0.762\;5,\;0.521\;7,\;0.145\;0,\; \\& 0.144\;4,\;0.170\;0,\;0.776\;9{]^{\rm T}},\; \end{aligned} $

$\begin{aligned} {\dot{ q}}_1^ - = & [ - 1.483\;0,\; - 0.618\;2,\; 0.203\;1,\;0.601\;9,\;0.999\;8,\; \\& 1.097\;3,\; - 0.558\;1,\;1.832\;4{]^{\rm T}} ,\; \end{aligned} $

${{{P}}_{{\rm{sw}},\;{\rm m}}} = {[0,\; - 0.15,\;0.08]^{\rm T}}.$

根据初始状态 ${{x}}_1^ + $和末状态 ${{x}}_1^ - $，可以求出机器人摆动脚的初始状态和末状态，进而可根据第4章的内容规划出机器人各关节的参考轨迹.

2）随机过程的样本生成. 为了开展不确性扰动下的动态步行仿真实验，生成了4组随机过程的样本，如图5所示为样本随时间t的变化. 其中，样本 ${\eta _{\rm s,1}}$、 ${\eta _{\rm s,2}}$、 ${\eta _{\rm s,3}}$、 ${\eta _{\rm s,4}}$对应的中心频率 $\varOmega $分别为 $2\pi $、1、 $5\pi $、 $6$，对应的谱宽参数 $\sigma $分别为10、2、30、42. 此外，对于每个样本， ${w_{\rm l}} = 0.002\pi $， ${w_{\rm r}} = 120\pi $， $N = 20\;000$.

为了验证所设计滑模控制器对外部扰动幅值估计误差的适应性，假设机器人在第 $k$步内受到的外部扰动为 ${{f}}{\eta _k}$， ${\eta _k} = {\eta _{{\rm{s}},3}}$， $k = 1,2,\cdots$，扰动强度 ${{{f}}_{\rm{x}}} = 5.39$ N，相当于 $10\% $的躯干重量. 此时，机器人系统为周期扰动系统，可基于常规庞卡莱映射方法进行稳定性分析. 为了证明当外部扰动的幅值估计在一定范围内偏大或偏小时，控制系统仍然收敛，本研究考虑估计幅值比实际幅值分别偏差−20%、+20%的情形. 在2种情形下，滑模控制器的控制参数均选取 $\phi $= $0.1$，仿真结果如图6所示. 可以看出，在2种情形下，机器人在第2个步行周期时都收敛到了稳定的极限环，说明控制器对外部扰动的幅值估计误差具有一定的适应性.

为了进一步分析控制器参数 $\phi $对控制器性能的影响，假设机器人在每一步内受到的外部扰动均为 ${{f}}{\eta _{{\rm{s}},3}}$，扰动强度 ${{{f}}_{\rm{x}}} = 5.39$ N. 为了说明不同控制器参数 $\phi $对控制器性能的影响，选取4组控制器参数进行仿真实验，分别为 $\phi = 0.005$、0.025、0.100、0.500. 在各组仿真实验中，外部扰动的估计幅值等于实际幅值，最后的仿真结果如图7所示. 可以看出，在各组仿真实验中，机器人系统在第2个步行周期时都收敛到了稳定的极限环，说明在不同控制器参数下，控制系统都具有较好的收敛速度. 另外，还可以看出，控制器参数 $\phi $越小，机器人最后的状态轨迹越接近理想轨迹，其相应的庞卡莱映射也与理想步态越接近，说明控制器参数 $\phi $越小，控制器的跟踪性能越好. 这是由于机器人的状态轨迹收敛到了某个更靠近理想步态的薄边界层 $B = \{ ({{ q}},{\dot{ q}})|\left\| {{s}} \right\| < \phi {{ I}_{8 \times 1}}\} $. 不过，在实际中，过小的控制器参数 $\phi $易导致颤振问题，因此须根据实际情况合理地选取控制器参数.

为了验证不确定性扰动下机器人动态步行的鲁棒性，根据5.2节的内容，关注机器人在有限步数 $n = 4$内的步行稳定性，并对不确定性扰动做周期化处理. 假设机器人在第 $k$步受到的外部扰动为 ${{f}}{\eta _{{\rm{s}},k}}$， $k = 1, \cdots ,4$，每一步的扰动强度 ${{{f}}_{\rm{x}}} = 10.78$ N，相当于 $20\% $的躯干重量. 机器人在每一步的控制器参数 $\phi $= $0.1$，外部扰动的幅值估计 ${\hat{ F}} = { f}$，仿真结果如图8所示. 可以看出，机器人的状态轨迹相图分布在有界区域，且特定庞卡莱映射收敛. 因此，根据有限时间稳定性，机器人在有限步数 $n = 4$内是稳定的，说明所设计控制器能有效实现机器人在不确定性扰动下的鲁棒控制. 值得注意的是，在图8（a）中，机器人的状态轨迹并没有收敛到某一个周期轨道，而是分布于4条状态轨道的周围. 实际上，根据文献[19]，当 $n$足够大时，机器人的状态轨迹将密集地分布于相空间，从而难以基于常规的庞卡莱映射方法或状态轨迹相图来判断系统的稳定性. 因此，对于不确定性扰动下的双足机器人，只能基于特定庞卡莱映射方法分析其在有限步数内的步行稳定性.

考虑机器人在非平整路面受到不确定性扰动作用的情形，同样关注机器人在 $n = 4$步内的步行稳定性. ${\eta _k} = {\eta _{{\rm{s}},k}}$， $k = 1, \cdots ,4$，扰动强度 ${{{f}}_{\rm{x}}} = 5.39$ N，相当于 $10\% $的躯干重量，路面相对高度在 $[ - 30,30]$ mm内随机变化，同样以4步为一个周期，并假设机器人每次预估的落地位置与地面实际高度存在 $10\% $的误差. 为了控制落地速度，设定下落加速度 ${{{a}}_{\rm s}} = {[{\rm{ - 6}},0,{\rm{ - }}6]^{\rm T}}$，下落速度 ${{{v}}_{\rm s}} = {[{\rm{ - 0}}{\rm{.1}},0,{\rm{ - }}0.1]^{\rm T}}$. 为了实现鲁棒控制，取控制器参数 $\phi $=0.1， ${\hat{ F}} = {{f}}$，仿真结果如图9所示. 可以看出，机器人的状态轨迹相图分布在有界区域，且特定庞卡莱映射收敛. 因此，机器人在有限步数步内是稳定的，说明所设计控制器能够有效实现机器人在非平整路面上的鲁棒控制.

（1）将随机振动系统的特定庞卡莱映射方法拓展到不确定性扰动下双足机器人的稳定性分析，并给出了双足机器人特定庞卡莱映射的构建方法.

（2）基于滑模控制方法，提出自适应滑模控制器. 仿真结果表明，该控制器对外部扰动的幅值估计误差具有较好的适应性；控制器的跟踪性能与控制器参数 $\phi $有关， $\phi $越小，控制器的跟踪性能越好.

（3）将自适应滑模控制器拓展到非平整路面的鲁棒控制. 通过在线规划控制机器人的落地速度，并基于自适应滑模控制器进行反馈控制，机器人在仿真实验中能实现相对高度在 $[ - 30,30]$ mm内随机变化的非平整路面上的鲁棒行走.

（4）后续工作将考虑将自适应滑模控制器拓展到具有不确定性动力学模型的双足机器人系统.