Mine roadway topology-oriented monorail global energy consumption optimal path planning method

Abstract

The present application relates to the technical field of mine auxiliary transportation automation control, and specifically discloses a monorail global energy consumption optimal path planning method for mine roadway topology, collects roadway measurement data to construct a coarse-grained topological map composed of nodes and arc segments, encapsulates high-density elevation point cloud subsets within each arc segment; takes the map as a search space and elevation difference as a basic cost to search for a preliminary path sequence; retrieves the elevation point cloud subsets of each arc segment along the preliminary path, reconstructs the elevation continuous change curve and calculates the slope change rate; differentiates the arc segment intervals according to the slope change rate to form a non-uniformly distributed grid along the whole path; takes each distribution point on the grid as a discretization reference, adds the square of the second derivative of the speed to the penalty term to the traction net energy consumption optimization target, solves the speed value of each distribution point under the constraints of traction force, braking deceleration and centrifugal force of the curve, and outputs the speed instruction sequence.

Description

Technical Field

[0001] This invention relates to the field of automated control technology for mine auxiliary transportation, specifically to a method for global energy consumption optimization path planning for monorail cranes based on mine roadway topology. Background Technology

[0002] Monorail transport systems are a crucial component of mine auxiliary transportation, responsible for the underground transport of materials, equipment, and personnel. Due to the complex network of mine roadways, monorail locomotives frequently encounter uphill and downhill sections, curves, and switchovers during operation. Their traction energy consumption is closely related to the selection of the travel path, the control of the operating speed, and the topographical features of the roadways. Existing monorail path planning methods in mines mostly focus on the path search algorithm itself, typically using roadway length, travel time, or simple geometric distance as optimization objectives, and employing graph search algorithms to find a feasible path from the starting point to the destination within the roadway network.

[0003] The existing technology has the following shortcomings:

[0004] Existing mine monorail path planning methods simplify roadways into uniformly sloped straight sections discretized at fixed intervals, completely filtering out micro-terrain undulations with wavelengths of only tens of meters and equivalent slope variations of 2 to 3 degrees. This results in the inability to perceive and express the frequent switching between traction and regenerative braking caused by such high-frequency undulations. Consequently, it is impossible to avoid high-energy-consuming terrain in advance during the path search stage, nor can it suppress it through reasonable inertial ramp-up strategies in the subsequent speed optimization. Ultimately, this leads to an uncontrollable systematic deviation between planned energy consumption and actual operating energy consumption. Summary of the Invention

[0005] The purpose of this invention is to provide a global energy-optimal path planning method for monorail cranes based on mine roadway topology, in order to solve the problems mentioned above.

[0006] The objective of this invention can be achieved through the following technical solutions:

[0007] A method for global energy-optimal path planning for monorail cranes based on mine roadway topology includes the following steps:

[0008] S1 collects the three-dimensional coordinates of the centerline of the mine roadway, the elevation of the slope change point, the radius of curvature of the curve, the position of the turnout and the division of the power supply section, obtains the elevation point cloud with an interval of no more than five meters between adjacent measuring points, constructs a coarse-grained topology map composed of nodes and arc segments, and encapsulates the corresponding elevation point cloud subset inside each arc segment;

[0009] S2 uses a coarse-grained topological map as the search space and the elevation difference between the start and end points of each arc segment as the basic cost to search for a preliminary path sequence from the starting node to the target node.

[0010] S3. Along the initial path sequence, retrieve the subset of elevation point clouds encapsulated within each arc segment, reconstruct the continuous elevation change curve of each arc segment along the mileage direction, and calculate the derivative of the slope angle relative to the mileage at each location as the slope change rate.

[0011] S4. Based on the slope change rate values ​​at each location, a threshold is set, and the arc segments with slope change rates exceeding the threshold are subdivided along the mileage direction. The segments with slope change rates not exceeding the threshold are kept at their original span, forming a non-uniform point grid along the entire route.

[0012] S5 uses each point on a non-uniform grid as the discretization benchmark and the elevation change curve as the terrain input. The sum of the squares of the second derivatives of the running speed is superimposed as a penalty term on the traction net energy consumption optimization objective. Under the constraints of traction force, braking deceleration and curve centrifugal force, the speed values ​​of each point are solved to minimize the optimization objective, and the speed command sequence arranged along the mileage is output.

[0013] As a further aspect of the present invention: the construction of a coarse-grained topological map composed of nodes and arc segments, and the encapsulation of corresponding elevation point cloud subsets within each arc segment, specifically includes:

[0014] Using the centerline of the tunnel as the baseline, sampling points are picked at equal intervals along the direction of the baseline to obtain the three-dimensional coordinates and corresponding elevation values ​​of each sampling point, forming an elevation sequence along the tunnel direction.

[0015] Perform a first-order difference operation on the high-order sequence to obtain the slope value sequence at each sampling point, and perform a second-order difference operation on the slope value sequence to obtain the curvature value at each sampling point;

[0016] Sampling points with curvature values ​​exceeding a preset curvature threshold are selected as terrain change points. The continuous path between adjacent terrain change points in each tunnel section is defined as an arc segment, and terrain change points and tunnel intersections are defined as nodes.

[0017] The elevation point cloud data of all sampling points within each arc segment is compressed and encapsulated into an elevation point cloud subset of the arc segment. The network structure consisting of nodes and arc segments is stored as a coarse-grained topology map, in which each arc segment is associated with the node number of its starting and ending points, the arc segment length, and the corresponding elevation point cloud subset.

[0018] As a further aspect of the present invention: the formation process of the preliminary path sequence is as follows:

[0019] Starting from the current node, traverse all arcs adjacent to the current node, retrieve the elevation point cloud subset corresponding to each arc, and calculate the fluctuation amplitude of the elevation within each arc.

[0020] The weighted sum of the elevation difference between the start and end points of each arc segment and the fluctuation amplitude is used as the passage cost of the arc segment, and the adjacent arc segment with the minimum passage cost is selected as the preferred arc segment.

[0021] Update the end node of the preferred arc segment to the current node, and record the preferred arc segment in the path sequence. Repeat the traversal and selection process until the current node is the target node.

[0022] Delete arcs in the path sequence whose travel cost exceeds a preset limit, replace them with the arcs with the lowest cost among the remaining arcs, and output the corrected path sequence.

