Heuristic sequencing methods for time optimal tracking of nested, open and closed paths

Tracking sequences of predefined open and closed paths is of crucial interest for applications like laser cutting and similar production processes. The disconnected paths are connected by non-productive, four times continuously differentiable trajectories, which also account for the overall process time. Heuristic methods are applied in order to find a proper sequencing of the open and closed path and thereby minimize the overall process time while respecting constraints given by the system limits. To this end, the exact traversing times of the non-productive linking trajectories are computed, which also have to be time optimal subject to the system limits. While problems with only closed paths present can be formulated as travelling salesman problem, handling open and nested paths introduce additional constraints. Two heuristic algorithms for travelling salesman problems are extended for open path processing and compared with respect to solution quality and calculation time using randomly generated problems. Finally, problems with nested cutting path structures are addressed by extending the heuristic algorithm with the best performance to account for the precedence constraints, associated with these nested structures.


Introduction
For laser cutting applications a cutting job consists of predefined open and closed paths.Since these paths are often not connected, the question arises, how these paths should be sorted in order to minimize overall process time.Numerous contributions, concerning this Fig. 1 Gantry laser cutting machine with examplary task topic, can be found in the literature, as indicated by the review paper [4], where 72 contributions are collected and compared with respect to the used methods and algorithms.
Although a more general approach would be desirable, this paper focuses on paths in a plane and gantry like robotic systems as shown in Fig. 1.These systems are subject to restrictions like maximum velocity, acceleration, and jerk for each axis, respectively.Furthermore there are process specific constraints.In case of the laser cutting process, the maximum velocity tangentially to the path depends on the material which has to be cut, and in order to ensure a clean cut, the velocity at the beginning and at the end of each path has to be zero.For each individual path, a time optimal path tracking solution satisfying these constraints can be found.These can be connected by time optimal connecting trajectories along straight lines.With an increasing number of paths to be tracked the impact of the nonproductive traversing time introduced by these links on the overall process time is getting dominant if the sequencing is not handled properly.The sequencing for problems with only closed cutting paths present can be formulated as travelling salesman problem.Therefore two heuristic algorithms for solving travelling salesman problems are presented, one wellknown construction heuristics, the Christofides algorithm [1], and one iteratively improving heuristics, the Lin-Kernighan-Helsgaun algorithm [7].These two algorithms are then extended for problems with open cutting paths and compared with respect to an integer linear programming algorithm.The optimal traversing time of the non-productive connecting trajectories depends non-linearly on the end points as well as on the system limits.In contrast to e.g.[3], where a piecewise linear function approximating this non-linear behaviour is introduced, in this paper analytic, non-linear expressions are derived and used to achieve the optimal traversing time.
If some of the cutting paths fully enclose other open or closed cutting paths, special care has to be taken to ensure that the enclosing paths are processed after all the inner paths are cut, resulting in precedence constraints for the sequencing algorithms.This is because once a closed path is cut, the interior is no longer connected to the raw material and it may be slightly shifted, introducing inaccuracies once the enclosed paths are cut.This type of problem is often refered as sequential ordering problem or as traveling salesman problem with precedence constraints and is addressed e.g. in [10] by a branch-and-bound algorithm for small problem sizes or in [6] by an ant colony system.In order to detect such nested cutting paths structures an algorithm is proposed, based on the winding number algorithm [8].Finally, the provided heuristics are examined, whether and how they can be extended to solve problems with precedence constraints.

