Globally Optimal Joint Search of Topology and Trajectory for Planar Linkages

TL;DR

This paper introduces a globally optimal method for jointly optimizing the topology and trajectory of planar linkages, enabling efficient design of low-cost robots with complex motion capabilities.

Contribution

It formulates the joint topology and trajectory search as a mixed-integer convex programming problem solved via branch-and-bound, which is novel in this context.

Findings

More efficient linkage structure discovery

More accurate end-effector trajectory generation

Global optimality achieved through MICP and BB

Abstract

We present a method to find globally optimal topology and trajectory jointly for planar linkages. Planar linkage structures can generate complex end-effector trajectories using only a single rotational actuator, which is very useful in building low-cost robots. We address the problem of searching for the optimal topology and geometry of these structures. However, since topology changes are non-smooth and non-differentiable, conventional gradient-based searches cannot be used. We formulate 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 experiments show that the proposed approach finds complex linkage structures more efficiently and generates end-effector trajectories more accurately.

Figures

Equations

$$\displaystyle\mathbf{n}_{1}(t)=\left(\begin{array}[]{cc}{\mathbf{sin}(t)R+X_{C}},&{\mathbf{cos}(t)R+Y_{C}}\end{array}\right),$$

$$\exists U_i,F_i\in \{0,1\}$$

$$1-F_i\le U_i$$

$$U_1=U_K=1$$

$$F_1=0$$

$$\exists C_{ji},C_{ji}^1,C_{ji}^2\in \{0,1\}$$

$$C_{ji}=C_{ji}^1+C_{ji}^2$$

$$C_{ji}^1\le U_j\wedge C_{ji}^2\le U_j$$

$$\sum _{j=1}^{i-1}{C}_{ji}=2-2F_i$$

$$\exists C_{0i}^d\in \{0,1\}$$

$$\sum _{j=0}^{i-1}{C}_{ji}^d=1$$

$$\displaystyle\exists C_{ji},C_{ji}^{1},C_{ji}^{2}\in{\color[rgb]{0.81640625,0.15625,0.140625}[0,1]}$$

$$C_{ji}=C_{ji}^1+C_{ji}^2$$

$$C_{ji}^1\le U_j\wedge C_{ji}^2\le U_j$$

$$\sum _{j=1}^{i-1}{C}_{ji}=2-2F_i$$

$$\displaystyle\exists C_{0i}^{d}\in{\color[rgb]{0.81640625,0.15625,0.140625}[0,1]}$$

$$\displaystyle{\color[rgb]{0.81640625,0.15625,0.140625}\{C_{ji}^{d}|j=0,\cdots,i-1\}\in\mathcal{SOS}_{1}}$$

$$\sum _{j=0}^{i-1}{C}_{ji}^d=1$$

$$\exists Q_{ji}\in [0,\infty ]$$

$$Q_{ji}\le C_{ji}K$$

$$U_i+\sum _{j=1}^{i-1}{Q}_{ji}=\sum _{k=i+1}^K{Q}_{ik}$$

$$\exists R_{ji}\in [0,\infty ]$$

$$R_{ji}\le C_{ji}K\wedge R_{ji}\le (1-F_j)K$$

$$\sum _{j=1}^{i-1}{R}_{ji}=1-F_i+\sum _{k=i+1}^K{R}_{ik}$$

$$\Vert {\mathbf{n}}_j^d-{\mathbf{n}}_i^d{\Vert }^2=\Vert {\mathbf{n}}_j^{(d\mathrm{mod}T)+1}-{\mathbf{n}}_i^{(d\mathrm{mod}T)+1}{\Vert }^2\forall 1\le d\le T,$$

$$\exists \{dx,dy{\}}_{ki}^d\forall k=1,2\wedge i=2,\cdots ,K\wedge d=1,\cdots ,T$$

$$|\{dx,dy{\}}_{ki}^d-\{x,y{\}}_j^d+\{x,y{\}}_i^d|\le 2B(1-C_{ji}^k)\forall j=1,\cdots ,i-1,$$

$$\{dx,dy{\}}_{11}^d=\{dx,dy{\}}_{21}^d=\{x_1^d-X_C,y_1^d-Y_C\},$$

$$\{{\alpha }_i|-B={\alpha }_1<{\alpha }_2<\cdots <{\alpha }_S=B\}.$$