[0023] As a further aspect of the present invention: the calculation process for the fluctuation amplitude is as follows:

[0024] The elevation values ​​of all sampling points in the current arc segment are read from the elevation point cloud sub-set. After removing outliers whose elevation difference with adjacent points exceeds a predetermined multiple, the effective elevation sequence of the arc segment is generated.

[0025] Perform differential operations on the elevation values ​​of each adjacent point in the effective elevation sequence to obtain multiple differential values ​​along the arc direction, and calculate the sum of the absolute values ​​of the differential values;

[0026] Calculate the ratio of the sum of the absolute values ​​of the differences to the total number of sampling points within the arc segment minus one, and use this ratio as the average absolute difference fluctuation amplitude of the arc segment;

[0027] Multiplying the average absolute difference fluctuation amplitude by the normalization coefficient of the arc length yields the final fluctuation amplitude of the arc. The final fluctuation amplitude increases monotonically with the increase of the short wavelength fluctuation density within the arc.

[0028] As a further aspect of the present invention: S3 specifically includes:

[0029] Extract the mileage values ​​and corresponding elevation values ​​of all sampling points from the current arc segment's elevation point cloud subset, forming a discrete elevation sequence arranged along the mileage increasing direction;

[0030] Local extremum detection is performed on the discrete elevation sequence to extract all local maxima and minima within the arc segment and their corresponding mileage positions. The mileage interval between adjacent local maxima and minima is calculated to form the undulation cycle sequence of the arc segment.

[0031] The frequency of occurrence of each mileage interval value in the fluctuation cycle sequence is statistically analyzed, and the mileage interval value with the highest frequency is taken as the main fluctuation cycle of the arc segment.

[0032] Using the main fluctuation period as the fundamental wavelength, the discrete elevation sequence is fitted by superimposing the fundamental wave and its integer multiple harmonics to obtain a continuous elevation function with mileage as the independent variable, which serves as the continuous elevation change curve of the arc segment.

[0033] As a further aspect of the present invention: the process of obtaining the continuously varying elevation curve further includes:

[0034] Using the reciprocal of the main fluctuation period as the fundamental frequency, a fundamental cosine function with mileage as the independent variable is constructed, and the fundamental amplitude is determined by minimizing the elevation fitting residual of the fundamental cosine function at discrete points.

[0035] The elevation values ​​of the fundamental cosine function at each discrete point are subtracted from the discrete high-order sequence to obtain the first residual sequence. Then, each order of harmonic function is constructed in integer multiples of the fundamental frequency. The residual sequence is fitted one order at a time. After each fitting, the elevation value of the current order harmonic function is subtracted from the current residual sequence. The residual sequence is updated and then the next order is entered. The process continues until the root mean square value of the updated residual sequence is less than the preset tolerance.

[0036] By summing the elevation values ​​of the fundamental cosine function and all the determined harmonic functions at each mileage value, the continuous elevation variation curve of the arc segment is obtained.

[0037] As a further aspect of the present invention: the formation process of the non-uniform point grid is as follows:

[0038] Mark the arc segment interval where the slope change rate exceeds the threshold as the interval to be subdivided, and read the slope change rate value at each mileage position in the interval to be subdivided, and record the peak positions where the slope change rate reaches the maximum value and their corresponding mileage coordinates.

[0039] Using each peak position as a subdivision breakpoint, the interval to be subdivided is divided into multiple initial sub-intervals along the mileage direction, so that each initial sub-interval contains only one peak position;

[0040] Calculate the average slope change rate within each initial sub-interval, and use the ratio of the average value to the threshold as the re-segmentation coefficient of the initial sub-interval. Based on the re-segmentation coefficient, re-segment each initial sub-interval equally, so that the initial sub-interval with the higher average slope change rate has more re-segmentation segments.

[0041] All sub-intervals obtained after further subdivision of each initial sub-interval are arranged sequentially to replace the original sub-intervals to be subdivided, forming an updated non-uniform point grid.

[0042] As a further aspect of the present invention: the step of dividing the interval to be subdivided into multiple initial sub-intervals along the mileage direction specifically includes:

[0043] Read the mileage coordinates corresponding to all peak positions in the subdivided interval in ascending order of mileage, and form a peak mileage sequence.

[0044] Using the starting mileage position of the interval to be subdivided as the first cutting starting point, and taking each mileage coordinate in the peak mileage sequence as the cutting breakpoint, the cutting operation is performed at each cutting breakpoint along the direction perpendicular to the mileage, and the interval to be subdivided is successively divided into multiple mileage continuous sub-segments.

[0045] Each segment is defined as an initial sub-interval, and each initial sub-interval contains only one peak position in the peak mileage sequence. At the same time, the remaining mileage segment between the last peak position in the peak mileage sequence and the end mileage position of the subdivided interval is defined as the final initial sub-interval.

[0046] As a further aspect of the present invention: S5 specifically includes:

[0047] The entire route is divided into multiple terrain feature segments. Each terrain feature segment is a continuous mileage interval in which the sign of the slope change rate remains unchanged. The set of collocation points contained in each terrain feature segment is the set of collocation points for solving the corresponding terrain feature segment.

[0048] The sum of the absolute values ​​of the velocity differences between adjacent points within each terrain feature segment is used as the total velocity fluctuation of the terrain feature segment. The total velocity fluctuation of each terrain feature segment is used as the allocation weight of the initial velocity value of each point within the terrain feature segment. The overall weight coefficient of the penalty term is allocated to each point according to the weight, forming the differentiated penalty intensity of each point.

[0049] Under the condition of satisfying all constraints, the velocity value of each collocation point is adjusted in an iterative manner. In each round of iteration, the marginal violation of each constraint under the current velocity sequence is used as the reverse guide for the corresponding round velocity correction amount. The process continues until the change in the velocity value of each collocation point is less than the preset convergence tolerance, and then the velocity value at each collocation point is output.

[0050] The beneficial effects of this invention are:

[0051] (1) By collecting elevation point clouds with an interval of no more than five meters between adjacent measuring points, a coarse-grained topological map composed of nodes and arc segments is constructed, and a corresponding subset of elevation point clouds is encapsulated within each arc segment, thus realizing hierarchical storage of terrain data. Specifically, the coarse-grained topological map only records the connection relationship between nodes and arc segments for path search, while the high-density elevation point cloud is only encapsulated within a limited number of arc segments selected by the initial path, avoiding the direct loading of all high-precision measuring point data into the path search space. Compared with the scheme of directly using all measuring points as graph nodes for search, the present invention reduces the amount of data processed in the search stage. At the same time, since the elevation point cloud data is completely preserved within the arc segments, the original elevation information with an interval of no more than five meters within each arc segment can still be restored in the terrain reconstruction stage, taking into account both the computational efficiency of global path search and the accurate expression of micro-terrain undulation features.

[0052] (2) The arc segment is subdivided differently based on the slope change rate. Dense points are formed in areas where the slope change rate exceeds the threshold, while sparse points are maintained in areas where the slope change rate does not exceed the threshold. The resulting non-uniform point grid matches the actual terrain undulation characteristics of the tunnel. This invention divides the entire process into multiple terrain feature segments with unchanged slope change rate signs. Differentiated penalty intensity is allocated to each point according to the total speed fluctuation within each terrain feature segment, so that the constraint adjustment intensity of the speed optimization process is adapted to the actual terrain undulation degree of each segment. Through the above-mentioned terrain feature-driven point distribution and penalty intensity allocation method, the speed change rate between adjacent points in the speed command sequence can be adaptively adjusted according to the severity of the actual terrain undulation, avoiding the speed curve mismatch problem caused by insufficient capture of local terrain details under the fixed-interval point method. Attached Figure Description

[0053] The invention will now be further described with reference to the accompanying drawings.

[0054] Figure 1 This is a flowchart of the method of the present invention;

[0055] Figure 2 This is a flowchart illustrating the formation process of the initial path sequence in this invention.

[0056] Figure 3 This is a flowchart of the process for obtaining the continuously changing elevation curve in this invention. Detailed Implementation

[0057] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0058] Please see Figure 1 As shown, this invention is a method for global energy consumption optimization path planning for monorail cranes based on mine roadway topology, including the following steps:

[0059] S1 collects the three-dimensional coordinates of the centerline of the mine roadway, the elevation of the slope change point, the radius of curvature of the curve, the position of the turnout and the division of the power supply section, obtains the elevation point cloud with an interval of no more than five meters between adjacent measuring points, constructs a coarse-grained topology map composed of nodes and arc segments, and encapsulates the corresponding elevation point cloud subset inside each arc segment;

[0060] S2 uses a coarse-grained topological map as the search space and the elevation difference between the start and end points of each arc segment as the basic cost to search for a preliminary path sequence from the starting node to the target node.

[0061] S3. Along the initial path sequence, retrieve the subset of elevation point clouds encapsulated within each arc segment, reconstruct the continuous elevation change curve of each arc segment along the mileage direction, and calculate the derivative of the slope angle relative to the mileage at each location as the slope change rate.

[0062] S4. Based on the slope change rate values ​​at each location, a threshold is set, and the arc segments with slope change rates exceeding the threshold are subdivided along the mileage direction. The segments with slope change rates not exceeding the threshold are kept at their original span, forming a non-uniform point grid along the entire route.

[0063] S5 uses each point on a non-uniform grid as the discretization benchmark and the elevation change curve as the terrain input. The sum of the squares of the second derivatives of the running speed is superimposed as a penalty term on the traction net energy consumption optimization objective. Under the constraints of traction force, braking deceleration and curve centrifugal force, the speed values ​​of each point are solved to minimize the optimization objective, and the speed command sequence arranged along the mileage is output.

[0064] In S1, the three-dimensional coordinates of the centerline of the mine roadway, the elevation of the slope change point, the radius of curvature of the curve, the location of the turnout, and the division of the power supply section are collected. Elevation point clouds with an interval of no more than five meters between adjacent measuring points are obtained. A coarse-grained topological map composed of nodes and arc segments is constructed, and a corresponding subset of elevation point clouds is encapsulated within each arc segment. Specifically, this includes:

[0065] A mine explosion-proof total station was used to conduct traverse surveys along the centerline of the roadway. A sampling point was collected every 5 meters, and the three-dimensional coordinates of each sampling point in the mine's independent coordinate system were recorded. The X and Y axes of the three-dimensional coordinates represent the horizontal position, and the Z axis represents the elevation value. For curved sections, the sampling interval was increased to 2 meters. All sampling points were arranged in chronological order along the roadway direction, forming the original measurement point sequence for each roadway. The Z-axis coordinates of all sampling points in the original measurement point sequence of each roadway were extracted in ascending order of mileage to form the elevation sequence for that roadway. Each element in the elevation sequence corresponds to a unique mileage position.

[0066] For each roadway's elevation sequence, starting from the second sampling point in the sequence, the difference between the elevation value of the current sampling point and the elevation value of the previous sampling point is calculated sequentially. This difference is then divided by the mileage interval between the two sampling points. The results are arranged sequentially to form the roadway's slope value sequence, where each slope value corresponds to a sampling point location. Similarly, starting from the second slope value in the slope value sequence, the difference between the current slope value and the previous slope value is calculated sequentially. This difference is then divided by the mileage interval between the two sampling points corresponding to the two slope values. The results are arranged sequentially to form the roadway's curvature value sequence, where each curvature value corresponds to a sampling point location.

[0067] The preset curvature threshold is 0.05 per meter, determined based on the maximum allowable rate of gradient change in the mine transport roadway design specifications. The curvature value sequence of each roadway is traversed, and sampling points with curvature values ​​greater than 0.05 per meter are marked as terrain abrupt change points. Along the roadway direction, the continuous path between two adjacent terrain abrupt change points is defined as an arc segment, and the intersections of terrain abrupt change points and two or more roadways are defined as nodes. Each arc segment is uniquely determined by its starting and ending nodes, and each node is associated with its three-dimensional coordinate values.

[0068] For each arc segment, the 3D coordinates and corresponding mileage values ​​of all sampling points between the start and end points of the arc segment are extracted to form the original elevation point cloud data of that arc segment. The original elevation point cloud data is then compressed. Specifically, starting from the first sampling point within the arc segment, the 3D coordinates and mileage values ​​of every other sampling point are retained in ascending order of mileage. This retained data is then used as a subset of the elevation point cloud for that arc segment. All nodes and arc segments are organized into a network structure based on the connection relationships of the tunnels. Each node in this network structure records its node number and 3D coordinates, and each arc segment records its start node number, end node number, arc length, and the corresponding subset of elevation point clouds. This network structure constitutes the coarse-grained topology map.

