## Abstract

Today most origami crease patterns used in technical applications are selected from a handful of well-known origami principles. Computational algorithms capable of generating novel crease patterns either target artistic origami, focus on quadrilateral creased paper, or do not incorporate direct knowledge for the purposeful design of crease patterns tailored to engineering applications. The lack of computational methods for the generative design of crease patterns for engineering applications arises from a multitude of geometric complexities intrinsic to origami, such as rigid foldability and rigid body modes (RBMs), many of which have been addressed by recent work of the authors. Based on these findings, in this paper we introduce a Computational Design Synthesis (CDS) method for the generative design of novel crease patterns to develop origami concepts for engineering applications. The proposed method first generates crease pattern graphs through a graph grammar that automatically builds the kinematic model of the underlying origami and introduces constraints for rigid foldability. Then, the method enumerates all design alternatives that arise from the assignment of different rigid body modes to the internal vertices. These design alternatives are then automatically optimized and checked for intersection to satisfy the given design task. The proposed method is generic and applied here to two design tasks that are a rigidly foldable gripper and a rigidly foldable robotic arm.

## 1 Introduction

Origami has received considerable attention in recent decades when mathematicians, engineers, and artists recognized the benefits of origami. These benefits are numerous: origami is scale-independent [1], allows for a compact or flat stowage in the unfolded state as well as complex three-dimensional motion during deployment [2]. The deployment can be achieved by a low number of degrees-of-freedom (DOF), which minimizes the amount of resources required to actuate the folding motion and enables a reliable control [3]. The facilitated actuation and the complex motion then enable programmable structures that can change mechanical properties, shape, and function on demand [4,5]. In addition, origami can be produced in the flat state using additive manufacturing [6], which enhances realization possibilities and reduces assembly time.

Due to its benefits, the principle of origami offers potential for various scientific fields that pose different engineering design tasks. The transformation of origami into engineering applications includes an entire spectrum of more abstract to more direct implementations. The abstract end of the spectrum corresponds to so-called origami-inspired products [7] that fold and exhibit origami-like geometry but otherwise show little resemblance to origami. More closely related to traditional origami are origami-adapted mechanisms [8] that are based on origami crease patterns but use nonpaperlike materials and accommodate for finite thickness. Finally, origami-applied systems [7] use paperlike material and exhibit little to no alteration of the underlying crease pattern. Together, these categories constitute the entire spectrum of origami-based design [9,10], which finds applications in mathematics [11], material science [1215], DNA [16] and biomedical research [1721], mechanical engineering [22], robotics [2326], consumer goods [27], architecture [28], and space [2931].

The widespread applicability of origami is accompanied by the adaptation and introduction of various computational tools [10,3238] that facilitate the application of origami in engineering design tasks. However, despite the trend toward computer-aided processes, today most crease patterns used in technical applications are selected from a handful of well-known origami principles [2]. Computational algorithms capable of generating novel crease patterns either target artistic origami [39,40] or shape matching [41] without guaranteeing rigid foldability, are limited to quadrilateral creased paper [42] using only one pattern [33] and generation of origami tessellations with discrete sector angles [43]. Furthermore, most of the methods do not embed direct design knowledge on how to generate rigidly foldable crease patterns. Two examples for the latter case are works by Fuchi et al. [44] and Gillman et al. [45] that utilize ground structures that only allow for sector angle configurations that include angles of 45 deg or multiples thereof, which offers limited design freedom. This only leads to regular patterns, and thus resembles an undirected search for feasible solutions rather than the purposeful design of novel crease patterns that satisfy given engineering design tasks. In both cases, a ground structure is used as a genotype, in which each element can be switched on or off while optimizing using genetic algorithms (GA). Further approaches that rely on GA and meta-heuristics to generate solutions are works [46] and [47], respectively. Although those approaches follow encodings similar to Refs. [44,45], they are not limited to a fixed ground structure, but the question of how to perform a search that guarantees rigid foldability of generated candidate solutions is still unresolved. In an effort to improve the performance of GA and topology optimization in Ref. [45], the subsequent work [48] uses Bayesian optimization to design origami for a specific task. Rather than focusing on generating patterns, the focus is set on an efficient response evaluation since nonlinear material properties are considered. Again, a predefined ground structure is used for which the optimization process tunes the locations of vertices to satisfy the objective function.

The lack of computational methods for the generative design of crease patterns for engineering applications arises from a multitude of geometric complexities intrinsic to origami. First, many engineering applications require an origami to fold rigidly for the incorporation of rigid materials or electronics. Second, each internal vertex in a crease pattern exhibits two rigid body modes (RBMs) that allow for both “upward” and “downward” motion [49,50], which results in an exponential growth of the number of possible modes as the number of internal vertices increases. Third, the relation between kinematic determinacy, DOF, and pattern symmetry is still largely unexplored in origami. Fourth, real-world applications do not allow for self-intersection for which there still exist no intrinsic conditions. Finally, converting infinitely thin patterns into realizations with finitely thick materials imposes problems on crease pattern and hinge design [51].

In a recent work [52], whose findings will be briefly explained in Sec. 2, the authors addressed some of the above complexities by presenting the Principle of Three Units (PTU) that leads to an analytical, intrinsic condition for the rigid foldability of single degree-n vertices, where n is the number of crease lines incident to a vertex. In addition to this condition, the PTU yields an analytical kinematic model for single vertices to which the RBMs are inputs rather than the results of simulation [53]. The PTU [52] further offers guidelines for the generation of crease patterns whose kinematic motion can be explicitly determined. The condition for the rigid foldability of a single vertex can then be extended to each vertex in the generated crease pattern to achieve global rigid foldability if the underlying crease pattern is acyclic.

Based on these findings, in this paper we introduce a Computational Design Synthesis (CDS) method [54] for the generative design of novel crease patterns to develop origami concepts for engineering applications. We focus on the generation of simple crease patterns that fit the scope of engineering applications, and use the term “concept” to clarify that the method addresses the abovementioned geometric complexities except the problem of finite thickness. This leads to a method that generates rigidly foldable, infinitely thin, and intersection-free crease patterns that satisfy a given design task.

Section 2 outlines the PTU and the corresponding findings, after which the proposed CDS method is presented in Sec. 3. The method for the generative design of novel origami concepts is then applied to two engineering tasks that involve the design of rigid origami grippers and robotic arms are given in Sec. 4. Subsequently, the results are presented, the method is discussed in Sec. 5, and the contributions and findings of the paper are summarized in Sec. 6.

## 2 Background

This section describes the PTU and the condition for the rigid foldability of degree-n vertices, after which the corresponding kinematic model is summarized. The section concludes with the implications for crease pattern generation.

### 2.1 The Principle of Three Units.

The PTU [52] predicates on the fact that every single degree-n vertex requires n − 3 inputs and 3 outputs to fold in a kinematically determinable manner [55]. The inputs are the dihedral angles that drive the motion of the vertex, here called the driving angles, which need to be prescribed in order to fold the vertex. The outputs are the three remaining dihedral angles, here called the unknown dihedral angles, whose behavior is determined by the driving angles and the size of the sector angles around the vertex. As an example, Fig. 1(a) depicts a degree-8 vertex with five driving angles ρ1−4 and ρ6 assigned to an arbitrary set of crease lines, and the corresponding unknown dihedral angles ρ5, ρ7, and ρ8.

Fig. 1
Fig. 1
Close modal

Independent of the chosen set of driving angles and the degree of the vertex (for n ≥ 4), virtually cutting the vertex at the locations of the unknown dihedral angles reveals three parts [55] called unitsu1, u2, and u3 (Fig. 1(b)) [52]. What remains within these units are the sector and driving angles determined by the user, allowing for the entire kinematic behavior of each individual unit to be expressed analytically. Let a unit ui be denoted with an ordered set comprising the sector and driving angles within the unit, $ui=(αui,ρui)$, where the set of driving angles can be empty if a unit consists of a single sector. For the example in Fig. 1, this notation results in u1 = ((α1−5), (ρ1−4)), u2 = ((α6, α7), (ρ6)), and u3 = (α8).

Each unit ui spans an angle between its first and its last crease line called a unit angleUi that can be expressed analytically by alternatingly multiplying the rotation matrices of sector and driving angles and calculating the angle between the first and the last vector obtained. Unit angles change their size while folding if their respective unit contains a driving angle, such as U1(t) and U2(t) in Fig. 1(c), and they are constant if their respective unit consists of only a single sector angle, such as U3 in Fig. 1(c). In a rigidly foldable state, the surface around a vertex is continuous, i.e., without cuts, and the unit angles constitute the sides of a single spherical triangle called the vertex triangle (shown in blue in Fig. 1(c)). Conversely, if such a vertex triangle exists, then the underlying vertex is rigidly foldable. By applying the triangle inequality [56], the PTU leads to the following condition for the rigid foldability of degree-n vertices:
$Umax≤Umed+Umin$
(1)
where Umax, Umed, and Umin are the maximum, the median, and the minimum value of the unit angles in the set {U1(t), U2(t), U3(t)}, respectively, for each t during the folding motion.

### 2.2 The Kinematic Model of the Principle of Three Units.