$$\displaystyle\left(\begin{array}[]{c}\alpha\\ \tilde{\alpha}\end{array}\right)=\sum_{s=1}^{S}\lambda_{s}\left(\begin{array}[]{c}\alpha_{s}\\ \alpha_{s}^{2}\end{array}\right)$$

$$\{{\lambda }_{1,\cdots ,S}\}\in {\mathcal{S}\mathcal{O}\mathcal{S}}_2$$

$$\sum _{s=1}^S{\lambda }_s=1$$

$$\forall i=1,\cdots ,K\wedge d=1,\cdots ,T$$

$$\Vert {\mathbf{n}}_i^d-{\mathbf{n}}_i^{(d\mathrm{mod}T)+1}{\Vert }^2\le (2\sqrt{2}B{)}^2(1-F_i)$$

$$\forall k=1,2\wedge i=1,\cdots ,K\wedge d=1,\cdots ,T$$

$$\Vert {\mathbf{d}}_{ki}^{(d\mathrm{mod}T)+1}{\Vert }^2\le {\widetilde{dx}}_{ki}^d+{\widetilde{dy}}_{ki}^d+(2\sqrt{2}B{)}^2{F}_i$$

$$\Vert {\mathbf{d}}_{ki}^d{\Vert }^2\le {\widetilde{dx}}_{ki}^{(d\mathrm{mod}T)+1}+{\widetilde{dy}}_{ki}^{(d\mathrm{mod}T)+1}+(2\sqrt{2}B{)}^2{F}_i$$

$$\forall k=1,2\wedge i=1,\cdots ,K\wedge d=1,\cdots ,T$$

$${\widetilde{dx}}_{ki}^d+{\widetilde{dy}}_{ki}^d\ge l_{min}^2-((2\sqrt{2}B{)}^2+l_{min}^2)F_i.$$

$${\mathbf{d}}_{1i}^d\times {\mathbf{d}}_{2i}^d\ge ϵ,$$

$$\measuredangle {\mathbf{d}}_{1i}^d,{\mathbf{d}}_{2i}^d\ge ϵ.$$

$$<{\mathbf{v}}_l^L,{\mathbf{d}}_{1i}^d>\ge 0$$

$$<{\mathbf{v}}_l^R,{\mathbf{d}}_{1i}^d>\le 0$$

$$<\mathbf{R}(ϵ){\mathbf{v}}_l^L,{\mathbf{d}}_{2i}^d>\le 0$$

Full text

1 Introduction

A planar linkage is a mechanical structure built with a set of rigid bodies connected by hinge joints. This structure typically has one effective degree-of-freedom actuated by a rotational motor. Since they impose a minimal burden on controller design, these structures are widely used as building blocks for low-cost toys 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 can fulfill various requirements of different types of locomotion, including walking and swimming [11, 22].

A challenging problem in mechanics design is to find the linkage structure with an end-effector that will trace out a given curve. This problem is challenging in that it searches over three coupled variables: topology, geometry, and trajectory. The linkage topology determines which rigid bodies are connected and the order of their connections. Clearly, the topology is a non-smooth and non-differentiable decision variable. The linkage geometry determines the shape of each rigid link. Finally, the trajectory determines the pose of the linkage structure at each time instance. The last two variables are smooth and differentiable, but directly optimizing them induces non-convex functions. Previous works [9, 29] have proposed various solutions to address problems of this kind. These methods rely on random searches, such as $A^{*}$ [9] and covariance matrix adaptation [29], to try different topologies. Then, for each topology, they perform non-linear programming (NLP) under the given topology to determine the geometry and trajectory. However, these methods are computationally expensive because a huge number of samples are needed for the random search to converge. Moreover, even after determining the topology, these methods can find only sub-optimal solutions due to the non-convex nature of NLP.

Main Results: Given the input of a target trajectory, we present a new method that can efficiently compute a planar linkage structure with globally optimal topology and geometry configurations and an accurate trajectory reproduction. Based on recent advances in mixed-integer modeling [24, 6, 23], we relax this joint search problem as an MICP problem, the global optimum of which is arbitrarily close to the global optimum of the original problem. The main benefit of MICP relaxation is that the search can be accomplished efficiently using the BB algorithm. BB is more strategic than random search, as used by [29], because it cuts impossible or sub-optimal search spaces at an early stage, leading to higher efficiency. We have compared MICP with prior methods using different examples. The results show that our proposed MICP approach finds solutions more efficiently and that the resulting structure matches the target trajectory more closely.

