Coupling interpolation function selection method in industrial robot trajectory planning

Abstract

The invention discloses a coupling interpolation function selection method in industrial robot trajectory planning, comprising the following steps: step one: constructing an interpolation function set, and constructing an objective function that meets the requirements of time-impact optimal trajectory planning based on the interpolation function set; step two: initializing a golden eagle population, obtaining an initial population and an external archive; step three: judging whether the current number of iterations t reaches a set maximum number of iterations: if so, executing step seven; if not, executing step four: step four: updating the position of each golden eagle and calculating a fitness value; step five: using differential evolution and reverse learning strategies to cause two disturbances to the golden eagle population in turn, causing individual golden eagles to mutate, and calculating the fitness value of the individual after mutation; step six: updating the external archive, and cyclically executing step three; step seven: after the iteration ends, mapping the position of each golden eagle in the external archive to the optimal solution of the objective function, and obtaining a Pareto frontier containing the optimal solution set of the objective function.

Description

Technical Field

[0001] The invention belongs to the technical field of industrial robots, and specifically relates to a coupling interpolation function selection method in trajectory planning of an industrial robot. Background Art

[0002] As the backbone of the fourth industrial revolution, industrial robots are now widely used in various fields of manufacturing. In order to realize the functions of robots and achieve task goals, people need to plan the trajectory of robots. Practice has proved that the robot trajectory that only considers the task goals and function realization is very imprecise. Trajectory quality evaluation factors must also be considered, among which the most representative are: work efficiency, energy consumption and running stability. The work efficiency and running stability of industrial robots are two contradictory factors. The former is quantified as the total execution time for the robot to complete the task goal, and the latter is quantified as the third-order derivative of the joint position, called impact. The joint jitter caused by impact is the direct cause of the robot's trajectory smoothness, trajectory tracking accuracy and mechanical structure stress distribution. Optimizing impact is of great significance to improving the quality of trajectory planning and extending the service life of the robot. In order to shorten the total execution time and improve the work efficiency of the robot, it is necessary to increase the torque and power of the joint motor, which will inevitably increase the kinematic parameters of the robot joint, including joint impact; on the contrary, to reduce the joint impact, it is necessary to sacrifice work efficiency and extend the total execution time of the robot. Therefore, in order to enable the robot to maintain high work efficiency and run smoothly, time-impact optimal trajectory planning is required.

[0003] The basic idea of ​​trajectory planning for industrial robots is to use interpolation functions to interpolate between given passing points to generate several trajectory points. According to the planning space, it can be divided into two categories: Cartesian space trajectory planning and joint space trajectory planning. The trajectory points generated by Cartesian space trajectory planning are expressed by a series of robot end effector postures, and inverse kinematics is required to transform them into joint space to achieve control. Although the final planned end trajectory is accurate and intuitive, since each path point needs to be transformed by inverse kinematics, it takes up a lot of computing resources and it is difficult to avoid singular positions. Cartesian space trajectory planning is often used in working conditions such as milling, arc welding, and spraying that require high accuracy of the robot end trajectory. In contrast, joint space trajectory planning only requires the given passing points to be transformed by inverse kinematics, which greatly reduces the amount of calculation; at the same time, the entire interpolation process of generating trajectory points is carried out in the joint space, so singular points can be effectively avoided. Due to the high nonlinearity of the robot's forward kinematics model, the trajectory of the end effector is usually difficult to predict and accurately control. Therefore, joint space trajectory planning is generally used in working conditions such as spot welding and palletizing that do not pay attention to the robot's motion history. In order to obtain continuous velocity, acceleration and jerk at each passing point to ensure the smoothness of the robot joint trajectory, the commonly used interpolation functions are: polynomial spline, Bezier spline, B-spline, non-uniform rational B-spline (NURBS), etc. At present, these interpolation functions have been widely used.

[0004] The optimal trajectory planning of industrial robots means that after using the interpolation function to construct the robot trajectory, the optimization algorithm must be used to achieve trajectory optimization. The core idea is to obtain a set of interpolation function coefficients under the robot kinematic constraints to minimize the objective function. The time-impact optimal trajectory planning mentioned above belongs to the hybrid optimal trajectory planning problem. Among them, the time objective function is expressed as the total execution time, and the impact objective function is expressed as the square of the acceleration integrated in the time interval. The currently commonly used optimization methods for optimal trajectory planning are: considering kinematic or dynamic constraints to find the best, using optimization algorithms to find the best, and converting the optimization problem into a general model to find the best. However, in the prior art, the user selects an interpolation function in advance to construct the robot trajectory, and then uses the optimization algorithm to optimize the constructed trajectory. This system completely ignores the influence of the interpolation function itself on the optimization target, and the optimization result is interfered by human subjective choice. Summary of the invention

[0005] In view of this, an object of the present invention is to provide a method for selecting a coupling interpolation function in trajectory planning of an industrial robot, which can screen the best one among the interpolation functions to minimize the optimization target.

[0006] In order to achieve the above object, the present invention provides the following technical solutions:

[0007] A method for selecting a coupling interpolation function in trajectory planning of an industrial robot comprises the following steps:

[0008] Step 1: Construct an interpolation function set, and construct an objective function that meets the requirements of time-impact optimal trajectory planning based on the interpolation function set;

[0009] Step 2: Use the multi-objective Golden Eagle algorithm to optimize the objective function, initialize the Golden Eagle population, and obtain the initial population and external archive;

[0010] Step 3: Determine whether the current number of iterations t reaches the set maximum number of iterations: If yes, execute step 7; if not, execute step 4:

[0011] Step 4: Update the position of each golden eagle and calculate the fitness value;

[0012] Step 5: Use differential evolution and reverse learning strategies to cause two disturbances to the golden eagle population in turn, select and retain the better individuals, and calculate the fitness value of the individuals after mutation;

[0013] Step 6: Update the external archive and execute step 3 repeatedly;

[0014] Step 7: At the end of the iteration, the position of each golden eagle in the external archive is mapped to the optimal solution of the objective function, and the Pareto frontier containing the optimal solution set of the objective function is obtained.

[0015] Further, in step 1, the objective function is:

[0016]

[0017]

[0018] f2=min{J1,J2,...,J P}

[0019] The constraints are:

[0020]