In addition to the condition in Eq. (1), the work in Ref. [52] presents a kinematic model that derives the analytic expression for the unknown dihedral angles by applying the law of spherical cosines. Especially useful about this model is that the RBMs can be handled as an input, allowing for the design space to be enumerated based on the decision whether a vertex should fold up or downward. Based on this kinematic model, here we introduce the function
$fφ(M,u1,u2,u3)=(φ1(t),φ2(t),φ3(t))$
(2)
that takes a RBM M and the three units as inputs and returns the functions of the three unknown dihedral angles φ1(t), φ2(t), and φ3(t). The three unknown dihedral angles are adjacent to their respective units in the counter-clockwise direction as shown in Fig. 1(d). For the example given in Fig. 1, applying fφ would result in ρ5 = φ1(t), ρ7 = φ2(t), and ρ8 = φ3(t).
To be able to assemble the overall kinematics of a crease pattern once all of its vertices are modeled individually, we introduce another function
$fR(M,u1,u2,u3)=(R1,R2,R3)$
(3)
that takes the same input as fφ and returns three rotation matrices, here called the local rotation matrices $R1−3$, which correspond to the rotations at the crease lines of the unknown dihedral angles, as shown in Fig. 1(d). This calculation always starts from the first crease line of the first unit that serve as references for the x-axis and the xy-plane, respectively. For the example given in Fig. 1, this would result in the following local rotation matrices $R1=Rz(α1)Rx(ρ1)Rz(α2)⋅…⋅Rx(ρ4)Rz(α5)$, $R2=R1$$Rx(φ1)Rz(α6)Rx(ρ6)Rz(α7)$, and $R3=R2Rx(φ2)Rz(α8)$.

### 2.3 Implications for Crease Pattern Generation.

Every crease pattern is perceived as a network of vertices (nodes) that transform n − 3 inputs into 3 outputs [55], turning an origami into a directed graph where the inputs and the outputs are represented as incoming and outgoing edges, respectively. For an internal vertex vi of degree ni, its indegree (number of incoming edges) $deg−(vi)=ni−3$ linearly scales with ni and is thus unrestricted as long as ni ≥ 4. However, its outdegree $deg+(vi)$ (number of outgoing edges) must always satisfy the condition $deg+(vi)=3$. Thus, any number of incoming edges can be added to a vertex as long as the vertex has three outgoing edges. This is the first condition that needs to be satisfied while generating crease pattern graphs. The second condition stems from the kinematic model of the PTU that requires a graph to be acyclic in order to model its motion.

## 3 Method

The proposed CDS method (Fig. 2) represents crease patterns as graphs and uses a graph grammar system to generate all possible crease pattern topologies (“Topology” in Fig. 2). The inputs to the method are an initial graph G0, which serves as a starting point for the graph grammar generation of origami concepts, and a maximum number of internal vertices Nmax contained within these concepts. Both G0 and Nmax are user defined. The generation uses two different rule types, rule r1: extend vertex, and rule r2: combine vertices, within a matching process to identify the vertices to which the two rules can be applied. Simultaneously, the graph grammar automatically introduces the variables and constraints of the generated graph topologies for the optimization in subsequent steps. The generated topologies are then checked for redundancy and eliminated based on isomorphism. Further elimination criteria can be supplied manually.

Fig. 2
Fig. 2
Close modal

Subsequently, the method steps into a loop (“Geometry” in Fig. 2) in which a graph is first embedded within three dimensions by associating it with the vertex coordinates and the relevant variables to represent its geometry. Then, each graph is associated with all different RBM assignments to its internal vertices, corresponding to the enumeration of all design alternatives. Using the formulation of a design task provided by the user that involves an objective function and other optimization-related inputs, the evaluation step first optimizes each alternative and, if an optimized alternative meets the design criteria supplied with the design task, subjects it to an intersection check. A design alternative that fails this intersection check is discarded, otherwise it is added to the list of feasible origami concepts. In the following subsections (Secs. 3.13.4.2), we detail the method according to the steps in Fig. 2.

### 3.1 Representation.

A crease pattern is represented as a directed, labeled graph G = (V, E, LV, LE), where V is a nonempty set of vertices vi and E is a nonempty set of directed edges ei,j from vertices vi to vj such that ij. Vertices vi are called the predecessors of vj and vj are called the successors of vi. LV is a nonempty set of vertex labels and LE is a nonempty set of edge labels. Let ΣV be a map ΣV:VLV that labels each vertex vi to its ordered predecessors Pi and a type Ti, so that each vertex label Li in LV is an ordered pair defined as Li = (Pi, Ti). If a vertex vi has no predecessors, $Pi=(∅)$, then vi is a source vertex. Otherwise, Pi is populated with all predecessors of vi ordered in the counter-clockwise direction with respect to the crease pattern graph in the neighborhood of vi. The type Ti of vertex vi adopts the symbol χ when the vertex can be extended with new outgoing edges, or it adopts the symbol $Ti=∅$ when the vertex is either a source or has already been extended. Let ΣE be a map ΣE: ELE that assigns an edge label Li,j in LE to each directed edge ei,j, where $Li,j=(FLi,j,FRi,j)$ is an ordered pair of the facets f located on the left and on the right side of ei,j, respectively. While the sets of vertices and edges define the topology of a graph, the labels are introduced to control the generation of new graphs.

To clarify the notation, Fig. 3 shows the simplest possible initial graph G0. This graph consists of just two vertices v1 and v2 connected by a directed edge e1,2. The vertex label $L1=((∅),∅)$ denotes that the vertex v1 is a source and cannot be extended, respectively. The vertex label L2 = ((v1), χ) denotes that v1 is the single predecessor of v2 and that v2 can be extended, respectively. Finally, the two components $FL1,2=f1$ and $FR1,2=f2$ of the edge label L1,2 signify that the edge e1,2 is adjacent to the (yet invisible) facet f1 on its left and f2 on its right side.

Fig. 3
Fig. 3
Close modal

In addition to the initialization associated with the representation of the topology, the user has to relate the initial graph to its engineering purpose and thus to its geometry. This relation is achieved by allocating in-plane coordinates $xi$ to all vertices vi and defining a driving angle ρi,j for every incoming edge incident to a vertex vj that can be extended. As an example, for the graph in Fig. 3 this could be defined as $x1=(0,0)$, $x2=(1,0)$, and ρ1,2 = t where t represents a linear driving angle that goes from zero to some value tmax ∈ [−π, π]. In short, Fig. 3 thus represents an initial graph with a single extendable vertex (v2), with a geometry (vertex coordinates), and with a single DOF (driving angle ρ1,2 = t).

In general, the definition of the initial graph is an input to the method and can involve different topologies, including multiple DOF. However, the two guidelines given above also restrict the definition of possible initial graphs. First, all initial graphs must be acyclic. Second, the outdegree of all initial vertices must be smaller than 3; the method models the kinematics of a vertex while extending a vertex, and it does not extend vertices that already satisfy $deg+(vi)=3$.

### 3.2 Generation.

The origami graph grammar GG is defined by the triple $GG=(G0,R,∅)$ where G0 is the initial graph, $R=(r1,r2)$ is the set of rules containing r1 and r2, and $∅$ is the terminal symbol with respect to the type Ti of a vertex vi that prevents rule r1 from being applied to a vertex with $Ti=∅$.

The rule set $R=(r1,r2)$ is designed to conform to the two guidelines for the generation of rigidly foldable crease pattern graphs. The first rule r1 extends a vertex by new outgoing edges and incident vertices and ensures that each extended vertex has three outgoing edges after the application. The second rule r2 combines two vertices and their incoming edges to enable the generation of higher order vertices (degree-5, 6, etc.) while guaranteeing that the generated graphs stay acyclic. In addition to the graph transformations, the rules r1 and r2 embed a graph G in the plane, model extended vertices kinematically, and build or adjust the sets of the optimization variables and constraints. To automate this process, the method initializes empty sets for the in-plane coordinates x, the units u, the dihedral angles ρ, the global rotation matrices R, the optimization variables Φ, and the optimization constraints ψ. Moreover, the method initializes a set of three-dimensional vertex coordinates $X={Xi}$ where $Xi$ are the in-plane coordinates of all vertices vi contained within the initial graph G0, which are appended with a zero that represents the z-coordinate. Then, these sets are automatically filled and adjusted when the rules r1 and r2 are applied.

#### 3.2.1 Rule r1: Extend Vertex.

Rule r1 is a production of the type r1: LHS1RHS1, where the LHS1 is a single vertex va to which the rule is applied and where the RHS1 is the same vertex va with a number of new outgoing edges and incident vertices. Depending on the initial graph G0, the outdegree $deg+(va)$ on the LHS1 can be 0, 1, or 2, and r1 should accordingly generate $3−deg+(va)$ new edges and vertices. Thus, rule r1 is parametric with respect to the number of successors of va on the LHS1, which results in six possible rule application scenarios r1,1r1,6. For clarity, only the scenario r1,1 (Fig. 4) is described in detail, whereas all other rule application scenarios r1,2r1,6 are listed in a more concise form in the Supplementary Material (see Supplemental Figs. 1 and 2 available in the Supplemental Materials on the ASME Digital Collection, and Supplemental Tables 1 and 2 available on the ASME Digital Collection). Figure 4 depicts the case r1,1 in which va has no successors on the LHS1,1, where LHS1,1 and RHS1,1 are highlighted in red and the graph surrounding va in gray as a reference. This surrounding graph consists of any number of predecessors $vpm$ for 1…m ordered in the counter-clockwise order, where the respective edges are divided by the sector angles $αp1$ to $αpm−1$.

Fig. 4
Fig. 4
Close modal

#### 3.2.2 LHS Matching of r1,1.

A match $M1,1$ of the LHS1,1 is found if va is extendable as defined by its type Ta = χ and if $deg+(va)=0$ (which corresponds to this specific scenario).

#### 3.2.3 Graph Transformation of r1,1.

