四足机器人因其出色的运动能力和较强的地形适应性成为应用前景非常广阔的特种机器人，常被用于探索、搜救和军事领域等，其在具有挑战性的障碍^[1]和自然地形中具有巨大的应用潜力，如Minitaur^[2]、MIT Cheetah^[3]、HYQ^[4]和ANYmal^[5]等，这些四足机器人在复杂的地形中具有良好的运动稳定性. Raibert于20世纪80年代实现了机器人的步态控制，提出的“弹跳高度-前进速度-机体位姿”三体解耦控制理论^[6]实现了腿足式机器人的平衡控制，此后高动态控制策略对四足机器人的发展和应用作出了巨大的贡献，如应用于Atlas^[7]的质心动力学(centroidal dynamics, CD)^[8]、线性倒立摆模型(linear inverted pendulum model , LIPM)和模型预测控制(model predictive control, MPC)^[9-10]等.

动物表现出多模式的行为，而受到生物学启发的机器人可以在具有挑战性的地形中更有效地保持平衡和协调. 朱雅光等^[11]的研究发现，基于响应外部信息并实现各种行为的神经架构的机器人^[12]在复杂地形上的运动具有卓越的协调控制能力，如四足机器人的鲁棒感知运动^[13]、多环境下机器人的自适应控制^[14-15]. 作为人与环境交互的工具，具有神经控制器的四足机器人非常必要，该控制器能够响应外部感觉和自我感知信息，实现多步态平滑切换^[16]和包括小跑在内的多种运动行为. 作为基本的神经组件，中央模式生成器(central pattern generator, CPG)主要负责为动物生成有节奏的运动，如爬行、跳跃、行走和跑步. 生物启发式的CPG控制已被证明在各种机器人^[17-18]中是成功的，包括利用人工神经网络穿越不规则地形的机器人^[19]、使用CPG实现滑行步态过渡控制的蛇形机器人^[20]及通过CPG控制两足动物运动^[21]等.

目前，大多数具有神经控制机制的机器人都将CPG信号与运动模块耦合，很难考虑腿部的独立控制，增大了模型的复杂性. 作为一类受生物启发的神经网络，能够自主、自我调节控制的CPG是四足机器人有节律运动必不可少的. 在复杂多变的环境中，具有简单模型和高调节能力的神经控制机制是必需的.

本文融合CPG和MPC，提出神经网络动力学控制策略，提高机器人运动的响应性和步态的多样性. CPG的作用体现在信号控制精度、行为模式多样性和平滑运动过渡方面，通过各种步态参数生成多种运动行为，使用MPC作为控制器，在不同步态参数下寻求保持机器人平衡的最佳反作用力，提高腿足式机器人的运动稳定性. 基于CPG-MPC控制架构的机器人可以实时实现多种行为，通过实验验证了基于CPG-MPC策略的不同步态模式的性能.

四足机器人的样机如图1所示，该机器人由机身龙骨和腿部结构组成，共有12个自由度. 简化腿部结构，将单腿设计为根关节、髋关节和膝关节，腿部关节配备直流无刷电机进行驱动，实现力矩及位置控制. 为了均衡腿部关节的质量，膝关节和髋关节的驱动电机放置在同轴侧大腿与机身的连接处，膝关节通过连杆完成运动传递. 为了获得关节的转动角度及角速度，关节驱动电机配备了14位编码器. 此外，在质心位置处搭载惯性测量单元，用来获取机身的加速度与姿态角信息. 为了缓解冲击带给机体结构的损坏，增大足端与地面接触的摩擦力，选用圆饼型橡胶作为足端，并且配备有六维力传感器. 四足机器人采用全肘式结构，平衡不均匀受力部位；材料选为铝合金，降低四足机器人在运动过程中产生的惯性力，结构的具体参数如表1所示.

合理的运动学模型是实现上层控制策略的基础，如图2所示，四足机器人的左后、左前、右后与右前腿依次表示为LH、LF、RH与RF. 定义机身坐标{B}位于机器人质心(center of mass, CoM)，设置x轴正方向为机器人的前进方向，z轴正方向垂直于机身平面竖直向上，y轴由右手定则确定. 世界坐标{W}位于机身坐标投影于地面的垂直方向，到机身坐标{B}的距离为H. 为了便于正逆运动学的求解，定义根关节坐标系{G_i}（i=1,2,3,4）位于髋关节与根关节轴线的相交处， 机身坐标系和根关节坐标系的x、y、z指向相同. 由于4条单腿结构的参数相同，建立任意条单腿的D-H模型，根据几何法求得4条腿的运动学方程. 通过计算确定足端相对根关节坐标系{G_i}下的位置为