Problem definition and combinatorial background
Before setting up the problem definition, some graph theoretical terminology should be introduced, which is used through out the subsequent sections.An undirected graph is defined by a set of nodes N = {n i | 0 ≤ i < N}, where N denotes the number of nodes, and a set of edges E ⊆ {e i,j | 0 ≤ i < N, i < j < N} connecting some or all of these nodes, where edge e i,j is equivalent to e j,i .Following the notion of e.g.[12], an undirected graph can be -weighted, if each edge e i,j is associated with a weight w i,j .
-simple, if each edge e i,j is unique in E.
-an undirected multigraph, if edges are allowed to be not unique in E.
-connected, if each pair of nodes can be connected by a sequence of edges in E -complete, if each pair of nodes is connected by exactly one edge.
A cycle in the graph is a sequence of edges, which starts and ends at the same node.Special cycles are Eulerian cycles, which contain every edge of a graph once, and Hamilton cycles, which start and end at the same node and traverse every other node of the graph exactly once.A Hamilton cycle is also called a tour of the graph.
A spanning tree of a connected graph is a subset of edges, which does not contain cycles and contains all nodes.Consequently the minimum spanning tree of a connected, weighted graph is the spanning tree with the minimum total edge weight.
Finding the tour with the minimum total edge weight in a complete, undirected, weighted graph is equivalent to solving the corresponding symmetric traveling salesman problem.This very well-known NP-hard combinatorial problem is not solvable in polynomial time.
The number of all tours in such a graph is N T = (N−1)! 2 , which grows extremely fast with an increasing number of nodes.Therefore a brute force approach can only be suitable for a very low number of nodes.Although there exist many exact algorithms to solve a traveling salesman problem, like branch-and-cut or branch-and-bound, an exact solution is getting more and more impractical with an increasing number of nodes.Actually, for most applications a good approximation of the optimal solution would suffice.This can be achieved efficiently by heuristic algorithms.
A cutting job, which has to be processed by the laser cutting machine, consists of several cutting paths, denoted c i , with 1 ≤ i ≤ N c .Each of these paths c i is defined by two points in the work plane, a start point r s,i and an end point r e,i .These two points are connected by a path, for which it is assumed that a time optimal trajectory with respect to the machine and process constraints is known.Additionally a job should start and end at a defined idle position denoted r 0 .Since the trajectories for each cutting path are already assumed to be optimal, the overall processing time can be reduced by finding a processing sequence of these paths, which minimizes the total non-productive traversing time between the consecutive paths.To this end, a complete, undirected, weighted graph of the start and end points of the cutting paths is constructed with edge weights according to the traversing time between these points.The optimal solution of the problem then corresponds to the solution of the symmetric traveling salesman problem on this graph.
The cutting paths c i can be closed, i.e. r s,i = r e,i , or open, i.e. r s,i = r e,i .Without loss of generality it is assumed that the first N cc paths are closed and the remaining N co paths are This can always be obtained by reordering indices.Since the start and end point coincide, closed paths can be represented by a single node in the graph.Based on this, the N = N cc + 2N co + 1 nodes of the graph can be associated with the start (and end) points of the cutting paths as well as the idle position according to the mapping which is shown for an examplary task in Fig. 2. With the definition of the edge weights, provided in Sect.3, the graph is fully defined and the traveling salesman problem can be solved.If open cutting paths are present, i.e.N co > 0, the algorithm used must ensure that these paths are actually traversed in the solution.Additionally, if the cutting job contains nested cutting paths, the according precedence constraints have to be taken into account.