[0021] Where f1 represents the execution time objective function; f2 represents the impact objective function; J k represents the square time integral of the trajectory constructed by the kth interpolation function; n represents the number of passing points; h i represents the time interval between the i-th passing point and the i+1-th passing point; F represents the number of robot joints; t f represents the number of robot joints; k represents the number of interpolation functions to be selected in the interpolation function set; Θ j represents the position of robot joint j; represents the velocity of robot joint j; represents the acceleration of robot joint j; represents the jerk of robot joint j; θ jmax represents the position constraint of robot joint j; V jmax represents the velocity constraint of robot joint j; A jmax represents the acceleration constraint of robot joint j; J jmax represents the jerk constraint of robot joint j.

[0022] Furthermore, in step 2, the method for initializing the golden eagle population is:

[0023] Use Tent mapping to generate a chaotic sequence in the range [0,1]:

[0024]

[0025] Where i = 1, 2, 3, ..., N represents the population index; j = 1, 2, 3, ..., D represents the dimension of the decision variable; α is the Tent mapping parameter; the first dimension in the chaotic sequence is a random number between [0, 1], ω i,j and ω i,j+1 represent the jth and j+1th dimension values ​​of the i-th golden eagle respectively; c n represents a chaotic sequence;

[0026] Then the j-th dimension of the i-th golden eagle in the chaotic population S1 is Initialized to:

[0027]

[0028] Among them, ub j and lb j They represent the upper and lower boundary conditions of the j-th dimension of the decision variable respectively;

[0029] Using the reverse learning strategy, we can obtain the jth dimension of the i-th golden eagle in the reverse population S2.

[0030]

[0031] The chaotic population and the reverse population S2 of S1 are merged to obtain the population S, and the fitness value of each golden eagle in the population S is calculated. The population S is graded in a gradually increasing manner using the fast non-dominated sorting method. The golden eagles with smaller non-dominated levels are preferentially selected to join the initial population. If the non-dominated levels of two golden eagles are equal, the golden eagle with a larger crowding distance is preferentially added to the initial population; until the initial population is filled; at the same time, the golden eagle with the smallest non-dominated level is added to the external archive.

[0032] Furthermore, the method of using the fast non-dominated sorting method to classify the population S is:

[0033] 21) Obtain a Pareto frontier by judging the dominance and non-dominance relationship between the golden eagles, and assign a level d to all the golden eagles in the Pareto frontier; d = 1;

[0034] 22) Determine whether all golden eagles in population S are assigned grades: if so, the population S is graded; if not, execute step 23);

[0035] 23) Eliminate the golden eagles of level d from the population S, d = d + 1;

[0036] 24) Obtain a Pareto frontier by judging the dominance and non-domination relationship among the remaining golden eagles, and assign grade d to all golden eagles in the Pareto frontier;

[0037] 25) Loop through step 22).

[0038] Furthermore, the crowding distance of each golden eagle in each level of Pareto frontier for:

[0039]

[0040] Where R represents the non-dominated level; j = 1, 2, 3, ..., nobj represents the objective function number; i represents the i-th non-dominated solution in the Pareto frontier with level R, and Respectively represent the maximum and minimum values ​​of the j-th objective function; and They represent the j-th objective function values ​​of the i-th and i-1-th non-dominated solutions in the Pareto front of level R, respectively.

[0041] Furthermore, in step 3, the method of updating the position of each golden eagle using the multi-objective golden eagle optimization algorithm is:

[0042]

[0043] in, represents the position of the i-th golden eagle at the t+1 iteration step; represents the position of the i-th golden eagle at the t-th iteration step; r1 and r2 are both random numbers between [0,1], p a and p c Respectively represent the attack coefficient and patrol coefficient of the current iteration step; A i represents the attack vector of the i-th golden eagle; C i represents the cruise vector of the i-th golden eagle; Weight represents the weight; and:

[0044]

[0045] Among them, t is the current iteration number, and T is the maximum iteration number;

[0046]

[0047] in, represents the best prey searched by the f-th golden eagle so far, and also represents the prey selected by the i-th golden eagle; X i Indicates the current position of the Golden Eagle;

[0048] C i =[rand(0,1),rand(0,1),...,C i,m ,...,rand(0,1)]

[0049]

[0050] Among them, C i,m represents the m-th dimension component of the patrol vector of the i-th golden eagle; d represents the n-dimensional hyperplane equation, which can be expressed as: Among them A i =[a1,a2,...,a n ] represents the attack vector of the i-th golden eagle, represents the prey selected by the i-th golden eagle; a j , a m Represents the attack vector A i The j-th and m-th components of , where m represents a random integer between [1,n] and rand(0,1) represents a random number between [0,1];

[0051]

[0052] in, and Represent the initial and final values ​​of the attack coefficient respectively; and Represent the initial value and final value of the patrol coefficient respectively.

[0053] Furthermore, the method of using the differential evolution algorithm to move the golden eagle population is as follows:

[0054] Use the DE / rand / 2 strategy to perform differential mutation on Golden Eagle:

[0055]

[0056] Where r1, r2, r3, r4, r5 represent random population indexes between [1, N], and r1≠r2≠r3≠r4≠r5≠i, N is the size of the golden eagle population; F is the scaling factor; V i is the mutation vector;

[0057] The crossover operation is:

[0058]

[0059] Among them, rand i,j represents a random number between [0,1]; k is a random integer between [1,D]; CR∈[0,1] is the crossover rate; U i is the test vector, also known as offspring; X i is called the target vector, also known as the parent; and:

[0060]

[0061] Among them, CR0 is the crossover constant; t is the current iteration number, and T is the maximum iteration number;

[0062] Select the action as:

[0063]

[0064] in, Represents the i-th golden eagle after mutation.

[0065] Furthermore, the principle of using the reverse learning method to disturb the golden eagle population is:

[0066]

[0067] Among them, Y i represents the reverse individual; ub and lb represent the upper and lower boundary conditions of the decision variable respectively;

[0068] Similarly, select the operation:

[0069]

[0070] Among them, X i Represents the i-th golden eagle after mutation.

[0071] Furthermore, in step 6, the external archive update method is:

[0072] If the updated new position of the Golden Eagle can dominate some members in the current external archive, these members will be deleted and the updated Golden Eagle will be added to the external archive; if the updated new position of the Golden Eagle is in a non-dominant relationship with all members in the current external archive, and it indicates that the number of members in the external archive is less than the preset external archive size M, the updated Golden Eagle will be added to the external archive; otherwise, the Golden Eagle with the smallest crowding distance in the external archive will be taken as the deletion object, and the updated Golden Eagle will replace the deleted members in the external archive.