(1)$ \left.\begin{gathered} {}^G{{{P}}_{i,x}} = - {a_3}{s_{23}} - {a_2}{s_2}, \\ {}^G{{{P}}_{i,y}} = \lambda {a_1}{c_1}+{a_2}{s_1}{c_2}+{a_3}{s_1}{c_{23}}, \\ {}^G{{{P}}_{i,z}} = \lambda {a_1}{s_1} - {a_2}{c_1}{c_2} - {a_3}{c_1}{c_{23}}. \\ \end{gathered} \right\}$

式中：s_i = sin θ_i, c_i = cos θ_i, s_ij = sin (θ_i+θ_j), c_ij = cos (θ_i +θ_j)；λ为左、右腿的符号变量，右腿为λ = −1，左腿为λ = 1.

将机身形心在世界坐标系{W}的位置表示为^WP_c=[^Wx, ^Wy, ^Wz]^T ，则任意腿足端相对于世界坐标系{W}下的位置表示为

(2)$ {}^W{{\boldsymbol{P}}_{\mathrm{f}i}} = \left[ {\begin{array}{*{20}{c}} {{}^W{{\boldsymbol{P}}_{\rm{c}}}} \\ 1 \end{array}} \right]+{}_B^W{\boldsymbol{R}} {}^B{{\boldsymbol{P}}_{{\mathrm{f}}i}}. $

式中：^BP_fi为足端基于坐标系{B}中的位置描述，$ {}_{B}^{W}{\boldsymbol{R}} $为旋转矩阵.

足式机器人作为欠驱系统，在某一运动方向成为欠驱方向前，为了保证机器人在该方向上具有稳定性，需要提前计算这个方向的足端作用力. 这需要控制器具有一定的预测性，可以提前调整各项状态参数，模型预测控制则成为足式机器人的主流控制策略.

除了提出的世界坐标系{W}和机身坐标系{B}，为了描述机器人的惯性张量、足端位置和足端反作用力等其他参数，引入定向机身坐标系{C}，坐标系的原点与坐标系{B}的原点重合，方向与坐标系{W}相同，如图3所示为简化后的单刚体模型(SRBD). 图中，足端圆点表示足端与地面的接触点，足端竖直向上的直线是基于坐标系{W}描述的足端反作用力，记作^Wf_i(i = 1, 2, 3, 4)，机器人本体的重力记作G(= mg)，虚线描述基于坐标系{B}的足端位置，记作^BP_fi(i=1, 2, 3, 4)，后面描述足端位置参数简写为P_fi，实线描述基于坐标系{W}的足端位置，记作^WP_fi(i=1, 2, 3, 4). ^WP_c是机器人质心基于坐标系{W}的描述，后面质心位置参数简写为P_c.

在刚体动力学模型(SRBD)中，足端反馈力U作为输出直接施加于机器人，控制机器人的位置和姿态. 根据牛顿-欧拉运动定律，可以得到机身平动和转动的关系式，即

(3)$ \left. \begin{gathered} {{\ddot {\boldsymbol{P}}}_{\mathrm{c}}} = m^{-1}{{\displaystyle \sum\nolimits_{i = 1}^4 {{{\boldsymbol{f}}_i}} }} - {{{\boldsymbol{g}}}}, \\ \sum\limits_{i = 1}^4 {{{\boldsymbol{P}}_{{\mathrm{f}}i}}} \times {{\boldsymbol{f}}_i} =\frac{{\mathrm{d}}}{{{\mathrm{d}}t}}({\boldsymbol{I}}{\boldsymbol{\omega }}) = {\boldsymbol{I}}\dot {\boldsymbol{\omega}} +{\boldsymbol{\omega}} \times ({\boldsymbol{I}}{\boldsymbol{\omega}} ) \approx {\boldsymbol{I}}\dot {\boldsymbol{\omega}} . \end{gathered}\right\}$

式中：I为基于坐标系{B}的惯性张量，ω为机器人的角速度，P_fi×定义为足端位置矢量的叉乘形式.将近似后的状态量之间的映射关系进行整合，可得机器人的状态空间表达式，即整理式（3），可得刚体动力学模型的连续状态空间表达式：

(4)$ \dot {\boldsymbol{X}}(t) = {{\boldsymbol{A}}_{\rm{c}}}({\boldsymbol{\varTheta}} ){\boldsymbol{X}}(t)+{{\boldsymbol{B}}_{\rm{c}}}({{\boldsymbol{P}}_{{\mathrm{f}}1}},{{\boldsymbol{P}}_{{\mathrm{f}}2}},{{\boldsymbol{P}}_{{\mathrm{f}}3}},{{\boldsymbol{P}}_{{\mathrm{f}}4}}){\boldsymbol{U}}(t). $

