Generating the Isocurve Representation for Configuration Space of Mechanisms

This paper proposes an approach for generating a set of isocurves as a representation for configuration space (CS) of a mechanism. An isocurve here is a curve in CS with some parameters fixed. Compared with conventional methods like box approximations, sampling points, and boundary tessellations, the isocurve-based representation has some advantages in space parameterization and data management. This approach directly formulates the joint loop equations in the form of kinematic matrices, which does not need any extra conversions and solves the equations with an isocurve-tracing method applying ODE solvers. Since the isocurves are connected in certain orders with the guide isocurves in a lower-dimensional space, tracing all the isocurves only needs one initial solution point for an isolated solution component. In addition, the proposed approach includes an interference-handling step, which trims off the collision portions of the isocurves by checking their feasibility according to the previously defined halfspace constraints, and a measure for identifying the forward direction at singular points where the first-order derivatives vanish. The approach is implemented through programming and the results for a few examples show its effectiveness.


Introduction
The configuration of a part in mechanism is its position and orientation in a world coordinate system and the configuration variables for a single part are simply its three-position coordinates and three-orientation angles, which usually serve as its configuration representation [1].Configuration space (CS) of a mechanism is a representation for all the feasible positions and orientations of all parts in the mechanism.Taking into account the motion constraints on closed chains and the collision-free requirements between different parts, the CS is often one or more regions of certain manifolds in a high-dimensional space [2] and it reflects overall kinematic characteristics of the mechanism.Usually, the CS coordinates are the kinematic parameters, its manifold dimensionality shows the degree of parameter coupling imposed by loop constraints, and its boundaries indicate the range of feasibility from interference-free requirements.Therefore, configuration space is frequently used in mechanical design for kinematical identification, function verification, and design optimization [1,3].Operation planning and control is another important application of configuration space after the design is finished [2,4].In this case, a subset of mechanism's CS can be utilized to describe the feasible motion ranges or the reachable space of a specific point of the mechanism after considering its interference with environmental obstacles and other actuation requirements, which is known as workspace in robotics.
Generally, there exist three kinds of representation schemes for mechanism's CS; they are box approximations, sampling points, and boundary tessellations.Box approximation methods, which describe the range of CS with a number of boxes with various sizes, are mostly adopted in the design research community because the workspace volume as an important performance metric can be calculated conveniently from this representation [5,6].Sampling points methods are often used in path planning area, in which the points randomly generated within CS together with the edges connecting them serve as a representation for the feasible configurations and their transition relations [7,8].The graph form of the sampling points representation is well suitable for the purpose of optimal path search and there is no special difficulty in handling higher-dimensional problems.Boundary tessellation methods that only focus on the positions of boundary points of CS are also preferred by some researchers [9,10] because of their simplified target in generating the representation and probably more precise results.Although the above CS representations can support the functions in design and planning in some ways, they cannot meet the requirements in more sophisticated design analyses and process planning that involves frequent queries of the CS points, especially for the CSs trimmed from lowerdimensional manifolds in higher-dimensional spaces because the identification of each CS point requires solving nonlinear equations with a large number of unknowns.For example, when one is planning the motion of a complex parallel mechanism, he may need to frequently retrieve the neighboring CS points in order to identify an optimal motion step.In this case, the well-prepared CS representation can greatly reduce the computation cost spent on the CS point retrievals.
This paper proposes to represent a CS with a set of trimmed isocurves in certain orders.An isocurve here is a curve in the CS manifolds with some parameters fixed.Compared with the conventional methods mentioned above, the isocurve-based representation has some advantages in space parameterization and data management, and thus it is convenient for querying individual points in CS.However, identifying these points in large number needs efficient and precise algorithms for solving constraint equations.
There exists plenty of research work on solving constraints in mechanisms and a good brief review could be found in [11].Generally, the approaches developed in the past can be classified into two categories: algebraic methods and numerical methods.The algebraic methods include the elimination method [12,13] and the Gröebner basis method [14].The methods have advantages in finding all the solutions with high precision, but they usually suffer from the difficulties emerging in implementation that involves the equation form conversion and symbolic computation.Most of the numerical methods reported in literature are the continuation method [15][16][17] and the interval-based methods [11,[18][19][20]; the latter is more frequently used in practical cases, which generate box approximations for a CS.To avoid missing a solution, the interval-based methods check all the boxes in a CS whether they include a solution point and then minimize the boxes to achieve an appropriate precision.Although the computation cost may not be low, the calculation process is relatively simple compared to approaches involving the inverse computations.Usually, the forward computations from CS to constraint function values are carried out with the interval arithmetics [11,18].But partially using the inverse computations, such as the extreme point identification with optimization approaches [19] and the region exclusion by checking proper conditions for branch-and-bound approaches [20], is helpful for improving the efficiency.Newton-Raphson method, as a popular numerical method for solving general nonlinear equations, was also used for solving constraints in mechanisms [21][22][23].While it usually produces the solutions with higher precision compared to the interval-based methods, it can only identify the solutions point by point.In this paper, an ODE-based approach is proposed to obtain the isocurves in the solution space of mechanism constraints.
This paper is organized into seven sections.In Section 2, the concept and representation of isocurves for a CS are presented.After this, the formulas for the ODE-based approach that generates isocurves are derived in Section 3. Section 4 focuses on the method for trimming the isocurves with interference-free constraints, while Section 5 gives a second-order derivative approach for identifying the forward direction at singular points.Section 6 presents the computational results of three examples obtained with the proposed methods.Finally, the paper is ended with some conclusions and outlooks in Section 7.

Isocurve Representation for CS of Mechanism
For a given mechanism, let  = ( 1 ,  2 , . . .,   ) be the vector of motion parameters consisting of all its joint variables and assume that  is constrained by equality constraints ℎ() = (ℎ 1 (), ℎ 2 (), . . ., ℎ  ()) = 0,  < , which are mainly composed of the loop equations on all the closed kinematic chains.Then its configuration space without considering the collision-free requirement between parts in the mechanism is Here,  is the real number set and   is the p-dimensional space.
Partition  is  = (, ), where  ∈   is the vector for active parameters while  ∈  − is the vector composed of passive parameters.Here, active parameters are the motion parameters whose values can be changed independently while passive parameters are those whose values are determined by solving the equality constraints for given active parameters.Then, the domain of active parameters  can be defined as follows: In addition, use () = { |  = (, ) ∈ } to denote the set of -points in  that are associated with a -point  in .For a given , () consists of all the -points, made of  = (, ), which satisfy the equality constraints.Generally, () includes m points for a mechanism with m solution components.For a fixed point  0 = ( 0 1 ,  0 2 , . . .,  0  ) ∈ , define the set   = {( 1 , . . .,   ,  0 +1 , . . .,  0  ) ∈ } (1 ≤  ≤ ) extended from  0 by making its first k coordinates varying and the CS 's subset   = ⋃ ∈  ().An isocurve of  is a curve which has all the active parameters fixed except one of them.Let   be the variable active parameter and then the   -V at  0 can be expressed with Here,  0 defines the starting point of an isocurve while   gives its variable.It is worth noting that the isocurve variable   may not vary monotonically on the curve like the parameter of an ordinary parametric curve.Furthermore, we define   = {(;  +1 ) |  ∈   } to be an isocurve cluster, where   is called the guide of the  cluster and it specifies all the starting points of the isocurves in   .Actually,  −1 forms an isocurve representation of configuration space  because  = ⋃ Therefore,   , that is, , is represented with the   -V cluster  −1 at guide  −1 .Recursively,   is represented with the   -V cluster  −1 at guide  −1 , and, finally,  1 is represented with a single  1 -V at guide  0 = { 0 }.Summarily, they are According to the above formulation, starting with a given feasible point  0 ∈  or a point in ( 0 ) ⊂ , a hierarchical isocurve representation can be created in the order implied in (5).The representation mainly has three characteristics (see Figure 1): (1) it is parametric because the points in  can be arranged in the hierarchical order; (2) the elemental entities of the representation are curves in   ; (3) the curves are grouped and connected in a specific way, which facilitates the generation and query of them.

ODE-Based Approach for Isocurve Generation
Suppose that ℎ() = 0 represents the equality constraints for the mechanism under investigation.In detail, the equality constraints ℎ() = 0 are further divided into loop constraints   () =  ( = 1, 2, . . ., ) and selection constraints   () = 0 ( = 1, 2, . . ., ).Here, m is the number of independent joint circuits in the mechanism,  is the identity matrix, and For an isocurve, there is a parametric expression  = () for the solution of ℎ() = 0 when appropriate selection constraints   () = 0 are imposed on some parameters.From ℎ(()) ≡ 0 for variable t, there are ∇  ⋅ Ẋ = 0 ( = 1, 2, . . ., ) and ∇  ⋅ Ẋ = 0 ( = 1, 2, . . ., ), where ∇ ⋅ Ẋ = (∇  ⋅ Ẋ) 4×4 for  = (  ) 4×4 and Ẋ = /.The former is equivalent to  −1  (∇  ⋅ Ẋ) = 0, which can be expressed in the form of screws: The purpose of using the screw form here is to reduce the number of equations from 12 to 6; this is true for all the individual loop constraints of a mechanism because the constraint equations in a homogeneous 4 × 4 matrix are converted into the 6×1 vector form for screws.The reduction of constraints is mainly due to the fact that the derivatives   /  of the 12 constraints in matrix equation   () =  only have 6 independent elements when they are expressed in the frame of   as  −1  (  /  ).Assembling the equations for Ẋ, we have Here J is a (6 + ) ×  matrix.For an appropriately defined isocurve that has only one DOF for expressing the point variation in the curve, the rank of J is  − 1 and its onedimensional null space can be obtained with the singular value decomposition.Usually, an appropriate isocurve is defined by specifying the right selection constraints   () = 0. Suppose that the unit basis vector of its null space is d, which satisfies  = 0. Then the Ordinary Differential Equations (ODE) for the isocurve is Here  is a sign chosen to make the right side () change continuously along the isocurve.Since the solutions obtained from a linear system solver for  = 0 are two vectors ±, the choice of the solution's sign is needed in order to ensure that the curve tracing is conducted in a consistent direction.In this research,  is chosen at each step of the numerical calculation so that the inner product of the current () with the one at the previous step is a positive value.

Interference-Free Constraint between Point and Convex Body.
Construction of interference-free constraints is a complicated work and so we simplify the problem with some approximations (see Figure 2).Here, the interferencechecking is made for a pair of parts   and   with relative motion  , , which could be an obstacle to avoid in space or a link of the mechanism.First, both of the parts are represented with a set of implicit functions   = { ,1 (),  ,2 (), . ..} that bound a region in the way   = { | () ≤ 0, () ∈   } ( = , ), where  is a point coordinate defined in the part's local coordinate frame.Second, both of the parts are approximated with a set of scattered points   = {V ,1 , V ,2 , . ..} ( = , ) where V , are the points' coordinates in   's local frame.Then, the interference-free constraints between   and   are defined as follows: where  , = −min Here,  , is the transformation from   's frame to   's and  , V gives the coordinates in the   's frame for point V.In addition, the coordinate transformation  , =  , () is composed of the joint matrices in the shortest joint chain connecting   and   .
For the example in Figure 2, constraint  , ≤ 0 is constructed from the implicit functions  , ( = 1, 2, . ..) of   , the boundary points V , ( = 1, 2, . ..) of   , and the motion matrix  , transforming   's coordinates to   's.To apply the above constraints to some objects like nonconvex parts, which is difficult to be expressed with the intersections of half-spaces defined with a set of implicit functions, they should be decomposed into multiple convex bodies first and then the inequality constraints are created for each pair of convex bodies, respectively, belonging to the two parts under consideration.The one overall inequality constraint can be expressed as follows: While the above proposed approach is generic in detecting collisions of parts in a mechanism, it suffers from some restrictions.For example, all the part surfaces need be expressed with implicit functions, concave parts should be decomposed into multiple convex bodies, and the checking points on part boundaries should be selected beforehand.

Interference-Free Constraint between Line Segment and Convex Body.
Here, let us consider the interference-free constraint between a line segment  = V 1 V 2 and a convex body  defined by  = { 1 (),  2 (), . ..}. Suppose that the points on  are expressed as and the constraint is formed according to (10) as follows: In the above formula, and  is the transformation from the 's local frame to that of the convex body.Since the part faces in industry are often planar or quadratic surfaces,   () are considered as linear or quadratic functions here, which can be expressed in the form   () =   +    +    2 , 0 ≤  ≤ 1.Therefore, the evaluation of the constraint function in (12) mainly includes the operations of identifying the uppermost curve segments among different   ().As shown in Figure 3, the upper-most curve segments in orange can be determined after calculating the intersection points (marked singular point with circles) between different   (), which is trivial for linear or quadratic functions.
Taking into account the interference-free requirement and motion parameter ranges, the end points of the open feasible segments of the curve should be a solution of (()) = 0, which is defined by (11).Therefore, the trimming points can be obtained by solving (()) = 0. Usually, the equation can be solved with Newton-Raphson method, but since it is hard to calculate the derivative of (()), which is required by the method, we adopt the bisection method here.Assume that two consecutive points  in and  out satisfying ( in ) < 0 and ( out ) > 0 on the discrete expression of the isocurve  = () are already detected.Then we search for  * ∈ (0, 1) that makes ((1 −  * ) in +  *  out ) = 0, by iteratively subdividing the interval into two at the midpoint.

Handling the Singularity Issue
In the ODE approach for generating isocurves presented in Section 3, the dimension of 's null space is supposed to be 1.But at singular points of a mechanism, this may not be true.One case is that additional infinitesimal motions in certain directions exist at a singular point, which makes the rank  of  less than  − 1.Another case is kinematotropic mechanisms [24] which have different finite motions in other directions at a singular point.So, the problem is to determine the forward direction at a point where  is less than  − 1 in solving (9).When  is less than  − 1, its null space has two or more basis vectors for expressing the solution of () = 0 in the way The forward direction   at   should satisfy the following conditions (see Figure 4): (1)   should belong to the null space; that is, , where V , are its basis vectors; (2)   should be continuous with its previous forward direction on the curve, which requires 1 −  < (  )   −1 < 1 +  for a small  > 0; (3)   should represent a direction for finite motion that satisfies constraints ℎ(()) = 0.
In the above three conditions, the first two are easy to handle, but the last one is not.Here, a second-order derivative method is proposed to identify the finite motion direction in the linear combinations of two or more basis vectors.For simplicity, only one loop constraint () −  = 0 in ℎ() = 0 is considered here.From (())− ≡ 0, there is ∑  =1 $  () ẋ  () ≡ 0, where $  = screw( −1 (/  )).Furthermore, performing the second-order differential operation on the identity, the following can be obtained: where  , =  −1  2      . ( where ] and the symmetric matrix at the right side is Let   R  = (  ) 3×3 ; there are Since the right side of ( 15) is completely determined by the first-order derivative, the number of independent variables in the 's second-order derivative is 3, and the total independent-variable number for the 's secondorder derivative is 6 after considering those in   ö  .In the way similar to the definition for the first-order screw,  , can be expressed with a skew symmetric matrix formed using the 6 independent second-order variables Φ = ( 1 ,  2 ,  3 ;  1 ,  2 ,  3 ) as follows: where   =       .In the above formulation, the second matrix is calculated from the first-order derivative while the first one is related to the second-order derivative in the way: After  , are computed with (( 16)-( 17)), the matrix  = (  ) 4×4 in ( 13) can be obtained for a given linear combination  = ( ẋ of the null space basis vectors.When   = 0 for all 1 ≤ ,  ≤ 4, then  is considered a finite motion direction in this research.

Examples
In this section, three examples are presented to show how the proposed approach can be used to generate the isocurve representation of configuration space of a mechanism, among which the first and third examples are planar mechanisms while the second one is a spatial mechanism.

Four-Stop Geneva Mechanism.
The first example is a four-stop Geneva mechanism, which is shown in Figure 5.The mechanism is formed with three parts  0 ,  1 , and  2 , and it has three different motion phases, respectively, illustrated in Figures 5(b)-5(d).In the figure, symbols  1 ,  2 , , , and  represent the five motion parameters of the mechanism;  1 and  2 , respectively, denote the rotational motions of  1 and  2 with respect to  0 , , and  which are the rotational angles from  1 to  2 , respectively, for the first and second phases, and  is the distance between the pin of  2 and the center of  1 at the first phase.The dimensions of the mechanism are  = 100,  1 = 80, ℎ = 65,  = 5,  = √ 2,  = cos −1 ((3 2 −  2 1 )/2 √ 2 2 ),  = sin −1 ( sin / 1 ) + /4, and  =  − /4.
In the first phase, the mechanism has one loop whose matrix equation is  01 ( 1 ) 12 (, ) −1 02 ( 2 ) −  = 0.The equation is solved with the ODE approach presented in Section 3, and its solution is presented in Figure 6.In the figure, all the three horizontal axes are the sequence numbers of points on the isocurve while the vertical axes are the values of its individual coordinate components.The figure shows that the approach can generate the isocurve with a high accuracy; the matrix equation errors in the Frobenius norm are less than 7 × 10 −6 and the larger errors occur at points where the curves have steeper slope.For example, the curve for parameter  changes rapidly near the 25th and 250th points and then the errors there fluctuate severely.The second phase shown in Figure 5(c) has a smooth motion and its precision reaches as high as 6 × 10 −13 .
The third phase shown in Figure 5(d) has no loop constraints but it needs to handle the interference-free constraints.Here, the methods presented in Section 4 are used to generate its configuration space and the result is displayed in Figure 7.In this example, the point set  1 for  1 includes the 12 corner points on its outer boundary and its surface set  1 has 9 half-space functions; they are the four  1 arcs, the four small  circles at outlets of the 4 slots and one big circle .Part  2 is represented with two separate bodies for interferencechecking; one is the C-shape area and the other is the circular pin with radius .So, it has two surface sets  2 =  2,1 ∪  2,2 and two point sets  2 =  2,1 ∪  2,2 whose points are sampled on their respective boundaries.In Figure 7, the interferencefree areas are represented with straight-line isocurves and the end points of the lines exactly give the locations where  1 and  2 contact with each other.
The overall configuration space of the mechanism is a union of multiple domains in the 5-dimensional space of parameters  1 ,  2 , , , and .Figures 8(a)-8(c) show the isocurve CS representation of phases 1 and 2, respectively, in the subspaces of ( 2 , , ), (,  1 , ), and (,  2 , ), where the red curves are for the CS of phase 1, the black are for phase 2, and the dashed lines show the CS domain adjacency relations.The two points connected with a dashed line represent the same position and orientation of the mechanism, belonging to different phases as a transitional state or different periods in cyclic motions separated by this point.In addition, Figures 8(d) and 8(e) give the CS for phases 1 and 3, respectively, in subspaces ( 1 ,  2 , ) and ( 1 ,  2 , ), where the blue lines are for phase 3. The final subfigure in Figure 8 presents the CS for phases 2 and 3 in ( 1 ,  2 , ), and it can be seen that the projection of the CS of phase 2 on subspace ( 1 ,  2 ) is completely included in the CS of phase 3, which means that phase 2 can be looked as constrained states of phase 3 and represented with one more parameter.Finally, the CS for all the three phases can be represented with isocurves in the same space; Figure 9 gives this representation in ( 1 ,  2 , ), which shows the complete feasible configurations and their transitional relations of the mechanism.
The mechanism has only one loop and its loop constraint  01  12  23  −1 03 −  = 0 is easy to construct after the local coordinate frames for each parts are determined.The difficulty comes from the change of the number for active parameters.
Starting with a given initial point where the Jacobian matrix of the equations has two zero-singular values, the loop equations can be solved with the ODE method developed in Section 3 after one selection constraint () =  01 =  0 01 is added, where  0 01 is the fixed value for one isocurve, and its isocurve representation of the solution can be obtained as the plane of  12 = 0 shown in Figure 10(b).However, if the singular values on the plane are checked exactly, it can be found that the number of zero-singular values at configuration ( 01 ,  01 ,  12 ) = (0, 100, 0) becomes three.This corresponds to the state shown in Figure 10(a), in which the mechanism has mobility with parameters  01 ,  01 , and  12 .Figure 11 shows the distributions of nonzero-singular values on the CS plane  01 - 01 .Since the parameter number is 6, the number of nonzero singular values should be 4 for two active parameters, but it can be seen in Figure 11(d) that the fourthsingular value becomes zero at point   , which indicates that the mechanism may have another DOF there.
Similarly, if we start with an initial point where the Jacobian matrix of the equations has one zero-singular value, a single line CS can be obtained as the vertical line at ( 01 ,  01 ) = (0, 100) shown in Figure 10 (see Figure 12).Please note that the green curve for the second-singular value does not reach zero at  12 = 0 though it has a small value there.So, the singular value distributions indicate that the mechanism may have two other DOFs at  12 = 0.
Although the singular value decomposition (SVD) of Jacobian matrix can identify the singular configuration and its number of additional DOFs for infinitesimal motions, it cannot give the directions of the additional finite motions.For example, when the singular point   is found during generating the isocurves on the plane  01 - 01 , we need to determine the additional finite motion direction; in fact, it is (0, 0, 1) in ( 01 ,  01 ,  12 ).Here, we use the secondderivative method presented in Section 5 for this purpose.For matrix  = (  ) 4×4 in (13) for  =  01  12  23  −1 03 , the value distribution of each element   on the surface of unit sphere at   is evaluated (see Figure 13), where the directions on the sphere are parameterized with  and  (see Figures 14(a) and 14(b)).Actually,   = (0, 100, 0, 0, 100, 0) for point   and the three basis vectors of the null space of the Jacobian matrix at   can be obtained as   ( = 1, 2, 3).Then, the vector on the sphere has the form (, ) = sin () cos () 1 + cos () 2 + sin () sin () 3 .
In order to show the direction  with () = 0, we use the following measure: directions of (, ) = (0.5, 0.5) and (, ) = (1.5,0.5).So, they represent the forward and backward directions for the additional finite motion DOF, respectively.

3-RPR Parallel Mechanism.
The third example is a 3-RPR parallel mechanism whose notations of links and joints are shown in Figure 15.The side length of the equilateral triangular platform is 100, the three slideway parts   ( = 1, 2, 3) have a length of 60, and the three legs have the maximum length 600 and the minimum length 60.The mechanism includes 8 link parts   ( = 0, 1, 2, . . ., 7), 9 joints with motion parameters   ,   , and   ( = 1, 2, 3), and two independent loops; here loops  0  1  4  7  6  3  0 and  0  2  5  7  6  3  0 are chosen as the independent loops to construct loop constraints.In addition to the 9 parameters, we also use three nonjoint parameters   ,   , and , which describe the location of platform  7 relative to the frame of the fixed part  0 .Based on the above notations, the parameter vector is defined as  = ( 1 ,  2 , . . .,   ) = ( 3 ,  2 ,  1 ,  1 ,  2 ,  3 ,  1 ,  2 ,  3 ,   ,   , ),  = 12, and the two loop constraints are where The equality constraints also include three additional equations that define the nonjoint parameters (  ,   , ): where  07 =  01  14  −1 74 ,  10 =   ,  11 =   , and  07 (, ) is the (, )-element of matrix  07 .The interference-free constraints include that the three legs should not intersect with each other, the triangular platform  7 is only allowed to move within the rectangular area of  0 , and the legs should not interfere with the platform.In order to create these constraints, the methods presented in Section 4 for point-body and line-body interference-free requirements are utilized here.
The results of the proposed approach include the detailed information about all the parameter values, and thus projections of the isocurves on various subspaces  Many spatial parallel mechanisms have a similar topology to the 3-RPR except they usually have more DOFs and the interferences should be checked in 3D space.For example, the spatial parallel mechanisms 3P-S-S/S and 3S-P-S/S, which are studied in [25], are also actuated via three prismatic joints and have three leg chains as well to support a platform.Actually, the isocurve generations for the mechanisms each are taken into account, they totally have 6 more DOFs.
Besides the DOF difference, the interference checking of the spatial mechanisms would take much more time for the evaluations of (( 11)-( 12)).Nevertheless, the constructions of the equations may not have special difficulties because the 3D shapes of their link parts are relatively simple, mainly involving cylinders, spheres, and cubes.Once the isocurve representations for their CSs are generated, the mappings between their three prismatic actuation parameters and other

Conclusions
Configuration space is a vital tool for reachability checking and path planning when its points can be retrieved conveniently.However, conventional CS representations like box approximations, sampling points, and boundary tessellations do not organize the points in a structured format and hence they have difficulty supporting frequent queries of the CS points.Especially, when CS is a trimmed region on a lower-dimensional manifold of high-dimensional space, the cost of the point computation is high but the conventional methods fail to provide accurate or detailed information of the points.To address this issue, this paper proposes an isocurve representation for configuration space and its generation algorithm.
In this research, the configuration points in a CS are represented with a cluster of curves in the motion parameter space of the mechanism.To generate the representation, three novel techniques are developed.First, the isocurve-tracing method for the homogeneous transformation matrix equations is presented on the basis of common ODE solvers.Second, a half-space approach is proposed to define interferencefree constraints.Third, a measure with the second-order derivatives of the matrix equations is given to identify the forward direction of isocurve at singular points.
The advantages of the isocurve method are mainly reflected in two aspects: (1) the accurate CS points on manifolds can be generated efficiently with the isocurve-tracing method, which not only identifies the points in one isocurve with low computation cost but also traces to the points of other isocurves in the same connected solution component; (2) the points in CS are parameterized by isocurves and the isocurves are further organized in the hierarchical order from their lower-dimensional guide spaces to higher ones, which facilitates the retrieval of CS points.
Although the proposed method is novel and has advantages in some aspects, it still needs to be improved in certain ways.First, we have not elaborated upon the issue of how to choose the stepsizes or distances between the isocurves for a better approximation to the whole CS.Second, the calculated isocurve endpoints can only give somewhat rough representation for CS boundary.Third, the issue of initial point generation, especially for multiple disconnected solution components, is not addressed here.All of these issues could be addressed in the future research along this direction.

Figure 2 :
Figure 2: Entities used for constructing the interference-free constraint  , .

Figure 3 :
Figure 3: Constraint function for line-body interference-free requirement.

Figure 4 :
Figure 4: Forward direction determination at a singular point.

Figure 6 :
Figure 6: The isocurve representation for the CS of the first phase.

2 Figure 7 :
Figure 7: The isocurve representation for the CS of the third phase.

Figure 8 :
Figure 8: The CS for two phases and their adjacency relations.

Figure 14 (Figure 9 :
Figure 14(c) shows the contours of () while Figure 14(d) gives its surface graph.It can be seen from the figure that () = 0 for the lines  =  and  =  where  ∈ Z.In addition, the figure also clearly indicates that () = 0 at the

Figure 13 :
Figure 13: The distributions of the nonzero second-order derivatives in different directions at the singular point.

MathematicalFigure 14 : 3 Figure 15 :Figure 16 :
Figure 14: The directions in which all the second-order derivatives are zero at the singular point.