On the RHS1,1, rule r1,1 generates three new successor vertices $vs1−3$ numbered in counter-clockwise order that are connected to va with directed edges $ea,s1−3$. All generated vertices have identical vertex labels $Ls1−3=((va),χ)$, making them applicable for a successive rule application of r1. While the edge labels adjacent to existing facets adopt the edge labels of the corresponding edges ($FLp1,a$ and $FRpm,a$ in Fig. 4 right), two new facets are generated and denoted with $FLa,s1$ and $FLa,s2$. Hence, the edge labels result in $La,s1=(FLa,s1,FRpm,a)La,s2=$$(FLa,s2,FLa,s1)$, and $La,s3=(FLp1,a,FLa,s2)$. All newly generated indices s1−3 are in accordance with the enumeration of the nodes in G prior to the rule application. The application of r1 to a graph G, $G⇒r1G′$, results in the graph $G′=(V′=V⊎{vs1−3},E′=E⊎{ea,s1−3}$, $L′V=LV⊎{Ls1−3},L′E=LE⊎{La,s1−3})$ where $⊎$ represents the disjoint union of sets [57]. In addition, the type of the vertex va is set to $Ta′=∅$ on the RHS1,1 to prevent any further extension of va.

#### 3.2.4 Coordinates and Optimization Variables of r1,1.

To embed the graph G′ within the plane, a set of coordinates is required to describe the locations of all generated successors. Hence, all new edges $ea,s1−3$ are located with their respective sector angles $αs1−3$ (Fig. 4 right) and assigned with an edge length $ls1−3$. Using the normalized direction vector $di,j$ between vertex locations $xi$ and $xj$ of vertices vi and vj:
$di,j=xj−xi‖xj−xi‖$
(4)
the in-plane locations of the generated successors with respect to the position $xa$ of va and the predecessors Pa can be expressed as
$xs1=xa+la,s1R(αs1)da,pmxs2=xa+la,s2R(αs1+αs2)da,pmxs3=xa+la,s3R(−αs3)da,p1$
(5)
where $R$ is the two-dimensional rotation matrix. The locations of the newly generated vertices are added to the set of in-plane coordinates, $x′=x⊎{xs1−3}$.

In addition to the sector angles and edge lengths, all scenarios of r1 introduce an optimization variable Ma that stands for the RBM of the vertex va. The set of optimization variables is then expressed as $Φ′=Φ⊎{αs1−3,la,s1−3,Ma}$.

#### 3.2.5 Kinematic Modeling of r1,1.

When r1 is applied, va becomes an internal vertex of the crease pattern and can be kinematically modeled using the PTU. According to the scenario of r1,1 in Fig. 4, the units uia of va are expressed in terms of their respective sector and dihedral angles as
$u1a=((αs3,αp1,…,αpm−1,αs1),(ρp1,a,…,ρpm,a))u2a=((αs2))u3a=((2π−∑i=13αsi−∑j=1mαpj))$
(6)

If va exhibits only a single predecessor, the sum over $αpj$ is neglected. Then, the units in Eq. (6) are added to the set of units, $u′=u⊎{u1a−3a}$. By setting these units and the RBM Ma into fφ, the method determines the dihedral angles $φ1a−3a$ and associates them with the generated edges, $ρa,s1−3=φ1a−3a$, which are then added to the set of dihedral angles $ρ′=ρ⊎{ρa,s1−3}$.

The expression of units and dihedral angles enables the method to kinematically model each extended vertex individually. However, the locations of vertices need to be expressed in three-dimensional space to represent an actual origami, which requires the method to transfer the individual (local) kinematics to a global kinematic behavior. To do so, the method first finds the local rotation matrices $R1a−3a$ of va by applying $fR$. Then, these local rotation matrices are turned into global rotation matrices, $Ra,si=Rp1,aRz(π)Ria$, which are added to the set of global rotation matrices $R′=R∪{Ra,s1−3}$. A special case arises when va belongs to the initial graph since the global rotation matrix $Rp1,a$ has not been determined by an application of r1. In this case, the method sets $Rp1,a=Rz(αp)$ where αp is the angle between the edge $ep1,a$ and the x-axis. As described in Sec. 2, the first crease line of a vertex is locally associated with the x-axis, with which the method expresses the three-dimensional coordinates $Xsi$ as $Xsi=Xa+la,si*Ra,si⋅(1,0,0)T$. These coordinates are subsequently added to the set of three-dimensional vertex coordinates, $X′=X∪{Xs1−3}$.

#### 3.2.6 Optimization Constraints of r1,1.

Independent of the scenario, two types of optimization constraints are produced by r1, one defining the boundaries of the optimization variables while the other pertains to the rigid foldability. The method introduces finite lower and upper bounds for the sector angles that are greater than zero and smaller than π, respectively. As a result of the subsequent optimization, the shapes of the facets then become less sheared and less thin than their counterparts without these constraints, which facilitates the practical realization of crease patterns using finitely thick materials in the future. Equation (7) lists the optimization constraints applied to the sector angles as
$0<αmin≤αs1−3≤αmax<π∑i=13αsi+∑j=1mαpj≤2π−αmin$
(7)
where the second row relates to the sector angle of the sector corresponding to $FLa,s2$ in Fig. 4. If va exhibits only a single predecessor, the sum over $αpj$ in Eq. (7) is neglected. The same reasoning is applied to the lengths of the generated edges:
$lmin≤ls1−3≤lmax$
(8)

The second type of constraint introduced by r1 corresponds to the condition for rigid foldability in Eq. (1) that states Umax(t) ≤ Umed(t) + Umin(t). The difficulty for the application of this condition arises from the fact that the unit angles change their size over time depending on the parameter t, which effectively means that Eq. (1) has to be satisfied for the whole folding motion.

The rule system presented here always adds vertices to “the outside of” graphs, which results in only one time-dependent unit u1, whereas u2 and u3 as well as their unit angles U2 and U3 are constant, analogous to Eq. (6). A set of optimization constraints that guarantees rigid foldability for the whole folding procedure thus involves the application of Eq. (1) with both min(U1) and max(U1). In general, however, it is reasonable to assume that max(U1) occurs at t = 0, in which case Eq. (1) is satisfied because of the initial flat state. Then, the unit angle U1 decreases monotonously once t deviates from zero toward t = tmax. This assumption is fully true for degree-4 vertices [58] and depends on the sector angle configuration and the RBMs of higher order vertices [52]. Thus, the CDS method examines only the case min(U1) that occurs at t = tmax:
$Umax(tmax)≤Umed(tmax)+Umin(tmax)$
(9)

The constraints in Eqs. (7)(9) are then added to the set of optimization constraints ψ′ after each application of r1.

#### 3.2.7 Rule r2: Combine Vertices.

The second rule r2 is defined as a production of the type r2:LHS2RHS2 in Fig. 5, where LHS2 and RHS2 are both shown in red, the surrounding graph is shown in gray, and adjustments to the constraint system are shown in blue.

Fig. 5
Fig. 5
Close modal

Rule r2 accepts two vertices $vc1$ and $vc2$ that share an edge label on the LHS2 and combines them into a single vertex $vc1$ on the RHS2. The incoming edges formerly incident to $vc2$ are reassigned to $vc1$ and a new constraint is introduced between the reassigned edges and the edges incident to $vc1$. Then, r2 adjusts the set of optimization variables and constraints by replacing all instances of the lengths and sector angles corresponding to the reassigned edges, such as $lp2,c2$ and $αc2$ in Fig. 5 left, respectively, with the parts colored in blue in Fig. 5 right. The vertex vb signifies any vertex that was previously coupled to the sector angles of the reassigned edge ($αc2$).

#### 3.2.8 LHS Matching of r2.

A match $M2$ of the LHS2 in G involves four conditions. Both vertices $vc1$ and $vc2$ have to be extendable, $Tc1=Tc2=χ$. This condition guarantees the generation of acyclic graphs since no vertex can be combined with any of its predecessors.

The second condition requires that both vertices $vc1$ and $vc2$ belong to the same facet. This can be checked by comparing the edge labels of the incoming edges incident to $vc1$ and $vc2$, which is shown in Fig. 5 on the LHS2 where $FRp1,c1=FLp2,c2$ has to be satisfied. Note that while Fig. 5 illustrates the case in which both vertices $vc1$ and $vc2$ each exhibit one predecessor, multiple predecessors to both vertices are possible. In this case, the method compares the edge labels of the first and last predecessors of both vertices. Identifying these first and last predecessors is straightforward since the predecessors of a vertex are always ordered in the counter-clockwise direction as explained for rule r1.

To prevent the combination of vertices that share the same predecessor, such as $vs1−3$ in Fig. 4, the third condition states that $vc1$ and $vc2$ cannot exhibit the same predecessors, $Pc1∩Pc2={∅}$.

Since r2 deletes a vertex, the fourth condition states that $vc1$ and $vc2$ cannot both be contained in the initial graph, which can be assessed by checking if their in-plane coordinates are dependent on any optimization variables. If one vertex is contained within the initial graph, it becomes the “dominant” vertex that is not deleted, and if neither belongs to the initial graph then either vertex can be eliminated. The following graph transformation describes the generic case in which $vc1$ is dominant.

#### 3.2.9 Graph Transformation of r2.

Once r2 is applied, $G⇒r2G″$, the RHS2 includes $vc1$ while $vc2$ is subtracted from the set of vertices, $V″=V∖{vc2}$. The edge formerly incident to $vc2$ is subtracted from the set of edges and a new edge is introduced from the predecessors of $vc2$ to $vc1$, $E″=(E∖{ep2,c2})⊎{ep2,c1}$. The predecessors of both vertices are united in the correct order such that the predecessors of the vertex contributing the right side of the edge label ($FRp1,c1$ in Fig. 5) go first. In addition, vertex $vc1$ on the RHS2 can still be extended, resulting in the vertex label $Lc1″=((vp1,vp2),χ)$. The edge label $Lc2$ is then subtracted from the set of vertex labels, $LV″=LV∖{Lc2}$. The edge label formerly corresponding to the edge incident to $vc2$ is transferred identically to the reassigned edge, $Lp2,c1″=Lp2,c2$. Finally, the edge label of the edge incident to $vc2$ is subtracted from the set of edge labels, $LE″=LE∖{Lp2,c2}$. If $vc2$ is incident to multiple incoming edges on the LHS2, the above procedure is applied to all edges and edge labels.

