关节机器人轨迹跟踪PID(proportional integral derivative，PID)控制实质上是调整每个电机驱动的关节所对应的PID控制器参数，实现关节位置或速度的反馈控制。但是，电机驱动的饱和非线性、关节摩擦等机器人所固有的制约因素会导致控制稳定性降低并加大轨迹跟踪误差。为此，一些学者提出了将模糊控制与PID控制相结合的模糊PID控制方法，两者的结合可达到更好的控制稳定性和鲁棒性^[1-5]。尽管如此，对于模糊控制的隶属函数中心顶点位置和模糊规则等模糊控制器设计参数以及PID控制参数等的选取和设计需要经验和专家知识，存在设计不合理时难以取得最佳控制效果的不足，于是有学者提出了运用遗传算法解决这一问题^{[4, 6]}。但是，采用遗传算法又可能引起染色体编码冗长，导致搜索空间过大、寻优效率降低，为此，Young等^[7]提出了一种用特征参数表征模糊规则、隶属函数等设计参数来简化模糊控制器设计过程的方法，随后有学者提出直接对模糊控制器的隶属度函数参数进行编码的矩阵式个体编码^[8]、模糊控制器模糊规则集和隶属函数联合编码并自适应选取交叉和变异概率^[9]等简化方法，Jahedi等^[10]则更进一步，将特征参数编码的遗传优化方法从应用于模糊控制器的设计过程推进到模糊PID控制器的设计过程，实现了动态供能系统能效增益的过程控制。

本文将模糊控制器设计参数直接编码和特征参数编码的混合编码策略应用于模糊PID控制器隶属函数和模糊规则的遗传算法优化，并应用优化后的模糊PID控制器实现关节机器人的轨迹跟踪控制。对比分析未经遗传算法优化的模糊PID控制器(FPID)、经直接编码的遗传算法优化后的模糊PID控制器(GAdirect-FPID)和经混合编码的遗传算法优化后的模糊PID控制器(GAhybrid-FPID)等3种控制器实现六关节机器人轨迹跟踪控制效果。

1 关节机器人动力学模型

考虑n关节机器人，忽略摩擦力及机器人末端负载，其动态性能可由二阶非线性动力学方程描述为：

$ \boldsymbol{\tau}=\boldsymbol{M}(\boldsymbol{q}) \ddot{\boldsymbol{q}}+\boldsymbol{C}(\boldsymbol{q}, \dot{\boldsymbol{q}}) \dot{\boldsymbol{q}}+\boldsymbol{G}(\boldsymbol{q}) $

(1)

式中：τ∈R^n×1是关节力矩向量; q∈ R^n×1是机器人关节角度向量; $\mathit{\boldsymbol{\dot q}}$∈ R^n×1是机器人关节角速度向量; $\mathit{\boldsymbol{\ddot q}}$∈R ^n×1是机器人关节角加速度向量；M(q)∈R^n×n是惯性对称矩阵; C(q, $\mathit{\boldsymbol{\dot q}}$)∈ R^n×n是哥氏力和离心力系数矩阵; G(q)∈R^n×1是重力向量。驱动关节运动所需要的力矩τ可通过牛顿-欧拉迭代动力学方程计算得到^[11]。

当求得M(q_t)和C(q_t, $\mathit{\boldsymbol{\dot q}}$_t)$\mathit{\boldsymbol{\dot q}}$_t+G(q_t)，并已知某t时刻机器人关节角度q_t、关节角速度$\mathit{\boldsymbol{\dot q}}$_t、轨迹跟踪控制器输出的关节力矩τ_t，由动力学方程式(1)，可通过仿真求出当前t时刻关节角加速度$\mathit{\boldsymbol{\ddot q}}$_t：

$ \ddot{\boldsymbol{q}}_{\mathrm{t}}=\boldsymbol{M}^{-1}\left(\boldsymbol{q}_{t}\right)\left[\boldsymbol{\tau}_{t}-\boldsymbol{C}\left(\boldsymbol{q}_{t}, \dot{\boldsymbol{q}}_{t}\right) \dot{\boldsymbol{q}}_{t}-\boldsymbol{G}\left(\boldsymbol{q}_{t}\right)\right] $