Objective function
Regardless of whether the problem is solved by an exact or by a heuristic algorithm, the objective is to find the tour with the minimum total edge weight.In order to compute the edge weights, a special metric is introduced.Since the overall goal is to achieve a time optimal solution, all the considered links connecting two points on the plane have to be time optimal on their own.Therefore the distance between two arbitrary points can be expressed by the minimum time needed to traverse a straight line between these two points satisfying the system limits.The resulting distance measure satisfies the triangle inequality as well as the remaining requirements of a metric.
The transition from one point r n (n i ) = r i = x i y i in the plane to another point r n (n j ) = r j = x j y j , corresponding to the nodes n i and n j , respectively, is performed according to a time optimal sin 2 -jerk trajectory.This trajectory has the property that it is four times continuously differentiable, which is beneficial since the excitation of vibrations is reduced.Another advantage of this trajectory is that the minimum traversing time between two points can be calculated analytically with respect to maximum velocity, acceleration and jerk.These maximum values can be defined for each axis individually as v max,k , a max,k and j max,k with k ∈ {x, y}.The maximal tangential velocity, acceleration, and jerk can be constrained by v max,t , a max,t and j max,t , respectively.In order to calculate the minimal traversing time, the distances x i,j = x j − x i and y i,j = y j − y i along each axis are needed, which can be used to calculate the Euclidean distance between the two points.With these distance measures the actual maximum values for a straight path between the points r i and r j are The trajectory over time t can then be defined as with the function ξ(t, d, v max , a max , j max ) satisfying the properties By partitioning the traversing time t E in 7 phases by interleaving 4 jerk phases of length t j with 3 jerk-free phases of lengths t a , t v , and t a corresponding to phases of constant acceleration respectively velocity, ξ(t) and its derivatives (omitting the last arguments for brevity) can be stated as ...
In order to get a time optimal trajectory, four cases have to be distinguished, which have the following solutions: - The terminal time of ξ(t, d, v max , a max , j max ) can then be stated as Finally the minimum traversing time between the points r i and r j , which defines the edge weight w i,j between the nodes n i and n j , is The objective function is then defined as for a tour according to the set of edges {e i 1 ,j 1 , e i 2 ,j 2 , . . ., e i N ,j N } forming the tour.

Nested cutting path detection
As stated in Sect. 2 a cutting path c i is defined by a start point, an end point and the actual cutting path.In this paper, it is assumed that c i can be represented as a series of connected cubic Bèzier curves c i,j (σ ) ∈ R 2 with j ∈ {1, . . ., N b,i }, the number of Bèzier curves N b,i for the i-th cutting path and the curve parameter σ ∈ [0, 1].If a path would not be directly representable as a series of cubic Bèzier curves, it could always be approximated as such, which would be sufficient for the following consideration, provided the appoximation error is small enough to be neglected.The algorithm used to detect, whether a point lies within a closed series of Bèzier curves, is based on the point in polygon test via winding numbers, shown in [8].As indicated by Fig. 3, the intersections of the individual Bèzier curves and a ray casted from the test point in positve x-direction are examined.Since cubic Bèzier curves are used, the intersection points can be derived by solving polynomials of degree 3. Every time a curve traverses the ray from below the winding number N w is incremented, otherwise the winding number is decremented.The direction of traversal can easily be identified by evaluating the derivatives of the polynomials for the y-coordinates of the curves.Thereby saddle points in the interior of a Bèzier curve can be ignored, but special care has to be taken, if the ray intersects a start or end point of a curve segment.If the winding number results to N w = 0 after all the intersections are processed, the test point is outside the closed series of Bèzier curves, otherwise, if N w = 0, the point is inside.To speed up calculation, only those Bèzier curves are tested against intersection, for which the respective bounding boxes are intersected by the ray.
Assuming that the individual cutting paths do not intersect each other, one cutting path is enclosed by another cutting path, if and only if an arbitrary point of the former is enclosed by the latter.Additionally, a fast test based on the bounding boxes of the entire paths can be performed, since a cutting path may only be enclosed by another cutting path, if the bounding box of the former is enclosed by the bounding box of the latter.If this assumption does not hold, i.e. two cutting paths intersect each other, the proposed algorithm including the bounding box test may return the false positive result that the path with the smaller bounding box is enclosed by the other.This would result in a constraint which is not strictly necessary and might slightly influence the overall solution quality, since for intersecting paths neither way of prioritizing one path against the other would be ideal in general.In Fig. 3 Winding number algorithm for a closed series of Bèzier curves order to detect such corner cases excessive mutual intersection tests would be necessary, which is omitted in favor of overall calcuation time.
By performing this test for all the cutting paths, a precedence graph can be constructed.This graph is a forest where the root of each tree represents a path not enclosed by any other path.The leaf nodes represent the inner most paths not containing any paths.Consequently every edge in this precedence graph corresponds to a precedence constraint for the sequencing algorithm, resulting in a total number of constraints N pc .A precedence constraint indicates that a cutting paths has to be sequenced before another.