#### 3.2.10 Coordinates and Optimization Variables of r2.

The in-plane location and thus the coordinates $xj$ of a vertex vj are always determined by a single sector angle αj and a single edge length $lpi,j$. Since r2 deletes the vertex $vc2$, the corresponding sector angle and edge lengths of $vc2$ are eliminated from the set of optimization variables, $Φ″=Φ∖{αc2,lp2,c2}$.

#### 3.2.11 Optimization Constraints of r2.

Independent of the number of predecessors in $Pc1⊎Pc2$, r2 introduces one new constraint (shown in Fig. 5 in red) on the sector angle between the edges incident to $vp1$ and $vp2$. Since this sector angle is dependent on the locations of the vertices involved, it needs to be expressed as an angle with respect to the in-plane coordinates, $∠xp1xc1xp2$. The constraint introduced by r2 applies to the upper and lower bounds equivalent to Eq. (7):
$αmin≤∠xp1xc1xp2≤αmax$
(10)

Equation (10) is then added to the set of optimization constraints ψ.

Rule r2 eliminates the sector angle and the edge lengths that formerly determined the coordinates of $vp2$, but the constraints that apply to these variables are still valid and present in the set of constraints ψ. Thus, r2 replaces all instances of these variables within ψ. Analogous to the constraint introduced above, these variables need to be expressed in terms of the in-plane coordinates (Fig. 5, blue). Variable $αc2$ is replaced with $∠xc1xp2xb$ and $lp2,c2$ is replaced with the Euclidean distance $||xc1−xp2||$ between $vp2$ and $vc1$. Hence
$ψ″=ψ(αc2←∠xc1xp2xb,lp2,c2←||xc1−xp2||)$
(11)

When $vc2$ has multiple predecessors on the LHS2, this replacement procedure needs to be analogously performed for all reassigned edges. This also applies when the sector angles and edge lengths are already expressed in terms of the in-plane coordinates, i.e., when r2 has already been applied to $vc2$.

#### 3.2.12 Automatic Graph Generation and Filtering.

The rule set $R=(r1,r2)$ is designed to comply with the guidelines for the generation of rigidly foldable and acyclic crease pattern graphs and is compact given the small number of rules required for the generation of such graphs [52]. Together, the rules r1 and r2 are able to generate existing crease patterns such as slender origami [59], origami strings [60], or the gripper in Ref. [61]. To automate the generation step (Fig. 2), the method applies the two rules whenever a respective LHS match $M1$ or $M2$ is found and enumerates all possible rule application sequences that arise from the initial graph G0 and the maximum number of internal vertices Nmax.

However, distinct rule application sequences do not guarantee distinct graphs [62], which is why the method checks the generated crease pattern graphs for redundancy and filters them before they are subjected to the guidance step. For generic design tasks, the method filters graphs based on isomorphism while filters specific to a design task have to be introduced manually if required. The filtering of isomorphic graphs is described here, while two additional filters specific to the design tasks of the gripper and the robotic arm are explained in Sec. 4.2.

For the rule system presented, two graphs are isomorphic if they exhibit an edge-preserving vertex bijection. If this condition is satisfied, the set of vertex labels LV is guaranteed to coincide, whereas the set of edge labels LE does not play a role in terms of graph morphology (i.e., there is no difference in kinematics when facets are numbered differently). Isomorphism checks are implemented in many standard programs, and the method uses the function IsomorphicGraphQ integrated in Mathematica 11.

Both the automated generation of crease pattern graphs and the isomorphism check are listed in the pseudocode in Supplemental Table 3 in the Supplementary Material available on the ASME Digital Collection. In short, the method automatically enumerates all possible crease pattern graphs that arise from the initial graph and Nmax by applying both rules r1 and r2 whenever a LHS match is found. As shown in Fig. 2, after the generation step all graphs G are forwarded to the guidance step.

### 3.3 Guidance.

The RBMs M contained in the set of the optimization variables Φ play a special role since every M only adopts the discrete states “up,” M, and “down,” M. Each internal (or extended) vertex in a graph G contributes one variable M, which results in 2N different design alternatives for a crease pattern with N internal vertices. Due to the lack of more knowledge about the search space associated with the RBMs, the method enumerates all possible RBM assignments and subjects each design alternative to the evaluation.

### 3.4 Evaluation.

As illustrated in Fig. 2, the evaluation includes both the optimization of a crease pattern as well as an intersection check. While all design alternatives are optimized, only the optimized alternatives that meet the design criteria are subjected to the intersection check.

#### 3.4.1 Optimization.

Much of what is required for the optimization is intrinsic to the rule system. By applying the rules, the method automatically generates the set of optimization variables Φ and constraints ψ of each graph G, both of which are expressed analytically. What remains to be defined by the user to run the optimization is a design task involving an objective function Ω and the numerical values for the variable boundaries lmin, lmax, αmin, and αmax. The objective function can include and should be dependent on the dihedral angles ρi,j or the three-dimensional vertex coordinates $Xi$. In practice, a user defines one or multiple states t = topt ∈ [ − π, π] at which the dihedral angles ρi,j(topt) or the vertex coordinates $Xi(topt)$ should equal certain values or locations, respectively. The definition of the objective function for the exemplary design tasks is given in Sec. 4.3.

Depending on the objective function Ω, the proposed method adjusts the set of optimization variables as a last step before the optimization. For an origami crease pattern, the lengths of the edges that are incident to degree-1 vertices on the border of the paper do not influence the kinematics of the crease pattern. Thus, all variables corresponding to the edge lengths li,j that are present neither in the objective function nor in the set of constraints ψ are discarded from the set of optimization variables Φ, resulting in an adjusted set Φopt. For the subsequent representation of the origami, the discarded edge lengths in the set of three-dimensional vertex coordinates are substituted with the value of the minimum edge length, X = X(li,jlmin).

This procedure results in the following optimization scheme to optimize the geometry of all design alternatives:
$minΦoptΩs.t.ψ$
(12)
where Ω is the objective function defined by the user, Φopt are the optimization variables containing sector angles α and crease line lengths l, and ψ are the optimization constraints consisting of variable boundaries and constraints for rigid foldability. Both Φopt and ψ are automatically generated and assigned to Eq. (12) by the method. Equation (12) is then solved by the function NMinimize integrated into Mathematica 11, which uses global numerical solvers to find a set of variables $Φopt*$ that optimizes the objective function such that $Ω(Φopt*)=Ω*$. NMinimize is applied with default settings in which the nonlinear optimization problem is first reformulated as an unconstrained problem using penalty functions and then solved with Differential Evolution [63].

#### 3.4.2 Intersection.

If an optimized design alternative satisfies $Ω*≤ϑ$, where $ϑ$ is the target criterion provided by the user in the design task, the design alternative is subjected to an intersection check that requires the method to first transform the set of optimized variables $Φopt*$ into an origami with facets. To do so, the optimized variables are first substituted into the set of three-dimensional vertex coordinates $X*=X(Φopt←Φopt*)$ that is then only dependent on t. To find the vertices that constitute a facet, the method iterates through all edges and for each edge assigns both incident vertices to the facets contained in the respective edge label. This procedure results in sets of vertices that represent the facets, and the method assigns to each of these sets a polygon. By associating each vertex vi with its three-dimensional coordinates $Xi*$, a design alternative is represented as a set of polygons that move with respect to t. Subsequently, the design alternative is subjected to the intersection check described in Ref. [64]. Since the intersection check is purely numerical, the method performs the check 16 times throughout the folding motion from t = 0 to t = tmax in even intervals. The specific number of checks is determined empirically to achieve an appropriate trade-off between a low time-consumption and a reliable intersection check.

## 4 Engineering Design Task: Origami Grippers

In this section, the proposed method is applied to two design tasks that include a gripper and a robotic arm. The principle of rigid origami has been used for gripping tasks before, such as in the Oriceps [23] or in Ref. [61], or trajectory following motions for surgical applications [20]. However, all of these existing origami-inspired designs were designed by hand. To display the usefulness of the proposed method, this work additionally introduces an obstacle that lies between the crease pattern and the point object that the final origami should be able to grip (Fig. 6(a)). In the robotic arm task, 11 points are approximated so that the tip of the arm follows a given trajectory (Fig. 6(b)), equivalent to the path generation synthesis of spatial mechanisms [65]. The remainder of the section first defines the input for both design tasks, and then in addition to the filter based on isomorphism, presents two additional filters to check for redundant graphs within the generation step. The section concludes with the optimization objectives definition, and the results produced by the method.

Fig. 6
Fig. 6
Close modal

### 4.1 Input.

The design tasks of the gripper and the robotic arm are illustrated in Figs. 6(a) and 6(b), respectively. The initial graph G0 (identical for both tasks) with vertices v1 and v2 end edge e1,2 is shown in black and corresponds to the graph described in Fig. 3 with the same vertex labels $L1=((∅),∅)$ and L2 = ((v1), χ), edge label L1,2 = (f1, f2), as well as in-plane coordinates $x1=(0,0)$ and $x2=(1,0)$. The driving angle of the edge e1,2 in both cases is linear, ρ1,2 = t, where t goes from zero to tmax = π/2, applied mirror-symmetrically with respect to the xz-plane. The actuation thus corresponds to a single DOF, and the maximum number of internal vertices for both design tasks is set to Nmax = 3. Note that graphs with one or two internal vertices are possible, Nmax is just the upper limit.

