Excavating Trajectory Planning of a Mining Rope Shovel Based on Material Surface Perception

The mining rope shovel (MRS) is one of the core pieces of equipment for open-pit mining, and is currently moving towards intelligent and unmanned transformation, replacing traditional manual operations with intelligent mining. Aiming at the demand for online planning of an intelligent shovel excavation trajectory, an MRS excavating trajectory planning method based on material surface perception is proposed here. First, point cloud data of the material stacking surface are obtained through laser radar to perceive the excavation environment and these point cloud data are horizontally calibrated and filtered to reconstruct the surface morphology of the material surface to provide a material surface model for calculation of the mining volume in the subsequent trajectory planning. Second, kinematics and dynamics analysis of the MRS excavation device are carried out using the Product of Exponentials (PoE) and Lagrange equation, providing a theoretical basis for calculating the excavation energy consumption in trajectory planning. Then, the trajectory model of the bucket tooth tip is established by the method of sixth-order polynomial interpolation. The unit mass excavation energy consumption and unit mass excavation time are taken as the objective function, and the motor performance and the geometric size of the MRS are taken as constraints. The grey wolf optimizer is used for iterative optimization to realize efficient and energy-saving excavation trajectory planning of the MRS. Finally, trajectory planning is carried out for material surfaces with four different shapes (typical, convex, concave, and convex–concave). The results of experimental validation show that the actual hoist and crowd forces are essentially consistent with the planned hoist and crowd forces in terms of the peak value and trend variations, verifying the accuracy of the calculation model and confirming that the full bucket rate and various parameters meet the constraints. Therefore, the trajectory planning method based on material surface perception are feasible for application to different excavation conditions.


Introduction
A mining rope shovel (MRS) is a kind of large mechanical excavator which is widely used in open-cast mining operation; it is one of the core pieces of equipment used in opencast mining [1]. Currently, MRS excavation work is mainly completed through manual operation of four key actions, namely, excavating, rotation, unloading, and returning, which are performed in a loop. However, this method has low digging efficiency, high energy consumption, and a high failure rate, which seriously hinders the application and development of MRS [2,3]. Therefore, there is an urgent need to replace traditional manual operation with intelligent excavation [4] in order to ensure that the MRS can autonomously and efficiently achieve the expected goals in any environment, as well as to establish an intelligent operation mode to ensure efficient, safe, and energy-saving excavation.
The reasonable planning of excavation trajectory is the basis of intelligent MRS. Currently, domestic and foreign scholars have conducted a great deal of research on excavation trajectory planning for excavators. Based on the theory of soil damage, Awuah-Offei et al. [5] constructed a digging resistance model of an electric shovel, analyzed the kinematics and dynamics of the digging process, and optimized the crowd and hoist speed with the unit energy consumption of digging materials as the objective function. Wang X et al. [6,7] proposed a point-to-point (PTP) trajectory planning method with the goal of minimizing the energy consumption per unit volume of mining, and designed an intelligent mining system accordingly. Wei B et al. [8] proposed a new three-degrees-of-freedom parallel excavating mechanism with higher flexibility for the MRS, and optimized the size and excavation trajectory of the new mechanism with the unit excavation energy consumption as the objective function. Bi Q et al. [9] used a multi-target genetic algorithm to optimize the excavation trajectory in stages for electric shovel automation and verified the effectiveness of the optimization method through field tests. Meng Y et al. [10] established a force model of the bucket during the excavation process based on Coulomb theory, optimized the excavation trajectory using the optimal energy consumption, and verified it by EDEM simulation. Zhang T et al. [11] proposed a multi-target trajectory planning framework for an MRS based on a pseudo-spectral method, which was used for autonomous excavation of the MRS and verified by simulation and experiment. Fan J et al. [12,13] proposed a new electro-hydraulic hybrid MRS and optimized its hoist and crowd speeds with the goal of minimizing energy consumption per unit excavation mass. It can be seen from the above research that most domestic and foreign scholars have focused on theoretical analysis for trajectory planning of excavators, and have not combined environmental perception with trajectory optimization. Moreover, the existing research often takes unit excavation energy consumption as the objective function, meaning that the optimization goal is relatively singular.
To address the aforementioned issues, this article takes a 1/30 scale model of the Chinamade WK series MRS as its research object and proposes an MRS excavating trajectory planning method based on material surface perception. In this method, the point cloud data of the material stacking surface is obtained by lidar in order to perceive the mining environment. The energy consumption per unit mass and the mining time per unit mass are taken as the optimization objectives, and the full bucket rate, speed, and force are taken as constraints. The gray wolf algorithm is used to plan the excavating trajectory. The first part of this paper introduces the research background, the second part analyzes the kinematics and dynamics of the MRS working device, and the third part introduces the method of geometric modeling the material surface using 3D point cloud data. In the fourth part, the trajectory planning model is established and the trajectory planning of the material surface is carried out under various working conditions. In the fifth part, a scale model testbed of the MRS is established and the trajectory planning results are verified by experiments. Finally, the full text is summarized and the paper is concluded.