Christofides algorithm
The Christofides algorithm [1] is a construction heuristics based on the minimum spanning tree of the complete, undirected, weighted graph.The main steps of the algorithm are as follows: -Build the minimum spanning tree -Select all nodes of the minimum spanning tree with an odd degree -Find a minimum-weight perfect matching of the complete subgraph of these nodes -Combine the edges from the perfect matching with the minimum spanning tree -Form an Eulerian cycle on the result of the previous step -Convert the Eulerian cycle to a Hamiltonian circle (tour) by skipping all repeating nodes The total sum of the edge weights (tour length) of the so obtained tour over all nodes of the graph is guaranteed to be within 1 and 1.5 times the length of the optimal tour if the edge weights satisfy the triangle inequality.
As shown in [9] this algorithm can be extended to obtain a near optimal solution for problems with N co > 0, i.e. open paths are present.In order to guarantee that an open path c j , j > N cc , with start and end nodes n i and n i+1 , respectively, is traversed, the weight w i,i+1 has to be modified according to This ensures that the edge corresponding to an open cutting path is selected, when building the minimum spanning tree.Additionally it has to be ensured that these edges are not skipped when converting the Eulerian cycle to the Hamilton cycle.As stated in [9], the shortcoming of this procedure is that the triangle inequality is no longer satisfied and the key feature of the Christofides algorithm, which is the upper bound of 1.5 times of the optimal tour length, does no longer apply.Since the early stages of the algorithm are only based on the minimum spanning tree and the perfect matching, sequencing is not present before building the Eulerian cycle.At this stage there is only a limited number of degrees of freedom how the sequencing has to be built.Therefore it is not generally possible, to satisfy all of the precedence constraints derived from nested cutting path structures.Although different ways of incorporating the constraints would be conceivable, like partitioning the problem according to the levels in the precedence graph, starting with the leafs and working bottom up, the solution quality strongly depends on the problem and can be arbitrary low.Thus extending the Christofides algorithm for nested cutting path structures is omitted here.
The algorithm used for solving the examples without nested cutting paths is implemented in Python based on the packages numpy and networkx.This is mainly because the networkx package contains an algorithm for finding the minimum-weight perfect matching based on [5].In order to speed up the calculation time of the minimum-weight perfect matching, a heuristic approach is chosen for problems with N ≥ 1000.This means that the matching is not performed on the complete subgraph of the odd degree nodes, but on a minimum-weight sparse subgraph thereof, with a minimum degree of the occurring nodes of five.

Lin-Kernighan-Helsgaun (LKH) algorithm
In contrast to the Christofides algorithm, which terminates when a feasible solution is found, the Lin-Kernighan-Helsgaun (LKH) algorithm [7], which is based on the algorithm of Lin-Kernighan [11], is an iteratively improving heuristics.After a preprocessing phase, essentially based on an extended minimum spanning tree, a suboptimal initial solution is generated with any suitable and fast construction heuristics.This initial solution is then improved by so called k-opt exchanges, which means that in every iteration step, k edges of the current tour are replaced by k other edges, in order to improve the tour.Special heuristic rules are applied to decide which edges should be removed and which edges should be used instead.This drastically reduces the corresponding search spaces and consequently also the calculation time.To decrease the calculation time even further, the k-opt exchanges are constructed sequentially from k = 2 to k = 5.Once an improvement is found, the exchange is applied immediately and the algorithm proceeds with the next iteration step.The algorithm terminates when no further improvement of the tour length with changes k ≤ 5 and the applied heuristic rules can be found.
This algorithm can be applied straightforwardly to a problem with only closed paths, i.e. problems with N cc > 0 and N co = 0, and without nested cutting paths present.In order to use this algorithm with open and closed paths, i.e. problems with N cc ≥ 0 and N co > 0, or nested cutting path structures some modifications in the problem setup and the algorithm are necessary.The algorithm used for solving the examples is implemented in C++ based on the standard library.

