R. Sachse, F. Geiger and M. BischoffConstrained motion design with distinct
actuators and motion stabilization

Renate Sachse, Institute for Structural Mechanics, University of Stuttgart,
Pfaffenwaldring 7, D-70550 Stuttgart, Germany, E-mail<EMAIL_ADDRESS>stuttgart.de

# Constrained motion design with distinct actuators and motion stabilization

Renate Sachse, 1  Florian Geiger and Manfred Bischoff  1  1
11affiliationmark:  Institute for Structural Mechanics, University of
Stuttgart,
Pfaffenwaldring 7, D-70550 Stuttgart, Germany

###### Abstract

The design of adaptive structures is one method to improve sustainability of
buildings. Adaptive structures are able to adapt to different loading and
environmental conditions or to changing requirements by either small or large
shape changes. In the latter case, also the mechanics and properties of the
deformation process play a role for the structure’s energy efficiency. The
method of variational motion design, previously developed in the group of the
authors, allows to identify deformation paths between two given geometrical
configurations that are optimal with respect to a defined quality function. In
a preliminary, academic setting this method assumes that every single degree
of freedom is accessible to arbitrary external actuation forces that realize
the optimized motion. These (nodal) forces can be recovered a posteriori. The
present contribution deals with an extension of the method of motion design by
the constraint that the motion is to be realized by a predefined set of
actuation forces. These can be either external forces or prescribed length
chances of discrete, internal actuator elements. As an additional constraint,
static stability of each intermediate configuration during the motion is taken
into account. It can be accomplished by enforcing a positive determinant of
the stiffness matrix.

###### keywords:

Adaptive structures, motion design, constraint enforcement, compliant
mechanism, actuation, stabilized motion

##  1 Introduction

Construction industry is responsible for a substantial part of global energy
consumption and requirement of material resources. Therefore, the improvement
of energy efficiency and sustainability in this field represents a major
challenge for architects and engineers, who are determined to design extremely
efficient structures. Adaptive structures are one technology to improve
structural efficiency by adapting to changing circumstances.

Two types of adaptive structures can be distinguished. The first type,
sometimes denoted as smart structure, adapts its internal forces or
deformation to varying loads. With the use of sensors and actuators, a light
weight design can be accomplished while maintaining structural performance. In
doing so, forces and deformations counteracted by actuation [9, 27, 28, 13].
Various types of actuation are used, e.g., piezoelectric elements [12] or
discrete actuator elements in truss systems, which are able to adjust their
length. The question of where to place these actuators as efficiently as
possible is answered by methods for actuator placement [1, 8, 18, 22]. In this
type of adaptive structures, actuation typically induces only small
displacements, such that linear structural analyses suffices.

The second type of adaptive structures adapts to environmental changes or to
varying service demands through major shape changes. Here, the individual
geometric configurations vary significantly. Examples are retractable roofs of
stadiums (e.g. the Commerzbank-Arena in Frankfurt, Germany [6]) or closing and
opening of facade elements (e.g. the biomimetic facade element Flectofin
[15]), which contribute to energy efficiency of a building. But also in
aviation, shape changes by morphing wings to improve efficiency of airplanes
are investigated [20, 3, 30, 32, 2]. The design of flexible or morphing
structures that allow for such large geometry changes poses a major challenge.
Classically, mechanisms such as folding or sliding with the help of joints are
used. For discrete structural typologies, technologies that use length-varying
actuator elements have been developed, for instance for truss structures [11,
25] and tensegrity structures [7, 33, 29]. In the case of continuous
structures, such as shells, optimization procedures can be used to increase
compliance [26, 14, 16, 4, 21, 18], thus facilitating continuous bending
deformation instead of discrete hinges.

However, not only the geometry of the morphing structure is important for its
efficiency. The deformation process has to meet certain requirements and the
actuators that control it consume energy. The task of efficient actuation is
tackled in control theory, especially in optimal control and robotics.
However, pertinent methods mostly focus on rigid systems with discrete
kinematics or flexible multibody systems. This is different in the field of
continuum robots and hyper-redundant manipulation [23], but the continuity and
thus a large number of degrees of freedom still poses a challenge.
Nevertheless, mechanics and structural analysis have already been combined
with optimal control procedures [10, 19, 31].

The topic of the present study is the design of deformation processes between
an initial geometry and a prescribed target geometry that are optimal with
respect to a chosen objective function. In a previous work by the authors
[24], the method of _variational motion design_ has been presented. It is
based on a variational formulation and was developed using an exemplary
objective function, namely the internal energy integrated along the
deformation path. In robotics, the quantity “cost of transport” allows to
quantitatively compare different types of transport between two locations,
e.g. flying or floating. This quantity can be reinterpreted for deformable
structures in motion design to a “cost of deformation”, which represents a
measure for the required energy to deform the structure into the desired
target geometry. The method was verified and validated by means of some
numerical benchmarks and experiments. A focus was laid on studies of rigid
body motions, kinematic structures and inextensional deformations of shells. A
brief description of the method is given in Section 2.

The method of motion design provides an optimized deformation trajectory that
minimizes a certain cost function (e.g. the integrated internal energy). A
major hypothesis in this context is that all degree of freedom are accessible
to actuation by prescribed forces. These can be computed a posteriori from
equilibrium, but of course technical realization of such an omnipotent
actuation is rather academic. It is one of the objectives of the present study
to remove this limitation. The restriction to a given set of potential
actuator forces represents a constraint within the variational framework to
compute the resulting motion.

The constraint that the motion is to be enabled with a prescribed set of
actuation forces is enforced by an extended functional with constraints using
the Lagrange multiplier method, as explained in Section 3. Also length changes
of discrete actuator elements can be considered as actuation modes. For this
purpose, Section 4 first introduces an actuator element formulation that is
suitable for implementation into the framework of motion design. Furthermore,
also other constraints may be introduced in this manner. In particular, static
stability of every intermediate state during deformation can be guaranteed by
enforcing a positive determinant of the stiffness matrix. This overall
procedure as well as a combination of both types of constraints (actuation and
stability) are presented in Section 5.

##  2 Motion design of structures

The bases of variational motion design is the functional $J$. As objective
function it represents the property that is to be assigned to the motion. The
procedure is described in detail in Sachse and Bischoff [24] and only the most
important basics are repeated herein. Geometrically non-linear, quasi-static
structural behavior is assumed. Following [24], the internal energy,
integrated along the motion path $s$ is considered. With the Green-Lagrange
strain ${\mathbf{E}}$ and a linear elastic St. Venant-Kirchhoff material law
this leads to

$\displaystyle J$
$\displaystyle=\int_{s}{\Pi}_{\mathrm{int}}\,{\mathrm{d}}s=\int_{s}\int_{\Omega}\frac{1}{2}{\mathbf{E}}^{\mathrm{T}}{\mathbf{C}}{\mathbf{E}}\,{\mathrm{d}}{\Omega}\,{\mathrm{d}}s=\min.$
(1)

Alternative objective functions can be used.

Figure 1: Motion design on the example of a two bar truss and a snap-through
problem.

A normalized arc length $\bar{s}$ for a clear specification of the path
integration bounds and a mapping parameter
$s_{\mathrm{u}}=\frac{\mathrm{d}s}{\mathrm{d}\bar{s}}$ are introduced. The
first variation of $J$ is then

$\displaystyle{\delta}J$
$\displaystyle=\int_{0}^{1}\bigg{[}\int_{\Omega}{\delta}{\mathbf{E}}^{\mathrm{T}}{\mathbf{C}}{\mathbf{E}}\,{\mathrm{d}}{\Omega}s_{\mathrm{u}}+\int_{\Omega}\frac{1}{2}{\mathbf{E}}^{\mathrm{T}}{\mathbf{C}}{\mathbf{E}}\,{\mathrm{d}}{\Omega}{\delta}s_{\mathrm{u}}\bigg{]}\,{\mathrm{d}}\bar{s}=0.$
(2)

With the usual spatial discretization of the domain ${\Omega}$ with finite
elements, the variation can be rewritten as

$\displaystyle{\delta}J$
$\displaystyle=\int_{0}^{1}\bigg{[}({\delta}{\mathbf{d}}^{s})^{\mathrm{T}}{\mathbf{F}}_{\mathrm{int}}^{s}s_{\mathrm{u}}^{s}+({\delta}{\mathbf{d}}_{,s}^{s})^{\mathrm{T}}{\mathbf{s}}_{\mathrm{u}}^{s}{\Pi}_{\mathrm{int}}^{s}\bigg{]}\,{\mathrm{d}}\bar{s}=0,$
(3)