Construction of the MRS Scale Model
The study focused on a certain model of WK series MRS made in China. On the basis of the model and based on similarity theory [14,15], an MRS scale model was established at a ratio of 1/30 and an experimental platform for the MRS prototype model was built, as shown in Figure 1.
The structural parameters corresponding to the front-end working device of the MRS scale model are defined by the diagram shown in Figure 2 [16]. The specific numerical values of each parameter and the geometric size of the dipper can be found in Table 1 [16]. The structural parameters corresponding to the front-end working device of the MRS scale model are defined by the diagram shown in Figure 2 [16]. The specific numerical values of each parameter and the geometric size of the dipper can be found in Table 1 [16].

Kinematic Analysis of the MRS Excavating Device
The digging cycle of the MRS consists of four stages: digging, rotating, unloading, and returning. This article mainly carries out trajectory planning for its excavating process. During the excavation process, the upper carriage rotation and the crawler walking mechanism are not working, and the boom remains stationary. Therefore, the MRS can be simplified to a 1R-1P system, that is, the rotary movement of the bucket handle around the 1  direction and the translational movement along the  The structural parameters corresponding to the front-end working device of the MRS scale model are defined by the diagram shown in Figure 2 [16]. The specific numerical values of each parameter and the geometric size of the dipper can be found in Table 1 [16].

Kinematic Analysis of the MRS Excavating Device
The digging cycle of the MRS consists of four stages: digging, rotating, unloading, and returning. This article mainly carries out trajectory planning for its excavating process. During the excavation process, the upper carriage rotation and the crawler walking mechanism are not working, and the boom remains stationary. Therefore, the MRS can be simplified to a 1R-1P system, that is, the rotary movement of the bucket handle around the 1  direction and the translational movement along the

Kinematic Analysis of the MRS Excavating Device
The digging cycle of the MRS consists of four stages: digging, rotating, unloading, and returning. This article mainly carries out trajectory planning for its excavating process. During the excavation process, the upper carriage rotation and the crawler walking mechanism are not working, and the boom remains stationary. Therefore, the MRS can be simplified to a 1R-1P system, that is, the rotary movement of the bucket handle around the ξ 1 direction and the translational movement along the ξ 2 direction, as shown in Figure 3a. Figure 3a,b shows the general pose and initial pose in the process of the MRS excavation, where l 1 and l 2 are the horizontal distance and perpendicular distance from the rotary center of saddle block to the coordinate origin, respectively; l 3 is the distance from the saddle block rotary center to the bucket tooth tip in the perpendicular direction of the bucket handle; θ 1 is the rotary angle of the bucket handle; d 2 is the elongation of the bucket handle; and d 20 is the initial value of the bucket handle extension length. In the article, the kinematics modeling and analysis of the working device of the front end of the shovel is carried out using the product of exponentials (PoE) [17][18][19]. The construction of the working device is shown in the above figure. First, the fixed base coordi In the above equation, the expression for the rotational axis spinor is Therefore, the specific expression of () BW Tq can be obtained by calculation as follows [16]: In the article, the kinematics modeling and analysis of the working device of the front end of the shovel is carried out using the product of exponentials (PoE) [17][18][19]. The construction of the working device is shown in the above figure. First, the fixed base coordinate system {B} and the end bucket tooth tip coordinate system {W} are established as shown. Second, the following process variables need to be introduced: T BW (q) denotes the pose transformation matrix of {W} relative to {B}; T BW (0) denotes the initial pose transformation matrix of {W} relative to {B}; and ξ 1 and ξ 2 represent the joint spinors of the rotating joint and prismatic joint in {B}, respectively. Therefore, the kinematic model from the base coordinate system {B} of the working device to the end bucket tooth tip coordinate system {W} can be expressed as [16].
In the above equation, the expression for the rotational axis spinor is ξ = (ω, ν) T , where ν = r × ω, r denotes an arbitrary point on the rotary axis, and ω is the unit vector representing the direction of the rotary axis. The spinor expression for the translational axis is ξ = (0, ν) T , while the kinematic matrix of the rotary axis is [16].
Therefore, the specific expression of T BW (q) can be obtained by calculation as follows [16]: where l 1 , l 2 , and l 3 are fixed structural parameters. Therefore, by measuring the values of θ 1 and d 2 we can obtain the pose of the bucket tooth tip relative to the body, thereby completing the solution of the forward kinematics problem. The solution to the inverse kinematics problem for the MRS scaled model is to calculate the joint variables θ 1 and d 2 based on the pose of the bucket tooth tip. Thus, assuming that the bucket pose transform matrix T is known [16], In the above formula, a is the rotary transformation matrix and b is the translational transformation matrix. By comparing (3) and (4), we obtained [16] By solving Equation (6), we can obtain the expressions for θ 1 and d 2 [16]: Therefore, when the pose transformation matrix from the machine body to bucket tooth tip is known, it is possible to calculate the values of θ 1 and d 2 from Formulas (7) and (8) to solve the inverse kinematics problem.