Problems with open and closed paths
In the preprocessing phase of the algorithm, constructive sets are associated to each node, which are used to decide which edge should be added in the iteration phase.With open Fig. 4 Edge selection for open paths with edges considered to be added to the current solution (dashed) and actually chosen edges (full) cutting paths present, it has to be ensured that the corresponding edges are part of these sets.Additionally when constructing the initial tour, it has to be ensured that these edges are part of the initial tour.In the iterative improving phase two modifications of the algorithm are applied, in order to ensure the validity of the solution and to reduce the calculation time.Firstly edges according to the open cutting paths are never removed from the current solution.Secondly each time an edge e i,j is considered to be added to the current solution, it is determined whether the start and end node, n i and n j , of this edge is either start or end node of an open cutting path.If so, it is checked whether flipping the sequential order of the two end points of the affected open path may further improve the current solution.Figure 4 shows the edges (dashed and fully drawn lines), which are taken into account, if any edge adjoint to one of the two open paths in the figure is considered to be added to the current solution.The fully drawn lines represent the edges, which are then actually chosen.

Problems with nested paths
Since each iteration step results in an intermediate solution tour, the ordinal numbers of the nodes in these tours can be used to check, whether all precedence constraints, associated with nested cutting paths, are satisfied after each iteration step.Furthermore the differences of the ordinal numbers give rise to a constraint violation measure by summing up the ordinal distance from an enclosed node to the node of the enclosing path, if the respective constraint is violated.This results in a measure, which is equal to zero, if all constraints are satisfied, or strictly positive otherwise.Using the total edge weights of a tour and the constraints violation measure as a combined cost function, the LKH heuristic can be extended to solve problems with nested cutting path structures.Thereby an improvement of the tour occurs either if the constraints violation decreases or if the constraint violation remains constant and the sum of the edge weights of the tour decreases.If the initial solution tour satisfies all the constraints, it is guaranteed that also the final solution satisfies all the constraints, since only improving edge exchanges are applied to the tour.Such an inital solution can always be constructed, e.g. by a nearest neighbor construction heuristic inserting nodes only after all its child nodes in the precedence graph have already been inserted.