The point object PG to grip in Fig. 6(a) is located at PG = (0, 0, 1) and the obstacle is represented by a cylinder with radius 0.25 whose axis is coincident with the line segment going from (0, − 1, 0.5) to (0, 1, 0.5). This setup allows for the generation of only one side of the gripper that emerges from the initial graph shown in black, after which the arm is rotated, once optimized, by $180deg$ around the z-axis, which is illustrated in Fig. 6(a) by the graph in gray corresponding to $v2′$. Figure 6(b) shows the trajectory to approximate and the 11 points $PTi$ located at (1, 0, 0) + (2cos ti, 0, 2 sin ti) where ti goes from 0 to π/2 in steps of π/20.

To facilitate the description of a crease pattern graph, here the following notation is introduced for the rule application sequences: when a vertex vi is extended by r1, i is added to the sequence, and when vertices vi and vj are combined by r2, (i, j) is added to the sequence. As an example, extending vertices v2 and v3 in the given order and then combining vertices v4 and v8 yields a sequence denoted as 2, 3, (4, 8).

Since both design tasks illustrated in Fig. 6 are mirror symmetric with respect to the xz-plane, a graph can be filtered if it is mirror symmetric to another graph with respect to the x-axis. An example for the two symmetric graphs 2, 3 and 2, 5 is given in Fig. 7(a). These two graphs will result in the same (symmetric) optimized configuration, which renders one of the graphs redundant. In addition, any rule applied to, e.g., 2, 5 will generate a graph that is symmetric to a graph that can be generated by applying the corresponding rule to 2, 3. Thus, in both design tasks, the method filters all symmetric graphs and their descendants that result from applying more rules to symmetric graphs.

Fig. 7
Fig. 7
Close modal

The second additional filter for both problems results from the task that the origami concept must perform. One of the generated vertices contained in the graphs has to approximate the point to grip PG in the gripper task or the trajectory points $PTi$ in the robotic arm task. For simplicity, such a vertex is called the gripping vertex for both tasks. Extending a vertex by r1 only adds value to the resulting graph, not its descendants, if one of its successors becomes the gripping vertex, which is why only the successors of the vertex extended last can perform the gripping task. Figure 7(b) shows the example of a crease pattern graph 2, 4, 3 that contains unnecessary vertices and edges, here called a semantically invalid graph [66]. Since v3 is the vertex extended last, v9, v10, or v11 are the gripping vertex candidates. In this case, extending v4 is semantically invalid since all of its successors v6, v7, and v8 can be cut from the graph without changing the performance of the origami gripper. In contrast to symmetric graphs, however, the descendants of semantically invalid graphs can lead to useful crease pattern graphs, which is why their descendants are not filtered.

### 4.3 Optimization.

As explained in Sec. 4.2, all successors of the vertex extended last are candidates to grip PG or approximate $PTi$, which is why each design alterative is optimized three successive times with different gripping vertices. These successors of the vertex extended last are denoted as vl,j and the respective three-dimensional coordinates are denoted as $Xl,j(t)$ with j = 1, 2, 3. Thus, the objective function for each optimization in the gripper design task is the Euclidean distance between PG and the location of the gripping vertex at topt = tmax = π/2:
$Ω=‖PG−Xl,j(π2)‖$
(13)
An optimization of the gripper is considered successful if the objective value $Ω*$ of a design alternative is smaller than or equal to $ϑ=10−5$. In the robotic arm task, each optimization minimizes the sum of the Euclidean distance between all trajectory points $PTi$ and the respective locations of the gripping vertex at topt = ti = (0, π/20, …, 9π/20, π/2):
$Ω=∑i=111‖PTi−Xl,j(ti)‖$
(14)

An optimization of the robotic arm is considered successful if the objective value $Ω*$ of a design alternative is smaller than or equal to $ϑ=3*10−1$. The variable boundaries for both tasks are lmin = 0.1, lmax = 1.5, where αmin = π/18 and αmin = 17π/18. The length boundaries prevent the generation of heavily distorted faces, whereas the sector angle boundaries contribute prevent sickle-like patterns or trivial folds for α = π, respectively.

### 4.4 Results.

With the initial graph G0 in Fig. 6, Nmax = 3, and rules r1 and r2, the graph grammar GG generates a total of 291 crease pattern graphs for both design tasks. After filtering isomorphic, symmetric, and semantically invalid graphs, 52 of the original 291 graphs remain. The design space of all possible crease pattern topologies and rule application sequences for both design tasks is visualized in the search tree in Fig. 8. The filled nodes illustrate the 52 meaningful graphs generated, whereas both the symmetric (sym.) as well as the semantically invalid (inv.) graphs are illustrated using white nodes. Each level of the search tree contains the graphs that exhibit the same number of applied rules starting from one application at the top to seven applications at the bottom.

Fig. 8
Fig. 8
Close modal

From the 52 distinct graphs emerge 1170 distinct design alternatives (and as many objective function evaluations) that comprised the different assignments of RBMs and gripping vertices. While the topology generation of graphs is instantaneous, the optimization of a single design and the intersection check of a successful design takes from a few seconds to minutes depending on the complexity of the graph. The graphs corresponding to the first three levels in Fig. 8 each take between 0.5 and 30 s to optimize, graphs generated by four rule applications take about 10–60 s, and graphs with more rule applications show an increasing tendency where the maximum time recorded is 3 min for a single optimization. The intersection check averages at around 10 s for each design alternative. The total runtime involving the enumeration of the design space, the optimization within it, as well as the intersection check, takes 35 h for the gripper task and 40 h for the robotic arm task on an Intel i7 processor with 8 GB RAM.

With respect to the gripper, 836 of the 1170 design alternatives are able to grip the target point PG, and thereof, 148 alternatives do not self-intersect. When the cylinder is present, the number of feasible origami concepts reduces to 36. The gripping motion of three such feasible origami concepts is depicted in Fig. 9 at discrete folding states from left to right. Figure 9(a) shows a gripper concept with a thin stem and big, arrow-like arms whose motion runs along the axis of the cylinder. The gripper in Fig. 9(b) occupies little space because most of its surface in the flat state lies beneath the cylinder, and its gripping motion is angled at approximately $45deg$ to the axis of the cylinder. The gripper in Fig. 9(c) is larger than the two other concepts, grips perpendicularly to the cylinder, and contains a degree-5 vertex. The details involving the optimized vertex locations for the designs in Fig. 9 can be found in the Supplementary Material (see Supplemental Fig. 3 available in the Supplemental Materials on the ASME Digital Collection and Supplemental Tables 4 and 5 available on the ASME Digital Collection).

Fig. 9
Fig. 9
Close modal

With respect to the robotic arm task, 56 of the 1170 design alternatives are able to successfully approximate the trajectory points $PTi$, and thereof, 29 alternatives do not self-intersect. The folding motion of three design alternatives is depicted in Fig. 10 from left to right together with the trajectory points $PTi$. The three design alternatives exhibit an objective value of $Ω*=2*10−2$ (a), $Ω*=3*10−2$ (b), and $Ω*=4*10−2$ (c). The details involving the optimized vertex locations for the designs in Fig. 10 can be found in the Supplementary Material (see Supplemental Fig. 3 available in the Supplemental Materials on the ASME Digital Collection and Supplemental Tables 4 and 5 available on the ASME Digital Collection).

Fig. 10
Fig. 10
Close modal

## 5 Discussion

In comparison to related numerical approaches [67] that take multiple hours for the optimization of a single design alternative, the graph-based representation offered in this work enables the optimization of a design alternative in seconds to minutes. This is primarily due to the fact that the presented method incorporates the design knowledge required to conceptualize origami-adapted structures, involving the condition for rigid foldability, the kinematic model that determines the motion of all vertices depending on their RBM assignments, as well as the corresponding rule sets. The formalized knowledge allows to formulate the search space of crease pattern topologies, and subsequently exhaustively searches that space to generate rigidly foldable solution candidates. Such a generic approach to search space generation and exploration is not found in related work, e.g., [37,38,45], even including methods [40,41] that reside to graph-theoretical approaches for resolving some of the complexity of origami crease pattern generation. Due to the generality of the presented method, the developed concepts of the origami-adapted structures can satisfy three-dimensional design tasks, fold rigidly with the predefined number of DOF, and are free of self-intersection.

Figure 11(a) plots the distribution of generated design alternatives over the number of rule applications for the gripper design task. The distribution of the 836 designs that are able to grip PG peaks at (A) with 311 designs and four rule applications, in contrast to the total distribution of the 1180 designs having a peak with 432 designs and five rule applications at (B). Figure 8 shows that at four rule applications, all graphs except one are generated by three applications of r1 and one application of r2, while only r2 can be applied thereafter. Rule r2 introduces constraints into the system without adding variables, which is why the number of designs that can grip PG steadily decreases from four applications onward. For the same reason, the percentage of successful gripping designs per total number of designs is highest at three rule applications (C divided by D), where five of seven graphs are generated purely by r1 (Fig. 8, third level of the search tree). The distributions for the 148 designs without self-intersection and the 36 designs without any intersection peak at points (E) and (F) with 79 and 14 designs, respectively, showing a behavior similar to the distribution of the 836 designs that are able to grip PG.

Fig. 11
Fig. 11
Close modal

When at first one could have argued not to apply r2 since it only introduces constraints, the increase in feasible designs from three to four rule applications show that a certain number of applications of r2 can drive the optimized designs out of intersecting configurations. Without specific constraints that prevent intersection within the optimization, r2 thus provides more variability in the crease pattern graphs that satisfy the design task. However, if r2 is applied too many times, the designs are either unable to grip PG or lead to intersection, which leads to zero feasible designs at six (G) and seven rule applications for the gripper task.

