Optimally swimming Stokesian robots

We study self propelled stokesian robots composed of assemblies of balls, in dimensions 2 and 3, and prove that they are able to control their position and orientation. This is a result of controllability, and its proof relies on applying Chow's theorem in an analytic framework, similarly to what has been done in [3] for an axisymmetric system swimming along the axis of symmetry. However, we simplify drastically the analyticity result given in [3] and apply it to a situation where more complex swimmers move either in a plane or in three-dimensional space, hence experiencing also rotations. We then focus our attention on energetically optimal strokes, which we are able to compute numerically. Some examples of computed optimal strokes are discussed in detail.


Introduction
Self propulsion at low Reynolds number is a problem of considerable biological and biomedical relevance which has also great appeal from the point of view of fundamental science. Starting from the pioneering work by Taylor [25] and Lighthill [18], it has received a lot of attention in recent years (see the encyclopedia article [4] for an elementary introduction and the review paper [16] for a comprehensive list of references).
Both relevance for applications and theoretical interest stem from the fact that one considers swimmers of small size. Reynolds number Re = LV /ν gives an estimate for the relative importance of intertial to viscous forces for an object of size L moving at speed V through a newtonian fluid with kinematic viscosity ν. Since in all applications V rarely exceeds a few body lengths per second, if one considers swimming in a given medium, say, water, then Re is entirely controlled by L. At small L, inertial forces are negligible and, in order to move, micro-swimmers can only exploit the viscous resistance of the surrounding fluid. The subtle consequences of this fact (which are rather paradoxical when compared to the intuition we can gain from our own swimming experience) are discussed in [23]. For example, the motion of microswimmers is geometric: the trajectory of a low Re swimmer is entirely determined by the sequence of shapes that the swimmer assumes. Doubling the rate of shape changes simply doubles the speed at which the same trajectory is traversed. As observed in [24], this suggests that there must be a natural, attractive mathematical framework for this problem (which the authors, indeed, unveil). On the other hand, bacteria and unicellular organisms are of micron size, while artificial robots to be used non-invasively inside human bodies for medical purposes must be small. Discovering the secrets of biological micro-swimmers and controlling engineered micro-robots requires a quantitative understanding of low Re self-propulsion.
The basic problem of swimming is easy to state: given a (periodic) time history of shapes of a swimmer (a sequence of strokes), determine the corresponding time history of positions and orientations in space. A natural, related question is the following: starting from a given position and orientation, can the swimmer achieve any prescribed position and orientation by performing a suitable sequence of strokes? This is a question of controllability. The peculiarity of low Re swimming is that, since inertia is negligible, reciprocal shape changes lead to no net motion, so the question of controllability may become non trivial for swimmers that have only a few degrees of freedom at their disposal to vary their shape. The well known scallop theorem [23] is precisely a result of non-controllability.
Once controllability is known, i.e., it is shown that it is possible to go from A to B, one can ask the question of how to go from A to B at minimal energetic cost. This is a question of optimal control.
In spite of the clear connections between low Re self propulsion and control theory, this viewpoint has started to emerge only recently, and mostly in the mathematical literature. Examples are [15], and the more recent contributions [14], [19], and [9]. The papers [3], [4], [2] study in detail both controllability and optimal control for axisymmetric swimmers whose varying shapes are described by few (in particular, two) scalar parameters.
In this paper we analyze the problem of low Reynolds number swimming from the point of view of geometric control theory, and focus on a special class of model systems: those obtained as assemblies of a small number of balls. These systems offer an interesting balance between complexity of the analysis and richness of observable behavior. In addition, they enable us to resolve hydrodynamic interactions numerically by using a fast implementation of the boundary integral method, in which integration points are chosen based on Fibonacci numbers.
Our approach is similar in spirit to the one in [3], [4], [2], but we extend it to non-axisymmetric systems such as three spheres moving in a plane, and systems of four spheres moving in space. The motion of these systems is described by both positional and orientational variables, leading to a much richer geometric structure of the state space. For all the model swimmers described above, controllability is proved by using two main ingredients. The first is Chow's theorem, leading to local controllability in a neighborhood of a point X in state space. We verify the full rank (3.7) hypothesis of this theorem by showing that the vector fields of the coefficients of the governing ODEs (a linear control system without drift) and their first order Lie brackets span the whole tangent space to the state space at X. The second is the analyticity of the coefficients, and the fact that our shape space is always connected: this allows us to use a classical result by Hermann and Nagano (see e.g. [13]). A similar analyticity result has been already proved in [3], but we exhibit here a much simpler proof, which is also easier to generalize to more complex situations. In the same paper [3], Chow's theorem was also used to prove local controllability. In this simpler axisymmetric case, however, position is described only by one scalar parameter (the position p of one distinguished point along the axis of symmetry). The non-degeneracy condition reduces then to the non-vanishing of a single scalar quantity, namely, the curl of the vector field V, giving the rate of change of position as a consequence of shape changes at rates ξ i according toṗ = V 1 (ξ)ξ 1 + V 2 (ξ)ξ 2 . In the richer context of this paper, which involves also rotational degrees of freedom of the swimmers, proving controllability requires an explicit computation of all the first order Lie Brackets. In fact, all the systems we analyze here satisfy the condition where M is the number of controls (rate of change of shape variables) and d is the dimension of the state space (position, orientation, and shape). This is the necessary condition that first order Lie brackets alone suffice to show that the Lie algebra of the coefficients has full rank, so that controllability follows.
For controllable systems, it makes sense to confront the issue of how to achieve the desired target (position and orientation) at minimal energy cost. We present a method to address this optimal control question numerically, and we then examine in detail several examples of optimal strokes in a concrete example (three balls swimming in a plane with a prescribed lateral displacement). Depending on whether the final orientation is also prescribed, and whether the initial shape is prescribed as well or rather one treats it as a parameter to be optimized, we obtain dramatically different answers. Their variety illustrates the suprisingly richness of behavior of low Re swimmers.
The rest of the paper is organized as follows. In Section 2 we recall some results on Stokes flows and in Section 3 we give a precise mathematical statement of the swimming problem. In this latter section we also discuss in detail the three concrete model swimmers to which our analysis applies. Section 4 is devoted to the proof of the controllability of these model systems. In Sections 5 and 6 we state the optimal control problem and a numerical strategy for its solution.
Examples of optimal strokes for three balls swimming in a plane are discussed in detail in Section 7.