where the internal energy ${\Pi}_{\mathrm{int}}^{s}$ as well as the global
vector of internal forces ${\mathbf{F}}_{\mathrm{int}}^{s}$ can be identified.
Both quantities are still continuous with respect to the path parameter $s$,
which is indicated by the superscript $(\;)^{s}$. In order to solve this
equation, a second discretization, namely a discretization of the motion path
from the initial configuration to a prescribed deformed target configuration
with $\bar{n}_{{\mathrm{ele}}}$ elements is introduced

$\displaystyle\int_{0}^{1}(\ldots)\,{\mathrm{d}}\bar{s}=\sum_{{\bar{e}}=1}^{\bar{n}_{{\mathrm{ele}}}}\int_{s^{\bar{e}}}(\ldots)\,{\mathrm{d}}\bar{s}^{\bar{e}}.$
(4)

There is a formal analogy to space-time finite elements. However, the
corresponding variable $s$ does not comply with time and the underlying
functional $J$ is not obtained from the equations of motion. Successive
linearization and discretization of the motion path yields a linearized
algebraic system of equations for motion design problems

$\displaystyle{\mathbf{K}}_{{\mathrm{m}}{\mathrm{d}}}{\Delta}\bar{{\mathbf{D}}}$
$\displaystyle=-{\mathbf{R}}_{{\mathrm{m}}{\mathrm{d}}}\,,$ (5)

with the definitions

$\displaystyle{\mathbf{K}}_{{\mathrm{m}}{\mathrm{d}}}$
$\displaystyle=\sum_{{\bar{e}}=1}^{\bar{n}_{{\mathrm{ele}}}}\int_{\bar{s}^{\bar{e}}}\Big{(}\bar{{\mathbf{N}}}^{\mathrm{T}}{\mathbf{K}}_{\mathrm{T}}s_{\mathrm{u}}\bar{{\mathbf{N}}}+\bar{{\mathbf{N}}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{int}}{\mathbf{s}}_{\mathrm{u}}\bar{{\mathbf{N}}}_{,s}+\bar{{\mathbf{N}}}_{,s}^{\mathrm{T}}{\mathbf{s}}_{\mathrm{u}}{\mathbf{F}}_{\mathrm{int}}\bar{{\mathbf{N}}}+\bar{{\mathbf{N}}}_{,s}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{u}}{\Pi}_{\mathrm{int}}\bar{{\mathbf{N}}}_{,s}\Big{)}\,{\mathrm{d}}\bar{s}^{\bar{e}},$
(6) $\displaystyle{\mathbf{R}}_{{\mathrm{m}}{\mathrm{d}}}$
$\displaystyle=\sum_{{\bar{e}}=1}^{\bar{n}_{{\mathrm{ele}}}}\int_{\bar{s}^{\bar{e}}}\Big{(}\bar{{\mathbf{N}}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{int}}s_{\mathrm{u}}+\bar{{\mathbf{N}}}_{,s}^{\mathrm{T}}{\mathbf{s}}_{\mathrm{u}}{\Pi}_{\mathrm{int}}\Big{)}\,{\mathrm{d}}\bar{s}^{\bar{e}}.$
(7)

$\bar{{\mathbf{N}}}$ represents the matrix of shape functions for path
discretization and can contain any type of function, i.e., Lagrange
polynomials as well as spline functions. The vector of the total displacements

$\displaystyle\bar{{\mathbf{D}}}$
$\displaystyle=\begin{bmatrix}{\mathbf{D}}^{1}\\\ {\mathbf{D}}^{2}\\\
\vdots\\\ {\mathbf{D}}^{\bar{k}}\\\ \vdots\\\
{\mathbf{D}}^{\bar{n}_{\text{node}}}\end{bmatrix}$ (8)

consists of $\bar{n}_{\text{node}}$ subvectors ${\mathbf{D}}^{s}$ and gathers
the displacements of every degree of freedom in every single configuration
throughout the deformation process, where $\bar{n}_{\text{node}}$ denotes the
number of nodes in the path discretization. Every intermediate configuration
corresponds to one node $\bar{k}$ of the path discretization. Thus, the shape
functions in $\bar{{\mathbf{N}}}$ serve as interpolation between the
configurations. The total number of degrees of freedom is therefore

$\displaystyle\bar{n}_{\mathrm{dof}}=\bar{n}_{\text{node}}\cdot
n_{\mathrm{dof}}.$ (9)

The entire problem is solved monolithically by an iterative solution of the
linearized system of equations (5). However, the problem might become ill-
posed, as particular displacement trajectories may be representable by various
different solution vectors. A simple example for this is a straight line,
which can be represented by any sequence of straight steps with different
size. To guarantee uniqueness of the solution, the evolution of at least one
displacement degree of freedom has to be controlled throughout the entire
deformation process. The entire derivation and further explanations are given
in [24].

The following example briefly demonstrates the application of motion design.
The initial geometry is a shallow, discrete two-bar truss, shown in black in
Figure 1a. It ought to be deformed to the target configuration shown in blue.
The deformation path between the two configurations is supposed to minimize
the integrated internal energy, defined as cost of deformation. During the
deformation process, the point $P$ follows a yet unknown trajectory
(generically shown in red) until it arrives at its end position $P^{\prime}$.
This trajectory and the corresponding internal energy can be identified on the
diagram in Figure 1a. The area of the spanned surface represents the value of
the functional and is to be minimized by the resulting trajectory from motion
design. The principle of path discretization is illustrated by four linear
path elements in Figure 1b.

In our numerical experiment, however, actually 14 linear elements are used to
achieve a better approximation. Iterative solution of the system of equations
in eq. (5) provides the entire deformation path in one go. The result is shown
in Figure 1b. It turns out that, for the given objective function $J$, it is
optimal to first enforce a vertical snap through, followed by a horizontal
movement towards the location of $P^{\prime}$.

In order to technically realize the obtained deformation, forces are required.
These external forces are obtained in a post-processing step after solving the
motion design problem. From the displacement vectors ${\mathbf{D}}^{s}$ of the
optimized motion, the discrete internal forces can be computed for every
configuration. The required external (actuation) forces are then obtained from
equilibrium.

##  3 Constraints for prescribed sets of actuation forces

###  3.1 General concept

The motion design method from [24] relies on the hypothesis that discrete
actuation forces can potentially be applied to every degree of freedom. This
is illustrated in the simple truss system in Figure 2a. The initial
configuration is shown in black. In this example, the target configuration is
not fully prescribed, but only the vertical displacement $D_{2}$ of the second
node is supposed to take a specific value. In an unconstrained motion design,
a motion path is found on the basis that all forces $F_{1}$, $F_{2}$ and
$F_{3}$ are potentially non-zero and thus available for actuation, as
indicated in Figure 2b.

In constrained motion design, these forces are subject to conditions, which
can be quite general. Practically, these conditions emanate from
considerations concerning technical feasibility. In the simplest case, such a
condition can be that a certain force is zero, i.e., there is no external
actuator available to apply a force in the corresponding direction. In Figure
2c, this is represented by a system, in which only two forces $F^{\prime}_{1}$
and $F^{\prime}_{2}$ can be applied to the structure. $F_{3}$ has to be zero
during the entire motion. This also affects the end geometry, which is further
discussed in the following.

Figure 2: Illustrating example for unconstrained and constrained motion
design.

###  3.2 Equality constraints for uncontrolled degrees of freedom

The fact that the desired motion has to be realized with a limited set of
actuation forces represents a constraint within the motion design problem. The
main functional and objective remain the same. Generally, the task of motion
design with equality constraints can be formulated as: Minimize the functional
$J$ subject to constraints ${\mathbf{g}}={\mbox{$\mathbf{0}$}}$.

The constraints are enforced by an extension of the functional. Here, the
Lagrange multiplier method is applied. However, other methods can be used,
e.g., the penalty method or the Augmented Lagrange method. Within the Lagrange
multiplier method, additional variables, the Lagrange multipliers
${\bm{\lambda}}$, are introduced for each constraint. It then follows for the
extended functional

$\displaystyle\bar{J}=J+{\bm{\lambda}}^{\mathrm{T}}{\mathbf{g}}.$ (10)

According to the constraint, external actuation forces can only act in the
direction of specific degrees of freedom. For all other degrees of freedom,
the corresponding internal forces have to be zero due to equilibrium. This
condition can be expressed as a set of equality constraints on the motion.

As the described constraints must hold for the entire motion, they are
introduced for each configuration during the deformation process. They
therefore represent pointwise constraints, which are not included in the
integral of the main functional and can be formally expressed as

$\displaystyle{\mathbf{g}}$ $\displaystyle=\begin{bmatrix}{\mathbf{F}}^{1}\\\
{\mathbf{F}}^{2}\\\ \vdots\\\ {\mathbf{F}}^{\bar{k}}\\\ \vdots\\\
{\mathbf{F}}^{\bar{n}_{\text{node}}}\end{bmatrix}={\mbox{$\mathbf{0}$}}.$ (11)