Figure 11(b) refers to the robotic arm task and shows the distribution of the 56 successful designs that are able to approximate the given trajectory and the 29 designs that in addition do not self-intersect. The peaks appear at four rule applications (H and J), which again shows the usefulness of r2. In contrast to the results of the gripper, for the robotic arm task there are feasible concepts generated by six rule applications (K), which stems from the problem formulation (Fig. 10). Although the given trajectory is a perfect quarter circle that seems straightforward to approximate, the maximum crease line length lmax = 1.5 prevents a simple solution, allowing for complex concepts to be generated. The target criteria for both the gripper and the robotic arm task listed in Sec. 4.3 are determined empirically and would slightly change the plots in Fig. 11 if chosen differently. However, since the method uses a global optimization, the only influence of the target criteria on the results are the number of feasible design alternatives, leaving a choice for the user between precision and variability of the design alternatives.

Applied to both the gripper and the robotic arm design task, the enumeration of the design space yields 52 distinct crease pattern graph topologies depicted in Fig. 8. While the search tree is narrow at the top where r1 predominates, the midsection of the search tree becomes broader since more extended vertices offer more possibilities for the application of r2. The number of possible rule applications then progressively decreases with each line because of the maximum number of extended vertices Nmax = 3, which is why the search tree narrows again at the bottom. However, even with a small increase in Nmax to 8 vertices, for example, the maximal number of RBMs for only one particular topology is 28 making the enumeration of all crease patterns for their respective RBMs impractical. By adding to the search space all variants of crease pattern graphs that can be generated for certain limiting Nmax, the enumeration influences the scalability of the method significantly and blocks more complex engineering applications. For example, for Nmax = 8 the design space yields 3123 distinct graphs and 147,078 distinct design alternatives that needed to be optimized in the subsequent step. To improve the scalability of the method in the future, a more informed rule application should be performed. In both design tasks, there is no statistical difference in the performance of design alternatives with respect to the assignment of RBMs; both modes M and M are equally represented across the range. Considering the search method, a transfer from the design space enumeration to more efficient search methods, such as branch-and-bound algorithms [68], is required to expand the application of the method to larger design spaces.

Another interesting path for future work involves the expansion of the presented graph grammar. The rule set $R=(r1,r2)$ enables the generative design of a variety of novel crease patterns such as the ones shown in Figs. 9 and 10, as well as existing crease patterns such as slender origami [59], origami strings [60], or the gripper in Ref. [61]. However, many known crease patterns make use of specific sector angle configurations and mode assignments that fold rigidly although they cannot be obtained by using the graph grammar in this work. To generate such crease patterns, the rule system presented could be expanded by an adjusted version of rule r2, here denoted as $r2′$. Rule $r2′$ would also combine two vertices, but the vertex $vc1′$ on the RHS2 in Fig. 5 would lie on the line segment between the predecessors $vp1$ and $vp2$ and would become a sink vertex ($deg+(vc1′)=0$) of the type $Tc1′=∅$. Figure 12 illustrates how a Miura-ori pattern could be generated by the expansion of the presented rule system with rule $r2′$.

Fig. 12
Fig. 12
Close modal

The crease pattern graph depicted in Fig. 12(a) can be generated by the initial graph G0 in Fig. 3 with the rule application sequence 2, 3, 4, 8. By applying $r2′$ to vertices v9 and v14, resulting in the rule application sequence 2, 3, 4, 8, (9, 14)′ and the crease pattern graph shown in Fig. 12(b) that depicts a Miura-ori unit cell. The generation of the pattern in Fig. 12(b) would require a condition for the dihedral angles, $ρ8,9′=ρ4,9′$, and thus demand a condition for the sector angle configuration of the entire crease pattern. Such a rigidly foldable crease pattern can only be achieved by introducing some sort of symmetry [69] or by complying with the fold angle multipliers [70]. Introducing such conditions within the constraint system generated by the graph grammar is complex and requires future work.

Finally, the concepts developed in this paper are of zero thickness and applying them in real-life scenarios will require adaptations toward finitely thick materials. Adapting zero thickness concepts works best for crease patterns that exhibit only degree-4 vertices and uniform sector angle configurations. To achieve the latter, the variable boundaries αmin and αmax in Eqs. (7) and (10) can be defined to lie within a small interval. However, introducing such intervals might deplete the number of feasible design concepts, and the user of the method has to estimate the trade-off between the number and the realizability of these solutions. To achieve the former, users can select the feasible designs generated only by applications of rule r1 or a minimum number of applications of rule r2. Nonetheless, adapting the feasible concepts to finite thickness still requires future work, especially if the adaptations were to be integrated into the method itself. A range of possible adaptation techniques is given in Ref. [51].

## 6 Conclusion

In this paper, we introduced a generative method for the generation of novel origami concepts for engineering design tasks. The method includes an origami graph grammar that generates rigidly foldable crease patterns as well as the constraints for both the boundaries of the optimization variables based on the conditions for rigid foldability. Although the method generates rather simple patterns, the enumeration of all possible design alternatives based on the RBMs to the internal vertices shows a surprising richness of the underlying search space resulting in a considerable number of design concepts that have been generated. Finally, the paper demonstrates that the PTU enables the efficient optimization of design alternatives even including the self-intersection check. The future work will involve a development of a more efficient PTU-based generative method that does not rely on the exhaustive search and is able to generate origami patterns of increased complexity, involving origami strips, tessellations, and nondevelopable patterns.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper.

## Nomenclature

• f =

facet

•
• l =

crease line length

•
• n =

degree of a vertex

•
• r =

graph grammar rule

•
• t =

time parameter

•
• u =

set of units

•
• v =

vertex

•
• x =

set of in-plane vertex coordinates

•
• $d$ =

normalized direction vector

•
• $x$ =

in-plane vertex coordinates

•
• E =

set of edges

•
• G =

graph

•
• M =

generic rigid body mode

•
• N =

number of internal vertices

•
• R =

set of global rotation matrices

•
• V =

set of vertices

•
• X =

set of three-dimensional vertex coordinates

•
• $X$ =

three-dimensional vertex coordinates

•
• $G$ =

set of graphs

•
• $L$ =

language of a grammar

•
• $M$ =

left-hand-side match for graph grammar rules

•
• $R$ =

set of graph grammar rules

•
• ei,j =

edge from vertex vi to vertex vj

•
• ui =

unit, i = 1, 2, 3

•
• $fR$ =

function determining local rotation matrices

•
• fφ =

function determining unknown dihedral angles

•
• Li =

vertex label of vertex vi

•
• Li,j =

edge label of edge ei,j

•
• LE =

set of edge labels

•
• LV =

set of vertex labels

•
• Pi =

predecessors of a vertex vi

•
• Ti =

type of a vertex

•
• Ui =

unit angle, i = 1, 2, 3

•
• $Ri$ =

local rotation matrix, i = 1, 2, 3

•
• $Ri,j$ =

global rotation matrix of the edge ei,j

•
• $FLi,j,FRi,j$ =

facets on the left and right of edge ei,j, respectively

•
• M, M =

rigid body modes up and down, respectively

•
• PG, PT =

point to grip and points of trajectory, respectively

•
• α =

sector angle

•
• ρ =

set of dihedral angles

•
• ρi,j =

dihedral angles

•
• $ϑ$ =

target criterion

•
• ΣE, ΣV =

maps for edge and vertex labels

•
• φi =

unknown dihedral angle, i = 1, 2, 3

•
• Φ =

set of optimization variables

•
• χ =

symbol for type T meaning vertex can be extended

•
• ψ =

set of optimization constraints

•
• Ω =

objective function

•
• $∅$ =

empty and terminal symbol

•
• $⊎$ =

disjoint union of sets

## References