[0073] Furthermore, the method also includes step eight: using a fuzzy comprehensive evaluation method to evaluate a potential optimal solution from the Pareto frontier, and using the interpolation function corresponding to the potential optimal solution as the interpolation function for trajectory planning of the industrial robot.

[0074] The beneficial effects of the present invention are:

[0075] The coupled interpolation function selection method in the trajectory planning of an industrial robot of the present invention constructs an interpolation function set including a plurality of interpolation functions and obtains an objective function that meets the requirements of the time-impact optimal trajectory planning based on the interpolation function set. The objective function is set as the total execution time of the robot and the square integral of the jerk. An improved multi-objective Golden Eagle algorithm (MOGEO) is used to improve the population diversity and the ability to jump out of the local optimum of the algorithm by introducing chaotic mapping, reverse learning, differential mutation and adaptive weight strategy, solve the time-impact optimal trajectory planning problem, and use fuzzy comprehensive evaluation to evaluate the Pareto optimal frontier to obtain a potential optimal solution. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] In order to make the purpose, technical solution and beneficial effects of the present invention clearer, the present invention provides the following drawings for illustration:

[0077] Figure 1 A flow chart of a method for selecting a coupling interpolation function in trajectory planning of an industrial robot according to the present invention;

[0078] Figure 2 Trajectory curve graph generated for joint space spline interpolation;

[0079] Figure 3 This is a schematic diagram of the circling trajectory of the golden eagle;

[0080] Figure 4 This is a model diagram of the BRTIRSE1506A robot and mixing operation;

[0081] Figure 5 A schematic diagram of the robot's stirring path and the collected path points;

[0082] Figure 6 The Pareto optimal frontier for the time-shock optimal trajectory planning problem for the selection of coupled interpolation functions;

[0083] Figure 7 is the joint position curve;

[0084] Figure 8 is the joint velocity curve;

[0085] Fig. 9 is the joint acceleration curve;

[0086] Fig.10 is the joint acceleration curve;

[0087] Fig.11 Comparison of the acceleration curves of cubic B-spline and quintic B-spline. DETAILED DESCRIPTION

[0088] The present invention is further described below in conjunction with the accompanying drawings and specific embodiments so that those skilled in the art can better understand the present invention and implement it, but the embodiments are not intended to limit the present invention.

[0089] like Figure 1 As shown, the coupling interpolation function selection method in the industrial robot trajectory planning of this embodiment includes the following steps.

[0090] Step 1: Construct an interpolation function set, and based on the interpolation function set, construct an objective function that meets the requirements of time-impact optimal trajectory planning.

[0091] Figure 2 Example of generating trajectories for robot joint space spline interpolation: For a certain joint J of the robot (J = 1, 2, ..., F), define the time series t = [t1, t2, ..., t n+1 ] and the joint position sequence θ=[θ1,θ2,...,θ n+1 ]; define the time interval h between adjacent passing points i =t i+1 -t i ; Define the spline function θ of the robot's i-th trajectory with respect to time i (t), where t∈[t i ,t i+1 ], satisfying Θ i (t i )=θ i , Θ i (t i+1 )=θ i+1 In the above definition, i = 1, 2, ..., n + 1. The robot joint trajectory can be expressed as:

[0092]

[0093] In the prior art, the trajectory construction in the time-impact optimal trajectory planning of industrial robots is all carried out according to a certain artificially specified interpolation function. However, if it is only for the purpose of optimization, the influence of the interpolation function itself on the optimization target cannot be excluded. For example, for a given series of passing points, there is an interpolation function A, and the impact index of the trajectory constructed by it is always greater than that of B. In particular, this embodiment couples the selection of the interpolation function with the above-mentioned dual-objective optimization problem, and does not specify the interpolation function to be used in advance, but constructs a set of interpolation functions to be selected. It is hoped that while obtaining the optimal value of the robot time-impact dual-objective optimization, it is also possible to determine a certain interpolation function belonging to the interpolation function set corresponding to the optimal value, that is, to select the optimal interpolation function while performing the robot time-impact optimization, so as to avoid incomplete optimization caused by artificial selection. Specifically, the corresponding total execution time of the trajectory constructed by each interpolation function in the set is the same, so it is only necessary to select the one with the smallest corresponding impact index as the optimal one.

[0094] Specifically, the objective function is:

[0095]

[0096]

[0097] f2=min{J1,J2,...,J P}

[0098] The constraints are:

[0099]

[0100] Where f1 represents the execution time objective function; f2 represents the impact objective function; J k represents the square time integral of the trajectory constructed by the kth interpolation function; n represents the number of passing points; h i represents the time interval between the i-th passing point and the i+1-th passing point; F represents the number of robot joints; t f represents the number of robot joints; k represents the number of interpolation functions to be selected in the interpolation function set; Θ j represents the position of robot joint j; represents the velocity of robot joint j; represents the acceleration of robot joint j; represents the jerk of robot joint j; θ jmax represents the position constraint of robot joint j; V jmax represents the velocity constraint of robot joint j; A jmax represents the acceleration constraint of robot joint j; J jmax represents the jerk constraint of robot joint j.

[0101] Regarding the interpolation function set, several widely used interpolation functions can be selected and added to the interpolation function set. This embodiment selects cubic spline, cubic B-spline and quintic B-spline.

[0102] The essence of the cubic spline curve is that several cubic polynomial curves are smoothly connected at the interpolation nodes, which requires the robot joint velocity and acceleration curves to be continuous at the interpolation nodes. Given that the second-order derivative of the cubic spline function is a linear function, the three-bending-moment method is used to express the robot joint acceleration as:

[0103]

[0104] Integrate the above equation twice and substitute the boundary condition Θ i (t i )=θ i and θ i (t i+1 )=θ i+1 The joint position expression is:

[0105]

[0106] Obviously, only and is unknown, so we only need to obtain the acceleration of the robot at each interpolation node to get the spline function of the joint position. The joint velocity spline function can be derived from the above formula as follows:

[0107]

[0108] Known velocity continuity condition get:

[0109]

[0110] in:

[0111]

[0112] The above equation can only solve n-2 unknowns, and two more equations need to be added. In addition, consider the robot's initial and final joint velocities as the first-class boundary conditions of cubic spline interpolation. Here, let Substituting in:

[0113]

[0114]

[0115] Rewritten as a matrix, we get:

[0116] AX=B

[0117] in:

[0118]

[0119] Among them, μ i , i and d i They are all process parameters used to simplify the equation and have no actual meaning.

[0120] Obviously, the matrix A is a strictly diagonally dominant matrix with a unique solution. The pursuit method is used to solve the problem and finally the cubic spline function of the robot joint position is obtained.

[0121] B-spline interpolation is a commonly used method for robot trajectory planning. After the interpolation is completed, any node can be changed without affecting the overall shape of the curve by relying on its local support. The robot joint trajectory interpolated using the B-spline function is as follows:

[0122]

[0123] Among them: k' represents the number of times; d j represents the control point; m+1 represents the number of control points; represents the basis function, and the basis function is expressed as follows based on the DeBoor–Cox recursion formula:

[0124]

[0125]

[0126] Where: e represents the node vector. When the degree k' is an odd number, the node vector is defined as follows:

[0127]

[0128] The node vector e is a time series with initial and final times repeated k'+1 times. For a given number of times k', and time t i Interpolation node θ i The key to the B-spline trajectory planning problem is to solve the control point d j , to ensure that the equation Θ(t i )=θ i holds, where i=1,2,...,n+1. This system of equations contains m+1 unknowns and n+1 equations. Solving it requires determining the relationship between m and n. Number of nodes n knot +1=n+2k'+1, considering the number of control points of the B-spline function m+1 and the number of nodes n knot +1 and the relationship between the number of times p knot=m+p+1, and the relationship between the number of interpolation nodes n+1 and the number of control points m+1 is obtained: m+1=n+1+(p-1). Therefore, for odd-order B-spline, k'-1 equations need to be added to solve the equation system. For example, for quintic B-spline interpolation, the robot joint velocity and acceleration boundary conditions need to be added as follows:

[0129]

[0130]

[0131] Where v1 represents the starting speed; v n+1 represents the terminal velocity; a1 represents the starting acceleration; a n+1 represents the terminal acceleration.

[0132] The cubic B-spline only needs to supplement the joint velocity boundary conditions. The u-order derivative of the robot joint trajectory is obtained:

[0133]

[0134] Basis Functions The u-order derivative of is solved by the following formula:

[0135]

[0136] in:

[0137]

[0138] Taking the quintic B-spline as an example, written in matrix form:

[0139] CX=D

[0140] in:

[0141]

[0142] Using the formula X=C -1 D is solved, and the required control points are substituted into the robot joint trajectory equation to obtain the robot joint position B-spline function.

[0143] By observing the coefficient matrices A and C in the equation and the interpolation function expression, it is easy to find that when the joint position sequence θ is given, the time series t is the only factor that determines the interpolation function coefficient. Therefore, it can be determined that the robot joint impact index value is also only affected by the time series distribution. i It is input as decision variables into the meta-heuristic optimization algorithm to solve the proposed optimal problem.

[0144] Step 2: Use the multi-objective Golden Eagle algorithm to optimize the objective function, initialize the Golden Eagle population, and obtain the initial population and external archive

[0145] Population initialization is closely related to the quality of the final solution of the algorithm. This embodiment refers to a chaotic reverse initialization scheme, which aims to combine chaotic mapping and reverse learning to produce a better quality initialization population, thereby further improving the diversity of the population and accelerating the convergence of the algorithm. Tent mapping has the characteristics of simple structure, uniform distribution, and good ergodicity. Tent mapping is used to generate a chaotic sequence in the range of [0,1]:

[0146]

[0147] Where i = 1, 2, 3, ..., N represents the population index; j = 1, 2, 3, ..., D represents the dimension of the decision variable; α is the Tent mapping parameter; the first dimension in the chaotic sequence is a random number between [0, 1], ω i,j and ω i,j+1 represent the jth and j+1th dimension values ​​of the i-th golden eagle respectively; c n represents a chaotic sequence;

[0148] Then the j-th dimension of the i-th golden eagle in the chaotic population S1 is Initialized to:

[0149]

[0150] Among them, ub j and lb j They represent the upper and lower boundary conditions of the j-th dimension of the decision variable respectively;

[0151] Using the reverse learning strategy, we can obtain the jth dimension of the i-th golden eagle in the reverse population S2.

[0152]

[0153] The chaotic population and the reverse population S2 of S1 are merged to obtain the population S, and the fitness value of each golden eagle in the population S is calculated. The population S is graded in a gradually increasing manner using the fast non-dominated sorting method. The golden eagles with smaller non-dominated levels are preferentially selected to join the initial population. If the non-dominated levels of two golden eagles are equal, the golden eagle with a larger crowding distance is preferentially added to the initial population; until the initial population is filled; at the same time, the golden eagle with the smallest non-dominated level is added to the external archive.

[0154] In this embodiment, the method of using the fast non-dominated sorting method to classify the population S is:

[0155] 21) Obtain a Pareto frontier by judging the dominance and non-dominance relationship between the golden eagles, and assign a level d to all the golden eagles in the Pareto frontier; d = 1;

[0156] 22) Determine whether all golden eagles in population S are assigned grades: if so, the population S is graded; if not, execute step 23);

[0157] 23) Eliminate the golden eagles of level d from the population S, d = d + 1;

[0158] 24) Obtain a Pareto frontier by judging the dominance and non-domination relationship among the remaining golden eagles, and assign grade d to all golden eagles in the Pareto frontier;

[0159] 25) Loop through step 22).

[0160] Specifically, the crowding distance C of each golden eagle in each level of Pareto frontier i R for:

[0161]

[0162] Where R represents the non-dominated level; j = 1, 2, 3, ..., nobj represents the objective function number; i represents the i-th non-dominated solution in the Pareto frontier with level R, and Respectively represent the maximum and minimum values ​​of the j-th objective function; and They represent the j-th objective function values ​​of the i-th and i-1-th non-dominated solutions in the Pareto front of level R, respectively.

[0163] The pseudo code for initializing the golden eagle population is given in Algorithm 1.

[0164]

[0165] Step 3: Determine whether the current number of iterations t reaches the set maximum number of iterations: if so, execute step 7; if not, execute step 4.

[0166] Step 4: Update the position of each golden eagle and calculate the fitness value