Integer linear programming
In this section, an integer linear programming algorithm for solving problems without nested cutting paths is stated for reference purposes regarding solution quality.In order to solve the unconstrained traveling salesman problem with integer linear programming, each edge e i,j is associated with a boolean variable êi,j ∈ {0, 1}, which indicates whether the respective edge is part of the solution tour.According to the work of [2] the optimization problem can be stated by min êi,j s.t.
The first constraint (33) enforces that each node is part of the solution and has degree two, which means it has one incoming and one outgoing edge.The second constraint (34) ensures that the edges associated with the open cutting paths are part of the solution edges.Solving this problem does not necessarily result in a valid tour, since also a number of cycles containing all nodes fulfil the constraints.Therefore additional constraints, which eliminate cycles have to be added.The number of this constraints grows exponentially with an increasing number of nodes.Therefore the optimization problem is solved iteratively and in every iteration the constraints for each occurring cycle C k = {(i

Examples
The examples shown in this section are created randomly based on the Voronoi cells of randomly distributed points in the plane.Each problem consist of N cc closed and N co open cutting paths.The constraining maximal values used for calculating the edge weights are stated in Table 1.

Examples without nested cutting paths
At first the extended heuristic algorithms are compared regarding their performance on problems without nested cutting paths.Table 2 shows the total edge weight of the solution tours, excluding the weights of the edges according to the open cutting paths, acquired by the different algorithms.Since the three algorithms are implemented in different frameworks, comparing the absolute calculation time is not really meaningful.Nonetheless, the calculation times are stated in Table 3 to indicate the dependency on the problem sizes.For the problem sizes of N c ≥ 300 the integer linear programming approach did not terminate in a reasonable time span, which is noted by a "−" in the tables.Graphical representations are provided for the problems of sizes N c ≤ 100, as can be seen in the Figs. 5 and 6.The figures for the problems with a higher number of cutting paths are omitted, since the graphical representation is getting less and less expressive.As can be seen in Table 3 the Christofides heuristic and LKH heuristic perform equally well for small problem sizes, regarding solution quality as well as calculation time.But with an increasing number of nodes the LKH outperforms the Christofides heuristics in both criteria, especially regarding the calculation time.Additionally it has to be noted that for the Christofides heuristics and problems with N c ≥ 1000, the minimum-weight perfect matching has to be performed on a sparse subgraph of the nodes with odd degree, in order to get results within reasonable time.Remarkably, for all the examples for which an optimal solution has been found by the integer linear programming algorithm, the LKH heuristic leads to exactly the same solution.Furthermore, it is noteworthy that the calculation of the exact traversion time of the non-productive linking trajectories does not contribute substentially to the overall calculation time.

Examples with nested cutting paths
In this section the extended LKH algorithm is applied to problems with precedence constraints.The solution as well as the calculation time is compared to the results, derived by the LKH algorithm with deactivated constraints processing.Table 4 shows the total edge weight of the solution tours, excluding the weights of the edges according to the open cutting paths, for different problem sizes.It can be seen that the sum of the edge weights is  only slightly larger with precedence constraints active (last column) than with deactivated constraints processing (second last column).Although this does not imply optimality, it indicates that the solution quality is compareable to the unconstraint case.This can also be verified by the graphical representations for the example with N c = 20, provided by Fig. 7.
The solution, shown in Fig. 7-d, has many edges in common with the unconstraint solution, shown in Fig. 7-c, and differs only as needed to satisfy the precedence constraints, indicated by Fig. 7-b.As can be seen in Table 5, the calculation time for detecting the nested cutting path structures may be neglected.But satisfying all the constraints substantially increases the calculation time of the solution, especially for larger problem sizes and higher numbers of constraints.This is because applying k-opt exchanges involves almost always the reversal of parts of the tour.Therefore an increasing number of constraints means that more and more k-opt exchanges, which would lead to an improvement of the total tour length, have to be discarded, because constraints would be violated.

Conclusion
It has been shown that the Christofides heuristics as well as the LKH heuristics can be applied for scheduling open and closed cutting paths, while minimizing the overall traversing time of the non-productive trajectories connecting these paths.This has been achieved by amending these algorithms and by introducing a special metric for calculating the edge weights, which models the traversing time of a time optimal sin 2 -jerk trajectory along a straight path in the plane.The extended versions of the Christofides heuristics, the LKH heuristics, and an integer linear programming algorithm are compared regarding solution quality and calculation time using randomly generated examples without nested cutting paths of different problem sizes.The results of these examples indicate that both the Christofides and the LKH heuristics are capable of finding good approximations of the optimal tour for small problem sizes, and that the LKH heuristics performs better, especially regarding calculation time, for problems with a higher number of nodes.
Finally, the LKH heuristics is applied to examples with nested cutting paths of different problem sizes, resulting in good solution tours, while satisfying all the precedence constraints.Furthermore, it has been shown that the calculation time for the detection of the nested structures, using the algorithm provided in Sect.4, is neglectable.Whereas, the calculation time for finding the solution tour is substantially increased compared to finding the unconstrained solution, especially for a higher numbers of constraints.

Fig. 2
Fig.2Examplary indices and notation of cutting paths and nodes for examplary task of Fig.1

Fig. 5 Fig. 6 Fig. 7
Fig. 5 Example N c = 30, without nested cutting paths t a ≥ 0 and t v ≥ 0. The maximal values of ξ(t) and ξ(t) as well as the final value ξ(t E ) are thus max t 1 , j 1 ), (i 2 , j 2 ), . . ., (i N k , j N k )} are added with N k = |C k | according to the set of N k edges {e i 1 ,j 1 , e i 2 ,j 2 , . . ., e i N k ,j N k } forming the cycle.The algorithm used for solving the examples is implemented in Matlab based on intlinprog(. . .).

Table 2
Examples without nested cutting paths: Total edge weight of solution tour

Table 3
Examples without nested cutting paths: Approximate calculation time

Table 4
Examples with nested cutting paths: Total edge weight of solution tour

Table 5
Examples with nested cutting paths: Approximate calculation time