Mathematical setting of the problem
We investigate the problem where the swimmer is composed of N identical and non intersecting balls B 1 , · · · , B N of radius a > 0 not necessarily aligned.
We call x i ∈ R 3 the center of the ball B i , and Ω = R 3 \ ∪ N i=1 B i the unbounded domain filled by the fluid. Since the balls do not intersect, the vector x = (x 1 , · · · , x N ) ∈ Ω is restricted to belong to We work at low Reynolds number, so that the fluid obeys Stokes equations where (u, p) are respectively the velocity and the pressure of the fluid, η its viscosity, σ := η(∇u + ∇u t ) − pId is the Cauchy stress tensor and n is the outer unit normal to ∂Ω (hence, n points from the fluid to the interior of the balls). Existence and uniqueness of a solution to (2.1) is classical in the Hilbert Assuming that the force field belongs to H −1/2 (∂Ω) 1 , the solution of (2.1) can be expressed in terms of the associated Green's function, namely the stokeslet G(r) := 1 8πη where f is a distribution of forces on ∂Ω. The Neumann-to-Dirichlet map is a one to one mapping onto while, by the open mapping Theorem, its inverse (the Dirichlet-to Neumann map) T −1 is continuous.
Calling B = B(O, a) the ball of radius a centered in the origin O, the special structure of our domain allows us to identify ∂Ω with (∂B) N in such a way that T can be seen as a map in which H s N stands for (H s (∂B)) N and the dependence on x ∈ S N has been emphasized, defined by (2.6) 1 Here, and in all the paper the symbol H s denotes the Sobolev space of order s.