In the rest of the paper, we first review related work in Section 2 and then formulate our joint search problem in Section 3. The MICP model and various constraints required for the integrity of the planar linkage are presented in Section 4. Results and the evaluation of the proposed approach are given in Section 5.

2 Related Work

In this section, we review related work in robot design optimization, mixed-integer modeling, and topology optimization.

Robot Design Optimization: Robot design optimization is a superset of conventional topology and truss optimization [16] where the decision variables are only topology or geometry. This is because the specification of a robot design is 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 a result, many prior methods do not work since they only optimize a subset of decision variables [10, 22, 1, 19, 21]. Recent works [29, 9, 20] search for all variables simultaneously. However, these methods are based on random search techniques, which usually require a large amount of trial and error and find sub-optimal solutions.

Mixed-Integer Modeling: The main benefit of mixed-integer modeling is the use of the well-studied BB algorithm [14]. BB allows us to find the global optimum of non-convex programming problems, while only visiting a small fraction of the search space. Mixed-integer models have been applied to a large variety of problems including motion planning [7], inverse kinematics [6], network flows [5], and mesh generations [3]. By applying the big-M method [23], McCormick envelopes and piecewise approximations [15], and general non-convex problems can be easily relaxed as MICP problems. Prior works [12, 17] have also formulated topology optimization problems as MICP. However, our work is the first to formulate the planar linkage problem as MICP and we employ MICP to concurrently 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 and use gradient-based method to search for locally optimal structures over a search space of millions of dimensions. This technique has been widely used in the design of soft robots [27, 26, 28]. However, the optimization of articulated robots is more challenging 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 variables.

3 Joint Search for Planar Linkages

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 we have 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$ distinct positions, denoted as node set: $\mathbf{n}_{1,\cdots,N}$. Of these nodes, $\mathbf{n}_{1}$ is the rotational motor and $\mathbf{n}_{N}$ is the end-effector. Within one limit cycle, $\mathbf{n}_{1}$ follows a circular curve centered at \left(\begin{array}[]{cc}{X_{C}},&{Y_{C}}\end{array}\right) with a radius $R$:

[TABLE]

which induces trajectories of other nodes $\mathbf{n}_{i}(t)$ via forward kinematics. The other $N-2$ nodes can be one of two kinds: fixed or movable. In addition, a rigid body may exist between each pair of nodes $\mathbf{n}_{i,j}$, in which case $\|\mathbf{n}_{i}(t)-\mathbf{n}_{j}(t)\|$ is a constant.

Given these definitions, the input to our problem is a target end-effector trajectory $\mathbf{n}_{N}^{*}(t)$. The output of our method is the following set of variables defining both the topology and geometry of a planar linkage:

• •

An integer vector of size $N$ (the number of nodes), which containing the type of each node: fixed or movable.

• •

An $N\times N$ symmetric binary matrix $C^{N\times N}$ where $C_{ij}=1$ means a rigid body connects $\mathbf{n}_{i,j}$.

• •

The position of $\mathbf{n}_{1,\cdots,N}(t)$ at a certain, arbitrary time instance $t$.

The goal of our method is to find the globally optimal set of variables that minimizes $\int\|\mathbf{n}_{N}(t)-\mathbf{n}_{N}^{*}(t)\|^{2}dt$.

4 MICP Formulation of Joint Search

In this section, we present a set of linear constraints and quadratic objective functions for relaxing the joint search as an MICP problem. We first introduce the set of topology constraints to ensure the well-posed nature of the structure in Section 4.1 and then present constraints and objective functions for geometric correctness in Section 4.2.

4.1 Topology Constraints

As illustrated in Figure 2b, our method is based on the symbolic representation presented in [13, 22], which assumes that each movable node is attached to two other 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 linkage structures with closed loops.

Since the number of nodes is unknown, we assume that the maximal number of nodes is $K>N$. For each node other than the first motor node $\mathbf{n}_{1}$, we need a binary variable $U_{i}$ such that $U_{i}=1$ indicates $\mathbf{n}_{i}$ is used as a part of the planar linkage structure. In addition, we need another binary variable $F_{i}$ such that $F_{i}=1$ indicates $\mathbf{n}_{i}$ is fixed and $F_{i}=0$ indicates $\mathbf{n}_{i}$ is movable. These two sets of variables are under the constraint that only a used node can be movable. In addition, we assume that the last node $\mathbf{n}_{K}$ is the end-effector that must be used. In summary, we introduce the following sets of variables and node-state constraints:

[TABLE]

Our next set of constraints ensures local topology correctness. It ensures that each 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 and the position of the movable node can then be determined via the Law of Cosine [10]. We introduce auxiliary variables $C_{ji}^{1}$ to indicate whether $\mathbf{n}_{j}$ is the first node to which $\mathbf{n}_{i}$ is connected. $C_{ji}^{2}$ indicates whether $\mathbf{n}_{j}$ is the second node to which $\mathbf{n}_{i}$ is connected. In addition, we introduce two verbose variables $C_{0i}^{1,2}=1$ to indicate that $\mathbf{n}_{i}$ is connected to nothing. The resulting constraint set is:

[TABLE]

When $\mathbf{n}_{i}$ is fixed in the above formulation, then $F_{i}=1$ in Equation 4 and all $C_{ji}$ are zero except for $C_{0i}^{1,2}=1$ due to the sum-to-one constraints. If $\mathbf{n}_{i}$ is movable, then $F_{i}=0$ and $C_{ji}$ sums to two. As a result, there must be $j_{1},j_{2}<i$ such that $C_{j_{1}i}^{1}=1$ and $C_{j_{2}i}^{2}=1$. Note that $j_{1}$ and $j_{2}$ must be different because otherwise the constraint that $C_{ji}\in[0,1]$ will be violated. In addition, since the first node $\mathbf{n}_{1}$ is the motor node, it is excluded from these connectivity constraints. However, this naive formulation will require binary variables for each pair of $\mathbf{n}_{j}$ and $\mathbf{n}_{i}$, which requires $\mathcal{O}(K^{2})$ binary variables all together. Instead, we adopt the idea of special ordered set of type 1 ($\mathcal{SOS}_{1}$) [25] and model these constraints using $\mathcal{O}(K\lceil\mathbf{log}K\rceil)$ binary variables. Intuitively, $\mathcal{SOS}_{1}$ constrains that only one variable in a set can take a non-zero value and it can be achieved by using a logarithm number of binary variables. The improved constraint set is:

[TABLE]