假设时间步长Δt内$\mathit{\boldsymbol{\ddot q}}$_t不变化，在初始条件为q₀=$\mathit{\boldsymbol{\dot q}}$₀=0的情况下，对当前t时刻关节信息q_t和$\mathit{\boldsymbol{\dot q}}$_t进行数值积分，计算出下一时刻t+Δt的角速度$\mathit{\boldsymbol{\dot q}}$_t+Δt和关节角q_t+Δt为：

$ \begin{gathered} \dot{\boldsymbol{q}}_{t+\Delta t}=\dot{\boldsymbol{q}}_{t}+\ddot{\boldsymbol{q}}_{t} \Delta t \\ \boldsymbol{q}_{t+\Delta t}=\boldsymbol{q}_{t}+\dot{\boldsymbol{q}}_{t} \Delta t+\frac{1}{2} \ddot{\boldsymbol{q}}_{t} \Delta t^{2} \end{gathered} $

2 轨迹跟踪模糊PID控制器设计

机器人轨迹跟踪模糊PID控制器结构如图 1所示，图中q_d表示关节机器人的期望角位移；q表示根据关节机器人动力学模型仿真获得的角位移；e表示q_d与q的偏差，e=q_d-q；e_c表示偏差变化率，e_c=de/dt；ΔK_P、ΔK_I、ΔK_D表示PID调整参数；τ表示模糊PID控制器作用于关节机器人的关节力矩。

机器人轨迹跟踪模糊PID控制器具体控制原理为：首先将当前偏差e和偏差变化率e_c作为模糊控制器的输入变量，经过模糊化后得到e和e_c隶属不同模糊集的隶属度；然后根据模糊规则进行模糊推理得到输出变量ΔK_P、ΔK_I、ΔK_D的模糊结果，通过解模糊化得到输出变量精确结果；接着利用ΔK_P、ΔK_I、ΔK_D对PID控制器参数进行：

$ \left\{\begin{array}{l} K_{\mathrm{P}}=K_{\mathrm{P} 0}+\Delta K_{\mathrm{P}} \\ K_{\mathrm{I}}=K_{\mathrm{I} 0}+\Delta K_{\mathrm{I}} \\ K_{\mathrm{D}}=K_{\mathrm{D} 0}+\Delta K_{\mathrm{D}} \end{array}\right. $

(2)

的在线自适应调整，最后得到模糊PID控制器作用于关节机器人的关节力矩τ：

$ \begin{gathered} \boldsymbol{\tau}=K_{\mathrm{P}} \boldsymbol{e}+K_{I} \int \boldsymbol{e} \mathrm{d} t+K_{\mathrm{D}} \frac{\mathrm{d} \boldsymbol{e}}{\mathrm{d} t}=\left(K_{\mathrm{P} 0}+\Delta K_{\mathrm{P}}\right) \boldsymbol{e}+ \\ \left(K_{\mathrm{I} 0}+\Delta K_{\mathrm{I}}\right) \int \boldsymbol{e} \mathrm{d} t+\left(K_{\mathrm{D} 0}+\Delta K_{\mathrm{D}}\right) \cdot \frac{\mathrm{d} \boldsymbol{e}}{\mathrm{d} t}= \\ \left(K_{\mathrm{P} 0} \boldsymbol{e}+K_{\mathrm{I} 0} \int \boldsymbol{e} \mathrm{d} t+K_{\mathrm{D} 0} \frac{\mathrm{d} \boldsymbol{e}}{\mathrm{d} t}\right)+ \\ \left(\Delta K_{\mathrm{P}} \boldsymbol{e}+\Delta K_{\mathrm{I}} \int \boldsymbol{e} \mathrm{d} t+\Delta K_{\mathrm{D}} \frac{\mathrm{d} \boldsymbol{e}}{\mathrm{d} t}\right) \end{gathered} $

式中：K_P0、K_I0、K_D0为PID参数ΔK_P、ΔK_I、ΔK_D的初值。模糊控制器的输入变量e和e_c以及输出变量ΔK_P、ΔK_I、ΔK_D的模糊集个数均设计为7个，取为{NL，NM，NS，ZE，PS，PM，PL}，即{负大，负中，负小，零，正小，正中，正大}。模糊集隶属函数选取最常见的连续三角形函数，对于模糊集个数为7的隶属函数，指定以零为中心的单侧3个中心顶点值即可，其余输入输出变量的隶属函数同。由2个输入变量的模糊集个数为7可得输出变量ΔK_P、ΔK_I、ΔK_D的每一个变量所对应的模糊规则均为49条，每条模糊规则形式为：

$ \begin{aligned} &\text { if e is} \cdots \text {and ec is } \cdots \text { , }\\ &\text { then } \Delta K_{\mathrm{P}} \text { is } \cdots, \Delta K_{\mathrm{I}}\ \mathrm{is} \cdots, \Delta K_{\mathrm{D}} \text { is } \cdots \end{aligned} $

本文设计的ΔK_P、ΔK_I、ΔK_D模糊规则如表 1所示。

3 模糊PID控制器的遗传算法优化

遗传算法优化关节机器人轨迹跟踪控制模糊PID控制器的算法主要步骤有：

1) 根据模糊控制器的特征参数和设计参数取值范围，采用混合编码方法，随机生成N个个体染色体[a₁, a₂, …, a₆₃]作为当前种群。