Dynamic Analysis
As can be seen from the previous kinematics analysis, the working mechanism of the MRS can be simplified to a 1R-1P system, which includes the rotary movement of the bucket handle and bucket around the saddle block rotary center as well as the translational motion along the direction of the bucket handle. Therefore, this article defines the angle θ 1 between the bucket handle and the vertical direction and the elongation d of the bucket handle as generalized coordinates, as shown in Figure 4. In considering the structural characteristics of the working device of the MRS and providing a better calculation model for the subsequent trajectory optimization, this paper adopts the Lagrangian equation to dynamically model the working device of the MRS. The Lagrangian equation of the working device at the front end of the MRS is obtained as follows:  The Lagrangian equation of the working device at the front end of the MRS is obtained as follows: In the above equation, L represents the Lagrangian function, which is defined as the difference between the total kinetic energy K and total potential energy P in the relative inertial system, that is, L = K − P, while F d is the resulting force acting on the bucket along the bucket handle direction and F θ is the resulting force along the tangential direction of the bucket tooth tip. Based on the two generalized coordinates, the dynamic model of the digging device at the front end of the MRS can be obtained as follows: where m d is the mass of the bucket handle; L d is the length of the bucket handle; m c is the mass of the bucket, including the mass of the excavated materials and the weight of the bucket itself; L c is the length of the bucket; F ti is the hoist force of the wire rope; θ 2 is the angle between the wire rope and the bucket handle direction; F τ is the tangential digging resistance; F tui represents the crowd force of the bucket handle; and F n represents the normal digging resistance.

Material Surface Scanning Based on Laser Radar
To achieve rapid modeling of the material stack surfaces, in this paper we propose a method for using three-dimensional point clouds to rapidly construct a geometric model of the material surfaces. The 3D laser radar emits laser beams towards the material stack and obtains a large amount of point cloud data by calculating the reflection time and propagation speed of the laser. Then, the ordered point cloud is obtained by a combined filtering process to complete the modeling of the material surface.
The original point cloud data of the material stack, as shown in Figure 5, indicate that there is a large amount of noise in the point cloud due to the influence of scanning speed and the surrounding environment, making it difficult to generate an accurate threedimensional surface model. Therefore, the combined filtering method shown in Figure 6, which was jointly developed based on C++ and the PCL point cloud library, was used to reconstruct the surface. In considering the structural characteristics of the working device of the MRS and providing a better calculation model for the subsequent trajectory optimization, this paper adopts the Lagrangian equation to dynamically model the working device of the MRS. The Lagrangian equation of the working device at the front end of the MRS is obtained as follows: (9) In the above equation, L represents the Lagrangian function, which is defined as the difference between the total kinetic energy K and total potential energy P in the relative inertial system, that is, L K P = − , while d F is the resulting force acting on the bucket along the bucket handle direction and F θ is the resulting force along the tangential direction of the bucket tooth tip. Based on the two generalized coordinates, the dynamic model of the digging device at the front end of the MRS can be obtained as follows:

Horizontal Calibration of Point Cloud Data
Due to the uneven ground of the mine, the point cloud data have an error in the roll angle  and pitch angle  relative to the horizontal plane. Therefore, horizontal calibration was performed before removing the ground points. In the horizontal plane, the Cartesian coordinate system was established perpendicular to the direction of the MRS crawler as the x-axis, parallel to the direction of the crawler as the y-axis, and vertical to the direction of the horizontal plane as the z-axis.
The coordinate transformation matrix is shown below: cos 0 sin sin sin cos sin cos cos sin sin cos cos The formula for correcting the original point cloud data is as follows:

Horizontal Calibration of Point Cloud Data
Due to the uneven ground of the mine, the point cloud data have an error in the roll angle α and pitch angle β relative to the horizontal plane. Therefore, horizontal calibration was performed before removing the ground points. In the horizontal plane, the Cartesian coordinate system was established perpendicular to the direction of the MRS crawler as the x-axis, parallel to the direction of the crawler as the y-axis, and vertical to the direction of the horizontal plane as the z-axis.
The coordinate transformation matrix is shown below: The formula for correcting the original point cloud data is as follows: In the formula, x s , y s , and z s are the original point cloud data.

Point Cloud Data Combination Filtering Process
In this paper, the ground noise was filtered by straight-through filtering. Ground points were removed by setting a height threshold. If the z value of a point in the point cloud data was less than the maximum height, it was considered a ground point and discarded. If it was greater than or equal to the maximum height, it was considered a surface point of the material and was retained. This allows for fast removal of outliers, achieving the goal of rough processing in the first step. Additionally, a filter method combining radius and statistical filtering was used to remove surface and outlier noise. The formula used for calculating the Euclidean distance between any two points is as follows: In the equation, x, y, and z represent the three-dimensional coordinates of the point clouds. Radius filtering is used to calculate the Euclidean distance between a point cloud and other point clouds in its spatial neighborhood, and the point cloud is removed or retained by comparing the distance between the two along with the set threshold distance K s .
Based on this, the mean Euclidean distance d i of all point cloud neighborhoods can be calculated according to the number of point clouds in the neighborhood and the global neighborhood mean d can be calculated based on the total number of point clouds. Then, the standard deviation of the point cloud neighborhood is calculated as follows: where V represents the variance of the neighborhood and n represents the total number of point clouds. Assuming that the mean Euclidean distance d i of the point cloud neighborhood is normally distributed, we have By setting a filtering threshold factor α, the statistical filtering algorithm removes the point clouds as noise points when the following formula holds: After filtering, the number of point cloud data items is on the order of millions, making it difficult to store and process the data. Therefore, the uniform voxel algorithm was used for "down-sampling" of the point cloud data for compression using the following formula to calculate the minimum length of the voxel grid: (17) where x max , y max , and z max represent the maximum values of the three-dimensional coordinates of the point cloud and x min , y min , and z min represent the minimum values. Then, we can use the following formula to calculate the voxel grid size: where r is the length of the set voxel small grid. The centroid point within the voxel grid is used to replace all points within the grid, achieving the purpose of down-sampling.
After point cloud sampling, the improved bilateral filtering based on normality is used to smooth the surface of the point cloud to eliminate sharp noise caused by individual point cloud fluctuations, improve the quality of the point cloud, and complete surface reconstruction. This algorithm considers the distance between points and neighboring points and uses the distance along the normal direction as a judgment basis. For a certain point P in the point cloud, the unit normal vector of point P is first calculated using the points within the neighborhood range of point P, then the position of P is updated using the following formula: where P is the updated position of the point, n P represents the normal vector of point P, and δ P represents the bilateral filtering factor, as shown in the following equation: In the above formula, N r (P) represents the neighborhood of a point P, w d and w n are given Gaussian weights, and q represents the points within the neighborhood P.
Due to the fact that the shape of the material surface is often complex and variable, and may not necessarily be an ideal 40-degree slope [20,21], this paper piled four different shapes of material surfaces and scanned and filtered the material stacking surfaces under different working conditions to obtain the material stacking surfaces model shown in Figure 7. The color bars in the figure represent the vertical height of the material stacking surface from 0-400 mm. The material properties of the material stacking surface are shown in Table 2 below.
In the above formula, () rP N represents the neighborhood of a point P , d w and n w are given Gaussian weights, and q represents the points within the neighborhood P .
Due to the fact that the shape of the material surface is often complex and variable, and may not necessarily be an ideal 40-degree slope [20,21], this paper piled four different shapes of material surfaces and scanned and filtered the material stacking surfaces under different working conditions to obtain the material stacking surfaces model shown in Figure 7. The color bars in the figure represent the vertical height of the material stacking surface from 0-400 mm. The material properties of the material stacking surface are shown in Table 2 below.

Trajectory Planning Model Based on Grey Wolf Algorithm
Before establishing an optimization model, it is often necessary to determine the three elements of optimization: optimization variables, objective functions, and constraints. In this paper, information on the material stacking surface was obtained through laser scanning and the material surface data were imported into the trajectory planning model. The optimal excavation trajectory under the current operating conditions was obtained by using the grey wolf algorithm to perform optimization calculations.