式中：${\boldsymbol{X}} = {\left[ {{{\boldsymbol{\varTheta}} ^{\mathrm{T}}},{{\boldsymbol{P}}_{\rm{c}}^{\mathrm{T}}},{^W}{{\boldsymbol{\omega}} ^{\mathrm{T}}},{\boldsymbol{v}}_{\rm{c}}^{\mathrm{T}},{^W}{{g}}} \right]^{\rm{T}}} $，其中Θ为机器人姿态角，v_c为机器人速度；$ {\boldsymbol{U}} = {\left[ {{\boldsymbol{f}}_1^{\mathrm{T}},{\boldsymbol{f}}_2^{\mathrm{T}},{\boldsymbol{f}}_3^{\mathrm{T}},{\boldsymbol{f}}_4^{\mathrm{T}}} \right]^{\rm{T}}} $；A_c为状态矩阵, A_c∈R^{(3n+1)×(3n+1)}，n = 4 ；B_c为控制矩阵, B_c∈R^(3n+1)×3n，n = 4.

为了简化研究对象，将刚体动力学模型的状态空间表达式离散化，即用△t表示模型预测控制的周期，采用前向欧拉法对连续状态空间表达式离散化，则

(5)$ {{\boldsymbol{X}}_{k+1}} = {{\boldsymbol{A}}_k}{{\boldsymbol{X}}_k}+{{\boldsymbol{B}}_k}{{\boldsymbol{U}}_k}. $

式中：A_k=I+△t·A_c, B_k=△t·B_c, k = 0, 1, 2,···.

根据模型预测控制原理，将求解反作用力U视为求解k个预测时域最优控制量的模型预测控制问题，表达式如下：

(6)$\left.\begin{aligned} &\min\; \boldsymbol{J}(\boldsymbol{U}) = \left\|\boldsymbol{X}_{(i+1)} - \boldsymbol{X}_{(i+1, \text { ref })}\right\| \boldsymbol{Q}_{(i+1)} + \left\|\boldsymbol{U}_{(i+1)}\right\| \boldsymbol{R}_{(i+1)} =\\ &\qquad \qquad \;\;\;\;\left(\boldsymbol{A}_{\mathrm{qp}} \boldsymbol{X}_0-\boldsymbol{X}_{(i+1, \text { ref })}+\boldsymbol{B}_{\mathrm{qp}} \boldsymbol{U}_{(i+1)}\right)^{\mathrm{T}} \boldsymbol{Q}_{(i+1)}\times \\&\qquad \qquad \;\;\;\;\left(\boldsymbol{A}_{\mathrm{qp}} \boldsymbol{X}_0-\boldsymbol{X}_{(i+1, \text { ref })}+\boldsymbol{B}_{\mathrm{qp}} \boldsymbol{U}_{(i+1)}\right)+\\&\qquad \qquad \;\;\;\;\boldsymbol{U}_{(i+1)}^{\mathrm{T}} \boldsymbol{R}_{(i+1)} \boldsymbol{U}_{(i+1)};\\&{\rm{s.t.}} \;{\boldsymbol{X}}_{(i+1)}={\boldsymbol{A}}_i{\boldsymbol{X}}_{(i)}+{\boldsymbol{B}}_i{\boldsymbol{U}}_{(i)},\\&\quad\;\; c_{\min}\leqslant U_{i,m}\leqslant c_{\max}, \;i=0,1, \cdots ,k-1.\end{aligned}\right\} $

式中：X_(i+1)为前i+1个预测时域的状态量；U_(i+1)为前i+1个预测时域的控制量，U_i,m为第i次预测的控制量的第m个元素；X₀和X_ref分别为初始状态量和参考状态量；Q_(i+1)和R_(i+1)分别为前i+1个预测时域的状态权重矩阵和能耗权重矩阵；

$\begin{split} \boldsymbol{A}_{\mathrm{qp}}=&\left[\begin{array}{c}\boldsymbol{A}_0 \\ \boldsymbol{A}_1 \boldsymbol{A}_0 \\ \vdots \\ {\displaystyle \prod}_{k=i}^0 \boldsymbol{A}_k\end{array}\right], \\\boldsymbol{B}_{\mathrm{qp}}=&\left[\begin{array}{cccc}\boldsymbol{B}_0 & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \boldsymbol{A}_1 \boldsymbol{B}_0 & \boldsymbol{B}_1 & \mathbf{0} & \mathbf{0} \\ \vdots & \vdots & & \vdots \\ \left({\displaystyle \prod}_{k=i}^1 \boldsymbol{A}_k\right) \boldsymbol{B}_0 & \left({\displaystyle \prod}_{{k}={i}}^2 \boldsymbol{A}_k\right) \boldsymbol{B}_1 & \cdots & \boldsymbol{B}_{k-1}\end{array}\right],\end{split}$