2) 更新当前种群为父群。

3) 读取当前种群第j(j=1, 2，…, N)个染色体各个基因位编码，根据编码a₁-a₁₀生成新的输入变量e、ec和输出变量ΔK_P、ΔK_I、ΔK_D的隶属函数，更新隶属函数；根据编码a₁₁-a₆₃产生新的输出变量ΔK_P、ΔK_I、ΔK_D模糊规则，更新模糊规则。

4) 运行模糊PID控制器完成轨迹跟踪控制仿真，获得t时刻机器人仿真轨迹与期望轨迹的偏差e(t)，计算目标函数值J。选择ITAE准则^[12]设计控制器的目标函数：

$ J=\int_{0}^{T} t|e(t)| \mathrm{d} t \rightarrow \min $

(3)

式中: e(t)为t时刻的机器人仿真轨迹与期望轨迹偏差; T为运行模糊PID控制器完成一次轨迹跟踪控制仿真的总时间。

5) 若个体染色体数量j不满足j>N，则j+1，进入步骤3；若满足j>N，则进入6)。

6) 判断目标函数值的变化是否小于设定的阈值或者遗传代数是否达到设定值，若满足则结束；若不满足则进入7)。

7) 排序当前种群中每个个体染色体对应的目标函数值，对当前种群进行选择、复制、变异、交叉等遗传操作，生成下一代子群，进入2)。

3.1 模糊PID控制器的特征参数

模糊PID控制器的设计参数具体是指输入、输出变量的隶属函数中心顶点位置和由每条模糊规则if条件获得的输出规则then结果。本文引入模糊控制器的特征参数的目的是简化模糊控制器设计过程，分别包括隶属函数和模糊规则特征参数^{[7, 10]}。

对于隶属函数的特征参数，根据隶属函数设计的特性可用表示单侧三角形隶属函数中心顶点位置与零的疏密程度X_mi：

$ X_{m i}=\left|\frac{i}{\left(\frac{k-1}{2}\right)}\right|^{P_{m}}, \quad i=1,2, \cdots, \frac{k-1}{2} $

(4)

式中：k是隶属函数模糊集个数，此处k=7；P_m是幂函数的指数，取值将影响隶属函数各个中心值的疏密分布情况。可进一步采用伸缩因子G描述三角形隶属函数顶点的疏密程度X_i：

$ X_{i}=G X_{m i}, i=1,2, \cdots, \frac{k-1}{2} $

将P_m和G定义为隶属函数的特征参数，则隶属函数的设计参数X_i可由P_m和G生成，即：设计隶属函数时由需要指定单侧3个三角形隶属函数中心顶点位置转换为指定2个特征参数，可见隶属函数的设计过程变得简化了。

对于模糊规则的特征参数，由表 1可看出，ΔK_P和ΔK_D的输出规则可划分成几个具有相同规则结果的区域，区域的数量取决于输出变量隶属函数的模糊集个数k。以ΔK_P的模糊规则为例说明采用相平面进行区域划分的方法。如图 2所示，首先分别根据ΔK_P对应的输入变量e、e_c的模糊集取值{NL，NM，NS，ZE，PS，PM，PL}建立相平面XOY，X、Y轴不同模糊集值的交点为网格点，可见网格点代表了模糊规则的if部分，网格点距离X、Y轴的距离表征了e、e_c的各自不同模糊集值所对应的三角形隶属函数中心顶点位置X_i值；然后定义图 2中的黑点为分割点，坐标(x_si, y_si)计算方法为：