Finally, we introduce a third set of constraints to ensure global topology correctness. This set of constraints ensures that the linkage structure contains no wasted structures. In other words, each node must have some influence on the trajectory of the end-effector node and the first motor node must be connected to others. We model these constraints using the MICP formulation of network flows [5]. Specifically, each node $\mathbf{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 nodes with capacity $Q_{ji}$. We require inward-outward flux balance for each node except for the end-effector node:

[TABLE]

where we adopt the big-M method [23] in the second constraint to ensure that only edges between connected nodes can have a capacity up to $K$. Using a similar idea, we also formulate a constraint that a movable node must be connected to at least one other movable node. We assume that each node $\mathbf{n}_{i}$ generates a reversed outward flux that equals to $1-F_{i}$, and we assume that there is a flow edge defined between each pair of nodes with capacity $R_{ji}$. We require inward-outward flux balance for each node except for the motor node:

[TABLE]

These three constraints ensure that the planar linkage structure is symbolically correct, independent of the concrete geometric shape.

4.2 Geometric Correctness

The main utility of geometric correctness constraints is to compute the exact positions \mathbf{n}_{i}=\left(\begin{array}[]{cc}{x_{i}},&{y_{i}}\end{array}\right) of each node in the 2D workspace. These positions are functions of time $t$ and we sample a set of $T$ discrete time instances $t^{1,\cdots,T}$. In this section, we will always use superscripts for timestep indices. For example, at time instance $t^{d}$, the position of $\mathbf{n}_{i}$ is $\mathbf{n}_{i}^{d}$. We want to find a common geometric specification such that all the end-effector positions $\mathbf{n}_{K}^{1,\cdots,T}$ can be achieved.

The most important geometric variable is the length of each rigid rod. We define these parameters implicitly using a set of constraints such that, if $\mathbf{n}_{i}$ and $\mathbf{n}_{j}$ are connected, then the distance between these two nodes is a constant for all time instances. In other words, we need the following set of constraints if $C_{ji}=1$:

[TABLE]

after which any distance $\|\mathbf{n}_{j}^{d}-\mathbf{n}_{i}^{d}\|^{2}$ can be used as the rigid rod length.

However, there are two challenging issues in modeling these constraints that can affect the performance of the MICP solver. A first challenge is to minimize the use of binary variables. Because any pair of nodes $\mathbf{n}_{j}$ and $\mathbf{n}_{i}$ might be connected, a naive formulation will require a number of binary variables proportional to $K^{2}$. Instead, we introduce auxiliary term \mathbf{d}_{1i}^{d}=\left(\begin{array}[]{cc}{dx_{1i}^{d}},&{dy_{1i}^{d}}\end{array}\right), which indicates the relative position between $\mathbf{n}_{i}$ and the first other node connected to it at time instance $t^{d}$. Similarly, \mathbf{d}_{2i}^{d}=\left(\begin{array}[]{cc}{dx_{2i}^{d}},&{dy_{2i}^{d}}\end{array}\right) indicates the relative position between $\mathbf{n}_{i}$ and the second other node connected to it. These definitions induce the following big-M constraints:

[TABLE]

where $B$ is the big-M parameter, implying that all the node positions lie in a bounded region $[-B,B]^{2}$. Note that the first motor node $\mathbf{n}_{1}$ follows a circular curve (Equation 2), which requires special definitions of $\mathbf{d}_{11}^{d},\mathbf{d}_{21}^{d}$ as follows:

[TABLE]

where the center of rotation \left(\begin{array}[]{cc}{X_{C}},&{Y_{C}}\end{array}\right) is used as an additional auxiliary variable. The second challenge is that these constraints are non-convex because they involve quadratic terms. Fortunately, efficient formulations have been developed to relax non-convex functions using piecewise linear approximation [15] and a special ordered set of type 2 ($\mathcal{SOS}_{2}$) [25]. $\mathcal{SOS}_{2}$ effects a constraint that at most two of the variables in an ordered set with consecutive indices can take non-zero values. To use these formulations, we decompose the range $[-B,B]$ evenly into $S-1$ pieces with $S$ nodes:

[TABLE]

As a result, for any $\alpha\in[-B,B]$, a piecewise linear upper bound of $\alpha^{2}$ is $\tilde{\alpha}$, which is defined in Equation 10.

As illustrated in Figure 3, $\alpha^{2}\leq\tilde{\alpha}$ and this upper bound can be arbitrarily tight as $S\to\infty$. This formulation has been used in [6] to discretize the space of unit vectors. In the rest of the paper, we use a tilde to denote such an upper bound. Using these upper bounds, the equidistant constraints can be approximated using the following conic constraints:

[TABLE]

where the last term on the right-hand sides is the big-M term that excludes fixed nodes. The idea is to require the length of two vectors to be smaller than the upper bound of one another. Note that Equation LABEL:eq:equidistant converges to Equation 7 as $S\to\infty$. This formulation will require an upper bound for all $\mathbf{d}_{ki}^{d}$ and each upper bound requires $\lceil\mathbf{log}S\rceil$ binary variables. As a result, our formulation will introduce $\mathcal{O}(4TK\lceil\mathbf{log}S\rceil)$ binary decision variables. We also introduce a last constraint to ensure that rigid rods are not degenerate by ensuring minimal rod length $l_{min}$:

[TABLE]

By ensuring a fixed rigid rod length across all time instances, we can make sure that all the end-effector positions $\mathbf{n}_{K}^{d}$ can be achieved using the same planar linkage structure. In practice, however, we can only change the end-effector position by moving the first motor node $\mathbf{n}_{1}$, so we still need to ensure that the mechanics system will not glitch or does not have singular configurations. The most intuitive classification of singular configuration is the rank-deficiency of the Jacobian matrix [2]. However, this classification cannot be used in an MICP formulation because it is non-convex and the Jacobian matrix cannot be computed under our implicit representation of rigid rods. Instead, we adopt a heuristic proposed by [22], which avoids singularities by ensuring that, for any movable node $\mathbf{n}_{i}$, the two vectors $\mathbf{d}_{1i}^{d}$ and $\mathbf{d}_{2i}^{d}$ are not colinear. In other words, the triangle area formed by these two vectors is positive. This constraint takes the following bilinear form:

[TABLE]

where $\epsilon$ is a small constant. Although this constraint is bilinear, we can use McCormick envelopes [15] to relax it as a conic constraint. If the range $[-B,B]$ is cut into $S-1$ segments, then this formulation will introduce $\mathcal{O}(4TK\lceil\mathbf{log}S\rceil)$. However, a critical flaw of this formulation is that a McCormick envelope is an outer-approximation. As a result, the exact linkage structure can still be singular, although its conic relaxation is non-singular. To ensure strict non-singular formulation, we propose a constraint whereby the angle between the two vectors is larger than $\epsilon$, which is equivalent to the positive area constraint when combined with Equation LABEL:eq:miniLength:

[TABLE]

In addition, we propose an inner approximation such that the exact linkage structure is also guaranteed to be non-singular. Concretely, we cut the space of $\mathcal{SO}(2)$ into $S$ sectors, as illustrated in Figure 4a, so that $\mathbf{d}_{1i}^{d}$ will only fall into one of the $S$ sectors. If $\mathbf{d}_{1i}^{d}$ falls in a particular sector, then we restrict $\mathbf{d}_{2i}^{d}$ to its left half-space that is at least $\epsilon$-apart, as shown in Figure 4b. If we use an $\mathcal{SOS}_{1}$ constraint to select the sector in which $\mathbf{d}_{1i}^{d}$ falls, then only $\mathcal{O}(TK\lceil\mathbf{log}S\rceil)$ binary decision variables are needed. A minor issue with this formulation is that the allowed region of $\mathbf{d}_{2i}^{d}$ jumps discontinuously as $\mathbf{d}_{1i}^{d}$ changes continuously. We can fix this problem by double-covering the region of $\mathcal{SO}(2)$ using $2S$ sectors, as shown in Figure 4c, which will introduce $\mathcal{O}(TK\lceil\mathbf{log}2S\rceil)$ binary decision variables.

To formulate these constraints, we assume that each sector of $\mathcal{SO}(2)$ is flagged by a selector variable $\gamma_{l}$, which is bounded by its left/right unit-length plane-normal vectors $\mathbf{v}_{l}^{L}$/$\mathbf{v}_{l}^{R}$. Then the following constraints must be satisfied if $\mathbf{d}_{1i}^{d}$ falls inside the sector:

[TABLE]

where $\mathbf{R}(\bullet)$ is the $2\times 2$ counter-clockwise rotation matrix by angle $\bullet$. Combined with the fact that Equation 14 should only be satisfied for one particular sector and that only movable nodes satisfy these constraints, we have the following formulation:

[TABLE]

These constraints will avoid singular configurations.

4.3 The Complete MICP Formulation

Combining all the constraints, we minimize two objective function terms. First, we want the end-effector trajectory to match the target trajectory specified by users. Second, to minimize manufacturing cost, we want to use as few rigid rods as possible. To formulate the first objective term, we need to replace a trajectory with a discrete number of samples. However, the order of these samples is discarded. In practice, we find that better solutions can be found by preserving the order between these samples. This requirement is formulated by making sure that $\mathbf{n}_{K}^{d}$ will be visited by the end-effector sequentially when the motor node rotates by $2\pi$ either clockwise or counter-clockwise. This requirement is formulated using the following MICP constraints:

[TABLE]

where $D$ is a binary variable to choose which direction the motor rotates. Putting everything together, we arrive at the following MICP problem:

[TABLE]

where $\mathbf{n}_{K}^{d*}$ are the sampled points on the target trajectory and $w$ the regularization weight of the cost-efficiency term. Since non-convexity is not accepted by MICP, the solution returned by MICP is only a piecewise linear approximation of the original nonlinear problem. To return a solution with exact constraint satisfaction, we refine the solution by solving an additional NLP locally using the following formulation:

[TABLE]

where we fix all the binary variables $U_{i},F_{i},D$. Note that Equation 18 is a mixed-integer NLP (MINLP) generalization of Equation 17 and we have the following lemma:

Lemma 4.1**.**

Equation 17 converges to Equation 18 as $S\to\infty$, and the BB algorithm can find the global optimum for Equation 17.

5 Results and Evaluations

We have implemented our method using Gurobi [8] as our MICP solver for Equation 17 and Knitro [4] as our NLP solver for Equation 18. All the experiments are performed on a cluster with 4 CPU cores per process (2.5GHz E5-2680 CPU). Compared with prior work [22], the main benefit of our formulation is that we can search for planar linkage structures from a target trajectory of the end-effector that requires trivial effort from users. In Figure 5, we show a list of different target trajectories and the optimized planar linkage structures.

The performance and accuracy of our algorithm heavily depend on the three parameters: the max number of rigid rods $K$, the number of pieces for approximating $S$, and the number of samples on the target trajectory $T$. Since the cost of solving MICP grows exponentially with the number of binary decision variables, which is proportional to $K$, $\mathbf{log}S$, and $T$, our method cannot scale to large problems, as illustrated in Figure 6. In practice, we find that, given a maximal computational time of $10$ hours, we can compute globally optimal solutions for most benchmarks with $K\leq 7$, $S\leq 9$, and $T\leq 20$. This is enough if we design robots part-by-part, as is done in the Theo Jansen’s strandbeest. For other benchmarks, the computational time is longer than $10$ hours, but a feasible solution has been found, although it is sub-optimal. In Figure 7, we plot the average convergence history of a typical optimization. Since we express all the topology and geometric requirements as hard mixed-integer constraints, feasible solutions are quite rare in the search space and the optimizer takes most of the computational time pruning infeasible solutions. Once the first feasible solution is found, it is usually very close to the optimal solution and the optimizer refines it for less than 10 times to reach the optimal solution.

We have also compared our method with conventional global search algorithms such as simulated annealing (SA). We implemented a similar algorithm as proposed in [29]. In this algorithm, we randomly generate $1000000$ samples by random moves and accept these samples according to the simulated annealing rule. Each random move can be of one of three kinds: geometric change, node addition, and node removal. In geometric change, the length of a rigid rod is randomly perturbed. In node addition, a new node is added and the length of the new rigid rods are randomly picked. In node removal, the end-effector node is removed and the last movable node is used as the new end-effector node. We enhance 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 until the modified planar linkage structure satisfies all the topological constraints and has no singular configurations. As illustrated in Figure 8a, SA algorithm can find satisfactory results for simple target curves, but SA usually fails for more complex curve shapes (Figure 8bc). In Figure 8d, we also show the objective function values after convergence. The solution of MICP is almost always better than the solution of SA. However, SA outperforms MICP in one example, which is probably due to the inexact constraint satisfaction of MICP.

Usually, the design of a planar linkage structure is not only subject to a target end-effector trajectory, but also to various other user constraints. For example, the user might require certain nodes to be fixed, which can be easily achieved using our MICP formulation. The user may also reserve certain parts of the robot for some functional units that cannot be occupied by the planar linkages. This type of constraint can be expressed as collision avoidance between a planar linkage structure and a specified convex region, which can be formulated as MICP constraints using a prior method [7]. In Figure 9, we show results taking these constraints into consideration.

6 Conclusion & Limitations

We present a globally optimal formulation to jointly search for both the topology and geometry of a planar linkage structure. Our formulation relaxes the problem into a MICP, for which optimal solutions can be found efficiently using BB algorithms. Our results show that our formulation can search for complex structures from trivial and intuitive user inputs, i.e. target end-effector trajectories. Additionally, various design constraints can be easily incorporated. For moderately complex structures, the solve time using these formulations falls in the range between minutes and hours.

As a major limitation, the solve time increases quickly with the number of possible rigid bodies in the planar structure ($K$) and the number of samples on the target trajectory ($T$) because the number of decision variables depends on a multiplication of these two parameters. A related issue is that MICP only satisfies the geometric constraints approximately. As illustrated in Figure 10c, a predicted target trajectory with approximate constraints satisfaction can be different from a predicted target trajectory with exact constraints satisfaction after solving Equation 18. To reduce the approximation error, we have to increase the approximation granularity by using a larger $S$, which in turn increases the number of binary decision variables. Finally, note that our formulation does not generate all possible planar linkage structures but only those that allow sequential forward kinematic processing. This problem is inherited from [13, 1] by using the same representation as these works. Allowing more general planar linkages is also possible under the MICP formulation by using a new formulation of topology constraints.

6.1 Future Work

Our future research will focus on a balance between global optimality and formulation efficiency. Such a balance could possibly be achieved by using MINLP formulations. In addition, we observe that different planar linkages, as shown in Figure 10ab, can generate very similar target trajectories. This indicates that there exist many local optima with objective function close to the global optimum. However, a BB algorithm will only return the single global optimum. In addition, we found that we need to wait until the BB algorithm finds its global optimum; the intermediary solutions are not usually usable, as illustrated in Figure 10d. A potential future direction is to use algorithms such as Bayesian optimization that can explore multiple local optima and return many solutions for users to make a choice.