[0167] Inspired by the hunting behavior of the golden eagle, Mohammadi et al. proposed a new meta-heuristic algorithm called the Golden Eagle Algorithm (GEO) in 2020. With its fewer parameters, excellent solution speed and convergence performance, it has been gradually applied to various engineering problems. The cruise vector and attack vector are introduced to simulate the current hunting behavior of the golden eagle, such as Figure 3As shown. In the early stage of hunting, the golden eagle shows a high cruising-low attack tendency and hovers over the prey, with the purpose of searching for as many prey as possible; at the end of hunting, it shows a low cruising-high attack tendency to capture the prey. Obviously, the above two stages represent the search and development stages in the metaheuristic algorithm. The transition of the golden eagle from high cruising-low attack tendency to low cruising-high attack tendency during the hunting process is equivalent to the transition of the algorithm from the search stage to the development stage. Through this process, the intelligent balance of GEO search and development can be achieved and it is suitable for solving various global optimization problems. Multi-objective Golden Eagle Optimization Algorithm (MOGEO) based on GEO is proposed.

[0168] The attack vector can be expressed as the line between the current position of the golden eagle and the current best prey. MOGEO's development capabilities are mainly reflected through attack vectors. The attack vector is calculated as follows:

[0169]

[0170] in, represents the best prey searched by the f-th golden eagle so far, and also represents the prey selected by the i-th golden eagle; X i Indicates the current Golden Eagle position.

[0171] The cruise vector is perpendicular to the attack vector and tangent to its circling trajectory. For an n-dimensional cruise vector, it is located in the hypertangent plane of the trajectory. MOGEO's search capability is mainly reflected by the cruise vector. The cruise vector is calculated as follows:

[0172]

[0173] C i =[rand(0,1),rand(0,1),...,C i,m ,...,rand(0,1)]

[0174] Among them, C i,m represents the m-th dimension component of the patrol vector of the i-th golden eagle; d represents the n-dimensional hyperplane equation, which can be expressed as: Among them A i =[a1,a2,...,a n ] represents the attack vector of the i-th golden eagle, represents the prey selected by the i-th golden eagle; a j , a m Represents the attack vector A i The j-th and m-th components of , where m represents a random integer between [1,n] and rand(0,1) represents a random number between [0,1].

[0175] The step vector is calculated as follows:

[0176]

[0177] Where: r1 and r2 are random numbers between [0,1], p a and p c Respectively represent the attack coefficient and patrol coefficient of the current iteration step, and their size directly affects the patrol and attack tendency of the Golden Eagle. a and p c Updated as follows:

[0178]

[0179] in, and Represent the initial and final values ​​of the attack coefficient respectively; and Represent the initial value and final value of the patrol coefficient respectively.

[0180] Last updated Golden Eagle position:

[0181]

[0182] By introducing inertia weight, the algorithm's global search and local optimization capabilities are further balanced. That is, a larger inertia weight is used in the early iteration to expand the search range of the population, and a smaller inertia weight is used in the later iteration to improve the optimization capability. Although MOGEO linearly changes the attack coefficient p during the iteration process a and the cruise coefficient p c The values ​​of and make them increase and decrease respectively, which reflects the process of the algorithm from search to development as the iteration proceeds. However, in reality, this transformation is not obvious. When the algorithm is close to convergence, the cruise coefficient p c Still has a large value, that is, the golden eagle finally only hovers near the best prey and cannot capture the prey, and the golden eagle does not want to search for other prey at this time. In order to solve the problem, this embodiment introduces a nonlinear adaptive weight coefficient to replace the above formula for updating the golden eagle position:

[0183]

[0184] in, represents the position of the i-th golden eagle at the t+1 iteration step; represents the position of the i-th golden eagle at the t-th iteration step; r1 and r2 are both random numbers between [0,1], p a and p c Respectively represent the attack coefficient and patrol coefficient of the current iteration step; A i represents the attack vector of the i-th golden eagle; C i represents the cruise vector of the i-th golden eagle; Weight represents the weight; and:

[0185]

[0186] Among them, t is the current iteration number, and T is the maximum iteration number.

[0187] It is easy to see that the weight decreases nonlinearly between [0,1] and the rate of decrease gradually slows down when it approaches 0, thus strengthening the algorithm's transformation from search to development and further achieving a balance between the two.

[0188] Step 5: Use differential evolution and reverse learning strategies to cause two disturbances to the golden eagle population in turn, select and retain the better individuals, and calculate the fitness value of the individuals after mutation.

[0189] The differential evolution algorithm was first proposed by Storn et al. The algorithm introduces disturbances in the population to force the parent generation to mutate and cross over to produce offspring, and then uses the greedy criterion to select and retain the better offspring to complete the iterative optimization process. In MOGEO, the golden eagle only updates its position based on its own memory, and has little correlation with other golden eagles in the population, which makes MOGEO prone to falling into local optimality. Therefore, this embodiment introduces a differential evolution strategy to disturb the golden eagle population after the golden eagle's position is updated, and then integrates reverse learning to cause secondary disturbances, that is, using a mixed individual mutation strategy to help the golden eagle jump out of the local optimality, while also playing a role in expanding the diversity of the population.

[0190] The method of using differential evolution algorithm to move the golden eagle population is:

[0191] In order to cause as large a population disturbance as possible, the DE / rand / 2 strategy was used to perform differential mutation on the golden eagle:

[0192]

[0193] Where r1, r2, r3, r4, r5 represent random population indexes between [1, N], and r1≠r2≠r3≠r4≠r5≠i, N is the size of the golden eagle population; F is the scaling factor, and F=0.5 is taken in this paper; V i is the mutation vector;

[0194] The crossover operation is:

[0195]

[0196] Among them, rand i,j represents a random number between [0,1]; k is a random integer between [1,D]; CR∈[0,1] is the crossover rate; U i is the test vector, also known as offspring; X i It is called the target vector, also known as the parent. The first term of this formula is used to ensure that the test vector U i At least one dimension comes from the mutation vector Vi In the early stage of iteration, a larger value of the crossover rate is assigned to cause greater population disturbance, and in the later stage of iteration, a smaller value is assigned to accelerate the convergence of the algorithm. And:

[0197]

[0198] Among them, CR0 is the crossover constant, and this embodiment takes CR0=0.2; t is the current iteration step number, and T is the maximum iteration step number. Different from the selection operation of DE, this embodiment selects whether the offspring can replace the parent by judging the dominance relationship between the offspring and the parent. First, the fitness values ​​of the offspring and the parent are calculated respectively, and then the selection operation is performed as shown below:

[0199]

[0200] in, Represents the i-th golden eagle after mutation.

[0201] If the offspring can dominate the parent, the offspring will replace the parent; if the offspring does not dominate the parent, the parent has a 50% chance of being replaced by the offspring; if the offspring is dominated by the parent, the parent will be retained. So far, the individual mutation operation based on the differential evolution algorithm has been completed. In the future, reverse learning will be introduced to generate secondary perturbations of the population and individual selection strategies.

[0202] Specifically, the principle of using the reverse learning method to disturb the golden eagle population is:

[0203]

[0204] Among them, Y i represents the reverse individual; ub and lb represent the upper and lower boundary conditions of the decision variable respectively;

[0205] Similarly, select the operation:

[0206]

[0207] Among them, X i Represents the i-th golden eagle after mutation.

[0208] Step 6: Update the external archive and repeat step 3.

[0209] When the external archive is full, the original algorithm handles it as follows: Calculate the crowding distance C of each golden eagle i , using S i =1-C iThe sparse distance is obtained, and then a roulette algorithm is used to randomly select a golden eagle with a larger sparse distance, that is, the golden eagle in the denser area is discarded. The problem with this method is that the discarded golden eagle may not be in the most dense area, that is, a golden eagle in the secondary dense area is deleted, and the cycle is repeated, which eventually leads to poor distribution uniformity of the Pareto optimal solution obtained. In this embodiment, the update method of the external archive is: if the new position of the golden eagle after the update can dominate some members in the current external archive, then these members are deleted and the updated golden eagle is added to the external archive; if the new position of the golden eagle after the update is in a non-dominated relationship with all members in the current external archive, and it indicates that the number of members in the external archive is less than the preset external archive size M, the updated golden eagle is added to the external archive; otherwise, the golden eagle with the smallest crowding distance in the external archive is taken as the deletion object, and the updated golden eagle replaces the deleted members in the external archive, thereby ensuring the uniformity of the algorithm solution to a certain extent.

[0210] Step 7: At the end of the iteration, the position of each golden eagle in the external archive is mapped to the optimal solution of the objective function, and the Pareto frontier containing the optimal solution set of the objective function is obtained.

[0211] Step 8: Use the fuzzy comprehensive evaluation method to evaluate a potential optimal solution from the Pareto frontier, and use the interpolation function corresponding to the potential optimal solution as the interpolation function for industrial robot trajectory planning.

[0212] In this embodiment, MOGEO is used to solve the trajectory optimization problem of industrial robots. The decision variable is the time interval of each path segment, which is expressed as H = [h1, h2, ..., h n ], and created the objective function and constraints. It is worth noting that every time the algorithm calculates the fitness, it needs to record the corresponding optimal interpolation function. Similar to the population index, this value represents an additional label, accompanying the golden eagle to participate in the Pareto solution update. Algorithm 2 gives the pseudo code of the coupling interpolation function selection method in the trajectory planning of the industrial robot in this embodiment.

[0213]

[0214]

[0215] The coupling interpolation function selection method in the industrial robot trajectory planning of this embodiment is verified with reference to a specific example.

[0216] 1. Robot mixing operation

[0217] This embodiment intends to use the six-axis robot BRTIRSE1506A produced by Bronte as a model to conduct simulation experiments. The simulation condition is the robot mixing operation in an automated production line, such as Figure 4As shown in the figure, the robot's end effector is a wooden stirring rod, which is placed in a barrel with a cylindrical inner cavity. The barrel contains solid material, which is poured in advance by the previous process and piled into a cone shape. The robot's task is to use the stirring rod to level the material in the barrel in preparation for the next process. In order to clearly show the stirring path, the solid material has been hidden and 1 / 4 of the volume has been cut off along the circumference of the barrel. The robot's tool coordinate system is set at the end of the wooden stick, and the stirring path is set as a closed curve in space. Eight path points are collected in Cartesian space as the robot's trajectory planning points. Figure 5 As shown in Table 1, its expression in the joint space is shown in Table 1. It can be seen that the robot has a very small range of motion under this condition, and the mixing requires high-speed and efficient operation, so the execution time between each path point must be reasonably allocated to prevent excessive acceleration. The joint speed and acceleration at the start and end positions are set to 0. The kinematic constraints of the robot are set as shown in Table 2. Regarding the simulation settings, the population size is set to 100, the external file capacity is set to 100, and the maximum number of iterations is set to 200.

[0218] Table 1 Joint passing points

[0219]

[0220] Table 2 Robot kinematic constraints

[0221]

[0222]

[0223] 2. Fuzzy comprehensive evaluation

[0224] In order to screen out potential optimal solutions from the Pareto frontier, decision makers need to formulate corresponding standards. Commonly used methods include weighted method and fuzzy comprehensive evaluation method. The weighted method converts the multi-objective problem into a single-objective problem by assigning weights to the objective function. The disadvantage of this is that it is highly artificial and has poor objectivity. Moreover, when the orders of magnitude between the objective functions differ greatly, the objective function with a smaller order of magnitude is easily weakened. Therefore, it is necessary to introduce appropriate elastic parameters to unify the order of magnitude. This embodiment adopts the fuzzy comprehensive evaluation method, which can make a comprehensive evaluation of optimization problems affected by multiple objective functions. First, all optimal solutions are normalized using the membership function, and then the fuzzy membership of each solution is calculated, and the one with the largest fuzzy membership is used as the final solution to be selected. The linear fuzzy membership function is defined as follows:

[0225]

[0226] Where: f j (i) represents the i-th solution of the j-th objective function in the Pareto frontier, f jmax and fjmin Represent their maximum and minimum values ​​respectively. This function is a monotonically decreasing function, that is, a smaller objective function value corresponds to a larger fuzzy membership. The fuzzy membership is calculated by the following formula:

[0227]

[0228] 3. Simulation results and discussion

[0229] For the time-impact optimal trajectory planning problem of robot mixing operation, the Pareto optimal frontier obtained by MOGEO is as follows: Figure 6 As shown in , this is a finite set of possible optimal solutions to the optimization objective. Among them, solution A can obtain the optimal time index and the worst impact index; solution B can obtain the optimal impact index and the worst time index. From solution A to solution B, the time index changes from good to bad, and the impact index changes from bad to good, which shows the fact that time and impact are contradictory.

