## 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 [12–15], DNA [16] and biomedical research [17–21], mechanical engineering [22], robotics [23–26], consumer goods [27], architecture [28], and space [29–31].

The widespread applicability of origami is accompanied by the adaptation and introduction of various computational tools [10,32–38] 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}.

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 *units**u*_{1}, *u*_{2}, and *u*_{3} (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 *u*_{i} be denoted with an ordered set comprising the sector and driving angles within the unit, $ui=(\alpha ui,\rho 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 *u*_{1} = ((*α*_{1−5}), (*ρ*_{1−4})), *u*_{2} = ((*α*_{6}, *α*_{7}), (*ρ*_{6})), and *u*_{3} = (*α*_{8}).

*u*

_{i}spans an angle between its first and its last crease line called a

*unit angle*

*U*

_{i}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

*U*

_{1}(

*t*) and

*U*

_{2}(

*t*) in Fig. 1(c), and they are constant if their respective unit consists of only a single sector angle, such as

*U*

_{3}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:

*U*

_{max},

*U*

_{med}, and

*U*

_{min}are the maximum, the median, and the minimum value of the unit angles in the set {

*U*

_{1}(

*t*),

*U*

_{2}(

*t*),

*U*

_{3}(

*t*)}, respectively, for each

*t*during the folding motion.

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

*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*).

*f*

_{φ}and returns three rotation matrices, here called the

*local*rotation matrices $R1\u22123$, 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(\alpha 1)Rx(\rho 1)Rz(\alpha 2)\u22c5\u2026\u22c5Rx(\rho 4)Rz(\alpha 5)$, $R2=R1$$Rx(\phi 1)Rz(\alpha 6)Rx(\rho 6)Rz(\alpha 7)$, and $R3=R2Rx(\phi 2)Rz(\alpha 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 *v*_{i} of degree *n*_{i}, its indegree (number of incoming edges) $deg\u2212(vi)=ni\u22123$ linearly scales with *n*_{i} and is thus unrestricted as long as *n*_{i} ≥ 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 *G*_{0}, which serves as a starting point for the graph grammar generation of origami concepts, and a maximum number of internal vertices *N*_{max} contained within these concepts. Both *G*_{0} and *N*_{max} are user defined. The generation uses two different rule types, rule *r*_{1}: extend vertex, and rule *r*_{2}: 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.

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.1–3.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*, *L*_{V}, *L*_{E}), where *V* is a nonempty set of vertices *v*_{i} and *E* is a nonempty set of directed edges *e*_{i,j} from vertices *v*_{i} to *v*_{j} such that *i* ≠ *j*. Vertices *v*_{i} are called the *predecessors* of *v*_{j} and *v*_{j} are called the *successors* of *v*_{i}. *L*_{V} is a nonempty set of vertex labels and *L*_{E} is a nonempty set of edge labels. Let *Σ*_{V} be a map *Σ*_{V}:*V* → *L*_{V} that labels each vertex *v*_{i} to its ordered predecessors *P*_{i} and a type *T*_{i}, so that each vertex label *L*_{i} in *L*_{V} is an ordered pair defined as *L*_{i} = (*P*_{i}, *T*_{i}). If a vertex *v*_{i} has no predecessors, $Pi=(\u2205)$, then *v*_{i} is a *source* vertex. Otherwise, *P*_{i} is populated with all predecessors of *v*_{i} ordered in the counter-clockwise direction with respect to the crease pattern graph in the neighborhood of *v*_{i}. The type *T*_{i} of vertex *v*_{i} adopts the symbol *χ* when the vertex can be extended with new outgoing edges, or it adopts the symbol $Ti=\u2205$ when the vertex is either a source or has already been extended. Let *Σ*_{E} be a map *Σ*_{E}: *E* → *L*_{E} that assigns an edge label *L*_{i,j} in *L*_{E} to each directed edge *e*_{i,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 *e*_{i,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 G_{0}. This graph consists of just two vertices *v*_{1} and *v*_{2} connected by a directed edge *e*_{1,2}. The vertex label $L1=((\u2205),\u2205)$ denotes that the vertex *v*_{1} is a source and cannot be extended, respectively. The vertex label *L*_{2} = ((*v*_{1}), *χ*) denotes that *v*_{1} is the single predecessor of *v*_{2} and that *v*_{2} can be extended, respectively. Finally, the two components $FL1,2=f1$ and $FR1,2=f2$ of the edge label *L*_{1,2} signify that the edge *e*_{1,2} is adjacent to the (yet invisible) facet *f*_{1} on its left and *f*_{2} on its right side.

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 *v*_{i} and defining a driving angle *ρ*_{i,j} for every incoming edge incident to a vertex *v*_{j} 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 *t*_{max} ∈ [−π, π]. In short, Fig. 3 thus represents an initial graph with a single extendable vertex (*v*_{2}), 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,\u2205)$ where *G*_{0} is the initial graph, $R=(r1,r2)$ is the set of rules containing *r*_{1} and *r*_{2}, and $\u2205$ is the terminal symbol with respect to the type *T*_{i} of a vertex *v*_{i} that prevents rule *r*_{1} from being applied to a vertex with $Ti=\u2205$.

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 *r*_{1} 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 *r*_{2} 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 *r*_{1} and *r*_{2} 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 *v*_{i} contained within the initial graph *G*_{0}, which are appended with a zero that represents the *z*-coordinate. Then, these sets are automatically filled and adjusted when the rules *r*_{1} and *r*_{2} are applied.

#### 3.2.1 Rule **r**_{1}: Extend Vertex.

Rule *r*_{1} is a production of the type *r*_{1}: *LHS*_{1} → *RHS*_{1}, where the *LHS*_{1} is a single vertex *v*_{a} to which the rule is applied and where the *RHS*_{1} is the same vertex *v*_{a} with a number of new outgoing edges and incident vertices. Depending on the initial graph *G*_{0}, the outdegree $deg+(va)$ on the *LHS*_{1} can be 0, 1, or 2, and *r*_{1} should accordingly generate $3\u2212deg+(va)$ new edges and vertices. Thus, rule *r*_{1} is parametric with respect to the number of successors of *v*_{a} on the *LHS*_{1}, which results in six possible rule application scenarios *r*_{1,1} − *r*_{1,6}. For clarity, only the scenario *r*_{1,1} (Fig. 4) is described in detail, whereas all other rule application scenarios *r*_{1,2} − *r*_{1,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 *r*_{1,1} in which *v*_{a} has no successors on the *LHS*_{1,1}, where *LHS*_{1,1} and *RHS*_{1,1} are highlighted in red and the graph surrounding *v*_{a} 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 $\alpha p1$ to $\alpha pm\u22121$.

#### 3.2.2 LHS Matching of **r**_{1,1}.

A match $M1,1$ of the *LHS*_{1,1} is found if *v*_{a} is extendable as defined by its type *T*_{a} = *χ* and if $deg+(va)=0$ (which corresponds to this specific scenario).

#### 3.2.3 Graph Transformation of **r**_{1,1}.

On the *RHS*_{1,1}, rule *r*_{1,1} generates three new successor vertices $vs1\u22123$ numbered in counter-clockwise order that are connected to *v*_{a} with directed edges $ea,s1\u22123$. All generated vertices have identical vertex labels $Ls1\u22123=((va),\chi )$, making them applicable for a successive rule application of *r*_{1}. 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 s_{1−3} are in accordance with the enumeration of the nodes in *G* prior to the rule application. The application of *r*_{1} to a graph *G*, $G\u21d2r1G\u2032$, results in the graph $G\u2032=(V\u2032=V\u228e{vs1\u22123},E\u2032=E\u228e{ea,s1\u22123}$, $L\u2032V=LV\u228e{Ls1\u22123},L\u2032E=LE\u228e{La,s1\u22123})$ where $\u228e$ represents the disjoint union of sets [57]. In addition, the type of the vertex *v*_{a} is set to $Ta\u2032=\u2205$ on the *RHS*_{1,1} to prevent any further extension of *v*_{a}.

#### 3.2.4 Coordinates and Optimization Variables of **r**_{1,1}.

*G*′ within the plane, a set of coordinates is required to describe the locations of all generated successors. Hence, all new edges $ea,s1\u22123$ are located with their respective sector angles $\alpha s1\u22123$ (Fig. 4 right) and assigned with an edge length $ls1\u22123$. Using the normalized direction vector $di,j$ between vertex locations $xi$ and $xj$ of vertices

*v*

_{i}and

*v*

_{j}:

*v*

_{a}and the predecessors

*P*

_{a}can be expressed as

In addition to the sector angles and edge lengths, all scenarios of *r*_{1} introduce an optimization variable *M*_{a} that stands for the RBM of the vertex *v*_{a}. The set of optimization variables is then expressed as $\Phi \u2032=\Phi \u228e{\alpha s1\u22123,la,s1\u22123,Ma}$.

#### 3.2.5 Kinematic Modeling of **r**_{1,1}.

*r*

_{1}is applied,

*v*

_{a}becomes an internal vertex of the crease pattern and can be kinematically modeled using the PTU. According to the scenario of

*r*

_{1,1}in Fig. 4, the units

*u*

_{ia}of

*v*

_{a}are expressed in terms of their respective sector and dihedral angles as

If *v*_{a} exhibits only a single predecessor, the sum over $\alpha pj$ is neglected. Then, the units in Eq. (6) are added to the set of units, $u\u2032=u\u228e{u1a\u22123a}$. By setting these units and the RBM *M*_{a} into *f*_{φ}, the method determines the dihedral angles $\phi 1a\u22123a$ and associates them with the generated edges, $\rho a,s1\u22123=\phi 1a\u22123a$, which are then added to the set of dihedral angles $\rho \u2032=\rho \u228e{\rho a,s1\u22123}$.

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\u22123a$ of *v*_{a} by applying $fR$. Then, these local rotation matrices are turned into global rotation matrices, $Ra,si=Rp1,aRz(\pi )Ria$, which are added to the set of global rotation matrices $R\u2032=R\u222a{Ra,s1\u22123}$. A special case arises when *v*_{a} belongs to the initial graph since the global rotation matrix $Rp1,a$ has not been determined by an application of *r*_{1}. In this case, the method sets $Rp1,a=Rz(\alpha 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\u22c5(1,0,0)T$. These coordinates are subsequently added to the set of three-dimensional vertex coordinates, $X\u2032=X\u222a{Xs1\u22123}$.

#### 3.2.6 Optimization Constraints of **r**_{1,1}.

*r*

_{1}, 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

*v*

_{a}exhibits only a single predecessor, the sum over $\alpha pj$ in Eq. (7) is neglected. The same reasoning is applied to the lengths of the generated edges:

The second type of constraint introduced by *r*_{1} corresponds to the condition for rigid foldability in Eq. (1) that states *U*_{max}(*t*) ≤ *U*_{med}(*t*) + *U*_{min}(*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.

*u*

_{1}, whereas

*u*

_{2}and

*u*

_{3}as well as their unit angles

*U*

_{2}and

*U*

_{3}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(

*U*

_{1}) and max(

*U*

_{1}). In general, however, it is reasonable to assume that max(

*U*

_{1}) occurs at

*t*= 0, in which case Eq. (1) is satisfied because of the initial flat state. Then, the unit angle

*U*

_{1}decreases monotonously once

*t*deviates from zero toward

*t*=

*t*

_{max}. 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(

*U*

_{1}) that occurs at

*t*=

*t*

_{max}:

#### 3.2.7 Rule **r**_{2}: Combine Vertices.

The second rule *r*_{2} is defined as a production of the type *r*_{2}:*LHS*_{2} → *RHS*_{2} in Fig. 5, where *LHS*_{2} and *RHS*_{2} are both shown in red, the surrounding graph is shown in gray, and adjustments to the constraint system are shown in blue.

Rule *r*_{2} accepts two vertices $vc1$ and $vc2$ that share an edge label on the *LHS*_{2} and combines them into a single vertex $vc1$ on the *RHS*_{2}. 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, *r*_{2} 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 $\alpha c2$ in Fig. 5 left, respectively, with the parts colored in blue in Fig. 5 right. The vertex *v*_{b} signifies any vertex that was previously coupled to the sector angles of the reassigned edge ($\alpha c2$).

#### 3.2.8 LHS Matching of **r**_{2}.

A match $M2$ of the *LHS*_{2} in *G* involves four conditions. Both vertices $vc1$ and $vc2$ have to be extendable, $Tc1=Tc2=\chi $. 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 *LHS*_{2} 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 *r*_{1}.

To prevent the combination of vertices that share the same predecessor, such as $vs1\u22123$ in Fig. 4, the third condition states that $vc1$ and $vc2$ cannot exhibit the same predecessors, $Pc1\u2229Pc2={\u2205}$.

Since *r*_{2} 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 **r**_{2}.

Once *r*_{2} is applied, $G\u21d2r2G\u2033$, the *RHS*_{2} includes $vc1$ while $vc2$ is subtracted from the set of vertices, $V\u2033=V\u2216{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\u2033=(E\u2216{ep2,c2})\u228e{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 *RHS*_{2} can still be extended, resulting in the vertex label $Lc1\u2033=((vp1,vp2),\chi )$. The edge label $Lc2$ is then subtracted from the set of vertex labels, $LV\u2033=LV\u2216{Lc2}$. The edge label formerly corresponding to the edge incident to $vc2$ is transferred identically to the reassigned edge, $Lp2,c1\u2033=Lp2,c2$. Finally, the edge label of the edge incident to $vc2$ is subtracted from the set of edge labels, $LE\u2033=LE\u2216{Lp2,c2}$. If $vc2$ is incident to multiple incoming edges on the *LHS*_{2}, the above procedure is applied to all edges and edge labels.

#### 3.2.10 Coordinates and Optimization Variables of **r**_{2}.

The in-plane location and thus the coordinates $xj$ of a vertex *v*_{j} are always determined by a single sector angle *α*_{j} and a single edge length $lpi,j$. Since *r*_{2} deletes the vertex $vc2$, the corresponding sector angle and edge lengths of $vc2$ are eliminated from the set of optimization variables, $\Phi \u2033=\Phi \u2216{\alpha c2,lp2,c2}$.

#### 3.2.11 Optimization Constraints of **r**_{2}.

*r*

_{2}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, $\u2220xp1xc1xp2$. The constraint introduced by

*r*

_{2}applies to the upper and lower bounds equivalent to Eq. (7):

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

*r*

_{2}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,

*r*

_{2}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 $\alpha c2$ is replaced with $\u2220xc1xp2xb$ and $lp2,c2$ is replaced with the Euclidean distance $||xc1\u2212xp2||$ between $vp2$ and $vc1$. Hence

When $vc2$ has multiple predecessors on the *LHS*_{2}, 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 *r*_{2} 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 *r*_{1} and *r*_{2} 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 *G*_{0} and the maximum number of internal vertices *N*_{max}.

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 *L*_{V} is guaranteed to coincide, whereas the set of edge labels *L*_{E} 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 *N*_{max} by applying both rules *r*_{1} and *r*_{2} 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 2^{N} 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 *l*_{min}, *l*_{max}, *α*_{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* = *t*_{opt} ∈ [ − π, π] at which the dihedral angles *ρ*_{i,j}(*t*_{opt}) 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 *l*_{i,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*(*l*_{i,j} ← *l*_{min}).

*Ω*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 $\Phi opt*$ that optimizes the objective function such that $\Omega (\Phi opt*)=\Omega *$.

*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 $\Omega *\u2264\u03d1$, where $\u03d1$ 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 $\Phi 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(\Phi opt\u2190\Phi 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 *v*_{i} 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* = *t*_{max} 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.

### 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 *G*_{0} (identical for both tasks) with vertices *v*_{1} and *v*_{2} end edge *e*_{1,2} is shown in black and corresponds to the graph described in Fig. 3 with the same vertex labels $L1=((\u2205),\u2205)$ and *L*_{2} = ((*v*_{1}), *χ*), edge label *L*_{1,2} = (*f*_{1}, *f*_{2}), as well as in-plane coordinates $x1=(0,0)$ and $x2=(1,0)$. The driving angle of the edge *e*_{1,2} in both cases is linear, *ρ*_{1,2} = *t*, where *t* goes from zero to *t*_{max} = π/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 *N*_{max} = 3. Note that graphs with one or two internal vertices are possible, *N*_{max} is just the upper limit.

The point object *P*_{G} to grip in Fig. 6(a) is located at *P*_{G} = (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\u2032$. Figure 6(b) shows the trajectory to approximate and the 11 points $PTi$ located at (1, 0, 0) + (2cos *t*_{i}, 0, 2 sin *t*_{i}) where *t*_{i} goes from 0 to π/2 in steps of π/20.

### 4.2 Generation: Additional Filters.

To facilitate the description of a crease pattern graph, here the following notation is introduced for the rule application sequences: when a vertex *v*_{i} is extended by *r*_{1}, *i* is added to the sequence, and when vertices *v*_{i} and *v*_{j} are combined by *r*_{2}, (*i*, *j*) is added to the sequence. As an example, extending vertices *v*_{2} and *v*_{3} in the given order and then combining vertices *v*_{4} and *v*_{8} 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.

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 *P*_{G} 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 *r*_{1} 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 *v*_{3} is the vertex extended last, *v*_{9}, *v*_{10}, or *v*_{11} are the gripping vertex candidates. In this case, extending *v*_{4} is semantically invalid since all of its successors *v*_{6}, *v*_{7}, and *v*_{8} 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.

*P*

_{G}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

*v*

_{l,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

*P*

_{G}and the location of the gripping vertex at

*t*

_{opt}=

*t*

_{max}= π/2:

*t*

_{opt}=

*t*

_{i}= (0, π/20, …, 9π/20, π/2):

An optimization of the robotic arm is considered successful if the objective value $\Omega *$ of a design alternative is smaller than or equal to $\u03d1=3*10\u22121$. The variable boundaries for both tasks are *l*_{min} = 0.1, *l*_{max} = 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 *G*_{0} in Fig. 6, *N*_{max} = 3, and rules *r*_{1} and *r*_{2}, 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.

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 *P*_{G}, 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).

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 $\Omega *=2*10\u22122$ (a), $\Omega *=3*10\u22122$ (b), and $\Omega *=4*10\u22122$ (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).

## 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 *P*_{G} 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 *r*_{1} and one application of *r*_{2}, while only *r*_{2} can be applied thereafter. Rule *r*_{2} introduces constraints into the system without adding variables, which is why the number of designs that can grip *P*_{G} 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 *r*_{1} (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 *P*_{G}.

When at first one could have argued not to apply *r*_{2} since it only introduces constraints, the increase in feasible designs from three to four rule applications show that a certain number of applications of *r*_{2} can drive the optimized designs out of intersecting configurations. Without specific constraints that prevent intersection within the optimization, *r*_{2} thus provides more variability in the crease pattern graphs that satisfy the design task. However, if *r*_{2} is applied too many times, the designs are either unable to grip *P*_{G} 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 *r*_{2}. 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 *l*_{max} = 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 *r*_{1} predominates, the midsection of the search tree becomes broader since more extended vertices offer more possibilities for the application of *r*_{2}. The number of possible rule applications then progressively decreases with each line because of the maximum number of extended vertices *N*_{max} = 3, which is why the search tree narrows again at the bottom. However, even with a small increase in *N*_{max} to 8 vertices, for example, the maximal number of RBMs for only one particular topology is 2^{8} 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 *N*_{max}, the enumeration influences the scalability of the method significantly and blocks more complex engineering applications. For example, for *N*_{max} = 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 *r*_{2}, here denoted as $r2\u2032$. Rule $r2\u2032$ would also combine two vertices, but the vertex $vc1\u2032$ on the *RHS*_{2} in Fig. 5 would lie on the line segment between the predecessors $vp1$ and $vp2$ and would become a sink vertex ($deg+(vc1\u2032)=0$) of the type $Tc1\u2032=\u2205$. Figure 12 illustrates how a Miura-ori pattern could be generated by the expansion of the presented rule system with rule $r2\u2032$.

The crease pattern graph depicted in Fig. 12(a) can be generated by the initial graph *G*_{0} in Fig. 3 with the rule application sequence 2, 3, 4, 8. By applying $r2\u2032$ to vertices *v*_{9} and *v*_{14}, 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, $\rho 8,9\u2032=\rho 4,9\u2032$, 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 *r*_{1} or a minimum number of applications of rule *r*_{2}. 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

*e*_{i,j}=edge from vertex

*v*_{i}to vertex*v*_{j}*u*_{i}=unit,

*i*= 1, 2, 3- $fR$ =
function determining local rotation matrices

*f*_{φ}=function determining unknown dihedral angles

*L*_{i}=vertex label of vertex

*v*_{i}*L*_{i,j}=edge label of edge

*e*_{i,j}*L*_{E}=set of edge labels

*L*_{V}=set of vertex labels

*P*_{i}=predecessors of a vertex

*v*_{i}*T*_{i}=type of a vertex

*U*_{i}=unit angle,

*i*= 1, 2, 3- $Ri$ =
local rotation matrix,

*i*= 1, 2, 3- $Ri,j$ =
global rotation matrix of the edge

*e*_{i,j}- $FLi,j,FRi,j$ =
facets on the left and right of edge

*e*_{i,j}, respectively*M*_{↑},*M*_{↓}=rigid body modes up and down, respectively

*P*_{G},*P*_{T}=point to grip and points of trajectory, respectively

*α*=sector angle

*ρ*=set of dihedral angles

*ρ*_{i,j}=dihedral angles

- $\u03d1$ =
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

- $\u2205$ =
empty and terminal symbol

- $\u228e$ =
disjoint union of sets