Globally Optimal Joint Search of Topology and Trajectory for Planar Linkages
GGlobally Optimal Joint Search of Topology andTrajectory for Planar Linkages
Zherong Pan , Min Liu , , Xifeng Gao , Kai Xu , and Dinesh Manocha Department of Computer Science, University of North Carolina, North CarolinaNC 27514, USA, [email protected] School of Computer, National University of Defense Technology, Hunan HN410073, China, [email protected], [email protected] Department of Computer Science, Florida State University, Florida FL 32306, USA, [email protected] Department of Computer Science and Electrical & Computer Engineering,University of Maryland at College Park, Maryland MD 20742, USA, [email protected]
Abstract.
We present a method to find globally optimal topology andtrajectory jointly for planar linkages. Planar linkage structures can gen-erate complex end-effector trajectories using only a single rotational ac-tuator, which is very useful in building low-cost robots. We address theproblem of searching for the optimal topology and geometry of thesestructures. However, since topology changes are non-smooth and non-differentiable, conventional gradient-based searches cannot be used. Weformulate this problem as a mixed-integer convex programming (MICP)problem, for which a global optimum can be found using the branch-and-bound (BB) algorithm. Compared to existing methods, our experimentsshow that the proposed approach finds complex linkage structures moreefficiently and generates end-effector trajectories more accurately.
Keywords: mixed integer optimization, topology optimization, trajec-tory optimization
A planar linkage is a mechanical structure built with a set of rigid bodiesconnected by hinge joints. This structure typically has one effective degree-of-freedom actuated by a rotational motor. Since they impose a minimal burden oncontroller design, these structures are widely used as building blocks for low-costtoys and robots, as illustrated in Figure 1. By combining a series of hinge joints,the end-effector of the planar linkage will trace out a complex curve that canfulfill various requirements of different types of locomotion, including walkingand swimming [11,22].A challenging problem in mechanics design is to find the linkage structurewith an end-effector that will trace out a given curve. This problem is chal-lenging in that it searches over three coupled variables: topology, geometry, and a r X i v : . [ c s . R O ] M a y Pan et al. trajectory. The linkage topology determines which rigid bodies are connectedand the order of their connections. Clearly, the topology is a non-smooth andnon-differentiable decision variable. The linkage geometry determines the shapeof each rigid link. Finally, the trajectory determines the pose of the linkagestructure at each time instance. The last two variables are smooth and differ-entiable, but directly optimizing them induces non-convex functions. Previousworks [9,29] have proposed various solutions to address problems of this kind.These methods rely on random searches, such as A ∗ [9] and covariance matrixadaptation [29], to try different topologies. Then, for each topology, they performnon-linear programming (NLP) under the given topology to determine the ge-ometry and trajectory. However, these methods are computationally expensivebecause a huge number of samples are needed for the random search to con-verge. Moreover, even after determining the topology, these methods can findonly sub-optimal solutions due to the non-convex nature of NLP.Fig. 1: An example of planar linkages,used in a strandbeest robot for 2D walk-ing. See [18] for more details. Main Results:
Given the inputof a target trajectory, we present anew method that can efficiently com-pute a planar linkage structure withglobally optimal topology and geom-etry configurations and an accuratetrajectory reproduction. Based on re-cent advances in mixed-integer model-ing [24,6,23], we relax this joint searchproblem as an MICP problem, theglobal optimum of which is arbitrar-ily close to the global optimum of theoriginal problem. The main benefit ofMICP relaxation is that the search canbe accomplished efficiently using the BB algorithm. BB is more strategic thanrandom search, as used by [29], because it cuts impossible or sub-optimal searchspaces at an early stage, leading to higher efficiency. We have compared MICPwith prior methods using different examples. The results show that our proposedMICP approach finds solutions more efficiently and that the resulting structurematches the target trajectory more closely.In the rest of the paper, we first review related work in Section 2 and thenformulate our joint search problem in Section 3. The MICP model and vari-ous constraints required for the integrity of the planar linkage are presented inSection 4. Results and the evaluation of the proposed approach are given inSection 5.
In this section, we review related work in robot design optimization, mixed-integer modeling, and topology optimization. oint Search of Topology and Trajectory for Planar Linkages 3
Robot Design Optimization:
Robot design optimization is a superset ofconventional topology and truss optimization [16] where the decision variablesare only topology or geometry. This is because the specification of a robot designis given as a movement pattern [10], leading to a joint search in the space-time domain. The joint search problem greatly expands the search space. As aresult, many prior methods do not work since they only optimize a subset ofdecision variables [10,22,1,19,21]. Recent works [29,9,20] search for all variablessimultaneously. However, these methods are based on random search techniques,which usually require a large amount of trial and error and find sub-optimalsolutions.
Mixed-Integer Modeling:
The main benefit of mixed-integer modeling isthe use of the well-studied BB algorithm [14]. BB allows us to find the global op-timum of non-convex programming problems, while only visiting a small fractionof the search space. Mixed-integer models have been applied to a large varietyof problems including motion planning [7], inverse kinematics [6], network flows[5], and mesh generations [3]. By applying the big-M method [23], McCormickenvelopes and piecewise approximations [15], and general non-convex problemscan be easily relaxed as MICP problems. Prior works [12,17] have also formu-lated topology optimization problems as MICP. However, our work is the firstto formulate the planar linkage problem as MICP and we employ MICP to con-currently find the optimal topology, geometry, and trajectory of a linkage.
Topology Optimization:
Topology optimization of a continuum is a well-studied problem [16]. An efficient algorithm can smoothen the problem anduse gradient-based method to search for locally optimal structures over a searchspace of millions of dimensions. This technique has been widely used in the designof soft robots [27,26,28]. However, the optimization of articulated robots is morechallenging because the optimized structure must satisfy the joint constraints,making the decision variable non-smooth. Existing techniques use mixed-integer[12,17] or random search techniques [29,20] to optimize over these decision vari-ables.
In this section, we introduce the problem of joint searches for planar linkages.Our problem is to search for a structure, as illustrated in Figure 2a, where wehave a set of rod-like rigid bodies connected with each other using hinge joints.As a result, the end points of these rigid bodies can take at most N distinctpositions, denoted as node set: n , ··· ,N . Of these nodes, n is the rotationalmotor and n N is the end-effector. Within one limit cycle, n follows a circularcurve centered at (cid:0) X C , Y C (cid:1) with a radius R : n ( t ) = (cid:0) sin ( t ) R + X C , cos ( t ) R + Y C (cid:1) , (1)which induces trajectories of other nodes n i ( t ) via forward kinematics. The other N − Pan et al. n n n n n n n n → n → n → n → n → Fig. 2: (a): The Jansen’s mechanics used in Figure 1 is a planar linkage structureinvolving 7 nodes. The motor node n is green, the fixed node n is red, themovable nodes n , , , are black, and the end-effector node n is blue. Our goalis to find the topology and geometry of the linkage such that the end-effectorcurve matches the blue target curve. (b): Our MICP formulation is based onthe prior symbolic representation [13,1]. This representation assumes that eachnode is connected to exactly two other nodes with lower indices: n → , n → , n → , n → , n → .may exist between each pair of nodes n i,j , in which case (cid:107) n i ( t ) − n j ( t ) (cid:107) is aconstant.Given these definitions, the input to our problem is a target end-effectortrajectory n ∗ N ( t ). The output of our method is the following set of variablesdefining both the topology and geometry of a planar linkage: – An integer vector of size N (the number of nodes), which containing thetype of each node: fixed or movable. – An N × N symmetric binary matrix C N × N where C ij = 1 means a rigidbody connects n i,j . – The position of n , ··· ,N ( t ) at a certain, arbitrary time instance t .The goal of our method is to find the globally optimal set of variables thatminimizes (cid:82) (cid:107) n N ( t ) − n ∗ N ( t ) (cid:107) dt . In this section, we present a set of linear constraints and quadratic objectivefunctions for relaxing the joint search as an MICP problem. We first introducethe set of topology constraints to ensure the well-posed nature of the structurein Section 4.1 and then present constraints and objective functions for geometriccorrectness in Section 4.2. oint Search of Topology and Trajectory for Planar Linkages 5
As illustrated in Figure 2b, our method is based on the symbolic representationpresented in [13,22], which assumes that each movable node is attached to twoother nodes. These nodes can be of any type but must have lower node indices.As a result, forward kinematics can be processed sequentially even on linkagestructures with closed loops.Since the number of nodes is unknown, we assume that the maximal numberof nodes is
K > N . For each node other than the first motor node n , we needa binary variable U i such that U i = 1 indicates n i is used as a part of theplanar linkage structure. In addition, we need another binary variable F i suchthat F i = 1 indicates n i is fixed and F i = 0 indicates n i is movable. These twosets of variables are under the constraint that only a used node can be movable.In addition, we assume that the last node n K is the end-effector that must beused. In summary, we introduce the following sets of variables and node-stateconstraints: ∃ U i , F i ∈ { , } ∀ i = 1 , · · · , K − F i ≤ U i U = U K = 1 F = 0 . (2)Our next set of constraints ensures local topology correctness. It ensures thateach movable node is connected to exactly two other nodes with lower indices.As a result, the movable node and the two other nodes will form a triangle andthe position of the movable node can then be determined via the Law of Cosine[10]. We introduce auxiliary variables C ji to indicate whether n j is the first nodeto which n i is connected. C ji indicates whether n j is the second node to which n i is connected. In addition, we introduce two verbose variables C , i = 1 toindicate that n i is connected to nothing. The resulting constraint set is: ∃ C ji , C ji , C ji ∈ { , } ∀ j, i = 1 , · · · , K ∧ j < iC ji = C ji + C ji C ji ≤ U j ∧ C ji ≤ U ji − (cid:88) j =1 C ji = 2 − F i ∀ i = 2 , · · · , K ∃ C d i ∈ { , } ∀ d = 1 , i − (cid:88) j =0 C dji = 1 . When n i is fixed in the above formulation, then F i = 1 in Equation 3 and all C ji are zero except for C , i = 1 due to the sum-to-one constraints. If n i ismovable, then F i = 0 and C ji sums to two. As a result, there must be j , j < i such that C j i = 1 and C j i = 1. Note that j and j must be different because Pan et al. otherwise the constraint that C ji ∈ [0 ,
1] will be violated. In addition, since thefirst node n is the motor node, it is excluded from these connectivity constraints.However, this naive formulation will require binary variables for each pair of n j and n i , which requires O ( K ) binary variables all together. Instead, we adoptthe idea of special ordered set of type 1 ( SOS ) [25] and model these constraintsusing O ( K (cid:100) log K (cid:101) ) binary variables. Intuitively, SOS constrains that only onevariable in a set can take a non-zero value and it can be achieved by using alogarithm number of binary variables. The improved constraint set is: ∃ C ji , C ji , C ji ∈ [0 , ∀ j, i = 1 , · · · , K ∧ j < iC ji = C ji + C ji C ji ≤ U j ∧ C ji ≤ U ji − (cid:88) j =1 C ji = 2 − F i ∀ i = 2 , · · · , K ∃ C d i ∈ [0 , ∀ d = 1 , { C dji | j = 0 , · · · , i − } ∈ SOS i − (cid:88) j =0 C dji = 1 . (3)Finally, we introduce a third set of constraints to ensure global topologycorrectness. This set of constraints ensures that the linkage structure containsno wasted structures. In other words, each node must have some influence on thetrajectory of the end-effector node and the first motor node must be connectedto others. We model these constraints using the MICP formulation of networkflows [5]. Specifically, each node n i will generate an outward flux that equals to U i , and we assume that there is a flow edge defined between each pair of nodeswith capacity Q ji . We require inward-outward flux balance for each node exceptfor the end-effector node: ∃ Q ji ∈ [0 , ∞ ] ∀ j, i = 1 , · · · , K ∧ j < iQ ji ≤ C ji KU i + i − (cid:88) j =1 Q ji = K (cid:88) k = i +1 Q ik ∀ i = 1 , · · · , K − , (4)where we adopt the big-M method [23] in the second constraint to ensure thatonly edges between connected nodes can have a capacity up to K . Using a similaridea, we also formulate a constraint that a movable node must be connected to atleast one other movable node. We assume that each node n i generates a reversedoutward flux that equals to 1 − F i , and we assume that there is a flow edge definedbetween each pair of nodes with capacity R ji . We require inward-outward flux oint Search of Topology and Trajectory for Planar Linkages 7 balance for each node except for the motor node: ∃ R ji ∈ [0 , ∞ ] ∀ j, i = 1 , · · · , K ∧ j < iR ji ≤ C ji K ∧ R ji ≤ (1 − F j ) K i − (cid:88) j =1 R ji = 1 − F i + K (cid:88) k = i +1 R ik ∀ i = 2 , · · · , K, (5)These three constraints ensure that the planar linkage structure is symbolicallycorrect, independent of the concrete geometric shape. The main utility of geometric correctness constraints is to compute the exactpositions n i = (cid:0) x i , y i (cid:1) of each node in the 2D workspace. These positions arefunctions of time t and we sample a set of T discrete time instances t , ··· ,T . Inthis section, we will always use superscripts for timestep indices. For example, attime instance t d , the position of n i is n di . We want to find a common geometricspecification such that all the end-effector positions n , ··· ,TK can be achieved.The most important geometric variable is the length of each rigid rod. Wedefine these parameters implicitly using a set of constraints such that, if n i and n j are connected, then the distance between these two nodes is a constant forall time instances. In other words, we need the following set of constraints if C ji = 1: (cid:107) n dj − n di (cid:107) = (cid:107) n ( d mod T )+1 j − n ( d mod T )+1 i (cid:107) ∀ ≤ d ≤ T, (6)after which any distance (cid:107) n dj − n di (cid:107) can be used as the rigid rod length.However, there are two challenging issues in modeling these constraints thatcan affect the performance of the MICP solver. A first challenge is to mini-mize the use of binary variables. Because any pair of nodes n j and n i might beconnected, a naive formulation will require a number of binary variables pro-portional to K . Instead, we introduce auxiliary term d d i = (cid:0) dx d i , dy d i (cid:1) , whichindicates the relative position between n i and the first other node connected toit at time instance t d . Similarly, d d i = (cid:0) dx d i , dy d i (cid:1) indicates the relative positionbetween n i and the second other node connected to it. These definitions inducethe following big-M constraints: ∃{ dx, dy } dki ∀ k = 1 , ∧ i = 2 , · · · , K ∧ d = 1 , · · · , T |{ dx, dy } dki − { x, y } dj + { x, y } di | ≤ B (1 − C kji ) ∀ j = 1 , · · · , i − , (7)where B is the big-M parameter, implying that all the node positions lie in abounded region [ − B, B ] . Note that the first motor node n follows a circularcurve (Equation 1), which requires special definitions of d d , d d as follows: { dx, dy } d = { dx, dy } d = { x d − X C , y d − Y C } , (8) Pan et al. (cid:18) α ˜ α (cid:19) = S (cid:88) s =1 λ s (cid:18) α s α s (cid:19) { λ , ··· ,S } ∈ SOS S (cid:88) s =1 λ s = 1 , (9) α α α α α α Fig. 3: An illustration of the piecewise linear upper bound (blue) of the quadraticcurve α (red) with S = 5.where the center of rotation (cid:0) X C , Y C (cid:1) is used as an additional auxiliary vari-able. The second challenge is that these constraints are non-convex because theyinvolve quadratic terms. Fortunately, efficient formulations have been developedto relax non-convex functions using piecewise linear approximation [15] and aspecial ordered set of type 2 ( SOS ) [25]. SOS effects a constraint that at mosttwo of the variables in an ordered set with consecutive indices can take non-zerovalues. To use these formulations, we decompose the range [ − B, B ] evenly into S − S nodes: { α i | − B = α < α < · · · < α S = B } . As a result, for any α ∈ [ − B, B ], a piecewise linear upper bound of α is ˜ α ,which is defined in Equation 9. As illustrated in Figure 3, α ≤ ˜ α and this upperbound can be arbitrarily tight as S → ∞ . This formulation has been used in[6] to discretize the space of unit vectors. In the rest of the paper, we use atilde to denote such an upper bound. Using these upper bounds, the equidistantconstraints can be approximated using the following conic constraints: ∀ i = 1 , · · · , K ∧ d = 1 , · · · , T (cid:107) n di − n ( d mod T )+1 i (cid:107) ≤ (2 √ B ) (1 − F i ) ∀ k = 1 , ∧ i = 1 , · · · , K ∧ d = 1 , · · · , T (cid:107) d ( d mod T )+1 ki (cid:107) ≤ ˜ dx dki + ˜ dy dki + (2 √ B ) F i (cid:107) d dki (cid:107) ≤ ˜ dx ( d mod T )+1 ki + ˜ dy ( d mod T )+1 ki + (2 √ B ) F i (10)where the last term on the right-hand sides is the big-M term that excludes fixednodes. The idea is to require the length of two vectors to be smaller than theupper bound of one another. Note that Equation 10 converges to Equation 6as S → ∞ . This formulation will require an upper bound for all d dki and eachupper bound requires (cid:100) log S (cid:101) binary variables. As a result, our formulation willintroduce O (4 T K (cid:100) log S (cid:101) ) binary decision variables. We also introduce a lastconstraint to ensure that rigid rods are not degenerate by ensuring minimal rodlength l min : ∀ k = 1 , ∧ i = 1 , · · · , K ∧ d = 1 , · · · , T ˜ dx dki + ˜ dy dki ≥ l min − ((2 √ B ) + l min ) F i . (11) oint Search of Topology and Trajectory for Planar Linkages 9 v L v R γ γ γ γ γ γ γ γ (a) (cid:15) d d i d d i (b) γ γ γ γ γ γ γ γ γ γ γ γ γ γ γ γ (c) Fig. 4: Linear relaxation of angle constraints. (a): We cut SO (2) into 8 sectors,each of which is selected by a γ -flag using SOS constraints. A sector, e.g.the sector selected by γ , is bounded by its left/right unit-length plane-normalvectors v L / v R . (b): If d d i falls in the red area, then we restrict d d i to its lefthalf-space (gray), which is at least (cid:15) -apart (blue). However, note that when d d i moves across sector boundaries, the gray area will jump discontinuously. (c): Toavoid discontinuous changes in the restricted region for d d i when d d i undergoescontinuous changes, we propose to double cover SO (2) using 2 S = 16 sectors.By ensuring a fixed rigid rod length across all time instances, we can makesure that all the end-effector positions n dK can be achieved using the same pla-nar linkage structure. In practice, however, we can only change the end-effectorposition by moving the first motor node n , so we still need to ensure that themechanics system will not glitch or does not have singular configurations. Themost intuitive classification of singular configuration is the rank-deficiency of theJacobian matrix [2]. However, this classification cannot be used in an MICP for-mulation because it is non-convex and the Jacobian matrix cannot be computedunder our implicit representation of rigid rods. Instead, we adopt a heuristicproposed by [22], which avoids singularities by ensuring that, for any movablenode n i , the two vectors d d i and d d i are not colinear. In other words, the trianglearea formed by these two vectors is positive. This constraint takes the followingbilinear form: d d i × d d i ≥ (cid:15), where (cid:15) is a small constant. Although this constraint is bilinear, we can useMcCormick envelopes [15] to relax it as a conic constraint. If the range [ − B, B ]is cut into S − O (4 T K (cid:100) log S (cid:101) ).However, a critical flaw of this formulation is that a McCormick envelope isan outer-approximation. As a result, the exact linkage structure can still besingular, although its conic relaxation is non-singular. To ensure strict non-singular formulation, we propose a constraint whereby the angle between thetwo vectors is larger than (cid:15) , which is equivalent to the positive area constraintwhen combined with Equation 11: (cid:93) d d i , d d i ≥ (cid:15). (12) In addition, we propose an inner approximation such that the exact linkagestructure is also guaranteed to be non-singular. Concretely, we cut the space of SO (2) into S sectors, as illustrated in Figure 4a, so that d d i will only fall intoone of the S sectors. If d d i falls in a particular sector, then we restrict d d i toits left half-space that is at least (cid:15) -apart, as shown in Figure 4b. If we use an SOS constraint to select the sector in which d d i falls, then only O ( T K (cid:100) log S (cid:101) )binary decision variables are needed. A minor issue with this formulation is thatthe allowed region of d d i jumps discontinuously as d d i changes continuously. Wecan fix this problem by double-covering the region of SO (2) using 2 S sectors,as shown in Figure 4c, which will introduce O ( T K (cid:100) log S (cid:101) ) binary decisionvariables.To formulate these constraints, we assume that each sector of SO (2) is flaggedby a selector variable γ l , which is bounded by its left/right unit-length plane-normal vectors v Ll / v Rl . Then the following constraints must be satisfied if d d i falls inside the sector: < v Ll , d d i > ≥ < v Rl , d d i > ≤ < R ( (cid:15) ) v Ll , d d i > ≤ < R ( π ) v Rl , d d i > ≥ , (13)where R ( • ) is the 2 × • . Combinedwith the fact that Equation 13 should only be satisfied for one particular sectorand that only movable nodes satisfy these constraints, we have the followingformulation: ∀ i = 2 , · · · , K d = 1 , · · · , T< v Ll , d d i > ≥ √ B ( γ dl,i − − √ BF i < v Rl , d d i > ≤ √ B (1 − γ dl,i ) + 2 √ BF i < R ( (cid:15) ) v Ll , d d i > ≤ √ B (1 − γ dl,i ) + 2 √ BF i < R ( π ) v Rl , d d i > ≥ √ B ( γ dl,i − − √ BF i { ( γ d , i , · · · , γ d S,i } ∈ SOS S (cid:88) l =1 γ dl,i = 1 . (14)These constraints will avoid singular configurations. Combining all the constraints, we minimize two objective function terms. First,we want the end-effector trajectory to match the target trajectory specified byusers. Second, to minimize manufacturing cost, we want to use as few rigidrods as possible. To formulate the first objective term, we need to replace atrajectory with a discrete number of samples. However, the order of these samples oint Search of Topology and Trajectory for Planar Linkages 11 is discarded. In practice, we find that better solutions can be found by preservingthe order between these samples. This requirement is formulated by makingsure that n dK will be visited by the end-effector sequentially when the motornode rotates by 2 π either clockwise or counter-clockwise. This requirement isformulated using the following MICP constraints: ∀ d = 1 , · · · , T − (cid:107) R ( 2 πT ) d d − d d +111 (cid:107) ≤ (2 √ B ) D (cid:107) R ( − πT ) d d − d d +111 (cid:107) ≤ (2 √ B ) (1 − D ) , (15)where D is a binary variable to choose which direction the motor rotates. Puttingeverything together, we arrive at the following MICP problem: argmin T (cid:88) d =1 (cid:107) n dK − n d ∗ K (cid:107) + w K (cid:88) i =1 U i s . t . Equation 2, 3, 4, 5, 7, 8, 10, 11, 14, 15, (16)where n d ∗ K are the sampled points on the target trajectory and w the regulariza-tion weight of the cost-efficiency term. Since non-convexity is not accepted byMICP, the solution returned by MICP is only a piecewise linear approximationof the original nonlinear problem. To return a solution with exact constraintsatisfaction, we refine the solution by solving an additional NLP locally usingthe following formulation: argmin T (cid:88) d =1 (cid:107) n dK − n d ∗ K (cid:107) s . t . Equation 2, 3, 4, 5, 7, 8, 6, 12, 15, (17)where we fix all the binary variables U i , F i , D . Note that Equation 17 is a mixed-integer NLP (MINLP) generalization of Equation 16 and we have the followinglemma: Lemma 1.
Equation 16 converges to Equation 17 as S → ∞ , and the BB algo-rithm can find the global optimum for Equation 16. We have implemented our method using Gurobi [8] as our MICP solver forEquation 16 and Knitro [4] as our NLP solver for Equation 17. All the experi-ments are performed on a cluster with 4 CPU cores per process (2.5GHz E5-2680CPU). Compared with prior work [22], the main benefit of our formulation isthat we can search for planar linkage structures from a target trajectory of theend-effector that requires trivial effort from users. In Figure 5, we show a list ofdifferent target trajectories and the optimized planar linkage structures.
Fig. 5: We show 10 different optimized planar linkage structures with the end-effector trajectory in blue and the user-specified target trajectory in yellow. Forall these examples, we choose K = 5 ∼ S = 9, and T = 10 ∼
20. Theend-effector trajectory matches closely with the target trajectory.
Instance C o m p u t a t i o n a l T i m e ( h r ) K = 5 K = 6 K = 7 Computational Time vs. K S C o m p u t a t i o n a l T i m e ( h r ) Computational Time vs. S
15 20 25 30 T C o m p u t a t i o n a l T i m e ( h r ) Computational Time vs. T Fig. 6: We plot the average computational time for solving MICP in 35 exampleproblems using different parameters in Figure 5. The computational time forsolving MICP grows exponentially with K , log S , and T .The performance and accuracy of our algorithm heavily depend on the threeparameters: the max number of rigid rods K , the number of pieces for approx-imating S , and the number of samples on the target trajectory T . Since thecost of solving MICP grows exponentially with the number of binary decisionvariables, which is proportional to K , log S , and T , our method cannot scaleto large problems, as illustrated in Figure 6. In practice, we find that, given amaximal computational time of 10 hours, we can compute globally optimal so-lutions for most benchmarks with K ≤ S ≤
9, and T ≤
20. This is enough ifwe design robots part-by-part, as is done in the Theo Jansen’s strandbeest. Forother benchmarks, the computational time is longer than 10 hours, but a feasi-ble solution has been found, although it is sub-optimal. In Figure 7, we plot theaverage convergence history of a typical optimization. Since we express all thetopology and geometric requirements as hard mixed-integer constraints, feasiblesolutions are quite rare in the search space and the optimizer takes most of the oint Search of Topology and Trajectory for Planar Linkages 13 O b j e c t i v e F un c t i o n Convergence History O b j e c t i v e F un c t i o n Convergence History O b j e c t i v e F un c t i o n Convergence History
Fig. 7: We plot the convergence history curve for 3 typical optimizations byshowing the objective function values plotted against the number of nodes ex-plored in the BB search tree. The BB algorithm spends most of its time exploringinfeasible nodes and the first identified feasible solution is usually very close tothe optimal solution, so that the optimizer will return the globally optimal so-lution after refining the solution for 5 −
10 times.computational time pruning infeasible solutions. Once the first feasible solutionis found, it is usually very close to the optimal solution and the optimizer refinesit for less than 10 times to reach the optimal solution. (a) (b) (c) (d)
Fig. 8: SA can find good enough solutions for simple target curves (a). However,for more complex curve shapes, SA failes (b) while MICP succeeds (c). We alsoplot the objective function values returned by SA and MICP in 10 computationalexamples in (d), where MICP outperforms SA in 9 instances.We have also compared our method with conventional global search algo-rithms such as simulated annealing (SA). We implemented a similar algorithmas proposed in [29]. In this algorithm, we randomly generate 1000000 samples byrandom moves and accept these samples according to the simulated annealingrule. Each random move can be of one of three kinds: geometric change, nodeaddition, and node removal. In geometric change, the length of a rigid rod israndomly perturbed. In node addition, a new node is added and the length ofthe new rigid rods are randomly picked. In node removal, the end-effector nodeis removed and the last movable node is used as the new end-effector node. Weenhance standard SA algorithm by making sure that each random move is valid.In other words, we introduce an inner loop and repeated try random moves untilthe modified planar linkage structure satisfies all the topological constraints andhas no singular configurations. As illustrated in Figure 8a, SA algorithm canfind satisfactory results for simple target curves, but SA usually fails for morecomplex curve shapes (Figure 8bc). In Figure 8d, we also show the objectivefunction values after convergence. The solution of MICP is almost always better than the solution of SA. However, SA outperforms MICP in one example, whichis probably due to the inexact constraint satisfaction of MICP.Usually, the design of a planar linkage structure is not only subject to a targetend-effector trajectory, but also to various other user constraints. For example,the user might require certain nodes to be fixed, which can be easily achievedusing our MICP formulation. The user may also reserve certain parts of therobot for some functional units that cannot be occupied by the planar linkages.This type of constraint can be expressed as collision avoidance between a planarlinkage structure and a specified convex region, which can be formulated asMICP constraints using a prior method [7]. In Figure 9, we show results takingthese constraints into consideration. (a) (b) (c)FixedToo Close Bounded Region
Fig. 9: We show results taking two different user constraints into consideration.(a): Results with no constraints. The optimizer is guided by the regularizationterm to use as few nodes as possible. (b): We fix the center of rotation and theoptimizer finds a more complex structure with 6 nodes. (c): If we do not wantthe structure to be too close to the target curve, we can add a bounded regionand create a constraint that any nodes (other than the end-effector node) shouldbe inside the bounded region.
We present a globally optimal formulation to jointly search for both the topol-ogy and geometry of a planar linkage structure. Our formulation relaxes theproblem into a MICP, for which optimal solutions can be found efficiently usingBB algorithms. Our results show that our formulation can search for complexstructures from trivial and intuitive user inputs, i.e. target end-effector trajec-tories. Additionally, various design constraints can be easily incorporated. Formoderately complex structures, the solve time using these formulations falls inthe range between minutes and hours. oint Search of Topology and Trajectory for Planar Linkages 15(a) (b) (c) (d)Ambiguity Colinear Sub-Optimal
Fig. 10: Failure cases and issues with our formulation. (a,b): MICP only returnsthe single global optimum. But similar target trajectories can lead to two differ-ent linkage structures. (c): We only satisfy geometric constraints approximately,so that the linkage structure might not satisfy these constraints exactly. In thisexample, we have two rigid rods being colinear. (d): Usually, the early feasiblesolutions found by MICP are of low-quality, and we have to wait for the MICPto find the global optimum.As a major limitation, the solve time increases quickly with the number ofpossible rigid bodies in the planar structure ( K ) and the number of sampleson the target trajectory ( T ) because the number of decision variables dependson a multiplication of these two parameters. A related issue is that MICP onlysatisfies the geometric constraints approximately. As illustrated in Figure 10c,a predicted target trajectory with approximate constraints satisfaction can bedifferent from a predicted target trajectory with exact constraints satisfactionafter solving Equation 17. To reduce the approximation error, we have to increasethe approximation granularity by using a larger S , which in turn increases thenumber of binary decision variables. Finally, note that our formulation doesnot generate all possible planar linkage structures but only those that allowsequential forward kinematic processing. This problem is inherited from [13,1]by using the same representation as these works. Allowing more general planarlinkages is also possible under the MICP formulation by using a new formulationof topology constraints. Our future research will focus on a balance between global optimality and for-mulation efficiency. Such a balance could possibly be achieved by using MINLPformulations. In addition, we observe that different planar linkages, as shown inFigure 10ab, can generate very similar target trajectories. This indicates thatthere exist many local optima with objective function close to the global opti-mum. However, a BB algorithm will only return the single global optimum. Inaddition, we found that we need to wait until the BB algorithm finds its globaloptimum; the intermediary solutions are not usually usable, as illustrated inFigure 10d. A potential future direction is to use algorithms such as Bayesian optimization that can explore multiple local optima and return many solutionsfor users to make a choice.