[0230] The 54th solution is calculated as the potential optimal solution using the linear fuzzy membership function and fuzzy membership. Figure 6 Solution C in the figure. The final planning results are shown in Table 3. Then, the optimal execution time is 4.8190s, the optimal impact index value is 108653.79, and the optimal interpolation function is the quintic B-spline, indicating that the quintic B-spline interpolation robot mixing trajectory should be selected in the end. On the other hand, the quintic B-spline interpolation has a moderate amount of calculation, can achieve continuous acceleration, and inserting two unknown intermediate points at the start and end positions can also avoid sudden acceleration changes. It is itself a reasonable and widely used trajectory constructor, so the optimization result is reasonable. Figure 7-10 The motion function curve of the joint space robot's optimal trajectory is shown, including position, velocity, acceleration and jerk. The interpolation nodes are marked with hollow circles. Substituting their corresponding moments into the joint trajectory function, the calculated values ​​are exactly equal to the interpolation nodes, indicating that the optimal trajectory of the robot has no position error. The interpolation curves are all high-order continuous and differentiable, where the velocity and acceleration at the initial and final positions of the robot joints are controlled manually. Fig.10 As shown in Table 4, the jerk curves of each joint of the robot are smooth and continuous without mutation, which ensures the stability of the robot in completing the mixing operation and the smoothness of the formed trajectory. The maximum absolute values ​​of the velocity, acceleration and jerk of each joint are calculated as shown in Table 4. Compared with the data in Table 2, the obtained time-impact optimal trajectory meets the robot kinematic constraints.

[0231] In order to verify the effectiveness of the method proposed in this embodiment, it is necessary not only to output the optimal interpolation function, but also to require that for a given robot path, no matter how many times the algorithm is calculated, the type of interpolation function finally obtained should be consistent. Therefore, this embodiment uses MOGEO to perform 100 independent solutions, and the first 30 groups of calculation results are recorded in Table 5. 82% of the calculation results of the 100 times are the best for the quintic B-spline, and the other 18% are for the cubic spline and cubic B-spline, resulting in inconsistent optimization results.

[0232] Table 6 shows the time interval sequence obtained from the first 30 operations. In the last column, the serial number 1 indicates that the optimal interpolation function is cubic spline, 2 indicates cubic B-spline, and 3 indicates quintic B-spline. The solutions with results 1 and 2 are bolded. After layer-by-layer comparison, it is found that the solutions h1 and h8 with results 1 and 2 are significantly smaller than, and h3 and h6 are significantly larger than, the solution with result 3. There is no obvious difference in other values. Fig.11 The following table shows the comparison of the acceleration curves of the optimal trajectory and the quintic B-spline trajectory in a certain calculation of the cubic B-spline optimality. Table 7 shows the impact indicators of the two trajectories. Fig.11 As shown in Table 7, when using the quintic B-spline interpolation trajectory, the increase in the jerk of the robot joints in these two intervals is significantly greater due to the reduction of the interpolation spacing between the beginning and the end (the joints pass through the first and last segments of the path faster) compared to the spline of degree 3. Although increasing the interpolation spacing in the middle position can reduce the jerk to a certain extent, it is still a drop in the bucket compared to the increase. Fig.11 ,The impact of robot joints 1, 2, 4 and 6 near the endpoints in the quintic B-spline trajectory is significantly greater than that of the cubic B-spline, especially joints 4 and 6. From Table 7, the impact index of the former is even more than twice that of the latter. Therefore, in this case, MOGEO will choose cubic spline or cubic B-spline as the optimal interpolation function.

[0233] The fundamental reason for the appearance of non-quintic B-spline results is that MOGEO is disturbed by many random factors, resulting in different Pareto optimal frontiers obtained by the algorithm each time, and the quintic B-spline optimal solution expected by the user may be missed, which is unavoidable. For the mixing working condition of this embodiment, Table 5 shows that the optimization result is quintic B-spline, which still accounts for the majority. Therefore, it is believed that the proposed method using MOGEO can obtain a relatively stable output, which verifies the feasibility of selecting the coupled interpolation function for the robot time-impact optimal trajectory planning.

[0234] Table 3 Optimization results of MOGEO

[0235]

[0236] Table 4 Extreme values ​​of robot joint kinematic parameters

[0237]

[0238] Table 5 Results of multiple solutions to the optimal interpolation function

[0239]

[0240]

[0241] Table 7 Comparison of joint impact indicators of the two trajectories

[0242]

[0243] The above-described embodiments are only preferred embodiments for fully illustrating the present invention, and the protection scope of the present invention is not limited thereto. Equivalent substitutions or changes made by those skilled in the art based on the present invention are within the protection scope of the present invention. The protection scope of the present invention shall be subject to the claims.

Claims

1. A method for selecting a coupling interpolation function in trajectory planning of an industrial robot, characterized in that: The steps include: Step 1: Construct an interpolation function set, and construct an objective function that meets the requirements of time-impact optimal trajectory planning based on the interpolation function set; Step 2: Use the multi-objective Golden Eagle algorithm to optimize the objective function, initialize the Golden Eagle population, and obtain the initial population and external archive; The method to initialize the golden eagle population is: Use Tent mapping to generate a chaotic sequence in the range [0,1]: Where i = 1, 2, 3, ..., N represents the population index; j = 1, 2, 3, ..., D represents the dimension of the decision variable; α is the Tent mapping parameter; the first dimension in the chaotic sequence is a random number between [0, 1], ω i,j and ω i,j+1 represent the jth and j+1th dimension values ​​of the i-th golden eagle respectively; c n represents a chaotic sequence; Then the j-th dimension of the i-th golden eagle in the chaotic population S1 is Initialized to: Among them, ub j and lb j They represent the upper and lower boundary conditions of the j-th dimension of the decision variable respectively; Using the reverse learning strategy, we can obtain the jth dimension of the i-th golden eagle in the reverse population S2. The chaotic population and the reverse population S2 of S1 are merged to obtain the population S, and the fitness value of each golden eagle in the population S is calculated. The population S is graded in a gradually increasing manner using the fast non-dominated sorting method; the golden eagles with smaller non-dominated levels are preferentially selected to join the initial population. If the non-dominated levels of two golden eagles are equal, the golden eagle with a larger crowding distance is preferentially added to the initial population; until the initial population is filled; at the same time, the golden eagle with the smallest non-dominated level is added to the external archive; Step 3: Determine whether the current number of iterations t reaches the set maximum number of iterations: If yes, execute step 7; if not, execute step 4: Step 4: Update the position of each golden eagle and calculate the fitness value; Step 5: Use differential evolution and reverse learning strategies to cause two disturbances to the golden eagle population in turn, select and retain the better individuals, and calculate the fitness value of the individuals after mutation; Step 6: Update the external archive and execute step 3 repeatedly; Step 7: After the iteration, the position of each golden eagle in the external archive is mapped to the optimal solution of the objective function, and the Pareto frontier containing the optimal solution set of the objective function is obtained; Step 8: Use the fuzzy comprehensive evaluation method to evaluate a potential optimal solution from the Pareto frontier, and use the interpolation function corresponding to the potential optimal solution as the optimal interpolation function for industrial robot trajectory planning.