Here, the vectors ${\mathbf{F}}^{\bar{k}}$ gather the internal forces of only
those degrees of freedom that are not subject to actuation forces. 111When
using a path discretization with an approximation by B-spline shape functions,
these constraints should not be imposed on the control points but on the
deformation path itself. For this reason, the constraints are enforced at
collocation points, which are located on the path. In this case, Greville
points are used.

For the system presented in Figure 2 and a motion path discretization with
linear path elements, the internal force at degree of freedom number 3 must be
zero in every and five path increment. Consequently, the vector of constraints
can be written as

$\displaystyle{\mathbf{g}}$
$\displaystyle=\begin{bmatrix}F^{2}_{3}&F^{3}_{3}&F^{4}_{3}&F^{5}_{3}\end{bmatrix}^{\mathrm{T}},$
(12)

i.e., the internal force at degree of freedom number 3 must be zero in every
path increment. The internal force $F_{3}^{1}$ of the initial configuration is
omitted as the initial configuration is always assumed to be load-free.

If distributed line loads or surface loads are applied as actuation forces,
the discrete forces of the individual degrees of freedom share the same
scaling factor. In such cases, a master degree of freedom is selected and the
associated forces are coupled to this load value via an equality constraint.

After linearization, the algebraic system of equations of the constrained
problem reads

$\displaystyle\begin{bmatrix}{\mathbf{K}}_{{\mathrm{m}}{\mathrm{d}}}+{\lambda}_{i}{\mathbf{H}}_{i}&{\mathbf{G}}^{\mathrm{T}}\\\
{\mathbf{G}}&{\mbox{$\mathbf{0}$}}\end{bmatrix}\begin{bmatrix}{\Delta}\bar{{\mathbf{D}}}\\\
{\Delta}{\bm{\lambda}}\end{bmatrix}=-\begin{bmatrix}{\mathbf{R}}_{{\mathrm{m}}{\mathrm{d}}}+{\bm{\lambda}}^{\mathrm{T}}{\mathbf{G}}\\\
{\mathbf{g}}\end{bmatrix}.$ (13)

Here, the matrix ${\mathbf{G}}$ contains the first derivatives of all
constraints with respect to the degrees of freedom $\bar{{\mathbf{D}}}$,

$\displaystyle{\mathbf{G}}=\frac{{\partial}{\mathbf{g}}}{{\partial}\bar{{\mathbf{D}}}},$
(14)

the Hessian ${\mathbf{H}}_{i}$ contains the second derivatives of the
$i^{\text{th}}$ constraint,

$\displaystyle{\mathbf{H}}_{i}=\frac{{\partial}g_{i}}{{\partial}\bar{{\mathbf{D}}}\partial\bar{{\mathbf{D}}}}.$
(15)

This means that the constraints and thus the corresponding internal forces
${\mathbf{F}}^{\bar{k}}$, need to be differentiated twice with respect to all
degrees of freedom. The internal force at a degree of freedom $j$ in the load
configuration $\bar{k}$ only depends on the degrees of freedom at the
associated path node. Moreover, the tangent stiffness matrix
${\mathbf{K}}_{\mathrm{T}}^{\bar{k}}$ at this path node $\bar{k}$ already
incorporates the first derivative of the internal forces with respect to the
relevant degrees of freedom. Accordingly, the derivative of one component
$F^{\bar{k}}_{j}$ of the total internal force vector ${\mathbf{F}}^{\bar{k}}$
is provided by the associated column ${\mathbf{K}}_{j}^{\bar{k}}$ of the
stiffness matrix. Consequently, the vector ${\mathbf{G}}_{i}$ of first
derivatives of the $i^{\text{th}}$ constraint $g_{i}$ follows as

$\displaystyle{\mathbf{G}}_{i}=\frac{{\partial}g_{i}}{{\partial}{\mathbf{D}}^{\bar{k}}}=\frac{{\partial}F_{j}^{\bar{k}}}{{\partial}{\mathbf{D}}^{\bar{k}}}={\mathbf{K}}_{j}^{\bar{k}},$
(16)

where $j$ is the degree of freedom that corresponds to the constraint $g_{i}$.
Additionally, the second derivative of the $i^{\text{th}}$ constraint is
required for the Hessian matrix ${\mathbf{H}}_{i}$. In this case, this
represents the derivative of the associated column of the stiffness matrix.
Usually, the relevant column of the stiffness matrix ${\mathbf{K}}_{j}$ has to
be differentiated with respect to each degree of freedom. This yields for the
Hessian matrix of the $i^{\text{th}}$ constraint

$\displaystyle{\mathbf{H}}_{i}=\frac{{\partial}{\mathbf{K}}_{j}^{\bar{k}}}{{\partial}{\mathbf{D}}^{\bar{k}}}\,.$
(17)

However, the analytical derivative of the tangent stiffness matrix can only be
computed in special cases, such that numerical differentiation has to be
applied. For classical forward, backward or central difference schemes the
choice of the perturbation in order to minimize both the truncation error and
the elimination error is not trivial. Therefore, here complex step
differentiation, as presented in [17] is used. It allows for arbitrarily small
perturbations because there is no elimination error and thus accuracy up to
machine precision is always guaranteed.

Nevertheless, numerical differentiation with respect to each individual degree
of freedom is numerically expensive, as it involves repeated computation of
the stiffness matrix for each direction of differentiation. However, because
derivatives commute and the tangent stiffness matrix is already calculated
through displacement derivatives, the derivative of the $j^{\text{th}}$ column
with respect to the $m^{\text{th}}$ degree of freedom equals the derivative of
the $m^{\text{th}}$ column with respect to the $j^{\text{th}}$ degree of
freedom, i.e.,
$\frac{{\partial}{\mathbf{K}}_{j}}{{\partial}D_{m}}=\frac{{\partial}{\mathbf{K}}_{m}}{{\partial}D_{j}}$.
Therefore, the required differentiation with respect to column $j$ of the
stiffness matrix can be extracted from the derivative of the total stiffness
matrix by a scalar product with a unit vector ${\mathbf{I}}_{j}$ in direction
of $j$

$\displaystyle{\mathbf{H}}_{i}=\frac{{\partial}{\mathbf{K}}_{j}^{\bar{k}}}{{\partial}{\mathbf{D}}^{\bar{k}}}=\frac{{\partial}{\mathbf{K}}^{\bar{k}}}{{\partial}{\mathbf{D}}^{\bar{k}}}\cdot{\mathbf{I}}_{j}.$
(18)

Thus, the derivative in eq. (17) is transformed to a directional derivative in
eq. (18) and the advantage of using numerical differentiation regarding
numerical efficiency becomes apparent. Instead of differentiating a single
column of the stiffness matrix with respect to all degrees of freedom and
consequently evaluating it as many times for numerical differentiation, the
whole stiffness matrix now only has to be calculated one additional time in
the direction of ${\mathbf{I}}_{j}$. This approach also has major advantages
regarding implementation into existing finite element software, as usually the
stiffness matrix is calculated in an entirety and not only distinct columns.

The derivatives of the system in Figure 2 for the first ($i=1$) constraint
$F^{2}_{3}=0$ are exemplarily presented in the following. In this case, there
is only one relevant degree of freedom, namely $j=3$. The first derivative can
be extracted from the total stiffness matrix of the current load configuration
$\bar{k}=2$, where it represents the third column

$\displaystyle{\mathbf{G}}_{1}=\frac{{\partial}F^{2}_{3}}{{\partial}{\mathbf{D}}^{2}}={\mathbf{K}}^{2}_{3}.$
(19)

The Hessian matrix contains the whole stiffness matrix of the second
configuration and can be calculated with the help of a directional derivative

$\displaystyle{\mathbf{H}}_{1}=\frac{{\partial}{\mathbf{K}}^{2}}{{\partial}{\mathbf{D}}^{2}}\cdot{\mathbf{I}}_{3}.$
(20)

With these derivatives, the extended system of equations from eq. (13) can be
generated and solved. Thus, an optimized motion can be found that is realized
exclusively by the defined actuation forces. The result is the deformation
path and the development of the actuation force amplitudes during the motion.

###  3.3 Restrictions for the prescribed target geometry

In the basic motion design method, either the entire target geometry or only
parts of it have to be prescribed. This is similar for constrained motion
design, but with one important restriction: As equilibrium has to be fulfilled
using only the available actuation forces throughout the entire motion, this
must also be the case for the prescribed target geometry. Otherwise, the
problem is not well-posed.

Two different cases have to be distinguished, depending on the number of
prescribed target displacement values:

  * •

If the entire target geometry is prescribed, an optimization is carried out
prior to the motion design process. The objective function to be minimized is
the difference between the target geometry and the current geometry with the
constraint that the internal forces in direction of non-actuated degrees of
freedom must be equal to zero.

  * •