[0069] Please see Figure 2 As shown, in S2, using a coarse-grained topology map as the search space and the elevation difference between the start and end points of each arc segment as the basic cost, a preliminary path sequence from the starting node to the target node is obtained, specifically including:

[0070] After obtaining the coarse-grained topology map, a preliminary path sequence search is performed. The node where the monorail is currently located is set as the starting node, and the node where the unloading point specified by the transportation task is located is set as the target node. Using the starting node as the current search node, all adjacent nodes with arc connection relationships with the current search node are read from the coarse-grained topology map, and the elevation point cloud subsets associated with each of these adjacent arc segments are retrieved.

[0071] For each adjacent arc segment, fluctuation amplitude calculation is performed, the specific process of which will be detailed later. The starting elevation value and ending elevation value of each adjacent arc segment are subtracted, and the absolute value is taken to obtain the elevation difference between the starting and ending points of that arc segment. This elevation difference is multiplied by a first weighting coefficient to obtain a weighted elevation difference value. Simultaneously, the fluctuation amplitude is multiplied by a second weighting coefficient to obtain a weighted fluctuation amplitude value. The weighted elevation difference value and the weighted fluctuation amplitude value are then added together, and the result is used as the passage cost for that arc segment.

[0072] The first weighting coefficient is 0.6, the second weighting coefficient is 0.4, and the sum of the two weighting coefficients is 1. The adjacent arc segment with the smallest passage cost is selected as the preferred arc segment. The end node of the preferred arc segment is updated to the current search node, and the arc segment number of the preferred arc segment is recorded in the path sequence in sequence.

[0073] Repeat the above operations of traversing adjacent arcs, calculating the passage cost, selecting the preferred arc, updating the current node and recording it in the path sequence until the current search node is the target node. At this time, the path sequence contains the numbers of each arc segment passed from the starting node to the target node.

[0074] Perform a correction operation on the path sequence: Traverse the travel cost of each arc segment in the path sequence, mark the arc segments whose travel cost is greater than a preset cost limit as arc segments to be removed, the preset cost limit is 1.5 times the average travel cost of all arc segments in the path sequence; delete the arc segments to be removed from the path sequence to form a discontinuous path, search for alternative arc segments connecting the endpoints of each segment of the discontinuous path in the coarse-grained topology map, select the arc segment with the smallest travel cost from all alternative arc segments and insert it into the corresponding position of the path sequence to restore the continuity of the path sequence; output the corrected path sequence.

[0075] The specific calculation process for the above fluctuation amplitude is as follows:

[0076] Read the elevation values ​​and corresponding mileage values ​​of all sampling points from the elevation point cloud subset associated with the current arc segment, and arrange them in ascending order of mileage to form the original elevation sequence. Traverse each sampling point in the original elevation sequence, starting from the second sampling point and ending at the second-to-last sampling point, calculate the first difference between the elevation value of the current sampling point and the elevation value of its previous adjacent sampling point, and simultaneously calculate the second difference between the elevation value of the current sampling point and the elevation value of its next adjacent sampling point. Take the larger absolute value of the first difference and the second difference as the neighborhood change value of the sampling point.

[0077] When the change value in the neighborhood is greater than three times the change value in the neighborhood of its previous adjacent sampling point, the sampling point is marked as an outlier and removed from the original height sequence. The remaining sampling points in the height sequence after removing outliers are rearranged according to their original mileage order to form the effective height sequence for that arc segment. Starting from the second sampling point in the effective height sequence, the difference between the elevation value of the current sampling point and the elevation value of the previous sampling point is calculated sequentially. Each difference corresponds to a difference between two adjacent sampling points on the arc segment. The absolute values ​​of all the difference values ​​obtained along the arc segment direction are taken and summed to obtain the cumulative sum of the absolute values ​​of the difference values.

[0078] The total number of sampling points in the valid elevation sequence is counted. Subtracting 1 from this total number of sampling points gives the denominator. The absolute values ​​of the differences are summed to obtain the numerator, and division is performed to obtain the average absolute difference fluctuation amplitude for that arc segment. The arc length is read, and dividing the arc length by 100 defines the normalization coefficient for that arc segment. The average absolute difference fluctuation amplitude is multiplied by this normalization coefficient; the product is the final fluctuation amplitude for that arc segment. The final fluctuation amplitude increases with the frequency of alternating elevation changes along the mileage direction within the arc segment.

[0079] Please see Figure 3 As shown, in S3, along the initial path sequence, the subset of elevation point clouds encapsulated within each arc segment is retrieved, the continuous elevation change curve of each arc segment along the mileage direction is reconstructed, and the derivative of the slope angle relative to the mileage at each location is calculated as the slope change rate. Specifically, this includes:

[0080] After obtaining the corrected preliminary path sequence, the reconstruction operation of the continuously varying elevation curves is performed sequentially along each arc segment contained in the path sequence. The arc segment to be processed is extracted from the path sequence, and the associated elevation point cloud subset is read. From this subset, the mileage values ​​and their corresponding elevation values ​​of all sampling points are extracted and arranged in ascending order of mileage, forming a discrete elevation sequence for that arc segment. Each data point in this discrete elevation sequence contains two values: the mileage coordinates and the elevation value. The total number of data points in the sequence is equal to the total number of sampling points in the elevation point cloud subset of that arc segment.

[0081] Local extremum detection is performed on the above discrete high-order sequence. The specific process of local extremum detection is as follows:

[0082] Starting from the second data point of the discrete high-order sequence, iterate sequentially to the second-to-last data point. For the currently traversed data point, compare the elevation value of the point with the elevation value of the previous data point and the elevation value of the next data point. If the elevation value of the point is greater than both the elevation value of the previous data point and the elevation value of the next data point, then mark the point as a local maximum point.

[0083] If the elevation value of a point is less than both the elevation values ​​of the preceding and following data points, then the point is marked as a local minimum.

[0084] After traversal, extract all data points marked as local maxima and local minima within the arc segment, and record the mileage position of each extreme point and its corresponding elevation value.