$ \begin{array}{c} \left\{\begin{array}{l} x_{s i}=L \operatorname{sign}(i)\left|\frac{i}{\left(\frac{k-1}{2}\right)}\right|^{P_{s}} \\ y_{s i}=x_{s i} \tan \theta_{s} \end{array}\right. \\ i=-\frac{k-1}{2}, \cdots, 0, \cdots, \frac{k-1}{2} \\ L=\left\{\begin{array}{l} \operatorname{sign}\left(\tan \theta_{s}\right), \quad\left|1 / \tan \theta_{s}\right| \geqslant 1 \\ 1 / \tan \theta_{s}, \quad\left|1 / \tan \theta_{s}\right|<1 \end{array}\right. \\ \operatorname{sign}(i)=\left\{\begin{array}{l} 1, \quad i \geqslant 0 \\ -1, \quad i<0 \end{array}\right. \end{array} $

式中：L是用来限制分割点落在相平面内的范围参数; sign是符号函数; θ_s是分割点连线以X正半轴逆时针绕原点旋转的角度; P_s是幂函数指数；最后通过相邻两分割点间的垂直平分线生成一组相互平行的直线，这些平行线对相平面实现了区域划分，同时指定在同一区域内的网格点所对应的模糊输出规则then结果相同。根据PID控制器参数K_P、K_I、K_D调整的一般性原则^[13]，可将ΔK_P的相平面区域规则划分结果从左到右指定为NL、NM、NS、ZE、PS、PM、PL。模糊规则由θ_s和P_s获得，可将其定义为模糊规则的特征参数，即：设计模糊规则时所需指定的2个输入变量e和e_c的7×7个if条件可转换成指定2个特征参数，可见，模糊规则的设计过程变得大为简化。如果运行遗传算法程序随机生成P_s=2和θ_s=30°，则ΔK_P模糊规则相平面划分结果如图 2所示，相应获得ΔK_P的模糊规则如表 2所示。同样亦可由特征参数获得ΔK_D的模糊规则表。

3.2 遗传算法混合编码方法

运用遗传算法通常采用将每一个设计参数直接用个体染色体基因位表示的直接编码方法，模糊PID控制器的设计参数直接编码共需要162位基因位。为缩短直接编码时个体染色体长度，本文给出一种将特征参数和设计参数混合编码到染色体基因位上的混合编码策略：2个输入变量e和e_c，3个输出变量ΔK_P、ΔK_I、ΔK_D的隶属函数采用特征参数P_m和G编码，对应基因位a₁-a₅和a₆-a₁₀；ΔK_D与ΔK_P 2个输出变量的模糊规则用特征参数P_s和θ_s编码，对应基因位a₁₁-a₁₂和a₁₃-a₁₄；ΔK_I的模糊规则采用直接编码，用数字1~7分别表示模糊规则表 1中的NL-PL，则表中ΔK_I的模糊规则设计参数按列可表示为5544477 3333436 1123456 1223456 1233455 6333455 4444477共49位数字，对应基因位a₁₅-a₆₃，混合编码方法最终形成个体染色体总计63个基因位，比直接编码的染色体长度大为缩短。2种编码方法的染色体基因位参数如表 3所示。

4 仿真与结果分析

运用上述设计的遗传算法优化模糊PID控制器方法进行关节机器人轨迹跟踪控制仿真时，需先获得关节机器人各关节期望关节角度q_d，然后运用Matlab/Simulink仿真环境搭建关节机器人模糊PID控制仿真系统，再设置遗传算法参数，以及根据表 3所示的个体染色体各基因位对应的参数取值范围，即可按照第3节所述的遗传算法优化关节机器人轨迹跟踪控制模糊PID控制器算法步骤完成关节机器人轨迹跟踪模糊PID控制器的遗传算法优化。

选取机器人末端笛卡尔空间期望轨迹如下：

$ \left\{\begin{array}{l} x=0.25+0.05 t^{2} \\ y=0.10+0.1 t \sin \left(t 3 {\rm{ \mathsf{ π} }} / T_{c}\right) \\ z=-0.10+0.02 t \end{array}\right. $

(5)

式中：t是时间；T_c是周期时长；x、y、z分别为笛卡尔空间的坐标值。选取已知D-H参数如表 4所示的六关节机器人，其中a_i是关节i的连杆长度，α_i是关节i和关节i+1间的连杆夹角，d_i是关节i-1和关节i间的连杆偏置距离，θ_i是关节i的关节转角。

根据式(5)的曲线轨迹方程和机器人D-H参数，通过机器人逆运动学方程求得关节空间的期望关节角。采用LM(Levenberg-Marquardt，LM)方法^[14]求得逆运动学方程解，获得机器人6个关节角期望轨迹如图 3所示。

设置步长Δt=0.01 s，总时长T=10 s。运用关节机器人模糊PID控制仿真系统可得到各个时刻各关节角的仿真值。由各个时刻各关节角仿真值，运用机器人运动学方程可求解各个时刻对应的笛卡尔空间末端位置坐标，计算其与式(5)所示的笛卡尔空间理想曲线轨迹在相同时刻的欧氏距离，作为式(3)中的e(t)，得到目标函数值J。完成关节机器人轨迹跟踪模糊PID控制器的遗传算法优化，设置遗传算法参数如下：种群大小30，轮盘赌选择方式，复制比例为0.2，变异概率为0.01，两点交叉方式，交叉概率0.8，迭代代数为50。

将混合编码和直接编码方法分别应用于关节机器人模糊PID控制器遗传算法优化。设置2种编码方法的三角形隶属函数中心顶点值最大为10，最小为0.01；每个输出变量的模糊规则均为49条；直接编码和混合编码染色体基因位对应参数如表 3所示。运行遗传算法优化模糊PID控制器程序，得到如图 4所示的遗传算法迭代图，由图可知，混合编码的遗传算法在寻优过程中比直接编码更快地接近最优解，说明采用混合编码具有较高的寻优效率；直接编码遗传算法在迭代30代后最小目标函数值基本不变化，而混合编码遗传算法在30代后最小目标函数值仍有下降，说明采用混合编码方法不易过早收敛；混合编码方法的最小目标函数值小于直接编码方法的最小目标函数值，说明采用混合编码的遗传算法寻优结果更接近最优解。

将未经遗传算法优化的模糊PID控制器(FPID)、经直接编码的遗传算法优化后的模糊PID控制器(GAdirect-FPID)和经混合编码的遗传算法优化后的模糊PID控制器(GAhybrid-FPID)分别应用于六关节机器人轨迹跟踪，获得各关节在关节空间的轨迹跟踪偏差如图 5所示，由图可见，3种控制器均具有较好的快速响应性，但采用GAhybrid-FPID控制器的超调更小，震荡次数更少，鲁棒性更好。进一步根据式(3)可计算得到机器人在关节空间和笛卡尔空间的目标函数值如表 5所示，该值表征了仿真轨迹与期望轨迹的偏差，可见采用GAhybrid-FPID控制器具有更小的轨迹跟踪偏差，笛卡尔空间轨迹跟踪偏差分别比GAdirect-FPID控制器、FPID控制器小9.3%和16.6%。为进一步验证3种控制器的抗干扰能力，在t=0 s时刻分别在关节1~6处增加随机扰动信号。随机扰动信号服从均值μ=5，方差σ²=5的正态分布，其概率密度函数为f(x)=$\frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}} \sigma }}{{\rm{e}}^{ - \frac{{{{\left( {x - \mu } \right)}^2}}}{{2{\sigma ^2}}}}}$。同样可获得6个关节在有扰动时的关节空间轨迹跟踪偏差图和笛卡尔空间机器人末端轨迹跟踪偏差图，依然获得了如下结论：采用GAhybrid-FPID控制器具有更好的抗干扰能力。

5 结论

1) 关节机器人轨迹跟踪模糊PID控制器的遗传算法仿真和优化结果均表明，当设置相同的遗传算法参数时，相比于直接编码方法，隶属函数和模糊规则的设计参数与特征参数混合编码方法不仅简化了控制器设计过程，而且具有更高的寻优效率和更好的寻优结果，且不易过早收敛。

2) 经混合编码的遗传算法优化后的模糊PID控制器具有更好的控制效果和更小的轨迹偏差，显示了基于混合编码的遗传算法优化的模糊PID控制器在机器人轨迹跟踪控制中的有效性和优越性。

目前尚未开展关节机器人轨迹跟踪的混合编码遗传优化模糊PID控制器的实验论证，今后将工作重点转向于实验验证工作。