一種工業機器人多目標軌跡優化算法

李 琴，賈英崎，黃玉峰，李 剛，葉 闖

（1.西南石油大學機電工程學院，四川 成都 610500；2.中國石油天然氣集團有限公司東方地球物理勘探有限責任公司，河北 涿州 072750）

近幾十年來，隨著自動化技術的飛速發展，工業機器人在生產線上的應用愈加廣泛。軌跡規劃是機器人運動控制的基礎，良好的軌跡規劃可以提高機器人的工作效率、降低機器人的能耗以及延長機器人的使用壽命。

目前，針對機器人軌跡規劃優化的研究主要集中于采用多種算法相結合或改進原算法來獲取最優的規劃軌跡。例如：居鶴華等［1］根據機械臂的運動學約束，采用遺傳算法對3-5-3多項式進行軌跡優化，但是3-5-3多項式存在加加速度突變的問題，且傳統遺傳算法容易使解陷入局部最優。Huang等［2］以5次非均勻有理B樣條（non-uniform rational B-splines，NURBS）曲線作為工業機器人的軌跡規劃曲線，并采用非支配排序遺傳算法-II（non-dominated sorting genetic algorithm-II，NSGA-II）對運動時間和沖擊磨損這2個目標進行優化，結果表明基于5次NURBS曲線構造的軌跡曲線較平滑，但未解決NSGA-II的不足，最終的結果仍易陷入局部最優。Rout等［3］提出了一種成就標量化函數與NSGA-II相結合的混合多目標算法，并對六軸工業機器人進行多目標軌跡優化，但采用的三次多項式并未實際解決機器人加速度、加加速度突變的問題。

基于此，筆者針對文獻［1，3］中機器人加加速度不連續的問題，以工業機器人為研究對象，采用5次NURBS曲線作為其軌跡規劃曲線，并結合布谷鳥搜索（cuckoo search，CS）算法中的Levy飛行策略與NSGA-II，提出了一種用于機器人多目標軌跡優化的混合算法（簡稱為CSNSGA-II），并對初始種群進行Tent混沌優化處理，最終得到分布均勻的Pareto最優解，旨在提高算法的收斂速度以及解決NSGA-II易陷入早熟的問題。

1 工業機器人多目標軌跡優化模型

在工業生產中，機器人通常按提前預設的軌跡進行移動，通過求解逆運動學方程得到機器人各關節的軌跡點位置矢量pi=[pi，1pi，2…pi，N]（i=0，1，…，n；3≤N≤6）。其中：i為關節的軌跡點，共n+1個；N為關節數，即自由度數。各關節的軌跡點構成了工業機器人的運行軌跡。

通常情況下，工業機器人在作業過程中大多會重復同一種動作。因此，為提高機器人的工作效率、降低機器人的能耗以及延長機器人的使用壽命，基于其運動學約束（速度、加速度和加加速度約束），定義目標函數［4］，可表示為：

其中，約束條件為：

式中：S1(T)為機器人的總運動時間；S2(A)為機器人各關節加速度的平均值，用于衡量關節的能耗損失；S3(J)為機器人各關節加加速度的平均值，用于衡量關節的沖擊損耗；ti為經過第i個軌跡點時所對應的時間；am、jm分別為機器人第m個關節的加速度和加加速度；vi，m、ai，m和ji，m分別為第i個軌跡點處機器人第m個關節的速度、加速度和加加速度；vmax、amax和jmax分別為機器人關節的速度、加速度和加加速度的最大值。

2 基于NURBS曲線的工業機器人軌跡規劃

2.1 NURBS曲線的數學表示

根據工業機器人某一關節經過的軌跡點位置pi（i=0，1，…，n），采用k次NURBS曲線來規劃其軌跡曲線。k次NURBS曲線可表示為分段有理多項式函數形式：

式中：P(u)為分段有理多項式函數；ωi為NURBS曲線的權因子；di為NURBS曲線的控制頂點位置；Ni，k(u)為k次NURBS曲線的B樣條有理基函數。

Ni，k(u)是由節點矢量U=[u0u1…un+2k]按Cox De Boor遞推公式推導得到的（規定0/0=0），可表示為：