[0085] All extracted local maxima and minima are arranged in ascending order of mileage, forming an alternating sequence of extreme points where adjacent extreme points are of different types. Starting from the first extreme point in this sequence, the mileage difference between each extreme point and its next adjacent extreme point is calculated sequentially. All calculated mileage differences form a fluctuation cycle sequence for the arc segment, where each value represents the mileage interval between an adjacent peak and trough within the arc segment. The frequency of each mileage interval value in the fluctuation cycle sequence is counted, and the mileage interval value with the most frequent occurrences is selected as the main fluctuation cycle of the arc segment.

[0086] After obtaining the main fluctuation period, a fitting operation is performed on the discrete height sequence. Using the reciprocal of the main fluctuation period as the fundamental frequency, a fundamental cosine function with mileage as the independent variable is constructed. The amplitude of this fundamental cosine function is determined by minimizing the fitting residuals. Specifically, the mileage coordinates of each data point in the discrete height sequence of the arc segment are substituted into the fundamental cosine function to obtain the fundamental fitted elevation value at each data point. The original elevation value at each data point is subtracted from the fundamental fitted elevation value to obtain the fitting residual at each data point. The sum of squares of the fitting residuals at all data points is calculated, and the fundamental amplitude is adjusted to minimize this sum of squares. The amplitude at this minimum is the fundamental amplitude. The fundamental cosine function is represented as follows: ;

[0087] in, Indicates mileage The fundamental cosine function value at that point, Indicates the fundamental frequency amplitude. This represents the mileage value along the arc segment. This indicates the main fluctuation cycle.

[0088] After determining the fundamental amplitude, the original elevation value of each data point in the discrete elevation sequence is subtracted from the calculated value of the fundamental cosine function at that data point. The difference is used as the first-order residual value for that data point. All the first-order residual values ​​are arranged in mileage order to form a first-order residual sequence. Subsequently, the harmonic functions are fitted sequentially. The specific process is as follows: For the nth harmonic, the frequency of the harmonic is taken as n times the fundamental frequency. A cosine function of the nth harmonic with mileage as the independent variable is constructed. The amplitude of this harmonic function is determined by minimizing the fitting residual of the current residual sequence. The calculated value of this harmonic function at each mileage point in the current residual sequence is subtracted from the corresponding data point in the current residual sequence to obtain the updated residual sequence. After fitting each harmonic, the sum of squares of the values ​​of each data point in the updated residual sequence is calculated, divided by the total number of data points, and the square root is taken to obtain the root mean square value of the updated residual sequence. The preset tolerance value is 0.05 meters, which is determined based on the nominal accuracy of the elevation measurement of the mining total station. When the root mean square value of the updated residual sequence is less than 0.05 meters, the stepwise fitting operation is terminated, and the amplitude and order of all determined harmonic functions are recorded.

[0089] The elevation values ​​at each mileage are summed by adding the fundamental cosine function to the elevation values ​​of all determined harmonic functions. The sum is then used as the elevation change curve value at that mileage. This reconstruction operation is performed on each arc segment in the path sequence to obtain the elevation change curve for each segment. After reconstructing the elevation change curve for each arc segment, the slope change rate is calculated for each elevation change curve.

[0090] The slope change rate is calculated as follows: Within the arc segment to be calculated, starting from the initial mileage position of the arc segment, points are sequentially selected along the mileage direction with a step size of 0.1 meters. For each selected point, the slope of the tangent line to the continuous elevation change curve at that point is calculated, and this slope value is taken as the slope angle value at that position. Starting from the second selected point, the slope angle value at the current selected point is subtracted from the slope angle value at the previous selected point to obtain the slope angle difference value. This slope angle difference value is divided by 0.1 meters, and the result is taken as the slope change rate at the current selected point position. The slope change rate values ​​at each selected point position are used for subsequent non-uniform point grid division.

[0091] In S4, based on the slope change rate values ​​at each location, a threshold is set. Arc segments with slope change rates exceeding the threshold are subdivided along the mileage direction, while segments with slope change rates below the threshold retain their original span, forming a non-uniform point grid along the entire route. Specifically, this includes:

[0092] Set a slope change rate threshold. The calibration basis of this threshold is the ratio of the DC bus voltage fluctuation tolerance value of the monorail inverter to the thermal accumulation time constant of the braking resistor. The specific calibration process is as follows: read the rated value of the DC bus voltage of the monorail inverter, multiply the rated value by 5% to obtain the voltage fluctuation tolerance value, read the time required for the braking resistor to rise from the ambient temperature to the maximum allowable operating temperature as the thermal accumulation time constant, and divide the voltage fluctuation tolerance value by the thermal accumulation time constant as the slope change rate threshold. The unit of this threshold is degrees per meter.

[0093] The slope change rate values ​​at each location along the entire route are compared with the threshold. All continuous mileage intervals with slope change rate values ​​greater than the threshold are marked as intervals to be subdivided, and continuous mileage intervals with slope change rate values ​​not greater than the threshold are marked as maintenance intervals. The grid span of the points within the maintenance intervals remains unchanged from the original arc length.

[0094] For each marked subdivision interval, perform an initial subdivision operation. Scan all locations within the subdivision interval in ascending order of mileage, recording the slope change rate at each location. Compare the slope change rate at each location with the slope change rates at its preceding and following locations. If the current location's value is greater than both its preceding and following values, mark that location as a peak position and extract the mileage coordinates corresponding to all peak positions. Arrange the mileage coordinates of all peak positions in ascending order to form the peak mileage sequence corresponding to the subdivision interval.

[0095] Using the starting mileage position of the interval to be subdivided as the first cutting starting point, and the first mileage coordinate in the peak mileage sequence as the first cutting breakpoint, a cut is performed at the first cutting breakpoint along a direction perpendicular to the mileage direction, dividing the interval to be subdivided into a first sub-segment, which extends from the starting mileage position to the first cutting breakpoint. Using the first cutting breakpoint as the second cutting starting point, and the second mileage coordinate in the peak mileage sequence as the second cutting breakpoint, a cut is performed to obtain a second sub-segment. This process continues until the last mileage coordinate in the peak mileage sequence is used as the cutting breakpoint. Each segment obtained from each cut is defined as an initial sub-interval, containing only one peak position. The remaining mileage segment between the last peak position in the peak mileage sequence and the ending mileage position of the interval to be subdivided is defined as the final initial sub-interval, which does not contain any peak positions.