Determination of Optimization Variables
In this paper, following the comparative experiments performed by Wang et al. [4], trajectory planning was carried out by using the sixth-order polynomial interpolation method and the excavation process of the MRS was divided into x direction and y direction motion, as shown in Figure 8. The design variables were determined as x = [  To improve the efficiency of surface fitting [22,23], eight points with equal intervals were selected in the X axial direction of the point cloud and another eight points with equal intervals were selected in the Y axial direction. Surface fitting of the material stacking surface was performed using 64 points, with the goodness-of-fit parameter R 2 used to constrain the fitting effect (R 2 > 0.95). Finally, the mathematical expression of the material surface in the {Oo} coordinate system was obtained:

Trajectory Planning Model Based on Grey Wolf Algorithm
Before establishing an optimization model, it is often necessary to determine the three elements of optimization: optimization variables, objective functions, and constraints. In this paper, information on the material stacking surface was obtained through laser scanning and the material surface data were imported into the trajectory planning model. The optimal excavation trajectory under the current operating conditions was obtained by using the grey wolf algorithm to perform optimization calculations.

Determination of Optimization Variables
In this paper, following the comparative experiments performed by Wang et al. [4], trajectory planning was carried out by using the sixth-order polynomial interpolation method and the excavation process of the MRS was divided into x direction and y direction motion, as shown in Figure 8. The design variables were determined as x =  [a x0 , a x1 , a x2 , a x3 , a x4 , a x5 , a x6 , a y0 , a y1 , a y2 , a y3 ,a y4 , a y5 , a y6 ] and the excavation trajectory of the MRS can be represented as g x (t) = a x6 t 6 + a x5 t 5 + a x4 t 4 + a x3 t 3 + a x2 t 2 + a x1 t + a x0 g y (t) = a y6 t 6 + a y5 t 5 + a y4 t 4 + a y3 t 3 + a y2 t 2 + a y1 t + a y0 (22) Sensors 2023, 23, x FOR PEER REVIEW 11 of 23 Therefore, based on the state parameters of the initial and end positions of the bucket tooth tip, the coefficients of the sixth-order polynomial can be processed:  The intersection point of the material surface and the horizontal plane, which is the initial position of the bucket tooth tip, was set as the excavation starting point (0,0). Here, the termination position of the tooth tip at the end of the excavation is set to g x , g y and the excavation termination time is t d . To guarantee the smoothness of the MRS excavating process, the velocities and accelerations of the starting and ending points were both set to 0, that is: Therefore, based on the state parameters of the initial and end positions of the bucket tooth tip, the coefficients of the sixth-order polynomial can be processed: Therefore, the sixth-order polynomial only needs to determine the size of a x6 , a y6 , t d , g x , g y to obtain the excavation trajectory of the bucket tooth tip at any time, meaning that the variables to be optimized are x = [a x6 , a y6 , t d , g x , g y ].

Objective Function Determination
In this paper, the excavation energy consumption and excavation time are comprehensively considered in the excavation trajectory planning and the objective function of trajectory planning is constructed by the unit mass excavation energy consumption and unit mass excavation time of the MRS.
(1) Minimum unit mass excavation energy consumption target Reducing the excavating energy consumption while meeting the excavating requirements is the main problem of MRS trajectory planning. The energy consumption of the MRS excavating process is mostly the energy consumed by the crowd and hoist motors to overcome the excavation resistance and the gravity of the bucket handle, bucket, and material. Therefore, this paper takes the minimum unit mass excavating energy consumption as the primary objective function of the excavation process; its mathematical expression is Here, ∧ E represents the total energy consumption during excavation, v h represents the hoist speed, F h represents the hoist force, v c represents the crowd speed of the bucket handle, F c represents the crowd force, and m dig represents the quality of the excavated material.
(2) Maximum excavating efficiency target To improve the efficiency of MRS excavation, this article takes the time required to excavate the quality of material per unit as another optimization objective. Obviously, the less time it takes to excavate quality material per unit, the higher the excavation efficiency of the MRS. The mathematical expression for this objective is where t dig is the total excavation time. (

3) Total objective function
To simplify the calculation model used for trajectory planning, a weighted method is adopted to combine the above two objective functions into one. In this way, multi-target optimization is simplified into single-target optimization. According to the role of each objective function in the actual excavating process, the index tolerance method was used to confirm the weight coefficient from the variation range of each objective function, as shown in Table 3. Finally, the weight coefficients of each objective function were set to k 1 = 0.7 and k 2 = 0.3. The overall objective function is as follows:

Calculation of Excavation Volume
From the previous processing of point cloud data, we can obtain the expression of the material surface in the {Oo} coordinate system as z = f (x, y). Assuming the optimized excavation surface of the bucket is g(x, y), the excavation volume can be estimated through double integration. The calculation method is shown in Figure 9. The integration domain S is divided into k closed spaces ∆s i ; when ∆s i is very small, a point (x i , y i ) can be taken randomly in ∆s i , then the product of [ f (x i , y i ) − g(x i , y i )] and ∆s i can be used to calculate the volume of each part. Finally, the volumes of all parts are summed to obtain the excavation volume. The calculation formula is as follows:

Determination of Constraint Conditions
In order to ensure the feasibility of the excavation trajectory, constraint conditions need to be introduced into the MRS trajectory planning model, which can be mainly categorized into excavation performance constraints, motor performance constraints, and geometric dimension constraints, as follows.
(1) Constraint on bucket filling rate To guarantee the efficiency of the MRS excavating process and to avoid overloading or underloading during excavation, the variation range of the bucket filling rate is limited to 90% to 110%, that is:

Determination of Constraint Conditions
In order to ensure the feasibility of the excavation trajectory, constraint conditions need to be introduced into the MRS trajectory planning model, which can be mainly categorized into excavation performance constraints, motor performance constraints, and geometric dimension constraints, as follows.
(1) Constraint on bucket filling rate To guarantee the efficiency of the MRS excavating process and to avoid overloading or underloading during excavation, the variation range of the bucket filling rate is limited to 90% to 110%, that is: In the above equation, V represents the volume of the materials excavated by the bucket and V cd is the capacity of the bucket itself.
(2) Constraint on excavation time To improve the excavating efficiency of the MRS, and based on the relevant literature and actual working conditions, the time range for a single excavation operation of the MRS is set to 9 s to 15 s [11,12], that is: (

3) Digging back angle constraint
To guarantee the stability of the excavation process, the variation range of the digging back angle is limited as follows: where β max and β min represent the upper and lower limits of the digging back angle, respectively.
(4) Velocity constraint The digging process of the MRS is mainly driven by the crowd and hoist motors. Limited by the performance of the motors, the speed at any time must be less than the maximum output speed of the motor. In the meantime, to ensure the stability of the excavation, the crowd speed and hoist speed at any time should be greater than zero, that is: In the above equation, v hmax and v cmax represent the maximum rotational speeds of the hoist motor and crowd motor.

(5) Driving force and power constraints
To ensure that at any moment during the digging process, neither the crowd and hoist forces nor the crowd and hoist power should exceed the rated value of the motors, that is: In the above equation, F hmax and F cmax represent the maximum driving forces of the hoist motor and crowd motor, respectively, while P hmax and P cmax are the respective rated powers of the hoist and crowd motors.

(6) Geometric dimension constraints
The digging cycle of the MRS consists of four stages: digging, rotating, unloading and returning. To guarantee the smooth progress of the subsequent rotating and unloading processes, it is necessary to make sure that the bucket leaves the material surface when the digging is finished and that the bottom of the bucket is higher than the height of the shovel loader hopper, that is, where h is the height of the bucket tooth tip, H dm is the height when the bucket tooth tip leaves the material surface, H kk is the height of the shovel loader hopper, and H cd is the height from the bucket tooth tip to the bottom of the bucket. The maximum stroke of the MRS crowd motion is limited by the size of the bucket handle, while the hoist wire rope is limited by the top sheave of the boom, that is, where d is the elongation of the bucket handle, L dg is the length dimension of the bucket handle, and H tl is the height of the top sheave of the boom.
Therefore, the optimized design model is as follows:

Trajectory Planning Method Based on Grey Wolf Optimizer
The Grey Wolf Optimizer (GWO) is a swarm intelligence optimization algorithm inspired by the social hierarchy and hunting behavior of wolves. It was proposed by Australian scholar Mirjalili et al. [24] in 2014. The algorithm has the advantages of simple structure, few parameters to be adjusted, easy implementation, good solution accuracy and fast convergence speed.
The grey wolf algorithm first constructs the social hierarchy of a wolf pack. The fitness of each individual in the population is calculated and the three wolves with the highest fitness are designated as the α wolf, β wolf, and δ wolf, while the rest of the individuals are labeled as ω wolves. The process of searching for prey (seeking the optimal solution) is mainly achieved through the guidance of the α, β, and δ wolves. During the iterative calculation process, the three wolves α β δ with the highest fitness in the current population are preserved and the positions of other candidate wolves ω are continuously adjusted based on the coordinates of these three wolves; hunting (optimization) is completed by searching for prey, encircling prey, and attacking prey, and finally a set of optimal solutions is obtained. The specific mathematical model is shown below: Because the trajectory planning is for the intelligent operation mode service, it needs to meet real-time requirements; thus, the initialization of the population is particularly important. This article adopts the good point set method to obtain uniformly distributed initial values and then calculates iteratively using the GWO, which improves the convergence speed of the algorithm and can avoid its falling into local optima. Figure 11 shows Because the trajectory planning is for the intelligent operation mode service, it needs to meet real-time requirements; thus, the initialization of the population is particularly important. This article adopts the good point set method to obtain uniformly distributed initial values and then calculates iteratively using the GWO, which improves the convergence speed of the algorithm and can avoid its falling into local optima. Figure 11 shows the grey wolf population generated by the good point set method and by the random method in a three-dimensional search space. It can be seen from the figure that the initial population generated by the random method cannot be uniformly distributed in the entire search space and has strong randomness, making it easy for the algorithm to fall into local optima. In contrast, the initial population generated by the good point set method can be uniformly distributed in the entire search space, enriching the population diversity and enhancing the global search capability of the algorithm. The improved GWO process is shown in Figure 12.

Results of Trajectory Planning
Combining the improved grey wolf algorithm process, the basic parameter as follows: initial wolf population size pop = 50, search coefficient c = 2, and m iterations tmax = 100. After trajectory planning, the trajectory planning results for d excavation conditions are shown in Figures 13-16. The planned results indicate

Results of Trajectory Planning
Combining the improved grey wolf algorithm process, the basic parameters are set as follows: initial wolf population size pop = 50, search coefficient c = 2, and maximum iterations t max = 100. After trajectory planning, the trajectory planning results for different excavation conditions are shown in Figures 13-16. The planned results indicate that the trajectory planning method based on the grey wolf algorithm can be applied to different excavation conditions. hoist speed for various material surface. Figure 13c shows the crowd force and hoist force for various material surface. Figure 13 displays the planning results for a typical material surface, while Figure 14 shows the planning result for a concave material surface. Compared with the typical material surface, with the concave surface the hoist speed is significantly reduced and the excavation trajectory is deeper. Figure 15 shows the planning results for a convex material surface. Compared with the typical material surface, the crowd speed decreases, the hoist speed increases, and the excavation trajectory is shallower. Finally, Figure 16 shows the planning results for a convex-concave material surface. Compared with the typical material surface, the crowd and hoist speed are reduced and the excavation trajectory is between that of the typical material surface and the convex material surface.         Tables 4 and 5 show the optimization results obtained from trajectory planning under different excavating conditions. It can be found that the shallower the excavation trajectory, the less excavation time required; in addition, the convex-concave material surface requires the least excavating time and energy consumption. All optimization variables satisfy the constraint conditions, and the bucket fill factor meets the requirements; therefore, the trajectory planning method based on the grey wolf algorithm is applicable to different excavating conditions and has certain feasibility and reliability.  Figure 13a shows the optimal excavating trajectory under different excavation conditions; the abscissa represents the displacement in the horizontal direction, the ordinate represents the displacement in the vertical direction, the blue line represents the material surface fitting curve, the red line represents the planned excavation trajectory, and the black line represents the ideal 40 • material surface. Figure 13b shows the crowd speed and hoist speed for various material surface. Figure 13c shows the crowd force and hoist force for various material surface. Figure 13 displays the planning results for a typical material surface, while Figure 14 shows the planning result for a concave material surface. Compared with the typical material surface, with the concave surface the hoist speed is significantly reduced and the excavation trajectory is deeper. Figure 15 shows the planning results for a convex material surface. Compared with the typical material surface, the crowd speed decreases, the hoist speed increases, and the excavation trajectory is shallower. Finally, Figure 16 shows the planning results for a convex-concave material surface. Compared with the typical material surface, the crowd and hoist speed are reduced and the excavation trajectory is between that of the typical material surface and the convex material surface. Tables 4 and 5 show the optimization results obtained from trajectory planning under different excavating conditions. It can be found that the shallower the excavation trajectory, the less excavation time required; in addition, the convex-concave material surface requires the least excavating time and energy consumption. All optimization variables satisfy the constraint conditions, and the bucket fill factor meets the requirements; therefore, the trajectory planning method based on the grey wolf algorithm is applicable to different excavating conditions and has certain feasibility and reliability.

Experimental Verification
In this study, a scale model testbed of the MRS (as shown in Figure 17) was constructed for excavation testing. The feasibility and credibility of the trajectory planning method were then validated by comparing the planned and tested results.

Experimental Verification
In this study, a scale model testbed of the MRS (as shown in Figure 17) was constructed for excavation testing. The feasibility and credibility of the trajectory planning method were then validated by comparing the planned and tested results. The hoist and crowd movements of the digging device in the prototype were achieved by individual motor driven systems. The motor control was implemented through an upper computer, LabVIEW 2020 software, and a microcontroller. Dynamic torque sensors were employed to measure the real-time crowd force of the bucket handle, tension sensors to measure the real-time hoist force of the wire rope, and rotation encoders to measure the extension of the bucket handle and the lifting distance of the wire rope. LabVIEW software in cooperates with a data collector was used to realize the acquisition, storage, and display of sensor measurement data. The installation position of the motors and sensors is shown in Figure 18. The structure of the testbed system is shown in Figure  19.
To guarantee the credibility and stability of the test results, it is necessary to repeatedly take the average of trajectory planning tests under different excavating conditions. The data obtained from the torque and tensile sensors were processed using Matlab R2018a, resulting in the curves of crowd force and hoist force over time. The hoist and crowd movements of the digging device in the prototype were achieved by individual motor driven systems. The motor control was implemented through an upper computer, LabVIEW 2020 software, and a microcontroller. Dynamic torque sensors were employed to measure the real-time crowd force of the bucket handle, tension sensors to measure the real-time hoist force of the wire rope, and rotation encoders to measure the extension of the bucket handle and the lifting distance of the wire rope. LabVIEW software in cooperates with a data collector was used to realize the acquisition, storage, and display of sensor measurement data. The installation position of the motors and sensors is shown in Figure 18. The structure of the testbed system is shown in Figure 19.
To guarantee the credibility and stability of the test results, it is necessary to repeatedly take the average of trajectory planning tests under different excavating conditions. The data obtained from the torque and tensile sensors were processed using Matlab R2018a, resulting in the curves of crowd force and hoist force over time. Figure 20 compares the planned results and test results of the crowd and hoist forces under different excavation conditions. Due to the flexibility of the hoist wire rope and the resistance fluctuations during the excavation process, the test results are presented as an obviously fluctuating curve. Compared with the planned results, it can be seen that the values of the two are not much different and that the overall change trend is essentially the same. At the same time, the excavating quality and excavating time of the planned and test were compared to validate the credibility and feasibility of the trajectory planning method. Relevant coefficients (R 2 ) were introduced to describe the degree of match between the planned and tested results, as shown in Table 6.  Figure 20 compares the planned results and test results of the crowd and hoist forces under different excavation conditions. Due to the flexibility of the hoist wire rope and the resistance fluctuations during the excavation process, the test results are presented as an obviously fluctuating curve. Compared with the planned results, it can be seen that the values of the two are not much different and that the overall change trend is essentially the same. At the same time, the excavating quality and excavating time of the planned and test were compared to validate the credibility and feasibility of the trajectory planning method. Relevant coefficients ( 2 R ) were introduced to describe the degree of match between the planned and tested results, as shown in Table 6.    As shown in the table, the relative deviations between the planned and tested excavation quality under different excavation conditions are about 5%. The relative deviations of the digging times are below 1%. Additionally, the R 2 values for the planned and tested results of the crowd force and the hoist force are greater than 0.85, indicating that the two have a high degree of agreement. The test results show that the trajectory planning method based on material stack surface perception has certain feasibility and reliability, that the established dynamic model can more accurately predict the crowd and hoist force, and that the planned results of various performance indicators are reliable.

Conclusions
(1) A laser radar was used to obtain the point cloud data of the material stack surface in order to perceive the excavating environment, and the point cloud data were horizontally calibrated and filtered to establish a prediction model of the material stack surface. Furthermore, kinematic and dynamic analyses of the MRS excavation device were conducted using the Product of Exponentials and Lagrange equation. (2) A trajectory planning method for the MRS excavation based on material surface perception and the Grey Wolf Algorithm was proposed, with the unit mass excavation energy consumption and unit mass excavation time as the target functions and the electric motor performance and MRS geometry size as constraints. Trajectory planning was conducted on four different shapes (typical, concave, convex, and convex-concave) of material stack surfaces. (3) An MRS scale model testbed was constructed and used for experimental verification.
The test results show that the planned results for hoist force and crowd force were generally consistent with the test results in terms of the values and change trend under different excavation conditions and had R 2 values greater than 0.85, validating the feasibility and reliability of the proposed trajectory planning method.
This method can provide a theoretical basis for the intelligent and unmanned development of excavating machinery such as electric shovels and excavators, and has certain feasibility and applicability. In follow-up research, the initial pose of the MRS could be incorporated into the trajectory planning scope. Meanwhile, the point cloud processing and trajectory planning algorithm should be further improved and optimized in order to reduce the optimization time and improve computational efficiency to better meet practical needs.