2. The coupling interpolation function selection method in industrial robot trajectory planning according to claim 1 is characterized in that: In step 1, the objective function is: <h2 style=";text-align:left;direction:ltr">f2 = min{J1,J2,...,J<h2 style=";text-align:left;direction:ltr"> P <h2 style=";text-align:left;direction:ltr">} The constraints are: Where f1 represents the execution time objective function; f2 represents the impact objective function; J k represents the square time integral of the trajectory constructed by the kth interpolation function; n represents the number of passing points; h i represents the time interval between the i-th passing point and the i+1-th passing point; F represents the number of robot joints; t f represents the number of robot joints; k represents the number of interpolation functions to be selected in the interpolation function set; Θ j represents the position of robot joint j; represents the velocity of robot joint j; represents the acceleration of robot joint j; represents the jerk of robot joint j; θ jmax represents the position constraint of robot joint j; V jmax represents the velocity constraint of robot joint j; A jmax represents the acceleration constraint of robot joint j; J jmax represents the jerk constraint of robot joint j.

3. The coupling interpolation function selection method in industrial robot trajectory planning according to claim 1 is characterized in that: The method of using the fast non-dominated sorting method to classify the population S is: 21) Obtain a Pareto frontier by judging the dominance and non-dominance relationship between the golden eagles, and assign a level d to all the golden eagles in the Pareto frontier; d = 1; 22) Determine whether all golden eagles in population S are assigned grades: if so, the population S is graded; if not, execute step 23); 23) Eliminate the golden eagles of level d from the population S, d = d + 1; 24) Obtain a Pareto frontier by judging the dominance and non-domination relationship among the remaining golden eagles, and assign grade d to all golden eagles in the Pareto frontier; 25) Loop through step 22).

4. The method for selecting a coupling interpolation function in trajectory planning of an industrial robot according to claim 1, characterized in that: The crowding distance of each golden eagle in each level of Pareto frontier for: Where R represents the non-dominated level; j = 1, 2, 3, ..., nobj represents the objective function number; i represents the i-th non-dominated solution in the Pareto frontier with level R, and Respectively represent the maximum and minimum values ​​of the j-th objective function; and They represent the j-th objective function values ​​of the i-th and i-1-th non-dominated solutions in the Pareto front of level R, respectively.

5. The method for selecting a coupling interpolation function in trajectory planning of an industrial robot according to claim 1, characterized in that: In step 4, the position of each golden eagle is updated: in, represents the position of the i-th golden eagle at the t+1 iteration step; represents the position of the i-th golden eagle at the t-th iteration step; r1 and r2 are both random numbers between [0,1], p a and p c Respectively represent the attack coefficient and patrol coefficient of the current iteration step; A i represents the attack vector of the i-th golden eagle; C i represents the cruise vector of the i-th golden eagle; Weight represents the weight; and: Among them, t is the current iteration number, and T is the maximum iteration number; in, represents the best prey searched by the f-th golden eagle so far, and also represents the prey selected by the i-th golden eagle; X i Indicates the current position of the Golden Eagle; C i =[row(0,1),row(0,1),...,C i,m ,...,row(0,1)] Among them, C i,m represents the m-th dimension component of the patrol vector of the i-th golden eagle; d represents the n-dimensional hyperplane equation, which can be expressed as: Among them A i =[a1,a2,...,a n ] represents the attack vector of the i-th golden eagle, represents the prey selected by the i-th golden eagle; a j , a m Represents the attack vector A i The j-th and m-th components of , where m represents a random integer between [1,n] and rand(0,1) represents a random number between [0,1]; in, and Represent the initial and final values ​​of the attack coefficient respectively; and Represent the initial value and final value of the patrol coefficient respectively.

6. The coupling interpolation function selection method in industrial robot trajectory planning according to claim 1, characterized in that: The method of using differential evolution to disturb the golden eagle population is: Use the DE / rand / 2 strategy to perform differential mutation on Golden Eagle: Where r1, r2, r3, r4, r5 represent random population indexes between [1, N], and r1≠r2≠r3≠r4≠r5≠i, N is the size of the golden eagle population; F is the scaling factor; V i is the mutation vector; The crossover operation is: Among them, rand i,j represents a random number between [0,1]; k is a random integer between [1,D]; CR∈[0,1] is the crossover rate; U i is the test vector, also known as offspring; X i is called the target vector, also known as the parent; and: Among them, CR0 is the crossover constant; t is the current iteration number, and T is the maximum iteration number; Select the action as: in, Represents the i-th golden eagle after mutation.

7. The method for selecting a coupling interpolation function in trajectory planning of an industrial robot according to claim 6, characterized in that: The principle of using the reverse learning strategy to disturb the golden eagle population is: Among them, Y i represents the reverse individual; ub and lb represent the upper and lower boundary conditions of the decision variable respectively; Similarly, select the operation: Among them, X i Represents the i-th golden eagle after mutation.

8. The method for selecting a coupling interpolation function in trajectory planning of an industrial robot according to claim 1, characterized in that: In step 6, the external archive update method is: If the updated new position of the Golden Eagle can dominate some members in the current external archive, these members will be deleted and the updated Golden Eagle will be added to the external archive; if the updated new position of the Golden Eagle is in a non-dominated relationship with all members in the current external archive, and it indicates that the number of members in the external archive is less than the preset external archive size M, the updated Golden Eagle will be added to the external archive; otherwise, the Golden Eagle with the smallest crowding distance in the external archive will be taken as the deletion object, and the updated Golden Eagle will replace the deleted members in the external archive.

Images