設給定的軌跡點共有n+1個，為使機器人各關節經過給定的軌跡點，使k次NURBS曲線各分段的連接點與軌跡點一一對應，則需反求NURBS曲線的控制頂點。k次NURBS曲線由n+k個控制頂點與n+2k+1個節點組成。對于節點矢量U=[u0u1…un+2k]，本文采用規范節點矢量，即令u0=u1=…=uk=0，un+k=un+k+1=…=un+2k=1，其他節點均采用積累弦長參數化法對時間間隔矢量進行歸一化求解獲得：

式中：Δtj為軌跡點間的時間間隔。

2.2 5次NURBS曲線的控制頂點反算

由于3次NURBS曲線存在加加速度不連續變化的情況，本文采用5次NURBS曲線作為工業機器人的軌跡規劃曲線。但是，5次NURBS曲線的插值數據量龐大，因此本文主要推導NURBS曲線的矩陣表達式，以簡化計算過程。

基于上文推導的節點矢量U=[u0u1…un+2k]和軌跡點位置pi(i=0，1，…，n)，根據節點矢量、軌跡點與控制頂點的關系，反求5次NURBS曲線的控制頂點。由B樣條曲線的局部支撐性可知，5次NURBS曲線在節點區間ui≤u≤ui+1內可表示為：

為方便計算，令第i（i≠0）段曲線為：

由Cox De Boor遞推公式推導得到Ni中第1行的元素，分別為：

規定：

根據5次NURBS曲線的插值要求，可得：

其中，對于Qi(t)，當u=ui時，t=0；u=ui+1時，t=1。聯立式（5）和式（6）可得：

令 ：ai+2=n11ωi，bi+2=n12ωi+1，ci+2=n13ωi+2；ei+2=n14ωi+3，fi+2=n15ωi+4；qi=(ai+2+bi+2+ci+2+ei+2+fi+2)pi，則式（7）可以簡化為：

由式（8）構成的矩陣可以得出n+1個線性方程組，而所求的控制頂點為n+5個，根據線性方程組有唯一解的條件，還須根據邊界條件確定4個附加方程。指定啟、停速度為vs、ve以及啟、停加速度為as、ae，并以啟、停位置處5次NURBS曲線的一階和二階導數作為邊界條件構建方程。

k次NURBS曲線的求導公式為：

對5次NURBS曲線進行一階、二階求導，可得啟、停時刻的速度vs、ve和加速度as、ae：

由于直接求導法的計算量較大且過程復雜，本文通過德布爾遞推算法計算各階導數的控制點。NURBS曲線上某一點的r階導數P(r)(u)為：

結合式（10）和式（11），對上文的n+5個線性方程進行求解，由此得到5次NURBS曲線控制頂點的位置矢量：

3 CSNSGA-II構建

NSGA-II是目前在多目標優化問題中應用最為廣泛的算法之一，其是在NSGA的基礎上引入快速非支配排序方法，并結合精英策略，利用其對某些測試函數進行優化可收斂到真正的最優Pareto前沿。但隨著多目標優化問題的復雜性增加，傳統的NSGA-II出現了不足之處，主要體現在以下2個方面［5］：1）易出現局部早熟的情況；2）搜索后期會出現局部搜索能力不足等。因此，須對NSGA-II進行改進，以獲取更加真實且均勻的Pareto前沿。

CS算法是一種新型的啟發式算法，由Yang和Deb于2009年提出［6］。研究表明，CS算法具有局部和全局搜索性能，基于其獲得的最優解大多優于遺傳算法和粒子群優化算法［7］，且已有較多的學者將CS算法應用于機器人軌跡優化［8-10］。基于此，提出了一種改進CS算法與NSGA-II相結合的CSNSGA-II，其中CS算法的Levy飛行策略可以幫助NSGA-II跳出局部最優，找到真實的Pareto前沿。

3.1 基于Tent混沌映射的種群初始化

為提高NSGA-II對工業機器人軌跡的優化性能，采用混沌算法對種群進行初始化。文獻［4］采用一維logistic映射對種群進行初始化，但是一維logistic映射并不具有最好的遍歷性，因此本文引入Tent混沌映射來產生初始種群。Tent混沌映射相比于logistic映射具有更好的均勻性與隨機性［11］，可以提高算法的全局尋優能力。Tent混沌映射的表達式如下：

式中：β為混沌系數，取β=0.5；γ=1，2，…，v，v為決策變量數。

將混沌種群映射到對應的決策變量空間（tγmin，tγmax）中，得到初始種群的變量：

3.2 約束不可行度處理