If only a subset of degrees of freedom is prescribed, the rest of the geometry
can adjust freely to meet the condition of equilibrium with the available
actuation forces. The maximum number of prescribed values of displacement
degrees of freedom depends on the number of independent actuation forces. It
is not possible to prescribe more displacement values than this value, whereas
fewer are generally possible. In cases, where it is not obvious, whether an
equilibrium configuration is possible, the same optimization procedure as
described above can be applied.

###  3.4 Solution and interpretation of the results

Taking into account the described restrictions for the target geometry, the
extended system of equations (13) can be solved. In doing so, an optimized
motion path is obtained that considers the given constraints of enforcing an
equilibrium state with only the available actuation forces. The results for
the example in Figure 2 with a varying number of actuation forces are
summarized in Figure 3.

Figure 3: Solutions of the illustrating example with different numbers and
types of load cases.

If forces can be applied on every degree of freedom, no constraints are
enforced and the solution is identical to the result from an unconstrained
motion design. This is illustrated in Figure 3a. It can be seen in the force-
displacement curves that all the three actuation forces take non-zero values
during the deformation process to realize the calculated optimal deformation
path. This motion yields a minimized functional value of $J=600$.

If one load is suspended (Figure 3b), this optimal deformation path cannot be
followed anymore since all three possible point loads are required for keeping
an equilibrium state in the deformed configurations. Therefore, a different
motion is found that is enabled with the remaining forces. Moreover, as only
the vertical displacement of the target geometry is prescribed, the horizontal
displacement adapts such that the additional constraint is met, eventually
resulting in a different final configuration. The minimum of the functional
increases to $J=607$. If only one force, or one dependent couple of forces
$F_{1}=F_{2}$ (see Figure 3c/d) is allowed for actuation, the result is
identical to an equilibrium path obtained by a non-linear analysis with either
the arc length method or a displacement controlled algorithm. The motion does
not represent a “design” in this case because there is no other possible
equilibrium path (neglecting here the theoretical possibility of bifurcation
points and corresponding secondary paths).

In order to verify the applicability of the presented method, numerical
experiments are presented in the next section. Among others, problems from the
[24] are reused to demonstrate the correlation between unconstrained and
constrained motion design with restrictions concerning the actuation forces.

###  3.5 Numerical experiments for motion design with constrained actuation

####  3.5.1 Motion design with multiple snap-through processes