其中A_i和B_i分别为第i次预测的状态矩阵和控制矩阵；c_min和c_max为足端力的上、下限. 整理式（6），可得二次规划的标准形式：

(7)$ \left.\begin{gathered} \begin{array}{*{20}{c}} {\mathop {\min }\limits_U }&\end{array}{{\boldsymbol{U}}^{\mathrm{T}}}{\boldsymbol{HU}}/2+{{\boldsymbol{U}}^{\mathrm{T}}}{\boldsymbol{G}}; \\ \begin{array}{*{20}{c}} {{\mathrm{s.t.}}}&{{c_{\min }} \leqslant {U_{i,m}} \leqslant {c_{\max }}} .\end{array} \\ \end{gathered}\right\} $

优化问题(7)是二次规划问题，可以通过优化求解器来解决，如qpOASES. 在求解方程后，通过转矩映射将控制矢量U施加于机器人，即τ₁=J^TU. SRBD模型使用了忽略腿部惯性的假设，该假设适用于准静态运动. 若腿部与机身的质量比非常大，则施加于腿部的力矩在动力学中起着重要作用. 为了尽可能减少机器人腿部惯性的影响，增加了惯性补偿力矩，即

(8)$\underbrace {\frac{{\mathop{\rm d}\nolimits} }{{{\mathop{\rm d}\nolimits} t}}\left[ {\begin{array}{*{20}{c}}{\boldsymbol{\varTheta}} \\{{{\boldsymbol{P}}_{\rm{c}}}}\\{{}^W{\boldsymbol{\omega}} }\\{{{\dot {\boldsymbol{P}}}_{\rm{c}}}}\\{{}^W{{g}}}\end{array}} \right]}_{\dot {\boldsymbol{X}}} = \underbrace {\left[ {\begin{array}{*{20}{c}}{{{\boldsymbol{0}}_{3 \times 3}}}&{{{\boldsymbol{0}}_{3 \times 3}}}&{{\boldsymbol{R}}_z^{\mathrm{T}}(\alpha )}&{{{\boldsymbol{0}}_{3 \times 3}}}&{{{\boldsymbol{0}}_{3 \times 1}}}\\{{{\boldsymbol{0}}_{3 \times 3}}}&{{{\boldsymbol{0}}_{3 \times 3}}}&{{{\boldsymbol{0}}_{3 \times 3}}}&{{{\boldsymbol{I}}_{3 \times 3}}}&{{{\boldsymbol{0}}_{3 \times 1}}}\\{{{\boldsymbol{0}}_{3 \times 3}}}&{{{\boldsymbol{0}}_{3 \times 3}}}&{{{\boldsymbol{0}}_{3 \times 3}}}&{{{\boldsymbol{0}}_{3 \times 3}}}&{{{\boldsymbol{0}}_{3 \times 1}}}\\{{{\boldsymbol{0}}_{3 \times 3}}}&{{{\boldsymbol{0}}_{3 \times 3}}}&{{{\boldsymbol{0}}_{3 \times 3}}}&{{{\boldsymbol{0}}_{3 \times 3}}}&{ - {{\boldsymbol{I}}_{3 \times 1}}}\\{{{\boldsymbol{0}}_{1 \times 3}}}&{{{\boldsymbol{0}}_{1 \times 3}}}&{{{\boldsymbol{0}}_{1 \times 3}}}&{{{\boldsymbol{0}}_{1 \times 3}}}&{\boldsymbol{0}}\end{array}} \right]}_{{{\boldsymbol{A}}_{\rm{c}}}}\underbrace {\left[ {\begin{array}{*{20}{c}}{\boldsymbol{\varTheta}} \\{{{\boldsymbol{P}}_{\rm{c}}}}\\{{}^W{\boldsymbol{\omega}} }\\{{{\dot {\boldsymbol{P}}}_{\rm{c}}}}\\{{}^W{{g}}}\end{array}} \right]}_{\boldsymbol{X}} + \underbrace {\left[ {\begin{array}{*{20}{l}}{{{\boldsymbol{0}}_{3 \times 3}}}&{...}&{{{\boldsymbol{0}}_{3 \times 3}}}\\{{{\boldsymbol{0}}_{3 \times 3}}}&{...}&{{{\boldsymbol{0}}_{3 \times 3}}}\\{{}^B{{\boldsymbol{I}}^{ - 1}}[{{\boldsymbol{P}}_{{\mathrm{f}}1}}] \times }&{...}&{{}^B{{\boldsymbol{I}}^{ - 1}}[{{\boldsymbol{P}}_{{\rm{f}}n}}] \times }\\{{{{{\boldsymbol{I}}_3}} \mathord{\left/ {\vphantom {{{{\boldsymbol{I}}_3}} m}} \right. \kern-\nulldelimiterspace} m}}&{...}&{{{{{\boldsymbol{I}}_3}} \mathord{\left/ {\vphantom {{{{\boldsymbol{I}}_3}} m}} \right. \kern-\nulldelimiterspace} m}}\\{{{\boldsymbol{0}}_{1 \times 3}}}&{...}&{{{\boldsymbol{0}}_{1 \times 3}}}\end{array}} \right]}_{{{\boldsymbol{B}}_{\rm{c}}}}\underbrace {\left[ {\begin{array}{*{20}{c}}{{}^W{{\boldsymbol{f}}_1}}\\{{}^W{{\boldsymbol{f}}_2}}\\{{}^W{{\boldsymbol{f}}_3}}\\{{}^W{{\boldsymbol{f}}_4}}\end{array}} \right]}_{\boldsymbol{U}}.$