在NSGA-II處理約束優化問題時，常用到對不可行度的處理，以保留與可行解接近的解，舍棄與可行解相差較大的解。為了使不等式約束轉化為等式約束［12］，定義不可行度Cx（x=1，2，…，G，G為當前種群中解的數量）：

式中：gy、lz分別為不等式約束和等式約束；L、I分別為多目標優化模型中不等式和等式約束的數量。

不可行度Cx表示超過約束條件的距離，當個體在約束條件內時，Cx=0。

同時定義約束違規閾值Cxmean：

超過約束違規閾值的解定義為不可行解，滿足約束違規閾值的解定義為可行解。

3.3 自適應Levy飛行策略

傳統CS算法中的Levy飛行策略存在種群多樣性不易調節、搜索方式不能合理控制等問題，導致該算法的精度不高及收斂速度較慢，因此本文提出了一種自適應Levy飛行策略，其表達式如下：

式中：Xs，δ為第s只布谷鳥的第δ代鳥巢的位置；τ為動態慣性權重；α為變步長因子；⊕表示矩陣之間的各元素對應相乘，Levy(λ)為Levy隨機搜索策略。

Levy隨機搜索策略可表示為：

式中：U、V服從標準正態分布，μ=1.5。

動態慣性權重τ越大，CS算法的全局搜索能力越好；相反，τ越小全局搜索能力越差。因此采用自動調整慣性權重的方法，以使算法初期即具有較好的搜索能力。τ的表達式為：

式中：τmax為慣性權重的最大值；τmin為慣性權重的最小值；N()為服從(0，1)上均勻分布的隨機數。

變步長因子α可表示為：

其中：

式中：Xs，best為隨機選取的當前Pareto前沿的最優解；δmax為最大迭代次數。

3.4 CSNSGA-II流程

本文所提出的CSNSGA-II的流程如圖1所示，其中遺傳操作中的交叉算子采用獨立點交叉算子，變異算子采用多項式變異算子。

圖1 CSNSGA-II流程Fig.1 Flow chart of CSNSGA-II

由圖1可知，CSNSGA-II的具體步驟為：

步驟1 設定交叉、變異系數。采用Tent混沌映射生成初始化時間種群，并計算種群的適應度函數，在進行非支配排序后對排序后的解集進行擁擠距離計算，并按照擁擠距離進行排序。

步驟2 對排序后的所有解進行不可行度計算，得到可行解集W1和不可行解集V1。

步驟3 對可行解集W1進行二聯錦標賽選擇，對篩選的解進行下一代的交叉、變異操作，最終得到可行解集W2。

步驟4 采用自適應Levy飛行策略對不可行解集V1進行更新，得到新的解集V2。

步驟5 將解集W2與解集V2合并，并對新解集進行適應度計算、非支配排序和擁擠距離計算。

步驟6 重復步驟2、3和4，直到滿足結束條件。

4 6R勃朗特機器人多目標軌跡優化仿真分析

4.1 仿真參數設置

在MATLAB軟件中對6R勃朗特機器人（見圖2）進行運動學建模，并采用本文提出的CSNSGA-II對其軌跡進行優化。設置6R勃朗特機器人啟、停時刻的位置（笛卡爾坐標系中），并采用5次NURBS曲線構建其關節的軌跡，同時定義啟、停時刻的速度、加速度均為0，指定5次NURBS曲線的各項權因子均為1，初始化種群數量為200，最大迭代次數為100次。根據經驗，遺傳算法的交叉概率設為0.6，變異概率設為0.1。6R勃朗特機器人各關節軌跡點的位置如表1所示，基于MATLAB軟件可得其啟、停時刻位置分別如圖3和圖4所示，各關節的運動學約束如表2所示。

表1 6R勃朗特機器人各關節軌跡點的位置Table 1 Position of each joint trajectory point of 6R Bronte robot 單位：（°）

表2 6R勃朗特機器人各關節的運動學約束Table 2 Kinematic constraints of each joint of 6R Bronte robot

圖2 6R勃朗特機器人實物Fig.2 Physical object of 6R Bronte robot

圖3 6R勃朗特機器人啟動時刻的位置Fig.3 Position of 6R Bronte robot at the moment of start

圖4 6R勃朗特機器人停止時刻的位置Fig.4 Position of 6R Bronte robot at the moment of stop

4.2 仿真結果分析