[0096] After the initial sub-intervals are divided, a re-segmentation coefficient is calculated for each initial sub-interval. For an initial sub-interval, the slope change rate values ​​at each location point within the sub-interval are read. All values ​​are summed and divided by the total number of locations within the sub-interval to obtain the average slope change rate within that initial sub-interval. This average value is then divided by a slope change rate threshold, and the result is used as the re-segmentation coefficient for that initial sub-interval. When the re-segmentation coefficient is less than 1, the initial sub-interval is marked as no longer to be divided, maintaining its original mileage span.

[0097] When the resegmentation coefficient is greater than or equal to 1, an equidistant resegmentation operation is performed on the initial sub-interval. Specifically, the integer part of the resegmentation coefficient is taken as the base segment number. The base segment number is incremented by 1 to obtain the actual number of resegmented segments for the initial sub-interval. The mileage span of the initial sub-interval is divided by the actual number of resegmented segments to obtain the mileage length of each resegmented segment. Starting from the initial mileage position of the initial sub-interval, points are sequentially selected according to the mileage length of each resegmented segment, uniformly dividing the initial sub-interval along the mileage direction into the actual number of resegmented sub-intervals. The mileage position of the endpoint of each resegmented sub-interval is used as a matching point. When taking the integer part of the resegmentation coefficient, if the decimal part is greater than or equal to 0.5, the integer part is rounded up by 1; if the decimal part is less than 0.5, the integer part remains unchanged. When the actual number of resegmented segments is greater than 4, the actual number of resegmented segments for the initial sub-interval is forcibly set to 4, meaning each initial sub-interval can be divided into a maximum of 4 resegmented sub-intervals.

[0098] The original endpoints of the arc segments within each maintenance interval are used as the point locations for that interval. The endpoints of all the subdivided sub-intervals and the endpoints of all the initial sub-intervals of the final segment are merged in ascending order of mileage. If the absolute difference between two mileage locations is less than 0.01 meters, they are merged into one, and the average of the two merged locations is taken. All the merged mileage locations are then arranged sequentially to form a non-uniform point grid along the entire length. In this non-uniform point grid, the mileage intervals between adjacent points are unequal; areas with a higher rate of slope change have denser point grids, while areas with a lower rate of slope change have sparser point grids.

[0099] In S5, each point on a non-uniform grid is used as the discretization reference, and the terrain input is a continuously varying elevation curve. The sum of squares of the second derivatives of the operating speed is added as a penalty term to the traction net energy consumption optimization objective. Under the constraints of traction force, braking deceleration, and centrifugal force in curves, the speed values ​​of each point are solved to minimize the optimization objective. The output is a speed command sequence arranged along the mileage, specifically including:

[0100] After obtaining the non-uniform point grid, the velocity command sequence is solved. First, the entire route is divided into multiple terrain feature segments. Each terrain feature segment is defined as a continuous mileage interval along the mileage direction where the sign of the slope change rate remains unchanged. That is, within a terrain feature segment, the slope change rate is either all positive or all negative. When the sign of the slope change rate changes from positive to negative or from negative to positive, the position where the sign changes is the boundary point between two adjacent terrain feature segments.

[0101] After dividing the entire route according to the above rules, multiple terrain feature segments are obtained continuously arranged along the mileage direction, with each feature segment connected end to end and not overlapping. For each terrain feature segment, all points between the starting and ending mileage positions of the feature segment are extracted from the non-uniform point grid. These points are used as the solution point set for the feature segment, and the velocity value at each solution point is the unknown quantity to be solved.

[0102] For each terrain feature segment, the total velocity fluctuation is calculated. The specific calculation method is as follows: all the collocation points in the terrain feature segment are arranged in ascending order of mileage. Starting from the second collocation point and ending at the last collocation point, the absolute value of the velocity difference between the current collocation point and the previous collocation point is calculated in turn. The absolute values ​​of each velocity difference are added together, and the sum is taken as the total velocity fluctuation of the terrain feature segment.

[0103] Read the total velocity fluctuation values ​​of all terrain feature segments throughout the entire route and find the maximum value. For any terrain feature segment, divide the total velocity fluctuation of that segment by the maximum value among all terrain feature segments throughout the entire route to obtain the normalized velocity fluctuation value of that segment. Multiply the normalized velocity fluctuation value of that segment by the overall weight coefficient of the penalty term, and the product is used as the segment penalty weight value of that terrain feature segment.

[0104] The overall weight coefficient of the penalty item is determined based on the product of the monorail crane's curb weight and the inverter's current loop response bandwidth. The specific calibration process is as follows: read the curb weight value of the monorail crane in tons, read the inverter's current loop response bandwidth value in Hertz, multiply the curb weight value and the response bandwidth value, and divide the product by 100 to obtain the overall weight coefficient of the penalty item.

[0105] Within each terrain feature segment, the segment penalty weight is evenly distributed to each point within that segment. Specifically, the number of points within the segment is used as the denominator, and the segment penalty weight is used as the numerator; division is then performed to obtain the differentiated penalty intensity for each point within that segment. The differentiated penalty intensity for each point is directly proportional to the total velocity fluctuation of the terrain feature segment to which that point belongs; terrain feature segments with larger total velocity fluctuations receive higher differentiated penalty intensities.

[0106] After completing the differentiated penalty intensity allocation for each point, an iterative solution operation for the speed value is performed. The initial speed value for each point is set to the rated operating speed of the monorail, which is 2 meters per second.

[0107] In each round of iterative solution, the net traction energy consumption value for the entire journey is first calculated under the current speed value of each point. This energy consumption value is calculated based on the elevation change curve value at each point, the speed value at that point, and the total mass of the monorail.

[0108] Simultaneously, the traction force constraint value, braking deceleration constraint value, and curve centrifugal force constraint value are calculated separately: The traction force constraint value is the traction force required at each set point by the slope resistance and acceleration resistance. This value must not exceed the maximum traction force output by the monorail motor. The maximum traction force value is determined based on the rated power and rated speed of the motor and is taken as 120 kN; The braking deceleration constraint value is the deceleration value corresponding to the rate of change of speed at each set point. The absolute value of this value must not exceed the upper limit of braking deceleration, which is taken as 1.5 m / s² according to the mine transportation safety regulations; The curve centrifugal force constraint value is the centrifugal force value calculated at each set point from the speed value and the radius of curvature of the curve. This value must not exceed the maximum centrifugal force allowed for the monorail when running on a curve. The maximum centrifugal force value is determined based on the structural strength of the monorail and the track clamping force and is taken as 80 kN.