1.
Wang
,
Z.
,
Jing
,
L.
,
Yao
,
K.
,
Yang
,
Y.
,
Zheng
,
B.
,
Soukoulis
,
C. M.
,
Chen
,
H.
, and
Liu
,
Y.
,
2017
, “
Origami-Based Reconfigurable Metamaterials for Tunable Chirality
,”
,
29
(
27
), pp.
1700412
.
2.
Morgan
,
J.
,
Magleby
,
S. P.
,
Lang
,
R. J.
, and
Howell
,
L. L.
,
2015
, “
A Preliminary Process for Origami-Adapted Design
,”
ASME 2015 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference
,
Boston, MA
,
Aug. 2–5
, p. V05BT08A053.
3.
Stern
,
M.
,
Pinson
,
M. B.
, and
Murugan
,
A.
,
2017
, “
The Complexity of Folding Self-Folding Origami
,”
Phys. Rev. X
,
7
(
4
), p.
041070
.
4.
Hawkes
,
E.
,
An
,
B.
,
Benbernou
,
N. M.
,
Tanaka
,
H.
,
Kim
,
S.
,
Demaine
,
E. D.
,
Rus
,
D.
, and
Wood
,
R. J.
,
2010
, “
Programmable Matter by Folding
,”
Proc. Natl. Acad. Sci. U. S. A.
,
107
(
28
), pp.
12441
12445
.
5.
Martinez
,
R. V.
,
Fish
,
C. R.
,
Chen
,
X.
, and
Whitesides
,
G. M.
,
2012
, “
Elastomeric Origami: Programmable Paper-Elastomer Composites as Pneumatic Actuators
,”
,
22
(
7
), pp.
1376
1384
.
6.
Chen
,
T.
,
Bilal
,
O. R.
,
Lang
,
R. J.
,
Daraio
,
C.
, and
Shea
,
K.
,
2019
, “
Autonomous Deployment of a Solar Panel Using Elastic Origami and Distributed Shape-Memory-Polymer Actuators
,”
Phys. Rev. Appl.
,
11
(
6
), p.
064069
.
7.
Francis
,
K. C.
,
Rupert
,
L. T.
,
Lang
,
R. J.
,
Morgan
,
D. C.
,
Magleby
,
S. P.
, and
Howell
,
L. L.
,
2015
, “
From Crease Pattern to Product: Considerations to Engineering Origami-Adapted Designs
,”
ASME 2014 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference
,
Buffalo, NY
,
Aug. 17–20
, p. V05BT08A030.
8.
Zirbel
,
S. A.
,
Lang
,
R. J.
,
Thomson
,
M. W.
,
Sigel
,
D. A.
,
Walkemeyer
,
P. E.
,
Trease
,
B. P.
,
Magleby
,
S. P.
, and
Howell
,
L. L.
,
2013
, “
Accommodating Thickness in Origami-Based Deployable Arrays
,”
ASME J. Mech. Des.
,
135
(
11
), pp.
111005
.
9.
Morgan
,
J.
,
Magleby
,
S. P.
, and
Howell
,
L. L.
,
2016
, “
An Approach to Designing Origami-Adapted Aerospace Mechanisms
,”
ASME J. Mech. Des.
,
138
(
5
), p.
052301
.
10.
Meloni
,
M.
,
Cai
,
J.
,
Zhang
,
Q.
,
Sang-Hoon Lee
,
D.
,
Li
,
M.
,
Ma
,
R.
,
Parashkevov
,
T. E.
, and
Feng
,
J.
,
2021
, “
Engineering Origami: A Comprehensive Review of Recent Applications, Design Methods, and Tools
,”
,
8
(
13
), pp.
2000636
.
11.
Demaine
,
E. D.
,
Demaine
,
M. L.
, and
Mitchell
,
J. S. B.
,
2000
, “
Folding Flat Silhouettes and Wrapping Polyhedral Packages: New Results in Computational Origami
,”
Comput. Geom.
,
16
(
1
), pp.
3
21
.
12.
Liu
,
B.
,
Silverberg
,
J. L.
,
Evans
,
A. A.
,
Santangelo
,
C. D.
,
Lang
,
R. J.
,
Hull
,
T. C.
, and
Cohen
,
I.
,
2018
, “
Topological Kinematics of Origami Metamaterials
,”
Nat. Phys.
,
14
(
8
), pp.
811
815
.
13.
Schenk
,
M.
, and
Guest
,
S. D.
,
2013
, “
Geometry of Miura-Folded Metamaterials
,”
Proc. Natl. Acad. Sci. U. S. A.
,
110
(
9
), pp.
3276
3281
.
14.
,
T.
,
Ma
,
J.
,
Feng
,
H.
,
Hou
,
D.
,
Gattas
,
J. M.
,
Chen
,
Y.
, and
You
,
Z.
,
2020
, “
Programmable Stiffness and Shape Modulation in Origami Materials: Emergence of a Distant Actuation Feature
,”
Appl. Mater. Today
,
19
, pp.
100537
.
15.
,
A. A.
,
2016
, “
Mechanical Meta-Materials
,”
Mater. Horiz.
,
3
(
5
), pp.
371
381
.
16.
Han
,
D.
,
Pal
,
S.
,
Nangreave
,
J.
,
Deng
,
Z.
,
Liu
,
Y.
, and
Yan
,
H.
,
2011
, “
DNA Origami With Complex Curvatures in Three-Dimensional Space
,”
Science
,
332
(
6027
), pp.
342
346
.
17.
Kuribayashi
,
K.
,
Tsuchiya
,
K.
,
You
,
Z.
,
Tomus
,
D.
,
Umemoto
,
M.
,
Ito
,
T.
, and
Sasaki
,
M.
,
2006
, “
Self-Deployable Origami Stent Grafts as a Biomedical Application of Ni-Rich TiNi Shape Memory Alloy Foil
,”
Mater. Sci. Eng. A
,
419
(
1–2
), pp.
131
137
.
18.
Bobbert
,
F. S. L.
,
Janbaz
,
S.
,
van Manen
,
T.
,
Li
,
Y.
, and
,
A. A.
,
2020
, “
Russian Doll Deployable Meta-Implants: Fusion of Kirigami, Origami, and Multi-Stability
,”
Mater. Des.
,
191
, p.
108624
.
19.
Sargent
,
B.
,
Butler
,
J.
,
Seymour
,
K.
,
Bailey
,
D.
,
Jensen
,
B.
,
Magleby
,
S.
, and
Howell
,
L.
,
2020
, “
An Origami-Based Medical Support System to Mitigate Flexible Shaft Buckling
,”
ASME J. Mech. Rob.
,
12
(
4
), p.
041005
.
20.
Suzuki
,
H.
, and
Wood
,
R. J.
,
2020
, “
Origami-Inspired Miniature Manipulator for Teleoperated Microsurgery
,”
Nat. Mach. Intell.
,
2
(
8
), pp.
437
446
.
21.
Banerjee
,
H.
,
Li
,
T. K.
,
Ponraj
,
G.
,
Kirthika
,
S. K.
,
Lim
,
C. M.
, and
Ren
,
H.
,
2020
, “
Origami-Layer-Jamming Deployable Surgical Retractor With Variable Stiffness and Tactile Sensing
,”
ASME J. Mech. Rob.
,
12
(
3
), p.
031010
.
22.
Peraza-Hernandez
,
E. A.
,
Hartl
,
D. J.
,
Malak
,
R. J.
, Jr.
, and
Lagoudas
,
D. C.
,
2014
, “
Origami-Inspired Active Structures: A Synthesis and Review
,”
Smart Mater. Struct.
,
23
(
9
), p.
094001
.
23.
Edmondson
,
B. J.
,
Bowen
,
L. A.
,
Grames
,
C. L.
,
Magleby
,
S. P.
,
Howell
,
L. L.
, and
Bateman
,
T. C.
,
2013
, “
Oriceps: Origami-Inspired Forceps
,”
ASME 2013 Conference on Smart Materials, Adaptive Structures and Intelligent Systems
,
Snowbird, UT
,
Sept. 16–18
, p. V001T01A027.
24.
Felton
,
S. M.
,
Tolley
,
M.
,
Demaine
,
E. D.
,
Rus
,
D.
, and
Wood
,
R. J.
,
2014
, “
A Method for Building Self-Folding Machines
,”
Science
,
345
(
6197
), pp.
644
646
.
25.
Mintchev
,
S.
,
Daler
,
L.
,
L'Eplattenier
,
G.
,
Saint-Raymond
,
L.
, and
Floreano
,
D.
,
2015
, “
Foldable and Self-Deployable Pocket Sized Quadrotor
,”
2015 IEEE International Conference on Robotics and Automation (ICRA)
,
Seattle, WA
,
May 25–30
, IEEE, pp.
2190
2195
.
26.
Vander Hoff
,
E.
,
Jeong
,
D.
, and
Lee
,
K.
,
2014
, “
OrigamiBot-I: A Thread-Actuated Origami Robot for Manipulation and Locomotion
,”
2014 IEEE/RSJ International Conference on Intelligent Robots and Systems
,
Chicago, IL
,
Sept. 14–18
, IEEE, pp.
1421
1426
.
27.
Morris
,
E.
,
,
D. A.
, and
Malak
,
R.
,
2016
, “
The State of the Art of Origami-Inspired Products: A Review
,”
ASME 2016 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference
,
Charlotte, NC
,
Aug. 21–24
, p. V05BT07A014.
28.
Tachi
,
T.
,
2010
, “
Geometric Considerations for the Design of Rigid Origami Structures
,”
Proceedings of the International Association for Shell and Spatial Structures (IASS) Symposium
,
Shanghai, China
,
Nov. 8–12
, Vol. 12, No. 10, pp.
458
460
.
29.
Hanaor
,
A.
, and
Levy
,
R.
,
2001
, “
Evaluation of Deployable Structures for Space Enclosures
,”
Int. J. Space Struct.
,
16
(
4
), pp.
211
229
.
30.
Natori
,
M. C.
,
Katsumata
,
N.
,
Yamakawa
,
H.
,
Sakamoto
,
H.
, and
Kishimoto
,
N.
,
2013
, “
Conceptual Model Study Using Origami for Membrane Space Structures
,”
International Design Engineering Technical Conferences and Computers and Information in Engineering Conference
,
Portland, OR
,
Aug. 4–7
, p. V06BT07A047.
31.
Wilson
,
L.
,
Pellegrino
,
S.
, and
Danner
,
R.
,
2013
, “
Origami Sunshield Concepts for Space Telescopes
,”
54th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference
,
Boston, MA
,
Apr. 8–11
, p.
1594
.
32.
An
,
B.
,
Benbernou
,
N.
,
Demaine
,
E. D.
, and
Rus
,
D.
,
2011
, “
Planning to Fold Multiple Objects From a Single Self-Folding Sheet
,”
Robotica
,
29
(
1
), pp.
87
102
.
33.
Dudte
,
L. H.
,
Vouga
,
E.
,
Tachi
,
T.
, and
,
L.
,
2016
, “
Programming Curvature Using Origami Tessellations
,”
Nat. Mater.
,
15
(
5
), p.
583
.
34.
Tachi
,
T.
,
2009
, “
Simulation of Rigid Origami
,”
The Origami 4: Fourth International Meeting of Origami Science, Mathematics, and Education
,
,
Sept. 8–10
, pp.
175
188
.
35.
Filipov
,
E. T.
,
Liu
,
K.
,
Tachi
,
T.
,
Schenk
,
M.
, and
Paulino
,
G. H.
,
2017
, “
Bar and Hinge Models for Scalable Analysis of Origami
,”
Int. J. Solids Struct.
,
124
, pp.
26
45
.
36.
Kwok
,
T.-H.
,
Wang
,
C. C.
,
Deng
,
D.
,
Zhang
,
Y.
, and
Chen
,
Y.
,
2015
, “
Four-Dimensional Printing for Freeform Surfaces: Design Optimization of Origami and Kirigami Structures
,”
ASME J. Mech. Des.
,
137
(
11
), p.
111413
.
37.
Schreck
,
C.
,
Rohmer
,
D.
,
Hahmann
,
S.
,
Cani
,
M.-P.
,
Jin
,
S.
,
Wang
,
C. C. L.
, and
Bloch
,
J.-F.
,
2015
, “
Nonsmooth Developable Geometry for Interactively Animating Paper Crumpling
,”
ACM Trans. Graph.
,
35
(
1
), pp.
1
18
.
38.
Woodruff
,
S. R.
, and
Filipov
,
E. T.
,
2020
, “
A Bar and Hinge Model Formulation for Structural Analysis of Curved-Crease Origami
,”
Int. J. Solids Struct.
,
204–205
, pp.
114
127
.
39.
Lang
,
R. J.
,
1996
, “
A Computational Algorithm for Origami Design
,”
Proceedings of the Twelfth Annual Symposium on Computational Geometry
,
ACM
, pp.
98
105
.
40.
Tachi
,
T.
,
2009
, “
Origamizing Polyhedral Surfaces
,”
IEEE Trans. Vis. Comput. Graph.
,
16
(
2
), pp.
298
311
.
41.
Demaine
,
E. D.
, and
Tachi
,
T.
,
2017
, “
Origamizer: A Practical Algorithm for Folding Any Polyhedron
,”
33rd International Symposium on Computational Geometry (SoCG 2017). Schloss Dagstuhl-Leibniz-Zentrum Fuer Informatik
,
Brisbane, Australia
,
July 4–7
, Vol. 77, pp.
1
34
.
42.
Lang
,
R. J.
, and
Howell
,
L. L.
,
2018
, “
Rigidly Foldable Quadrilateral Meshes From Angle Arrays
,”
ASME J. Mech. Rob.
,
10
(
2
), p.
021004
.
43.
Dieleman
,
P.
,
Vasmel
,
N.
,
Waitukaitis
,
S.
, and
van Hecke
,
M.
,
2019
, “
Jigsaw Puzzle Design of Pluripotent Origami
,”
Nat. Phys.
,
16
(
1
), pp.
63
68
.
44.
Fuchi
,
K.
,
Buskohl
,
P. R.
,
Bazzan
,
G.
,
Durstock
,
M. F.
,
Reich
,
G. W.
,
Vaia
,
R. A.
, and
Joo
,
J. J.
,
2015
, “
Origami Actuator Design and Networking Through Crease Topology Optimization
,”
ASME J. Mech. Des.
,
137
(
9
), p.
091401
.
45.
Gillman
,
A. S.
,
Fuchi
,
K.
, and
Buskohl
,
P. R.
,
2019
, “
Discovering Sequenced Origami Folding Through Nonlinear Mechanics and Topology Optimization
,”
ASME J. Mech. Des.
,
141
(
4
), p.
041401
.
46.
,
D. A.
, and
Li
,
W.
,
2014
, “
A Novel Method to Design and Optimize Flat-Foldable Origami Structures Through a Genetic Algorithm
,”
ASME J. Comput. Inf. Sci. Eng.
,
14
(
3
), p.
031008
.
47.
Li
,
W.
, and
,
D. A.
,
2015
, “
Designing Optimal Origami Structures by Computational Evolutionary Embryogeny
,”
ASME J. Comput. Inf. Sci. Eng.
,
15
(
1
), p.
011010
.
48.
Shende
,
S.
,
Gillman
,
A.
,
Yoo
,
D.
,
Buskohl
,
P.
, and
Vemaganti
,
K.
,
2021
, “
Bayesian Topology Optimization for Efficient Design of Origami Folding Structures
,”
Struct. Multidiscipl. Optim.
,
63
(
4
), pp.
1907
1926
.
49.
Abel
,
Z.
,
Cantarella
,
J.
,
Demaine
,
E. D.
,
Eppstein
,
D.
,
Hull
,
T. C.
,
Ku
,
J. S.
,
Lang
,
J. L.
, and
Tachi
,
T.
,
2015
, “
Rigid Origami Vertices: Conditions and Forcing Sets
,”
Journal of Computational Geometry
,
7
(
1
), pp.
171
184
.
50.
Waitukaitis
,
S.
,
Menaut
,
R.
,
Chen
,
B. G.
, and
van Hecke
,
M.
,
2015
, “
Origami Multistability: From Single Vertices to Metasheets
,”
Phys. Rev. Lett.
,
114
(
5
), p.
055503
.
51.
Lang
,
R. J.
,
Tolman
,
K. A.
,
Crampton
,
E. B.
,
Magleby
,
S. P.
, and
Howell
,
L. L.
,
2018
, “
A Review of Thickness-Accommodation Techniques in Origami-Inspired Engineering
,”
ASME Appl. Mech. Rev.
,
70
(
1
), p.
010805
.
52.
Zimmermann
,
L.
,
Shea
,
K.
, and
Stanković
,
T.
,
2020
, “
Conditions for Rigid and Flat Foldability of Degree-n Vertices in Origami
,”
ASME J. Mech. Rob.
,
12
(
1
), p.
011020
.
53.
Zimmermann
,
L.
,
Stanković
,
T.
, and
Shea
,
K.
,
2017
, “
Finding Rigid Body Modes of Rigid-Foldable Origami Through the Simulation of Vertex Motion
,”
ASME 2017 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference
,
Cleveland, OH
,
Aug. 6–9
, p. V05BT08A046.
54.
Cagan
,
J.
,
Campbell
,
M. I.
,
Finger
,
S.
, and
Tomiyama
,
T.
,
2005
, “
A Framework for Computational Design Synthesis: Model and Applications
,”
ASME J. Comput. Inf. Sci. Eng.
,
5
(
3
), pp.
171
181
.
55.
Balkcom
,
D. J.
, and
Mason
,
M. T.
,
2008
, “
Robotic Origami Folding
,”
The International Journal of Robotics Research
,
27
(
5
), pp.
613
627
.
56.
Ramsay
,
A.
, and
Richtmyer
,
R. D.
,
1995
,
Introduction to Hyperbolic Geometry
,
Springer-Verlag
,
New York
.
57.
Kreowski
,
H.-J.
,
Klempien-Hinrichs
,
R.
, and
Kuske
,
S.
,
2006
, “
Some Essentials of Graph Transformation
,”
,
25
, pp.
229
254
.
58.
Zimmermann
,
L.
, and
Stanković
,
T.
,
2019
, “
Rigid and Flat Foldability of a Degree-Four Vertex in Origami
,”
ASME J. Mech. Rob.
,
12
(
1
), p.
011004
.
59.
Kamrava
,
S.
,
Ghosh
,
R.
,
Yang
,
Y.
, and
Vaziri
,
A.
,
2018
, “
Slender Origami With Complex 3D Folding Shapes
,”
Europhys. Lett.
,
124
(
5
), p.
58001
.
60.
Kamrava
,
S.
,
,
D.
,
Felton
,
S. M.
, and
Vaziri
,
A.
,
2018
, “
Programmable Origami Strings
,”
,
3
(
3
), p.
1700276
.
61.
Faber
,
J. A.
,
Arrieta
,
A. F.
, and
Studart
,
A. R.
,
2018
, “
Bioinspired Spring Origami
,”
Science
,
359
(
6382
), pp.
1386
1391
.
62.
Königseder
,
C.
,
Stanković
,
T.
, and
Shea
,
K.
,
2016
, “
Improving Design Grammar Development and Application Through Network-Based Analysis of Transition Graphs
,”
Des. Sci.
,
2
, pp.
1
34
.
63.
Storn
,
R.
,
1996
, “
On the Usage of Differential Evolution for Function Optimization
,”
Proceedings of North American Fuzzy Information Processing
,
Berkeley, CA
,
June 19–22
, IEEE, pp.
519
523
.
64.
Möller
,
T.
,
1997
, “
A Fast Triangle-Triangle Intersection Test
,”
J. Graph. Tools
,
2
(
2
), pp.
25
30
.
65.
Jiménez
,
J.
,
Alvarez
,
G.
,
Cardenal
,
J.
, and
,
J.
,
1997
, “
A Simple and General Method for Kinematic Synthesis of Spatial Mechanisms
,”
Mech. Mach. Theory
,
32
(
3
), pp.
323
341
.
66.
Ma
,
T.
,
Chen
,
J.
, and
Xiao
,
C.
,
2018
, “
Constrained Generation of Semantically Valid Graphs via Regularizing Variational Autoencoders
,”
Advances in Neural Information Processing Systems
,
,
Dec. 3–8
, pp.
7113
7124
.
67.
Zimmermann
,
L.
,
Shea
,
K.
, and
Stanković
,
T.
,
2018
, “
Origami Sensitivity—On the Influence of Vertex Geometry
,”
7th International Meeting on Origami in Science, Mathematics and Education (7OSME)
,
Oxford, UK
,
Sept. 4–7
, pp.
1087
1102
.
68.
Lawler
,
E. L.
, and
Wood
,
D. E.
,
1966
, “
Branch-and-Bound Methods: A Survey
,”
Oper. Res.
,
14
(
4
), pp.
699
719
.
69.
Tachi
,
T.
,
2010
, “
Freeform Rigid-Foldable Structure Using Bidirectionally Flat-Foldable Planar Quadrilateral Mesh
,”
Proceedings of The Conference Advances in Architectural Geometry
,
Vienna, Austria
,
Sept. 18–21
.
70.
Evans
,
T. A.
,
Lang
,
R. J.
,
Magleby
,
S. P.
, and
Howell
,
L. L.
,
2015
, “
Rigidly Foldable Origami Twists
,”
Origami 6: Sixth International Meeting of Origami Science, Mathematics, and Education
,
Tokyo, Japan
,
Aug. 10–13
, No. 1, pp.
119
130
.