(9)$ {{\boldsymbol{\tau}} _2} = {{\boldsymbol{M}}_{\mathrm{l}}}{\ddot {\boldsymbol{q}}}. $

式中：M_l为腿的惯性矩阵， $ \ddot{{\boldsymbol{q}}} $为角加速度. 实际的角加速度可以通过编码器获得，但获得的补偿力矩具有较高的频率，这导致机器人状态不稳定，因此使用期望的角加速度.

运动神经控制系统主要由中枢神经网络和落足点规划模块组成，是控制系统的规划层. 腿足式机器人的步态生成和步态模式的切换主要依赖可靠的节律信号，该信号由中枢神经网络产生. 落足点模块根据机器人的运动状态和环境信息求解适当的足端位置，使得机器人在复杂环境中不仅可以实现多种稳定行为，还能够及时调整位置和姿态，以躲避和跨越各种障碍物.

基于神经架构的机器人能够对外部信息作出响应并实现各种复杂的运动行为^[22]，因此CPG被广泛应用于机器人运动控制. 考虑到强耦合、控制复杂及传统CPG模型难以控制机器人足端轨迹的问题，通过改进的σ-Hopf谐振子来控制四足机器人运动时的占地系数、运动频率与相位差，其公式表示为

(10)$ \left. {\begin{array}{*{20}{l}} {\dot O_{\mathrm{X}} = \alpha (\mu - {r^2})(O_{\mathrm{X}} - {\mu _1}) - \omega (O_{\mathrm{Y}} - {\mu _2}),} \\ {\dot O_{\mathrm{Y}} = \alpha (\mu - {r^2})(O_{\mathrm{Y}} - {\mu _1}) - \omega (O_{\mathrm{X}} - {\mu _2}),} \\ {\omega = \dfrac{\text{π} }{{k ({{\mathrm{exp}}\;{( - aO_{\mathrm{Y}})}}+1) T}}+\dfrac{\text{π}}{{(1 - k) ({{\mathrm{exp}}\;{(aO_{\mathrm{Y}})}}+1) T}}.} \end{array}} \right\} $

式中：O_X和O_Y为步态信号，r为极限环半径，μ为振幅，μ₁和μ₂为外反馈，占地因子k用来控制单腿摆动相和支撑相之间的时间比率，相位调节因子γ实现机器人的步态切换，周期T控制机器人的运动周期.

仿生机器人由多个肢体或关节组成，需要通过具有一定相位差的控制信号完成运动协调，其运动协调神经调节模块如图4(a)所示. 左后腿(LH)和左前腿(LF)的相位移动分别为τ_L= kτ和τ+τ_L，其中τ = (1—γ) T，如图4(b)所示为不同γ、k和T对应的CPG信号. 四足机器人在自然环境中穿越不同的地形，需要通过调整步态模式来适应这些环境，引入基于CPG的控制系统可以极大地提高机器人对环境的适应性.

腿的摆动主要是为了足端达到期望的落足点位置，间接地控制机器人躯干的位姿和速度，保证机身的稳定性. 良好的足端轨迹可以减小足端落地时对地面的冲击力，提高腿部适应复杂地形的能力. 在设计足端的摆动轨迹时，需要重点考虑期望落足点、步长和抬腿高度等. 机器人期望的线速度和角速度参数由用户输入，其中横向速度$ {v}_{{\mathrm{c}},x}^{{\mathrm{d}}} $ 和纵向速度$ {v}_{{\mathrm{c}},y}^{{\mathrm{d}}} $转换为足端位移ΔL_xy0，如图5所示，表达式如下：