[0109] The amount by which the traction force constraint value exceeds the maximum traction force value is defined as the traction force marginal violation. The amount by which the absolute value of the braking deceleration constraint value exceeds the upper limit of the braking deceleration value is defined as the braking marginal violation. The amount by which the cornering centrifugal force constraint value exceeds the maximum centrifugal force value is defined as the centrifugal force marginal violation. These three marginal violations are collectively referred to as the marginal violations of each constraint in this iteration. The absolute values ​​of the marginal violations of each constraint are summed to obtain the total marginal violation. The total marginal violation is multiplied by the correction step size coefficient to obtain the total speed correction value for the current iteration. The correction step size coefficient is set to 0.01. This total speed correction value is distributed to each point according to the inverse relationship of the differentiated penalty intensity of each point, that is, the point with the greater the differentiated penalty intensity receives the smaller speed correction value.

[0110] Each point subtracts its assigned speed correction from its current speed value to obtain its updated speed value. After each iteration, the difference between the speed values ​​of each point in the current iteration and those in the previous iteration is calculated. The maximum absolute value of all differences is taken, and the iteration terminates when this maximum absolute value is less than a preset convergence tolerance, which is 0.01 meters per second. After the iteration terminates, the mileage position of each point is associated with its corresponding speed value and stored, arranged in ascending order of mileage to form a speed command sequence arranged sequentially along the mileage direction, and then output.

[0111] The working principle of this invention is as follows: By collecting the three-dimensional coordinates of the centerline of the mine roadway, the elevation of the slope change point, the radius of curvature of the curve, the position of the turnout, and the division of the power supply section, and obtaining the elevation point cloud with an interval of no more than five meters between adjacent measuring points, a coarse-grained topological map composed of nodes and arc segments is constructed, and the corresponding elevation point cloud subset is encapsulated within each arc segment; using this coarse-grained topological map as the search space and the elevation difference between the start and end points of each arc segment as the basic cost, a preliminary path sequence from the starting node to the target node is obtained; along this preliminary path sequence, the elevation point cloud subsets encapsulated within each arc segment are retrieved, the continuous elevation change curve of each arc segment along the mileage direction is reconstructed, and the elevation change curves of each arc segment are calculated. The derivative of the slope angle relative to the mileage at each location is used as the slope change rate. A threshold is set based on the slope change rate values ​​at each location. The arc segments with slope change rates exceeding the threshold are subdivided along the mileage direction, while the segments with slope change rates not exceeding the threshold retain their original span, forming a non-uniform grid of points along the entire route. Each point on this non-uniform grid is used as the discretization reference, and the continuous elevation change curve is used as the terrain input. The sum of the squares of the second derivatives of the running speed is added as a penalty term to the traction net energy consumption optimization objective. Under the constraints of traction force, braking deceleration, and centrifugal force in curves, the speed values ​​of each point are solved to minimize the optimization objective, and the speed command sequence arranged along the mileage is output.

[0112] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A method for global energy-optimal path planning for monorail cranes based on mine roadway topology, characterized in that, Includes the following steps: S1 collects the three-dimensional coordinates of the centerline of the mine roadway, the elevation of the slope change point, the radius of curvature of the curve, the position of the turnout and the division of the power supply section, obtains the elevation point cloud with an interval of no more than five meters between adjacent measuring points, constructs a coarse-grained topology map composed of nodes and arc segments, and encapsulates the corresponding elevation point cloud subset inside each arc segment; S2 uses a coarse-grained topological map as the search space and the elevation difference between the start and end points of each arc segment as the basic cost to search for a preliminary path sequence from the starting node to the target node. S3. Along the initial path sequence, retrieve the subset of elevation point clouds encapsulated within each arc segment, reconstruct the continuous elevation change curve of each arc segment along the mileage direction, and calculate the derivative of the slope angle relative to the mileage at each location as the slope change rate. S4. Based on the slope change rate values ​​at each location, a threshold is set, and the arc segments with slope change rates exceeding the threshold are subdivided along the mileage direction. The segments with slope change rates not exceeding the threshold are kept at their original span, forming a non-uniform point grid along the entire route. S5 uses each point on a non-uniform grid as the discretization benchmark and the elevation change curve as the terrain input. The sum of the squares of the second derivatives of the running speed is superimposed as a penalty term on the traction net energy consumption optimization objective. Under the constraints of traction force, braking deceleration and curve centrifugal force, the speed values ​​of each point are solved to minimize the optimization objective, and the speed command sequence arranged along the mileage is output.

2. The method for global energy consumption optimization path planning of monorail cranes based on mine roadway topology according to claim 1, characterized in that, The construction of a coarse-grained topology map composed of nodes and arc segments, and the encapsulation of corresponding elevation point cloud subsets within each arc segment, specifically includes: Using the centerline of the tunnel as the baseline, sampling points are picked at equal intervals along the direction of the baseline to obtain the three-dimensional coordinates and corresponding elevation values ​​of each sampling point, forming an elevation sequence along the tunnel direction. Perform a first-order difference operation on the high-order sequence to obtain the slope value sequence at each sampling point, and perform a second-order difference operation on the slope value sequence to obtain the curvature value at each sampling point; Sampling points with curvature values ​​exceeding a preset curvature threshold are selected as terrain change points. The continuous path between adjacent terrain change points in each tunnel section is defined as an arc segment, and terrain change points and tunnel intersections are defined as nodes. The elevation point cloud data of all sampling points within each arc segment is compressed and encapsulated into an elevation point cloud subset of the arc segment. The network structure consisting of nodes and arc segments is stored as a coarse-grained topology map, in which each arc segment is associated with the node number of its starting and ending points, the arc segment length, and the corresponding elevation point cloud subset.

3. The method for global energy consumption optimization path planning of monorail cranes based on mine roadway topology according to claim 1, characterized in that, The formation process of the preliminary path sequence is as follows: Starting from the current node, traverse all arcs adjacent to the current node, retrieve the elevation point cloud subset corresponding to each arc, and calculate the fluctuation amplitude of the elevation within each arc. The weighted sum of the elevation difference between the start and end points of each arc segment and the fluctuation amplitude is used as the passage cost of the arc segment, and the adjacent arc segment with the minimum passage cost is selected as the preferred arc segment. Update the end node of the preferred arc segment to the current node, and record the preferred arc segment in the path sequence. Repeat the traversal and selection process until the current node is the target node. Delete arcs in the path sequence whose travel cost exceeds a preset limit, replace them with the arcs with the lowest cost among the remaining arcs, and output the corrected path sequence.