test
First, an example is presented that was already previously studied in [24], in
order to demonstrate the difference between unconstrained and constrained
motion design. Three two-bar trusses are combined into one truss structure.
The members of the upper two-bar truss have a larger cross-sectional area and
therefore higher stiffness than the lower ones, by which it is supported (see
Figure 4a. The vertical displacement of the upper node is prescribed in the
target configuration. The result of an unconstrained motion design procedure
with a path discretization with 16 linear elements is shown in Figure 4b. As
already described in [24], a naive linear interpolation of all displacement
values between the initial and target values as “predictor motion” yields a
high value of $J=440$ for the quality functional. With unconstrained motion
design this reduces to $J=11.1$. Figure 4b visualizes the corresponding motion
along with the required actuation forces. To realize this optimal motion,
actuation forces need to be applied on all six degrees of freedom, which can
be seen in the load-displacement curves (green) in Figure 4d.

Figure 4: Combination of multiple snap-throughs with unconstrained and
constrained motion design.

For a constrained motion design, only two actuation forces, namely a vertical
and a horizontal one at the upper node are allowed. This results in four
constraints per load configuration, which enforce the remaining internal
forces at the uncontrolled degrees of freedom to be zero. Solving the extended
system of equations leads to the modified motion illustrated in Figure 4c. The
end geometries of the unconstrained and constrained motion design variants
differ slightly from each other and from the mirrored geometry, but the
vertical displacement of the upper node coincides in both cases, as this
represents the prescribed target value. Furthermore, it can be seen that
retaining the last equilibrium state requires lower forces in constrained
motion design. Since the members of the upper two-bar truss has a much larger
cross-sectional area than those of the two lower two-bar trusses, the small
deviation from the mirrored geometry has a significant influence on the
magnitude of the forces, recovered from the displacements for the case of
unconstrained motion design. In the constrained motion, however, forces are
allowed to only act on one single node. This reduces the described effect on
the loads attributable to the not perfectly matched mirrored end geometry. Due
to the additional constraints, the value of the functional increases to
$J=12.5$. The two resulting motions from unconstrained and constrained motion
design resemble each other, especially regarding the overall motion pattern
with the consecutive lateral snap-through processes. However, the magnitude of
the required forces differ significantly; the forces required to follow the
constrained motion are much lower (see Figure 4d).

####  3.5.2 Snap-through of a shallow arc

test
The next problem is a two-dimensional plane stress problem in the shape of a
shallow circular arc. It is discretized with displacement-based quadrilateral
finite elements. The target geometry, shown in blue in Figure 5a is defined to
be a circular arc that resembles – but is not – an equilibrium configuration
after snap-through. Only three vertical actuation forces are applied, shown in
red.

Figure 5: Constrained motion design of a shallow arc with three point loads.

Since the target geometry is not in equilibrium, point forces have to be
applied at every degree of freedom for its realization, see Figure 5b. This
contradicts the constraint that the motion ought to be realized with only
three actuation forces. Therefore, the target geometry has to be modified as
explained in Section 3.3. An optimization procedure is carried out, in which
the target displacements are the design variables and the displacement
difference between the prescribed and the current geometry is minimized. The
optimization problem is subject to the constraint that equilibrium is obtained
with the three available point forces only. Figure 5c illustrates that the
resulting geometry differs noticeably from the initially prescribed geometry.

This newly obtained geometry is the target configuration for the motion design
problems solved next. The deformation path is discretized by 24 linear
elements. An unconstrained motion design yields the motion in Figure 5d with a
functional value of $J=7.14$. It turns out that a non-symmetric snap-through
of the arc is beneficial for a low value of $J$. Although the target geometry
can be equilibrated with the three prescribed forces only, non-zero forces are
applied at the other degrees of freedom throughout the motion to realize the
optimal deformation path in unconstrained motion design.

This result is obviously different than the one from constrained motion
design, which is shown in Figure 5e. Here, only the three admissible actuation
forces are applied during the entire deformation process. Due to this
constraint, the functional value increases slightly to $J=7.20$. This
resulting optimal motion could not be followed within a non-linear analyses,
in which the relation between the three different actuation forces is locked.
Thus, even though the constraints increase the cost of deformation, it is
still lower than for the conventional approach of a non-linear analysis. This
is the consequence of enabling an independent evolution of the different
actuation forces throughout the deformation process.

##  4 Discrete actuator elements for internal actuation

###  4.1 Actuator element formulation

As an alternative to applying external forces, actuation of structures can
also be accomplished via actuator elements. Simple, linear actuators in truss
structures can contract or expand and thus enforce a deformation of the entire
structure or parts of it. To include such elements into the motion design
method, a new actuator element formulation is introduced. The basis of the
actuator element is a standard truss element with a vertical and horizontal
displacement degree of freedom at each node (cf. Figure 6 (left)). The
actuator element, however, also allows for a pure and independent elongation
or contraction. This is mapped by an additional parameter ${\alpha}$, which
represents a factor for the targeted actuator element elongation, as can be
seen in Figure 6 (right).

Figure 6: Parameter in the actuator element formulation.

In an actuator element, the total strain $E_{{\mathrm{d}}}$ consists of two
parts: the elastic strain $E_{{\mathrm{el}}}$ and the strain from the
actuation $E_{{\alpha}}$

$\displaystyle E_{{\mathrm{d}}}=E_{{\mathrm{el}}}+E_{\alpha}.$ (21)

The actuation strain $E_{\alpha}$ is determined with the additional elongation
parameter ${\alpha}$ and the total strain $E_{\mathrm{d}}$ is obtained from
the nodal degrees of freedom via the kinematic equations. This allows
evaluation of the elastic strain and eventually the normal stress from the
material law. Applying the linear elastic St. Venant-Kirchhoff material law,
the internal energy is obtained as

$\displaystyle{\Pi}_{{\mathrm{int}},{\mathrm{el}}}=\int_{\Omega}E_{{\mathrm{el}}}^{\mathrm{T}}{\mathbf{C}}E_{{\mathrm{el}}}\,{\mathrm{d}}{\Omega}=\int_{\Omega}(E_{{\mathrm{d}}}-E_{\alpha})^{\mathrm{T}}{\mathbf{C}}(E_{{\mathrm{d}}}-E_{\alpha})\,{\mathrm{d}}{\Omega}.$
(22)

From this, the internal forces can be derived, yielding for a single actuator
element

$\displaystyle{\mathbf{f}}_{{\mathrm{int}}}=\int_{{\Omega}_{e}}{\mathbf{B}}^{\mathrm{T}}{\mathbf{C}}(E_{{\mathrm{d}}}-E_{\alpha})\,{\mathrm{d}}{\Omega}.$
(23)

The strain-displacement matrix ${\mathbf{B}}$ is defined as the derivative of
the elastic strain with respect to the unknown parameters and is divided into
two parts: the derivative with respect to the local displacement degrees of
freedom ${\mathbf{d}}$ and the derivative with respect to the elongation
parameter ${\alpha}$. This results in

$\displaystyle{\mathbf{B}}=\begin{bmatrix}{\mathbf{B}}_{\mathrm{d}}\\\
B_{\alpha}\end{bmatrix}\qquad\text{with}\qquad{\mathbf{B}}_{\mathrm{d}}=\frac{{\partial}E_{\mathrm{d}}}{{\partial}{\mathbf{d}}}\qquad\text{and}\qquad
B_{\alpha}$ $\displaystyle=-\frac{{\partial}E_{\alpha}}{{\partial}{\alpha}}.$
(24)

Consistent linearization of the global residual equation for equilibrium
${\mathbf{R}}={\mathbf{F}}_{\mathrm{int}}-{\mathbf{F}}_{\mathrm{ext}}$ with
the global vector of internal forces ${\mathbf{F}}_{\mathrm{int}}$ and
external forces ${\mathbf{F}}_{\mathrm{ext}}$ yields the stiffness matrix. The
local stiffness matrix for the actuator element follows as

$\displaystyle{\mathbf{k}}_{{\mathrm{act}}}=\begin{bmatrix}{\mathbf{B}}_{{\mathrm{d}},{\mathbf{d}}}^{\mathrm{T}}{\mathbf{C}}E_{\mathrm{el}}+{\mathbf{B}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{C}}{\mathbf{B}}_{\mathrm{d}}&B_{\alpha}{\mathbf{C}}{\mathbf{B}}_{\mathrm{d}}\\\
{\mathbf{B}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{C}}B_{\alpha}&B_{{\alpha},{\alpha}}^{\mathrm{T}}{\mathbf{C}}E_{\mathrm{el}}+B_{\alpha}^{\mathrm{T}}{\mathbf{C}}B_{\alpha}\end{bmatrix}=\begin{bmatrix}{\mathbf{k}}_{{\mathrm{d}}{\mathrm{d}}}&{\mathbf{k}}_{{\mathrm{d}}{\alpha}}\\\
{\mathbf{k}}_{{\mathrm{d}}{\alpha}}^{\mathrm{T}}&{\mathbf{k}}_{{\alpha}{\alpha}}\end{bmatrix}.$
(25)

It includes the usual tangent stiffness matrix
${\mathbf{k}}_{{\mathrm{d}}{\mathrm{d}}}$ of a truss element in the upper left
corner. The global stiffness matrix can be obtained by standard assembly
operations for $n_{\mathrm{ele}}$ elements. These operations can be applied
separately for the different stiffness components

$\displaystyle{\mathbf{K}}_{{\mathrm{d}}{\mathrm{d}}}=\bigcup_{e=1}^{n_{\mathrm{ele}}}{\mathbf{k}}_{{\mathrm{d}}{\mathrm{d}}}\qquad\qquad{\mathbf{K}}_{{\mathrm{d}}{\alpha}}=\bigcup_{e=1}^{n_{\mathrm{ele}}}{\mathbf{k}}_{{\mathrm{d}}{\alpha}}\qquad\qquad{\mathbf{K}}_{{\alpha}{\alpha}}=\bigcup_{e=1}^{n_{\mathrm{ele}}}{\mathbf{k}}_{{\alpha}{\alpha}},$
(26)

which results in the global linearized system of equations for equilibrium

$\displaystyle\begin{bmatrix}{\mathbf{K}}_{{\mathrm{d}}{\mathrm{d}}}&{\mathbf{K}}_{{\mathrm{d}}{\alpha}}\\\
{\mathbf{K}}_{{\mathrm{d}}{\alpha}}^{\mathrm{T}}&{\mathbf{K}}_{{\alpha}{\alpha}}\end{bmatrix}\begin{bmatrix}{\Delta}{\mathbf{d}}\\\
{\Delta}{\bm{\alpha}}\end{bmatrix}=+{\mathbf{F}}_{\mathrm{int}}-{\mathbf{F}}_{\mathrm{ext}}.$
(27)

A prescribed elongation or contraction of an actuator can be treated as
inhomogeneous Dirichlet boundary condition. The procedure can also be
interpreted as an actuation with external forces, extracted from eq. (23) as

$\displaystyle{\mathbf{f}}_{{\mathrm{int}},{\mathrm{act}}}=-\int_{{\Omega}_{e}}{\mathbf{B}}^{\mathrm{T}}{\mathbf{C}}E_{\alpha}\,{\mathrm{d}}{\Omega}.$
(28)

The procedure is demonstrated in the following example, a statically
determinate structure, which was already introduced in the previous section.
Figure 7a shows the actuator element in red. A shortening by 50% is
prescribed. With the system of equations in eq. (27) and a load-controlled
geometrically non-linear analysis, the resulting nodal displacements can be
computed. The solution is shown in Figure 7a. Because the structure is
statically determinate, a stress-free, purely kinematic motion is obtained.

Figure 7b shows the same system with an additional bar, such that it is now
statically indeterminate by degree $n_{\mathrm{s}}=1$. The same actuator
element and the same target length change as before are used. The statical
indeterminacy of the structure in combination with the position of the single
actuator obstructs a contraint-free motion. Stresses are implied to the other
truss elements and the actuator does not reach the targeted length change of
50% when specifying ${\alpha}=-0.5$. The resulting deformed structure is shown
in Figure 7b.

When a second actuator element replaces a regular truss element, a constraint-
free length change as well as a kinematic motion are potentially possible.
However, it is not possible with the given actuator combination and the
prescribed elongation parameters of both ${\alpha}=-0.5$ as displayed in
Figure 7c. The truss elements still perform nearly the kinematic motion.
However, the length changes of both actuators do not exactly match this
specific motion, creating constraints in the actuators and leading to
different length changes than prescribed. This brief study shows that the
number and location of actuators and the degree of static indeterminacy are
crucial for the feasibility of a constraint-free kinematic motion.

Figure 7: Actuator length change as load case and the resulting deformation of
a statically determinate and statically indeterminate structure with one and
two actuator elements.

###  4.2 Motion design with actuator elements

The introduced actuator element formulation allows an easy implementation into
the motion design method. For this purpose, the elongation parameters
${\alpha}_{i}$ for every actuator element $i$ are discretized along the motion
path, as it has been done with the displacement degrees of freedom. Thus, the
elongation parameters ${\alpha}_{i}$, as part of the solution vector, are
directly solved for while optimizing the motion. Again, this is the same
process as for the displacement degrees of freedom. Therefore, the elongation
parameters vary according to the optimized motion and their evolution
throughout the deformation process represents an output of the motion design
method.

This is again illustrated with the truss example from Figure 2, prescribing
the same vertical displacement $D_{2}$ of the second node, using a path
discretization with twenty linear path elements. With the presented actuator
element formulation, the basic motion design method already includes the
actuation load case, i.e., without applying any external forces. In order to
ensure that only the actuator is used to realize the optimized motion without
any discrete point loads, external forces are constrained to be zero for every
degree of freedom (see Chapter 3).

First, the statically determinate structure with one actuator element is
investigated, cf. Figure 8a. Here, a purely kinematic motion of the rest of
the structure is obtained, which can be identified by the functional value
being $J=0.0$. The most important output of motion design is the evolution of
the elongation parameter throughout the deformation process. The actuator
adapts to the optimized motion, as can be seen in the plot of the elongation
parameter versus the displacement $D_{2}$ in Figure 8a (bottom). This can also
be observed in the statically indeterminate structure (degree
$n_{\mathrm{s}}=1$) with two actuator elements in Figure 8b. In this case, a
suitable combination of the elongation parameters of both actuator elements is
found by the motion design method, such that a purely kinematic motion is
obtained. Again, this results in a functional value of $J=0$. The motion is
different from the one in Figure 7, because now the elongation parameters of
both actuators may evolve independently.

Figure 8: Motion design with actuator elements.

Finally, a structure with a second additional truss element is studied (cf.
Figure 8c). This modification raises the degree of static indeterminacy to
$n_{\mathrm{s}}=2$ and leads to a slightly different motion with different
evolutions of both elongation parameters. A significant difference, however,
results from the fact the other bars cannot perform a pure rigid body motion
anymore, but build up internal stresses throughout the motion. This is the
main difference compared to the two previous motion designs and leads to a
functional value of $J=379$.

###  4.3 Bridge structure with a traveling load

To demonstrate the potential of motion design with actuator elements, a more
complex example is presented in the following. A bridge is modeled as a truss
structure, as shown in Figure 9a. Actuator elements are installed in the
entire lower chord, thus resulting in a total of ten actuators. In this
example, a parallel actuation mechanism is employed. This means that in
addition to the actuator element, a regular, passive truss element is
installed at the same place. As a result, this passive element is
automatically stressed when the actuator changes its length. This is supposed
to penalize the actuator extension in order to approximately represent the
corresponding actuation costs. The top chord and the vertical struts have a
cross sectional area of $A_{1}=0.1$, while the diagonals as well as the
actuator elements and the passive truss elements in the lower chord are built
of elements with half the cross section area, i.e., $A_{2}=0.05$.

Figure 9: Bridge structure and deformation due to static loading.

The bridge is exposed to a load $F$ traveling from left to right. The idea is
that the point of load application remains at the same height throughout the
entire travel, as indicated by the horizontal line in Figure 9a. In order to
realize this, the structure has to adapt and undergo large deformations. The
required motion should be solely enabled by the actuators in the lower chord.

The initial configuration of the bridge is the symmetric, unloaded state. To
find a suitable motion for the traveling load and to make the deformation
process as efficient as possible, motion design is carried out in two stages.
The two motion stages are not solved separately, but by one single motion
design procedure. However, they differ in the controlled and prescribed
displacement values, as described in the following.

Stage 1 of the motion represents the procedure of “picking up” of the point
load on the left side. Therefore, the structure needs to deform sideways until
the fourth node of the top chord (marked in blue) arrives at the starting
point of the traveling load. For motion design, the vertical and horizontal
target position of this node are prescribed and its horizontal displacement
value is controlled throughout the motion of the first stage (Configuration 0
to Configuration 5 in Figure 10a).

In Stage 2, the load travels from the left side to the central node of the top
chord. To guarantee that it always remains at the same vertical position, the
vertical displacements of the nodes of the top chord are prescribed in a
sequential way throughout the configuration: First, the vertical displacement
of the fourth node is prescribed, then the vertical displacement of the fourth
and fifth node are prescribed, then the vertical displacement of the fifth
node is prescribed and so forth (cf. Figure 10b, marked in blue). In the
motion design problems described up to now, at least one target displacement
value is prescribed for the motion and one degree of freedom is controlled
throughout the entire motion. This is not the case in Stage 2, but the nodes
with specified vertical displacement values change in every configuration.
However, these specifications still lead to a stable motion design procedure
with a unique result. Furthermore, the traveling load is applied to every node
reaching this specific vertical position. If two nodes are placed at the same
height, the load is split and distributed on both nodes. The application of
these point loads is also included as a constraint on the motion as described
in Section 3.

An alternative approach to design a motion is to find the geometry for every
single configuration separately with the specific constraints, i.e.,
displacement values and loading, by an optimization that minimizes the
internal energy. These optimized configurations can then be connected, thus
generating an entire motion. These connected optimized geometries are then
used as predictor motion for motion design. Figure 10a and Figure 10b show the
optimized motion that results from the solution of the motion design problem
including Stage 1 and Stage 2. Because of the applied constraints on the
admissible point loads, the deformation is solely realized by the actuators.
Furthermore, potential displacements due to the loading are compensated such
that the point load can be kept at exactly the same height during the entire
travel. The motions required to maneuver the load to the right side of the
bridge is obtained by symmetry. Using the motion design method, the functional
value reduces from $J=72.8$ in the predictor (connection of energy-minimal
configurations) to $J=51.6$, i.e., by 29%. This is due to the fact that the
entire motion is considered in the objective function and not only separate
configurations.

Figure 10: Optimized motion of the bridge structure with traveling load.

An output of the method are the targeted actuator length changes in the form
of the elongation parameters, required to realize the optimal motion. Their
evolution is exemplarily illustrated in Figure 10c for the first actuator
(Actuator 1) and the fourth actuator (Actuator 2) of the lower chord (marked
in Configuration 0). It can be seen that the elongation parameters do not
increase monotonously, but develop independently and can also decrease during
the process. This confirms that the method of motion design allows to design
efficient motions that are realized solely by actuator elements.

##  5 Motion Stabilization

###  5.1 Detection of critical points

In addition to the presented constraints regarding admissible actuation forces
for the optimized motion, other constraints can be applied to motion design
problems. One possibility is to enforce a stabilized motion, where no
structural instabilities occur during the deformation. Instability phenomena,
such as bifurcation or snap-through, are usually initiated at critical points,
namely bifurcation points and limit points. Often, but not always, reaching a
critical point is followed by entering an unstable equilibrium path. In quasi-
static problems, the determinant of the tangent stiffness matrix is zero when
a critical point is reached. This simple criterion is used in the following
for identification of critical points.

###  5.2 Inequality constraint for the stiffness determinant

The basic method for motion design again marks the starting point for
constrained motion design, however, now with different constraint types. To
avoid that the structure enters an unstable equilibrium path, an inequality
constraint can be introduced within the framework indicated in eq. (10). In
particular, an inequality constraint can be applied onto the determinant of
the tangent stiffness matrix ${\mathbf{K}}_{\mathrm{T}}$. The constraints are
summarized in a vector for all load configurations $\bar{k}$

$\displaystyle{\mathbf{g}}$
$\displaystyle=\begin{bmatrix}\det{{\mathbf{K}}^{2}}&\det{{\mathbf{K}}^{3}}&\ldots&\det{{\mathbf{K}}^{\bar{k}}}&\ldots&\det{{\mathbf{K}}^{\bar{n}_{\text{node}}}}\end{bmatrix}^{\mathrm{T}}.$
(29)

Thus, each entry $i$ of the vector ${\mathbf{g}}$, i.e., each constraint
$g_{i}$, needs to fulfill the inequality222For a B-spline discretization, the
nodes, where the constraints are enforced, are again the Greville collocation
points on the deformation path itself.

$\displaystyle g_{i}=\det{{\mathbf{K}}_{i}}\geq 0.$ (30)

Applying inequality constraints is realized by activating a corresponding
equality constrained, once the inequality constraint is violated. Hence, an
active set strategy must be used. This is done by the Karush-Kuhn-Tucker
conditions

$\displaystyle g_{i}$ $\displaystyle\geq 0,$ $\displaystyle{\lambda}_{i}$
$\displaystyle\geq 0,$ $\displaystyle g_{i}{\lambda}_{i}=0$ (31)

with the Lagrange multipliers ${\lambda}_{i}$. These conditions can also be
expressed in a semi-smooth way in the complementarity function for every
constraint $i$

$\displaystyle C_{i}={\lambda}_{i}-\max{(0,{\lambda}_{i}-cg_{i})}=0.$ (32)

The complementarity parameter $c$ has to be positive and stabilizes the system
of equations.

Like in Section 3, the derivatives of the constraints have to be computed to
set up the extended system of equations 13. Clearly, this is a non-trivial
task, but it can also be accomplished numerically. Also here, the complex step
differentiation from [17] is used to get the exact derivatives. However, only
the first derivative can be calculated exactly in this way. A possibility to
numerically compute second derivates is to make use of hyper-dual numbers, as
proposed by Fike and Alonso [5]. Even though these numbers enable the
computation of exact second derivates, they also suffer from the drawback that
they are not included in the standard C++ library and are therefore not
accessible for a variety of finite element research software. As an
alternative, the calculation of the Hessian and the associated second
derivatives of the determinant can simple be omitted. This leads to a not
consistently linearized problem and the Newton-Raphson method is transformed
into a modified Newton-Raphson solution scheme. It follows for the resulting
system of equations

$\displaystyle\begin{bmatrix}{\mathbf{K}}_{{\mathrm{m}}{\mathrm{d}}}+\cancel{{\lambda}_{i}{\mathbf{H}}_{i}}&{\mathbf{G}}^{\mathrm{T}}\\\
{\mathbf{G}}&{\mbox{$\mathbf{0}$}}\end{bmatrix}\begin{bmatrix}{\Delta}\bar{{\mathbf{D}}}\\\
{\Delta}{\bm{\lambda}}\end{bmatrix}=-\begin{bmatrix}{\mathbf{R}}_{{\mathrm{m}}{\mathrm{d}}}+{\bm{\lambda}}^{\mathrm{T}}{\mathbf{G}}\\\
{\mathbf{g}}\end{bmatrix}.$ (33)

The incomplete linearization may affect the convergence behavior of the non-
linear problem. It is expected to have an inferior, no longer quadratic
convergence behavior, compared with the complete system of equations including
the Hessian matrix. Nevertheless, the convergence behavior is expected to
remain acceptable, since ${\mathbf{K}}_{{\mathrm{m}}{\mathrm{d}}}$ is still
updated in each iteration. Due to the correct residual, the solution still
converges to the correct minimal solution, but usually with an increased
number of iterations.

To demonstrate the stabilization of motions, two simple two-bar structures,
incorporating snap-through phenomena as well as bifurcation points are
presented in the following.

###  5.3 Snap-through problem in a shallow two-bar truss

The method of stabilized motion design is first investigated on the example of
an unsymmetric two-bar truss, as shown in Figure 11a. Here, no constraints are
enforced regarding the actuation forces. Thus, a horizontal and a vertical
force can be applied onto the free node. This corresponds to an unconstrained
motion design, the result of which is shown in Figure 11b for a path
discretization with ten elements and the vertical displacement of the central
node being controlled. The sign of $\det{\mathbf{K}}_{\mathrm{T}}$ is
indicated in the background of the illustrated motion, where the red shaded
zone marks a region in which $\det{\mathbf{K}}_{\mathrm{T}}<0$. In the
unconstrained motion design solution, the node traverses this unstable region.
This is also visible in the plot of the determinant of the stiffness matrix
versus the displacement $D_{2}$ in Figure 11b. Between vertical displacement
values of approximately $D_{2}\approx 2$ and $D_{2}\approx 8$, the determinant
takes a negative value. For the given load-controlled process, this would
practically mean that an uncontrolled dynamic snap-through occurs instead of
the desired (designed) motion.

Applying constrained motion design, a stabilized motion is obtained (see
Figure 11c). The trajectory of the central node only touches the zone where
the determinant would become negative. In practical applications, one would
probably avoid exactly reaching this critical region, but keep a safe distance
by enforcing a (small) positive value $\det{\mathbf{K}}_{\mathrm{T}}>0$ in the
inequality constraint.

It can also be seen that the magnitude of the forces, especially the
horizontal point load at the midnode, changes in order to follow the
alternative deformation path. The motion design functional itself, i.e., the
minimum of the cost of deformation, stays the same, but its value increases
drastically from $J=1492$ to $J=6211$ due to the necessary detour. However, a
stabilization and a positive determinant are ensured throughout the entire
motion.

Figure 11: Unsymmetric two bar truss with an unconstrained and constrained
stabilized motion design.

###  5.4 Stabilization of a two-bar truss with bifurcation

In order to also investigate the stabilization in bifurcation problems, a high
two-bar truss is studied next (see Figure 12a. When the free node of a
_shallow_ two-bar truss is moved downwards, the first critical point reached
is a limit point, initiating a snap-through process. However, in this case, an
additional bifurcation point occurs before the limit point is reached. This
happens when the height to span ratio exceeds a value of $\sqrt{2}/2$. The
objective of motion design is again to reach a “mirrored” configuration.

This problem is already described in [24], where different solutions and
equilibrium paths are compared. Figure 12b shows the solution from
unconstrained motion design (with twenty linear path elements and the vertical
displacement controlled). Again, the red zone indicates a negative sign of the
determinant of the tangent stiffness matrix $\det{{\mathbf{K}}_{\mathrm{T}}}$.
Again, the unstable region is traversed by the deformation path, which can
also be seen in the progress of the stiffness determinant in the corresponding
diagram underneath. Thus, this result would be practically infeasible.

Applying the inequality constraint $\det{{\mathbf{K}}_{\mathrm{T}}\geq 0}$
results in the motion that can be seen in Figure 12c, where the red zone is
again only touched by the free node, but not traversed. The constraint of
enforcing a positive stiffness determinant can therefore also be successfully
applied to problems where bifurcation phenomena can occur.

Figure 12: High two-bar truss with a bifurcation point and an unstabilized and
stabilized motion.

###  5.5 Stabilization of a motion with restricted actuation

So far, the two constraint types of enabling an optimized motion by a
restricted set of actuation forces (equality constraints) and of stabilizing
the motion by enforcing a positive stiffness determinant (inequality
constraints) were treated separately. However, they can also be combined. The
vector ${\mathbf{g}}$ that gathers all constraints then contains entries with
the internal forces ${\mathbf{g}}_{F}$ as well as entries with the determinant
of the stiffness matrix ${\mathbf{g}}_{\det}$,

$\displaystyle{\mathbf{g}}=\begin{bmatrix}{\mathbf{g}}_{F}\\\
{\mathbf{g}}_{\det}\end{bmatrix}\,.$ (34)

As a consequence, the resulting system of equations is further expanded. This
is shown in the following with the system from Figure 2. The deformation path
is discretized by 16 linear elements (only eight are illustrated in Figure 3)
and the vertical displacement is controlled. It is already known at this point
that unconstrained motion design yields the functional value $J=600$, as shown
in Figure 3a. This represents the reference value for the following
modifications.

First, the horizontal force at the free node is suspended and a constrained
motion design is carried out. This leads to the motion in Figure 13a. Due to
the constraint, the value of the functional increases to $J=616$. Secondly,
all possible forces may again be applied, but a positive determinant of the
stiffness matrix is enforced to find a stabilized motion. Again, a different
motion is found, shown in Figure 13b, and the functional value increases
significantly to $J=1542$. In this structure, the horizontal bar is activated
much more in the form of a high compressive force to reduce the snap-through
effect of the two other bars. Finally, both constraints can be combined. The
result is a stable motion that always keeps a positive determinant of the
stiffness matrix, but is also realizable with only two point forces, as
illustrated in Figure 13c.

Figure 13: Constrained motion design on an exemplary system with different
types of constraints as well as their combination.

##  6 Conclusions

In this paper, the work presented in Sachse and Bischoff [24] is continued and
extended to a more practical approach. While in the previous study, only
unconstrained motions were investigated, this work presents the introduction
of constraints into the motion design process. This allows finding an optimal
motion between two prescribed geometries that is either realized by only a
specific number of load cases, actuator elements, or is stabilized by
enforcing a positive structure’s stiffness determinant.

After the basics of the unconstrained motion design method were summarized,
first, the introduction of load cases to motion design was presented. This was
introduced by enforcing a zero internal force at every unloaded degree of
freedom and can exemplarily be conducted with the Lagrange multiplier method.
Here, the second derivatives of the constraints need to be calculated to build
the resulting extended system of equations. The first derivative is
represented by the stiffness matrix, whereas the second derivative, the
derivative of the stiffness, can be obtained by a directional derivative. Due
to the particular (symmetric) character of the stiffness matrix, this
directional derivative allows for using only one more evaluation of the
function when applying numerical differentiation, increasing the numerical
efficiency of the calculation. Furthermore, additional considerations have to
be made for the definition of the desired end geometry. This approach is
further demonstrated in some numerical experiments with different types of
finite elements. They showed that the functional value of constrained
optimized motions increases compared to unconstrained optimized motions.
However, the constraints can be satisfied.

The length change of incorporated actuator elements in the structure can also
be regarded as a load case. Therefore, a new actuator element formulation with
an additional parameter, the targeted actuator length change ${\alpha}$, was
developed. This element formulation allows a straightforward implementation
into the presented motion design method because it is simply solved for this
additional parameter as well. By application of the actuator element,
kinematic motion of the remaining passive structure and the corresponding
actuator length changes can be determined. In addition, the potential of this
method with actuator elements is shown on the example of a bridge structure
with a traveling load.

Next, the constraint of a stabilized motion is introduced. In this case, the
fact that the determinant of the stiffness matrix of the structure becomes
negative when passing a critical point is utilized. The inequality condition
of always keeping a positive determinant is also enforced with the Lagrange
multiplier method while following an active set strategy. However, the
required derivatives are difficult to obtain. This is why, again, numerical
differentiation is used. To improve efficiency, only the first derivatives of
the constraints are calculated in this way leading to a modified Newton
scheme. Due to the type of constraint, where a determinant needs to be
calculated and differentiated, the solution process might be inefficient in
terms of the calculation time for large structures. Furthermore, the solution
of the resulting extended system of equations might exhibit numerical
instabilities as the entries are characterized by a significant difference in
magnitude, which is caused by the determinant’s value. However, in this study,
only the principle of the introduction of such a constraint is shown and its
functionality proven with simple structures. When the presented approach is
applied in the presented numerical examples, unstable deformation states are
avoided in snap-through problems as well as bifurcation problems. The
described numerical problems still need to be tackled in further research if
following this idea.

Finally, both type of constraints, enabling the designed motion only by a
prescribed number of loads as well as the motion stabilization, is combined
and presented in an exemplary structure. Similarly, other inequality
constraints on stress measures are thinkable. This way, it can be ensured that
they do not exceed any strength values or limits regarding local instability
behavior, i.e., buckling of single truss elements, during the entire
deformation process. For such constraints, usually only existing quantities in
a nonlinear finite element code, such as the first and second derivatives of
the stress tensor, are required and analytical sensitivities without the need
for numerical differentiation can be obtained.

The basic concept and structure of unconstrained and constrained motion design
is presented with an exemplary functional, the integrated internal energy over
the deformation path, and exemplary constraints in the previous publication by
the authors [24] as well as in this paper. The herein described constraints
should represent an approach for technical implementation and application of
the motion design method for real structures. However, also other functionals
and constraints may be introduced to adapt the method for the specific
requirements of the underlying task.

Besides, the type and number of loads were always defined prior to motion
design. One possibility for further research is to include an optimization
procedure also to find the optimal loads to trace the optimized deformation
path by motion design. This resembles actor placement procedures or, in this
case, load case placement procedures, that already exist as described in the
introductory section 1. Furthermore, the focus was laid on the application of
point forces up to now, but also other types of actuation are thinkable such
as discrete actuator elements or piezoelectric structural parts and may be
implemented.

## Acknowledgements

This research work by Renate Sachse was funded and supported by the State
Ministry of Baden-Wuerttemberg for Sciences, Research and Arts. Furthermore,
the contribution of Florian Geiger was conducted in the framework of the
Collaborative Research Centre 1244 “Adaptive Skins and Structures for the
Built Environment of Tomorrow”/project B01 funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) - project number
279064222. The authors are grateful for the generous support.

## References

  * [1] Makola M. Abdullah, Andy Richardson, and Jameel Hanif.  Placement of sensors/actuators on civil structures using genetic algorithms.  Earthquake Engineering & Structural Dynamics, 30(8):1167–1184, 2001.
  * [2] Rafic M. Ajaj, Christopher S. Beaverstock, and Michael I. Friswell.  Morphing aircraft: The need for a new design philosophy.  Aerospace Science and Technology, 49:154–166, February 2016.
  * [3] L.F. Campanile.  Shape-adaptive wings—the unfulfilled dream of flight.  In R. Liebe, editor, WIT Transactions on State of the Art in Science and Engineering, volume 2, pages 400–419. WIT Press, 1 edition, November 2006.
  * [4] L.F. Campanile.  Modal Synthesis of Flexible Mechanisms for Airfoil Shape Control.  Journal of Intelligent Material Systems and Structures, 19(7):779–789, July 2008.
  * [5] Jeffrey Fike and Juan Alonso.  The Development of Hyper-Dual Numbers for Exact Second-Derivative Calculations.  In 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Orlando, Florida, January 2011. American Institute of Aeronautics and Astronautics.
  * [6] Knut Göppert and Michael Stein.  A Spoked Wheel Structure for the World’s largest Convertible Roof – The New Commerzbank Arena in Frankfurt, Germany.  Structural Engineering International, 17(4):282–287, 2007.
  * [7] Albert Graells Rovira and Josep M. Mirats Tur.  Control and simulation of a tensegrity-based mobile robot.  Robotics and Autonomous Systems, 57(5):526–535, 2009.
  * [8] Vivek Gupta, Manu Sharma, and Nagesh Thakur.  Optimization Criteria for Optimal Placement of Piezoelectric Sensors and Actuators on a Smart Structure: A Technical Review.  Journal of Intelligent Material Systems and Structures, 21(12):1227–1243, August 2010.
  * [9] G. W. Housner, L. A. Bergman, T. K. Caughey, A. G. Chassiakos, R. O. Claus, S. F. Masri, R. E. Skelton, T. T. Soong, B. F. Spencer, and J. T. P. Yao.  Structural Control: Past, Present, and Future.  Journal of Engineering Mechanics, 123(9):897–971, September 1997\.
  * [10] A. Ibrahimbegovic, C. Knopf‐Lenoir, A. Kučerová, and P. Villon.  Optimal design and optimal control of structures undergoing finite rotations and elastic deformations.  International Journal for Numerical Methods in Engineering, 61(14):2428–2460, 2004.
  * [11] Fumihiro Inoue.  Development of Adaptive Construction Structure by Variable Geometry Truss.  In Carlos Balaguer and Mohamed Abderrahim, editors, Robotics and Automation in Construction. InTech, October 2008.
  * [12] H. Irschik.  A review on static and dynamic shape control of structures by piezoelectric actuation.  Engineering Structures, 24(1):5–11, January 2002.
  * [13] Sinan Korkmaz.  A review of active structural control: challenges for engineering informatics.  Computers & Structures, 89(23):2113–2132, 2011.
  * [14] Sridhar Kota, Jinyong Joo, Zhe Li, Steven M. Rodgers, and Jeff Sniegowski.  Design of Compliant Mechanisms: Applications to MEMS.  Analog Integrated Circuits and Signal Processing, 29(1):7–15, October 2001.
  * [15] J Lienhard, S Schleicher, S Poppinga, T Masselter, M Milwich, T Speck, and J Knippers.  Flectofin: a hingeless flapping mechanism inspired by nature.  Bioinspiration & Biomimetics, 6(4):045001, 2011.
  * [16] Kerr-Jia Lu and Sridhar Kota.  An Effective Method of Synthesizing Compliant Adaptive Structures using Load Path Representation.  Journal of Intelligent Material Systems and Structures, 16(4):307–317, April 2005.
  * [17] Joaquim R. R. A. Martins, Peter Sturdza, and Juan J. Alonso.  The Complex-step Derivative Approximation.  ACM Trans. Math. Softw., 29(3):245–262, September 2003.
  * [18] Helmut Masching and Kai-Uwe Bletzinger.  Parameter free structural optimization applied to the shape optimization of smart structures.  Finite Elements in Analysis and Design, 111:33–45, April 2016.
  * [19] Milenko Masic and Robert E. Skelton.  Path Planning and Open-Loop Shape Control of Modular Tensegrity Structures.  Journal of Guidance, Control, and Dynamics, 28(3):421–430, May 2005\.
  * [20] Kurt K. Maute and Gregory W. Reich.  Integrated Multidisciplinary Topology Optimization Approach to Adaptive Wing Design.  Journal of Aircraft, 43(1):253–263, January 2006.
  * [21] M Pagitz and J Bold.  Shape-changing shell-like structures.  Bioinspiration & Biomimetics, 8(1):016010, February 2013.
  * [22] Arka P. Reksowardojo, Gennaro Senatore, and Ian F. C. Smith.  Actuator Layout Optimization for Adaptive Structures Performing Large Shape Changes.  In Ian F. C. Smith and Bernd Domer, editors, Advanced Computing Strategies for Engineering, Lecture Notes in Computer Science, pages 111–129. Springer International Publishing, 2018.
  * [23] Daniela Rus and Michael T. Tolley.  Design, fabrication and control of soft robots.  Nature, 521(7553):467–475, May 2015.
  * [24] Renate Sachse and Manfred Bischoff.  A variational formulation for motion design of adaptive compliant structures.  International Journal for Numerical Methods in Engineering, Available online at: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6570 (Accessed December 5,2020), 2020.
  * [25] Gennaro Senatore, Philippe Duffour, and Peter Winslow.  Synthesis of minimum energy adaptive structures.  Structural and Multidisciplinary Optimization, March 2019.
  * [26] Ole Sigmund.  On the Design of Compliant Mechanisms Using Topology Optimization.  Mechanics of Structures and Machines, 25(4):493–524, January 1997\.
  * [27] Werner Sobek and Patrick Teuffel.  Adaptive systems in architecture and structural engineering.  In Smart Structures and Materials 2001: Smart Systems for Bridges, Structures, and Highways, volume 4330, pages 36–46. International Society for Optics and Photonics, 2001.
  * [28] Spencer B. F. and Nagarajaiah S.  State of the Art of Structural Control.  Journal of Structural Engineering, 129(7):845–856, July 2003.
  * [29] Ann C. Sychterz and Ian F. C. Smith.  Deployment and Shape Change of a Tensegrity Structure Using Path-Planning and Feedback Control.  Frontiers in Built Environment, 4, 2018.
  * [30] Srinivas Vasista, Liyong Tong, and K. C. Wong.  Realization of Morphing Wings: A Multidisciplinary Challenge.  Journal of Aircraft, 49(1):11–28, January 2012.
  * [31] Nicolas Veuve, Ann C. Sychterz, and Ian F.C. Smith.  Adaptive control of a deployable tensegrity structure.  Engineering Structures, 152:14–23, December 2017.
  * [32] Terrence A. Weisshaar.  Morphing Aircraft Systems: Historical Perspectives and Future Challenges.  Journal of Aircraft, 50(2):337–353, March 2013.
  * [33] J. van de Wijdeven and B. de Jager.  Shape change of tensegrity structures: design and control.  In Proceedings of the 2005, American Control Conference, 2005., pages 2522–2527 vol. 4, June 2005.