(11)$ \Delta {{\boldsymbol{L}}_{xy0}} = \frac{1}{{{f_{{\mathrm{sw}}}}}}\left[ {\begin{array}{*{20}{c}} {{{{v}}_{{\mathrm{c}},x}^{\mathrm{d}}}} \\ {{{{v}}_{{\mathrm{c}},y}^{\mathrm{d}}}} \\ 0 \end{array}} \right]. $

式中：f_sw为摆动相的频率. 航向角速度$ {\dot{\alpha }}^{{\mathrm{d}}} $也会产生足端位移偏差：

(12)$ \Delta {{\boldsymbol{L}}_{{\mathrm{h}}0}} = \frac{1}{{{f_{{\mathrm{sw}}}}}}\left[ {\begin{array}{*{20}{c}} 0 \\ 0 \\ {{{\dot \alpha }^{\mathrm{d}}}} \end{array}} \right] \times {{\boldsymbol{x}}_{{\mathrm{hip}}}}. $

式中：x_hip为坐标系{G}在坐标系{B}中的位置描述.

对“线速度-位移”向量ΔL_xy0和“角速度-位移”向量ΔL_h0矢量求和，可以得到“速度-位移”矢量ΔL_xy，即

(13)$ \Delta {{\boldsymbol{L}}_{xy}} = \Delta {{\boldsymbol{L}}_{{\rm{h}}0}}+\Delta {{\boldsymbol{L}}_{xy0}}. $

为了足端摆动相的起始点尽可能接近期望位置并减少位移误差累积，计算虚拟起始点与当前实际起始点的差值，即

(14)$ {{\boldsymbol{f}}_{{\mathrm{v}}xy}} = {\boldsymbol{x}}_{\mathrm{f}}^0 - {{\boldsymbol{x}}_{\mathrm{f}}}. $

式中：$ {{\boldsymbol{x}}}_{{\mathrm{f}}}^{0} $为虚拟足端起始点，x_f为当前的实际足端位置. 位移矢量ΔL_xy和f_vxy矢量叠加为足端的期望位移，即

(15)$ \Delta {\bar {\boldsymbol{L}}_{xy}} = \Delta {{\boldsymbol{L}}_{xy}}+{{\boldsymbol{f}}_{{\rm{v}}xy}}. $

通过$ {\Delta }{\bar{\boldsymbol{L}}}_{xy} $计算摆动相步长L_s和摆动轨迹所在平面的倾角φ_s, 即

(16)$ {{{L}}_{\mathrm{s}}} = \sqrt {{{\Delta {{\bar L }_x}}^2}+{{\Delta {{\bar L}_y}}^2}} . $

(17)$ {\varphi _{\mathrm{s}}} = {\mathrm{atan}}2 \;(\Delta {\bar L _y},\Delta {\bar L _x}). $

式中：$\Delta {\bar{{ L}}_x}、\Delta {\bar {{L}}_y} $为$\Delta {\bar{{ {\boldsymbol{L}}}}_{xy}} $的分量.

这种参考相关研究^[23]的落足点(见图6)的选取策略在面对崎岖地形时可以提高机器人的鲁棒性，避免了机器人在之前的足端位置采取步骤时可能出现的错误累积，提高了机器人的运动精度.

通过CPG规划足端的理论轨迹，如图7所示，耦合CPG控制信号和图5所示的落脚点规划，根据用户输入的运动参数(运动周期T、相移系数γ和抬腿高度h等)，CPG模块可以生成步态信号O_X和O_Y，通过摆动相步长、倾角和运动信号等运动特征参数对足端的摆动轨迹进行设计. 采用多项式函数进行轨迹插值，确定摆动腿足端整个摆动持续时间内的足端期望位置P_f,d. 根据弹簧阻尼模型，计算添加至摆动腿的虚拟力，即

(18)$ {\tau _3} = {K_{\mathrm{p}}}({\theta ^{\mathrm{d}}} - \theta )+{K_{\mathrm{d}}}({\dot \theta ^{\mathrm{d}}} - \dot \theta ). $

式中：K_p和K_d分别为摆动腿的刚度和阻尼系数，θ和 $ \dot{\theta } $分别为关节位置和速度.