圖5所示為分別采用CSNSGA-II、NSGA-II和多目標粒子群優化（multi-objective particle swarm optimization，MOPSO）算法對6R勃朗特機器人進行軌跡優化后得到的Pareto前沿分布結果。從圖中可以看出：相比于NSGA-II、MOPSO算法，基于CSNSGA-II得到的Pareto前沿分布得更加均勻，且更加靠近真實的Pareto前沿。

圖5 基于不同算法的6R勃朗特機器人Pareto前沿分布對比Fig.5 Comparison of Pareto front distribution of 6R Bronte robot based on different algorithms

為了進一步驗證CSNSGA-II的優勢，在基于3種算法優化后獲得的6R勃朗特機器人Pareto前沿中分別取3組解，結果如表3所示。其中：A1、B1、C1為CSNSGA-II優化后得到的3組解；A2、B2、C2為NSGA-II優化后得到的3組解；A3、B3、C3為MOPSO算法優化后得到的3組解。通過對比分析可得出，解A1支配解A2、A3，解B1支配解B2、B3，解C1、C2、C3互為非支配解。結果表明，相比于NSGA-II、MOPSO算法，基于CSNSGA-II能得到更優的Pareto前沿。

表3 基于不同算法的6R勃朗特機器人軌跡優化結果對比Table 3 Comparison of trajectory optimization results of 6R Bronte robot based on different algorithms

為了進一步分析CSNSGA-II的收斂性能，給定一個優化目標綜合函數S：

式中：ψ1、ψ2、ψ3為加權因子，在實際生產應用中，根據用戶的需求可以確定各加權因子的系數。在本文中，為了對比算法收斂的速度和收斂性能，取ψ1=ψ2=ψ3=1，η1=1，η2=100，η3=100。

利用式（21）計算得到基于不同算法的6R勃朗特機器人軌跡優化目標的歸一化加權迭代最優值，對比結果如圖6所示。圖6結果顯示，基于CSNSGA-II得到的機器人軌跡優化目標的最優值優于基于NSGA-II和MOPSO算法的，由此說明CSNSGA-II有助于跳出局部最優，且CSNSGA-II達到穩定收斂時的迭代次數比NSGA-II和MOPSO算法分別減少了63.16%和59.62%，有效提高了優化效率。

圖6 基于不同算法的6R勃朗特機器人軌跡優化目標歸一化加權迭代最優值對比Fig.6 Comparison of normalized weighted iterative optimal values of trajectory optimization objectives of 6R Bronte robot based on different algorithms

取解B1對應的時間間隔矢量Δt=［0 5.059 7 3.231 5 1.609 9 1.309 8 2.495 2 1.988 1 3.686 2］，其對應5次NURBS曲線的節點矢量U=［0000 0 0 0.261 0 0.427 8 0.510 9 0.578 4 0.707 2 0.809 7 1 1 1 1 1 1］，其對應的各關節的位置、速度、加速度和加加速度隨時間的變化曲線如圖7至圖10所示。從圖中可以看出，利用CSNSGA-II優化后，6R勃朗特機器人各關節的運動軌跡平滑且無沖擊，相比于文獻［13-15］所得結果更具有平穩性。

圖7 6R勃朗特機器人各關節的位置—時間曲線Fig.7 Position-time curve of each joint of 6R Bronte robot

圖8 6R勃朗特機器人各關節的速度—時間曲線Fig.8 Velocity-time curve of each joint of 6R Bronte robot

圖9 6R勃朗特機器人各關節的加速度—時間曲線Fig.9 Acceleration-time curve of each joint of 6R Bronte robot

圖10 6R勃朗特機器人各關節的加加速度—時間曲線Fig.10 Jerk-time curve of each joint of 6R Bronte robot

5 總 結

本文采用5次NURBS曲線作為工業機器人的軌跡規劃曲線，以運動時間、能耗、沖擊磨損這3個指標作為優化目標，并在運動學約束下進行優化求解。提出了一種基于改進CS算法與NSGA-II結合的混合算法——CSNSGA-II。通過對6R勃朗特機器人的軌跡進行優化后可知，CSNSGA-II能成功地跳出局部最優，解決了NSGA-II易陷入“早熟”的問題，同時基于該算法得到的Pareto前沿分布得更加均勻，有效提高了收斂速度。結果表明，所提出的CSNSGA-II具有較高的實用價值，可為工業機器人高效、持久、可靠地作業提供指導。