4. The method for global energy consumption optimization path planning of monorail cranes based on mine roadway topology according to claim 3, characterized in that, The calculation process for the fluctuation amplitude is as follows: The elevation values ​​of all sampling points in the current arc segment are read from the elevation point cloud sub-set. After removing outliers whose elevation difference with adjacent points exceeds a predetermined multiple, the effective elevation sequence of the arc segment is generated. Perform differential operations on the elevation values ​​of each adjacent point in the effective elevation sequence to obtain multiple differential values ​​along the arc direction, and calculate the sum of the absolute values ​​of the differential values; Calculate the ratio of the sum of the absolute values ​​of the differences to the total number of sampling points within the arc segment minus one, and use this ratio as the average absolute difference fluctuation amplitude of the arc segment; Multiplying the average absolute difference fluctuation amplitude by the normalization coefficient of the arc length yields the final fluctuation amplitude of the arc. The final fluctuation amplitude increases monotonically with the increase of the short wavelength fluctuation density within the arc.

5. The method for global energy consumption optimization path planning of monorail cranes based on mine roadway topology according to claim 1, characterized in that, S3 specifically includes: Extract the mileage values ​​and corresponding elevation values ​​of all sampling points from the current arc segment's elevation point cloud subset, forming a discrete elevation sequence arranged along the mileage increasing direction; Local extremum detection is performed on the discrete elevation sequence to extract all local maxima and minima within the arc segment and their corresponding mileage positions. The mileage interval between adjacent local maxima and minima is calculated to form the undulation cycle sequence of the arc segment. The frequency of occurrence of each mileage interval value in the fluctuation cycle sequence is statistically analyzed, and the mileage interval value with the highest frequency is taken as the main fluctuation cycle of the arc segment. Using the main fluctuation period as the fundamental wavelength, the discrete elevation sequence is fitted by superimposing the fundamental wave and its integer multiple harmonics to obtain a continuous elevation function with mileage as the independent variable, which serves as the continuous elevation change curve of the arc segment.

6. The method for global energy consumption optimization path planning of monorail cranes based on mine roadway topology according to claim 5, characterized in that, The process of obtaining the continuously varying elevation curve further includes: Using the reciprocal of the main fluctuation period as the fundamental frequency, a fundamental cosine function with mileage as the independent variable is constructed, and the fundamental amplitude is determined by minimizing the elevation fitting residual of the fundamental cosine function at discrete points. The elevation values ​​of the fundamental cosine function at each discrete point are subtracted from the discrete high-order sequence to obtain the first residual sequence. Then, each order of harmonic function is constructed in integer multiples of the fundamental frequency. The residual sequence is fitted one order at a time. After each fitting, the elevation value of the current order harmonic function is subtracted from the current residual sequence. The residual sequence is updated and then the next order is entered. The process continues until the root mean square value of the updated residual sequence is less than the preset tolerance. By summing the elevation values ​​of the fundamental cosine function and all the determined harmonic functions at each mileage value, the continuous elevation variation curve of the arc segment is obtained.

7. The method for global energy consumption optimization path planning of monorail cranes based on mine roadway topology according to claim 1, characterized in that, The formation process of the non-uniform point grid is as follows: Mark the arc segment interval where the slope change rate exceeds the threshold as the interval to be subdivided, and read the slope change rate value at each mileage position in the interval to be subdivided, and record the peak positions where the slope change rate reaches the maximum value and their corresponding mileage coordinates. Using each peak position as a subdivision breakpoint, the interval to be subdivided is divided into multiple initial sub-intervals along the mileage direction, so that each initial sub-interval contains only one peak position; Calculate the average slope change rate within each initial sub-interval, and use the ratio of the average value to the threshold as the re-segmentation coefficient of the initial sub-interval. Based on the re-segmentation coefficient, re-segment each initial sub-interval equally, so that the initial sub-interval with the higher average slope change rate has more re-segmentation segments. All sub-intervals obtained after further subdivision of each initial sub-interval are arranged sequentially to replace the original sub-intervals to be subdivided, forming an updated non-uniform point grid.

8. The method for global energy consumption optimization path planning of monorail cranes based on mine roadway topology according to claim 7, characterized in that, The process of dividing the interval to be subdivided into multiple initial sub-intervals along the mileage direction specifically includes: Read the mileage coordinates corresponding to all peak positions in the subdivided interval in ascending order of mileage, and form a peak mileage sequence. Using the starting mileage position of the interval to be subdivided as the first cutting starting point, and taking each mileage coordinate in the peak mileage sequence as the cutting breakpoint, the cutting operation is performed at each cutting breakpoint along the direction perpendicular to the mileage, and the interval to be subdivided is successively divided into multiple mileage continuous sub-segments. Each segment is defined as an initial sub-interval, and each initial sub-interval contains only one peak position in the peak mileage sequence. At the same time, the remaining mileage segment between the last peak position in the peak mileage sequence and the end mileage position of the subdivided interval is defined as the final initial sub-interval.

9. The method for global energy consumption optimization path planning of monorail cranes based on mine roadway topology according to claim 1, characterized in that, S5 specifically includes: The entire route is divided into multiple terrain feature segments. Each terrain feature segment is a continuous mileage interval in which the sign of the slope change rate remains unchanged. The set of collocation points contained in each terrain feature segment is the set of collocation points for solving the corresponding terrain feature segment. The sum of the absolute values ​​of the velocity differences between adjacent points within each terrain feature segment is used as the total velocity fluctuation of the terrain feature segment. The total velocity fluctuation of each terrain feature segment is used as the allocation weight of the initial velocity value of each point within the terrain feature segment. The overall weight coefficient of the penalty term is allocated to each point according to the weight, forming the differentiated penalty intensity of each point. Under the condition of satisfying all constraints, the velocity value of each collocation point is adjusted in an iterative manner. In each round of iteration, the marginal violation of each constraint under the current velocity sequence is used as the reverse guide for the corresponding round velocity correction amount. The process continues until the change in the velocity value of each collocation point is less than the preset convergence tolerance, and then the velocity value at each collocation point is output.

Images