机器人的控制系统可以简化为运动规划和运动控制2部分. 完善控制系统使机器人拥有智能性和对复杂环境的适应性已成为四足机器人研究的重要组成部分，而融合神经网络的四足机器人逐渐拥有了感知地形和路线自主规划的能力. 设计的四足机器人的神经控制器如图8所示. 控制框架主要由4个部分组成：步态机(CPG)、状态估计器、模型预测控制器(MPC)和参考轨迹生成器. 步态机采用σ-Hopf振荡器，主要的控制参数是相移系数γ、振荡频率ε(=1/T)和占空比k，参考轨迹生成器用于生成预测时域的机身状态信息X_ref. 落足点规划模块主要根据机器人当前的位置信息和期望信息生成步长L_s和轨迹平面倾角φ_s，结合CPG信号规划足端的期望位置P_f,d. 运动学模块将期望足端位置P_f,d转换成期望关节转角θ_d，与实际的关节转角一起组成基于关节空间的PD控制器. 模型预测控制器主要负责根据建立的单刚体模型完成控制量的优化，控制器根据约束目标、约束方程及对应的权重系数计算得到合理的足端反力U；除此之外，控制器内还包含腿惯性补偿模块，提高控制器的精度，使机器人更加精确地完成期望动作. 状态估计器是根据卡尔曼滤波原理并结合惯性测量单元向上位机反馈信息，估计机器人实际的位置和角度.

为了验证神经控制器的有效性，通过仿真和实验平台进行多种实验. 其中仿真平台完成神经控制器的对比实验，机器人位置跟踪、全向运动和多步态运动实验在实验平台上完成.

为了突出神经控制器的优势，开展2次对比试验. 如图9所示为根据用户输入速度参数积分的四足机器人期望的位移，如图10所示为未添加神经控制器的位置误差曲线. 如图11所示为添加神经控制器的位置误差曲线. 图中，e_x、e_y、e_z分别为x、y、z方向上的位置误差.

对比图10、11可知，基于神经控制器的质心位置误差小于PD控制器. 其中x和y方向最显著，误差为−0.005~0.005 m，z方向的振幅约为0.02 m.

除了关注机器人的位置误差外，本文记录了线速度误差和姿态角误差. 如图12(a)、(b)所示分别为添加神经控制器前、后的线速度误差曲线. 图中，e_vx、e_vy、e_vz分别为x、y、z方向上的线速度误差. 从图12可以看出，神经控制器具有较高的稳定性，沿x和y方向的线速度误差小于PD控制器，但沿z方向的速度误差基本一致，为−0.5~0.5 m/s，反映了图10、11位置跟踪误差的变化趋势.

如图13(a)、(b)所示分别为添加神经控制器前、后的姿态角误差曲线. 图中，e_λ、e_β、e_α分别为翻滚角λ、俯仰角β、偏航角α的误差. 从图13可知，神经控制器对姿态角控制稍弱，其中β和α的数值变化范围略大，俯仰角是−0.04~0.02 rad，偏航角为−0.005~0.005 rad.

仿真结果表明，与单一的PD控制器相比，在 MPC-CPG 控制器的调节下，机器人能够迅速地消除位置误差和姿态角误差并达到稳定状态. 为了进一步验证神经控制器的响应性和有效性，通过图1的机器人进行位置调节和多种步态行进实验，硬件模块如图14所示.

机器人的硬件控制模块包括电机、驱动器和传感器等. 位姿传感器负责获取机器人的位置和姿态信息，力传感器检测足端力，编码器检测关节电机转动的角度，用户通过遥控器向上位机NI-sbRIO发送命令，其板载有RS-232接口、RS-485接口、CAN接口、USB接口和以太网网口等，实现了上位机和位姿传感器、遥控器、驱动器和编码器之前的数据通信. NI-sbRIO主控制器处理器型号为Xilinx Zynq-7xxx，为双核异构处理器，其中ARM核为Cortex-A9架构，主频高达667 MHz，用于准确、高效地运行CPG-MPC控制算法.

四足机器人静态站立时利用手柄的摇杆改变期望的位置参数，控制器检测并接收来自手柄的模拟信号并将其转化为数字信号，优化求解器求解得到足端力. 在位置调节过程中，将期望质心位置(x_CoM和y_CoM)设置为0 m，期望的翻滚角、俯仰角和偏航角都是0，如图15所示为机器人的质心位置和高度的实际变化曲线和效果图. 从图15(a)可知，x和y在高度调节时有轻微抖动，抖动幅度小于0.01 m. 从图15(b)可知，控制器可以在短时间内响应期望位置参数，满足响应性要求.