Analyticity of the Dirichlet-to-Neumann map
This section is devoted to the following lemma which will be essential in the proof of the controllability theorem. We denote by L(E, F ) the Banach space of linear maps from E to F endowed with its usual norm.
The notion of analyticity in a Banach space is classical: it means that at all pointx x 0 ∈ S N , T x is equal to its Taylor series which converges in L(H Since the first term does not depend on x, we only need to prove the analyticity of Ψ : which follows from the analyticity of G in R 3 \ {0} since this formula does not involve the singular part of the Green's kernel G. 2 We have the following consequence.

Corollary 2
Since T x is an isomorphism for every x ∈ S N , the mapping x −→ T −1 x is also analytic from S N to L(H 1/2 Proof For x ∈ S N and for |y| small enough, writing the Taylor series of T x+y as where T x,q (y) denotes the q−homogeneous term in y of T x+y , we have where the last expression may be rearranged as a converging Taylor series.

Approximation for large distances
We will also use the asymptotic behavior of T −1 x as δ(x) := min i<j |x i − x j | tends to ∞. We introduce the notation for i, j ∈ {1, · · · , N } and, for every Proposition 3 For every f ∈ H . From (2.6) we have for i = 1, · · · , N and r ∈ ∂B, Using the expansion (valid uniformly for bounded z) Substituting into (2.9), and using (2.2), we obtain (2.7). For (2.8) we remark that (2.7) is the asymptotic expansion in x of the analytic operator T x at infinity. Since T is invertible, we obtain the expansion (2.8) by inverting (2.7). 2 In the sequel, we will use furthermore the fact that the spheres are non deformable and may thus only move following a rigid body motion. In this case, the velocity of each point r of the i−th sphere is given by (2.10) Using (2.8) we find the expansion of the total force that the system applies to the fluid where S ij is the matrix defined by while the total torque (with respect to the origin O) is given by The swimmers We will focus our attention on some special swimmers. Namely, we assume that the swimmer is composed of N balls, and we restrict ourselves to configurations which can be described by two sets of variables: • the shape variables, denoted by ξ ∈ S, where S is an open connected subset of R M , from which relative distances (x ij ) 1≤i,j≤N between the balls (B i ) 1≤i≤N are obtained, • the position variables, denoted by p ∈ P, which describe the global position and orientation in space of the swimmer. In our examples, the set P is typically a manifold of dimension less than or equal to six. The six-dimensional case is given by p ∈ P = R 3 × SO 3 and p consists of a translation and a rotation in the three-dimensional space.
We also assume that the orientation of the balls (B i ) 1≤i≤N and the distances (x ij ) 1≤i,j≤N depend analytically on (ξ, p), therefore the state of the system is analytically and uniquely determined by the variables (ξ, p), so that there exist N (analytic) functions which give the position of the current point of the i−th sphere of the swimmer in the state (ξ, p) in S × P. The non-slip boundary condition on ∂B i imposes that the velocity of the fluid is given by

Self propulsion
The question we want to address is whether it is possible to control the state of the system (i.e. both ξ and p) using as controls only the rate of shape changeṡ ξ. For this aim, we need to understand the way p varies when one changes ξ. This is done assuming that the swimmer is self propelled, and that the swimmer's inertia is negligible, which implies that the total viscous force and torque exerted by the surrounding fluid on the swimmer must vanish. As we shall see, using these conditions,ṗ is uniquely determined byξ. The condition that total viscous force and torque vanish is written as which lead, using (3.1), to a set of linear equations which linkṗ toξ. In all the examples we have in mind (see below), this allows us to solveṗ uniquely and linearly in terms ofξṗ In the preceding equation (V i (ξ, p)) 1≤i≤M are M vector fields defined on T P, the tangent bundle of P.
The state of the system (ξ, p) thus follows the ODE

Geometric control theory
When we use (ξ i ) 1≤i≤M as control variables, systems like (3.5) belong to the class of linear control systems with analytic coefficients of the form where X is a point of a d−dimensional manifold M, and (F i ) 1≤i≤M are p vector fields defined on T M, the tangent bundle of M. For such systems, the tools of geometric control theory can be applied. In our case, M and the vector fields F i are furthermore analytic, which enables us to use stronger results.
The theory for such systems may be found in [1,13], for instance.
In (3.6), the control which governs the system is the vector (α 1 , · · · , α p ), and the question of controllability of the system can be stated as follows: A classical result by Chow states that the system (3.6) is locally controllable near X 0 ∈ M (which means that question 4 possesses a positive answer for any X 1 in a suitable neighborhood of X 0 ) if the Lie algebra generated by Here, the Lie algebra Lie(F 1 , · · · , F M )(X 0 ) is the algebra obtained from the vector fields (F 1 , · · · , F M ) by successive applications of the Lie bracket operation defined as The rigorous proof of Chow's theorem may be found in [1,13] but we would like to give here a hint of why such a result holds. Given a vector field F : The main idea is that one can reach from X 0 points in the direction g = This formula shows that the global displacement of the solution of the dynamical system is proportional to time for these directions.
A more subtle move allows to reach points in the direction [F i , F j ]. Namely, we now need to consider which also shows that in time t, one can reach a displacement in the desired direction which is proportional to t 2 . Lie bracket directions of higher order are also attainable at the price of a higher order expression in t, leading to smaller displacements. If the Lie algebra has a full rank, every direction in M is attainable and the system is locally controllable.
As a special case in our applications, we will focus on systems for which the full rank condition (3.7) is satisfied with first order Lie brackets only, and a necessary condition for Our three systems are specially designed in such a way that they satisfy M + M (M − 1) 2 = d, and local controllability near X 0 follows from verifying that We then pass from local to global controllability by using the fact that our systems are, in addition, analytic. In this case, the following stronger result holds. In view of the preceding result, since in our case M is always connected, we simply need to prove local controllability at one point (ξ 0 , p 0 ) in order to obtain global controllability. The remaining of this section consists in considering systems of increasing complexity which all fit into the preceding framework.
For each case, we describe the system, explain how the model can be written in the form (3.5), with analytic vector fields F i , and show that (3.9) holds at one point.

The three sphere swimmer of Najafi and Golestanian (3S)
This swimmer, initially proposed in [20] has been studied thoroughly in [3] and [5]. It is composed of 3 spheres of radius a > 0 aligned along the x−axis, as depicted in Fig. 1.
We call ξ = (ξ 1 , ξ 2 ) the length of the arms, and p the position of the central sphere which leads us to consider the shape set S = (2a, +∞) 2 and the position set P = R. Using e 1 = (1, 0, 0) T , we can write and Due to axial symmetry, the only nontrivial equation in (3.2) is the one concerning the component of the total viscous force along the axis of symmetry. This reads and, similarly, while the coefficient ofṗ is simply given by Positivity of the viscous drag force λ 0 arising from a rigid translation at unit speed is a consequence of the positive-definiteness of the Oseen resistance matrix, see [11]. The coefficients above are independent of p because of translational invariance of the problem. From the analyticity of the coefficients and T x in x, we have thatṗ where α 1 (t) and α 2 (t) are the control functions.
The controllability condition (3.9) becomes which is enough to verify for one configuration. For this purpose, and in order to further simplify the computation, we make the observation that, due to symmetry Therefore, if we consider a symmetric shape ξ = (ζ, ζ), showing that condition (3.10) holds at ξ is equivalent to proving that Using the expansion for large distances of the total force (2.11) with u 1 = (ṗ −ξ 1 )e 1 , u 2 =ṗe 1 , and u 3 = (ṗ +ξ 2 )e 1 , leads to which does not vanish for ζ sufficiently large. This proves the global controllability of the system. This swimmer has been proposed and studied in [17]. It is composed of three balls B 1 , B 2 , B 3 of equal radii a > 0. We assume that the three balls can move along three horizontal axes that mutually meet at c with angle 2π 3 (see Fig. 2). We further assume that the balls do not rotate around their axes so that the shape of the swimmer is characterized by the three lengths ξ 1 , ξ 2 , ξ 3 of its arms, measured from the origin to the center of each ball. However, the swimmer may rotate around c in the horizontal plane.
We thus consider (S 1 , S 2 , S 3 ) a reference equilateral triangle with center O ∈ R 3 in the horizontal plane (O, x, y) such that dist(c, S i ) = 1 and we call t i = cS i . Position and orientation in the horizontal plane are described by the coordinates of the center c ∈ R 3 (but c stays confined to the horizontal plane) and the horizontal angle α that one arm -say arm number 1 -makes with a fixed direction -say (O, x) -, in such a way that d = 3. Therefore, we place the center of the ball B i at where R θ stands for the horizontal rotation of angle θ given for instance by the matrix: The swimmer is then fully described by the parameters X = (ξ, c, θ) ∈ S × P,

14
where S := (0, +∞) 3 and P = R 2 × R, and the functions X i are now defined as X i (ξ, c, α, r) = c + R θ (ξ i t i + r) ∀i ∈ {1, 2, 3} , Notice that the functions X i are still analytic in (ξ, c, θ), and we use them to compute the instantaneous velocity on the sphere B where e 3 is the vertical unit vector.
Due to the symmetries of the system, the swimmer stays in the horizontal plane. Therefore, the third component of the total force F 3 vanishes, while only the third component of the total torque T 3 might not vanish. The constraints of self-propulsion (3.2) and (3.3) -more precisely, the vanishing of the components F 1 , and F 2 of the force, and of the third component T 3 of the torque -lead to a linear system that can be written as where neither R, nor V i depend on c or θ because of translational and rotational invariance of the Stokes problem.
The complete system reads where V(ξ) is the 3 × 3 matrix whose columns are V i (ξ), and α(t) ∈ R 3 are the controls.
Calling F i = −R −1 V i , the special form of our system (in particular the fact that the coefficients do not depend on the position variables) allows us to rewrite the sufficient condition (3.9) for controllability at (ξ, c, θ) as We now turn to the computation of ∂F k ∂ξ l (ξ) for k = l or at least its asymptotic expansion in terms of negative powers of ζ.
From the definition of F k , one has and we need to compute all the factors of this expression atξ.
The first step consists in computing an expansion for R(ξ). We get from (2.11) and (2.13), and taking N = 3, and (3.14) Differentiating R with respect to ξ l gives (notice that from Section 2, R is analytic in ξ and we may therefore differentiate its expansion term by term) where S ij , originally given by (2.12), is now seen as its restriction to the horizontal plane, and the vector t l × e 3 as a horizontal vector. Similarly (here we have used the fact that, for a symmetric configuration, the displacement of one sphere induces no torque) and (3.18) From these expressions, and using (3.12), and the identity ∂S il ∂ξ l = 1 r 2 il (−(e il · t l )Id + t l ⊗ e il + e il ⊗ t l − 3(e il · t l )e il ⊗ e il ) one gets (after some calculations) The determinant ∆ given by (3.11) is thus equal to which does not vanish for ζ large enough. 2 We emphasize again that the global controllability result proven here means that the swimmer, by using suitable controls, is capable of moving anywhere in the horizontal plane and of rotating around the vertical axis by any angle.

The four sphere swimmer moving in space (4S)
We now turn to the more difficult situation of a swimmer able to move in the whole three dimensional space and rotate in any direction. In this case, we fix N = 4 and we consider a regular reference tetrahedron (S 1 , S 2 , S 3 , S 4 ) with center O ∈ R 3 such that dist(O, S i ) = 1 and as before, we call t i = OS i for i = 1, 2, 3, 4.
The position and orientation in the three dimensional space of the tetrahedron are described by the coordinates of the center c ∈ R 3 and a rotation R ∈ SO 3 , in such a way that d = 6. The four ball cluster is now completely described by the list of parameters X = (ξ, c, θ) ∈ S × P, where S := (0, +∞) 4 and P = R 3 × SO 3 .
We place the center of the ball B i at x i = c+ξ i Rt i with ξ i > 0 for i = 1, 2, 3, 4 as depicted in Fig. 3 and forbid possible rotation of the spheres around the axes. A global rotation (R = Id) of the swimmer is however allowed.
Furthermore, the function X i are now defined as which is still analytic in (ξ, c, R), from which we compute the instantaneous velocity on the sphere B i v i = ∂X i ∂t (ξ, c, R, r) =ċ + ω × (ξ i t i + r) + Rt iξi . where ω is the angular velocity. Constraints of self propulsion (3.2) and (3.3) now become where the 6 × 6 matrix known as the resistance matrix [11] is symmetric positive definite, and depends analytically on (ξ, c, R) as the vectors V i (ξ, c, R).
In this case, the complete system reads where V(ξ, c, R) is the 6×4 matrix whose columns are V i (ξ, c, R), and α(t) ∈ R 4 are the controls.
We call F i = −R −1 V i . Notice that for the symmetric configurationξ, the vectors F i have no torque components, and that F i does not depend on c because of the translational invariance of the problem. Taking into account these facts, the sufficient condition (3.9) for controllability at (ξ, 0, Id) simplifies to We now turn to the computation of ∂F k ∂ξ l (ξ, 0, Id) for k = l, or at least of its asymptotic expansion in terms of negative powers of ζ.
The definition of F k still gives and we need to compute all the factors of this expression at (ξ, 0, Id).
The first step consists in computing an expansion of R(ξ, 0, Id). We have from (2.11) and (2.13) Id for a regular tetrahedron) and (3.24) Differentiating R with respect to ξ l gives (notice that from section 2, R is analytic in ξ and we may therefore differentiate its expansion term by term) (3.25) Similarly (here we have used the fact that for a symmetric configuration, the displacement of one sphere induces no torque) and (3.28) From these expressions, and using (3.22), one gets after a tedious but straightforward computation The determinant ∆ given by (3.21) is thus equal to which does not vanish for ζ large enough. 2

Optimal swimming
A possible notion of optimal swimming consists in minimizing the energy dissipated while trying to reach a given position p T at time T from an initial position p 0 at time 0. In this sense, the optimal stroke is the one that minimizes where v satisfies the self propulsion conditions expressed by Equations (3.2) and (3.3), subject to the constraint that p(T ) is equal to p T .
The special structure of our problem allows us to express the dissipated energy solely in terms of the shape variables ξ, by virtue of equation (3.4). We start by expressing the velocity of the m−th ball due to a unit shape change in the i-th component ofξ as This allows us to simplify the energy of the system to an expression depending on the rate of shape and position changesξ: where the metric G is defined by Since dissipation is invariant with respect to the current position p, it is possible to show that G ij (ξ, p) = G ij (ξ, 0) := G ij (ξ) for any position p.
Optimal swimming is then given by the solution (ξ, p) ∈ (S × P) × [0, T ] of min (ξ,p)∈(S×P) T 0ξ · G(ξ) ·ξ dt =: E(ξ) (4.5a) where the contraints C(ξ, p, t) = 0 are given explicitly by Some of the constraints in (4.6) can be relaxed. For example, it is possible to fix only the periodicity of the shape, but not its actual initial and final values (so that the last constraint in (4.6) disappears), or to fix only some of the components of p(T ), leaving the others free to be optimized.

Numerical Approximation
Our strategy for the solution to problem (4.5) is to discretize in time both ξ(t) and p(t) using cubic splines, and rewrite everything as a finite dimensional constrained optimization problem, defined only on the coefficients of the spline approximation of ξ and p.
A custom C++ code for boundary integral methods, based on the deal.II library (see [7] and [6]), solves the Dirichlet to Neumann map T −1 , which is then used to construct G(ξ) and V(ξ, p).
The finite element method in the spline space is then used to produce a finite dimensional approximation to (4.5), which is solved using the reduced space successive quadratic programming algorithm (rSQP) provided with the Moocho package of the Trilinos C++ library [12].
In this section we touch only briefly on the numerical implementation of the discrete Dirichlet to Neumann map for a collection of spheres, and we refer to [2] for more details on the rSQP technique used for the search of constrained minima.
In principle, standard surface meshing techniques could be employed to integrate on the surface of our swimmer, but they lack the speed and accuracy that specialized methods can offer when integrating on such special domains like collections of spheres.
The method we use for the integration on the spheres is presented in [10] for use with smooth functions and experimental data. It uses an oblique array of integration sampling points based on a chosen pair of successive Fibonacci numbers.
The pattern which is obtained on each sphere has a familiar appearance of intersecting spirals, avoiding the local anisotropy of a conventional latitudelongitude array. We modified the original method to take explicitly into account multiple spheres and the singularity of the Stokeslet (see [21] for a more detailed introduction on boundary integral methods for viscous flow).
The discrete versions u h and f h of the continuous velocity and force densities u and f are obtained by sampling the original functions at the quadrature points, translated and replicated on each sphere according to ξ and p, to obtain N V scalars that represent N V /3 vectors applied to the quadrature points sampled on ∂Ω.
Integration on ∂Ω of the function u h (x) is then effectively approximated by where the weights ω i and quadrature points q i are derived according to [10].
The construction of a discrete Dirichlet to Neumann map for our swimmers follows the usual collocation method (see, for example, [22]), where the matrix A is obtained by writing explicitly the various components of Equation (5.1), and we have used the notation u h and f h to indicate the N V dimensional vector that samples the functions u(x) and f (x) at the quadrature points q.
The construction of the discrete Dirichlet to Neumann map allows us to compute explicitly, given a pair of shape and position variables (ξ, p), the numerical approximations of G h (ξ) and V h (ξ, p). These can in turn be used to solve the constrained minimization problem as in [2], i.e., by using the implementation of the rSQP method included in the Moocho package of the Trilinos library [8].

Numerical Results
In this section we analyze some of the optimal strokes that can be obtained with our software for the case of the three-sphere swimmer in the plane [17] (3SP), depicted in Figure 4. While this is one of simplest possible extension of the three-sphere swimmer by Najafi and Golestanian [20], it adds already a high degree of complexity to the numerical approximation scheme.
The added complexity comes both from the presence of three additional time dependent variables (one shape and two positions), and from the fact that the vector field that relates positional changes to shape changes (ṗ = V(ξ, p) ·ξ) depends also (though in an explicit way) on the orientation θ = p 3 of the swimmer. In the software used to compute the optimal strokes for axisymmetric swimmers [2], there is only one position variable p 1 (t) = c(t), which can be obtained by post-processing the vector field V(ξ) (which, in the axisymmetric case, is independent on the position p) and ξ itself.
In the present case of a planar swimmer, this is no longer possible, and the optimization software has to go from optimizing over the two shape variables ξ 1 (t) and ξ 2 (t) to optimizing over the entire shape-position space, composed of the six time dependent variables ξ(t) = (ξ 1 (t), ξ 2 (t), ξ 3 (t)) and p(t) = (c x (t), c y (t), θ(t)).
In the experiments that follow, we use water at room temperature (20 o C) as the surrounding fluid, and we express lengths in millimeters (mm), time in seconds (s) and weight in milligrams (mg). Using these units, the viscosity of water is approximately one (1mP a s = 1mg mm −1 s −1 ), and the energy is expressed in pico-Joules (1pJ = 1mg mm 2 s −2 = 10 −12 J).
The radius of each ball of the swimmer is .05mm, and additional constraints are added to the optimization software to ensure that each component of the shape variable ξ stays in the interval [.1mm, .7mm], so that the balls never touch each other and that they don't separate too much.
We present some of the computed optimal strokes that illustrate the richness of the problem. In all the examples that follow, we ask the swimmer to go from p(0) = (0mm, 0mm, θ 0 ) to p(T ) = (.01mm, 0mm, θ T ), for various combinations of θ 0 , θ T and ξ(0) = ξ(T ).
Since the swimmer is 2π/3 periodic, it is sufficient to study the optimal strokes varying θ 0 between 0 and 2π/3 while keeping c x and c y constant to explore all possible optimal swimming strokes. All other combinations can be obtained by rotations and permutations between the components of the optimal solution (ξ, p).
Moreover, since rotating the swimmer of π is equivalent to asking it to swim in the opposite direction, reversibility of Stokes flow allows us to further reduce the study of the optimal strokes to the interval θ 0 ∈ [0, π/3].
In particular, if the pair ξ(t) and p(t) = (c(t), θ(t)) is an optimal solution of (4.5), then we have that all the following permutations, symmetries and rotations of (ξ, p) are also optimal swimming strokes, but with different initial and final positions: Given a family of optimal solutions with different starting angles θ 0 ∈ [0, π/3], then properties (6.1) can be used to derive the optimal solutions in all of the interval [0, 2π].
The optimal strokes we present have been formatted to show the trajectory of the center of the swimmer amplified of a factor 60 with respect to the dimensions of the swimmer itself, and the orientation θ of the swimmer amplified of a factor 10. The full configuration of the swimmer is shown at ten equally spaced instants of time between t = 0 and t = T = 1.
Our first optimal stroke ( Figure 5) is obtained by fixing the initial and final rotation θ 0 = θ T to be zero and the initial shape of the swimmer to be ξ = (.4, .4, .4)mm. The total energy consumption to swim by .01mm in the x direction with a single stroke is 0.511pJ.
The first important thing to notice is that, although the initial configuration is symmetric with respect to the translation axis, the stroke is not. This implies that the optimal solution is not unique: exchanging the off axis shape variables ξ 2 (t) and ξ 3 (t) yields to the same final displacement with the same total energy consumption.
A configuration which is symmetric with respect the y−axis can be obtained by setting θ 0 = π/6, and also in this case the solution is non symmetric and non unique.
Given the permutation properties (6.1), we deduce that two distinct global solutions with the same energy dissipation exist for θ 0 = N π/6, for any integer N . From the existence of these different global minima we can deduce that at least two branches of local minima exist outside those particular points, starting off as continuous perturbations of the two global solutions.
Numerical evidence show that at least three stable branches exist, as shown in Figure 6, where we plot the energy dissipation of all the (local) optimal strokes with respect to the starting initial angle p 3 (0) = θ 0 . It is possible to steer the software towards one or the other local branch by feeding it with initial guesses which are close to the desired branch.
preferred swimming directions. The directions perpendicular to one of the arms (θ 0 = π/6+N π/3) are the most expensive ones (although they differ in energy only by approximately 2%).
From Figure 5 it is evident that a lot of energy is wasted in forcing the final rotation to be equal to the initial one. If we relax this requirement and allow the software to optimize the stroke without forcing the rotation to be periodic in time, then the swimmer does not follow an "eight-shaped" path anymore, and the actual target displacement is reached with a smaller stroke (Figure 7), and by dissipating less than half the energy (.209pJ instead of .511pJ). The final rotation is approximately .058 radians (about 3.3 degrees).
The main defect associated with this strategy is that all successive strokes will be rotated with respect to the initial swimming direction. If the stroke were to be repeated, we would observe a rotational drift and an average circular motion along a trajectory of radius r = .35mm. To ensure that N successive strokes produce an average trajectory along a straight line, one cannot use N copies of the same stroke. Instead, N different optimization problems need to be solved, in each of which the initial orientation is the one reached at the end of the previous stroke.
A complementary strategy to reduce the energy of a single stroke without generating a rotated final configuration consists in optimizing also on the starting shape, rather than on the final angle. Figure 8 shows this strategy, where the final angle is fixed to be equal to the initial angle. From the figure it is clear that the more complex hydrodynamic interaction due to the closeness of the spheres can be effectively exploited by the optimization software to reduce the dissipated energy.
The optimal initial shape is found to be (.15, .15, .19)mm, a non symmetric shape. The total energy dissipation with this stroke strategy is about .221pJ. While this is already less than half of the energy consumed by our first optimal stroke ( Figure 5), it is still about 10% more expensive than the second attempt (Figure 7). The combination of the two strategies, in which both the initial shape and the final angle are left free for the software to optimize, is shown in Figure 9. In this case the total energy dissipation drops down to .100pJ. The initial shape is optimized to the value (.23, .35, .2)mm and the final angle is 0.057 radians, or about 3.28 degrees. In this case the optimal stroke hits the boundary of the region of allowed shapes.