如图16所示为线速度的变化曲线，沿x轴和y轴方向的速度都有高频次的波动，但整体数据在−0.01~＋0.01 m/s变化，z方向的速度在进行高度调节时变化比较剧烈，其余时刻沿z轴速度基本在0附近波动. 相比于速度，姿态角变化的频次不是很高，但是波动的范围较大. 如图17所示，翻滚角在0~0.01 rad波动，俯仰角在−0.005~0.005 rad变化，偏航角在−0.027~−0.020 rad变化，偏航角和俯仰角后期都趋于稳定值.

在平坦地面上，对机器人全方向运动能力及运动状态受控性进行验证. 从零时刻开始，向机器人发送由纵向速度v_y、横向速度v_x、转向速度ω_z组成的混合指令，如图18、19所示分别为机器人线速度和质心位置曲线. 从v_x曲线图可以看出，在5~7 s时，实际速度曲线波动剧烈，为−0.2~0.2 m/s，除此之外，速度切换的时间点横向速度的曲线变化也很大，其余时刻的实际速度曲线一直在期望速度曲线附近稳定振荡. 速度决定实际位移，位移曲线斜率反映速度曲线的变化程度，如图19所示，z向期望速度为0，但实际上z向波动最剧烈，在0~0.02 m内变化.

全向运动中角速度和姿态角曲线变化如图20、21所示. 其中绕x轴和y轴的角速度为−0.5~0.5 rad/s，绕z轴的角速度部分区域变化剧烈，在8~12 s机器人以横向速度0.1 m/s和转向速度±0.05 rad/s进行复合运动，此时机器人的最大转向速度是整个运动过程中的最大值，约为0.18 rad/s，其余时刻绕z轴的角速度曲线变化平缓. 翻滚角和俯仰角在0 rad附近波动，最大航向角约为0.1 rad.

在进行全向运动的过程中，神经控制器的输出足端力f_x、f_y、f_z如图22所示，设定x和y方向的足端力无限定范围，z向足端力取值为[0, 300] N. 从图22可知，在以小跑步态行进过程中，X方向的足端力为0~100 N，结合图18可知，第12~14 s时，v_x达到最大值，对应的足端力也达到峰值. y方向的足端力为−20~20 N，z向足端力为0~150 N，满足机器人的质量要求.

为了测试和总结基于模型预测原理的神经控制器对多种常规和特殊步态的响应性以及运动参数对机器人稳定性的影响，开展30组实验，其中γ取0.50、0.55、0.60、0.65、0.70、0.75，k取0.75、0.76、0.78、0.80、0.82.

如图23所示为不同γ和k所对应的机器人质心速度和姿态角示意图. 图中，不同的颜色区域表示不同的数值，在图23的区域颜色变浅的过程中，区域代表的数值在递增. 机身高度方向速度v_z相对于v_x和v_y比较稳定，最大值出现在k =0.80~0.82和0.76~0.78处，约为±0.000 6 m/s. 翻滚角为0.001~0.006 rad，最大值出现在k = 0.8附近. 偏航角的最大值出现在右下角，约为−0.002 rad，当k = 0.75~0.78时，机器人的实际偏航角接近0. 从图23可以看出，随着γ的增加和k的减小，v_x接近期望速度0.1 m/s，v_y和v_z接近0，姿态角的实际值也接近0，即当足端下落时，对地面的冲击减小，同时有更多足端接触地面，机器人的速度和角度参数更接近期望值，这意味着机器人平地前行时的状态更稳定.

为了更加直观地探究机器人的运动状态对稳定性的影响，引入能耗参数（cost of transport，COT）评价指标，如图24所示为k和γ的COT. 可以看出，COT值在3.0~4.0的范围内波动，随着k的增加，足端触地时加剧了对地面的冲击，导致机器人的COT参数增大. 当k = 0.82，γ = 0.5时，COT最大，约为3.878. 当k = 0.8且γ = 0.75时，结合图23可知，此时质心沿z方向的速度和偏航角偏离期望值较大，所以COT增长至3.64. 总体分析可知，当机器人的k从0.5变化到0.75且γ从0.82下降到0.75，机器人以接近0.1 m/s的速度前进时，COT从3.878逐渐下降到3.372，机器人以更低的能耗率完成工作，这意味着机器人的状态更加稳定.

为了提高四足机器人的环境自适应性和行为多样性，本文设计融合CPG和MPC的神经控制器. 通过实验验证了利用提出方法可以实现机器人位置的跟踪调整，完成全向复合运动，同时可以根据工况及时调整机器人运动步态，达到降低机器人的COT，实现稳定运动的目的. 在所提方法的基础上，可以针对融合机器人复杂时变非线性动力学模型与中枢模式发生器开展进一步的研究，提高控制精度和强扰动情况下的机器人自适应调整能力.