Control and navigation problems for model bio-inspired microswimmers

Navigation problems for a model bio-inspired micro-swimmer, consisting of a cargo head and propelled by multiple rotating flagella or propellers and swimming at low Reynolds numbers, are formulated and solved. We consider both the direct problem, namely, predicting velocity and trajectories of the swimmer as a consequence of prescribed rotation rates of the propellers, and inverse problems, namely, find the rotation rates to best approximate desired translational and rotational velocities and, ultimately, target trajectories. The equations of motion of the swimmer express the balance of the forces and torques acting on the swimmer, and relate translational and rotational velocities of the cargo head to rotation rates of the propellers. The coefficients of these equations, representing hydrodynamic resistance coefficients, are evaluated numerically through a custom-built finite-element code to simulate the (Stokes) fluid flows generated by the movement of the swimmer and of its parts. Several designs of the propulsive rotors are considered: from helical flagella with different chirality to marine propellers, and their relative performance is assessed.

problem for an artificial microswimmer is to predict the swimming speed on the basis of the geometry of the swimmer and of the characteristics of the motor. This problem has been formulated by Purcell in his seminal papers [10,11], where the problem is also elegantly solved for a swimmer made of a rigid head and a rigid rotating helical tail (bio-inspired by the body architecture of the bacterium Escherichia coli), assuming that hydrodynamic interactions between head and tail are negligible. More recent extensions to include explicitly the effect of these hydrodynamic interactions, the presence of surrounding walls, or the effect of non-negligible inertial forces are discussed, for example, in [12][13][14].
Investigating biology with the functionalist viewpoint of the engineer often leads to new understanding of biological mechanisms and functions, while modeling the mechanics of life-like miniature machines can lead to the design of novel bio-inspired artificial constructs with superior performance. In this two-way interaction between biology and engineering, mathematical modeling plays an important role in that it allows to distill the essence of the success of some specific mechanism, phrase in the portable language of mathematical formalisms, and make it available to applications at different length scales, or based on different material systems, that can be quite different from the one that initially prompted the research. In this paper, we focus on robotic microswimmers whose body architecture is based on rotating tails attached to a rigid head and mimics the one of bacteria. We consider navigation and control problems inspired by those typically posed in the context of unmanned aerial vehicles (drones) see, e.g., [15,16]. We extend and adapt the techniques used in that field to the context of microswimming. In doing so, we follow Purcell [10] and consider, for simplicity, the case in which hydrodynamic interactions among the main body (hull) and the propulsive units (rotors) can be neglected. Through the analysis of some specific model systems, we show how to formulate and solve some navigation problems. The direct swimming problem consists in solving for the dynamics of the swimmer, namely, predicting velocity and trajectories of the swimmer as a consequence of prescribed rotation rates of the rotors. The inverse swimming problem (namely, find the rotation rates of the rotors to best approximate desired translational and rotational velocities and, ultimately, target trajectories) is a classical control problem. We anticipate and hope that this viewpoint will be useful also in future studies of microswimmers of biological interest, and in the design of new bio-inspired robotic constructs that can find interesting applications as biomedical devices.

Equations of motion
We consider a model micro-swimmer composed of a rotationally symmetric head (the cargo or resistive unit, which can be used as a hull to contain a payload, and that we model as a prolate ellipsoid of revolution along the longitudinal axis with unit vector 3 ) and four rotors (propulsive units shaped as either helical tails or marine propellers) with axes of rotation parallel to the longitudinal axis of the hull and rigidly attached to the hull with a fixed eccentricity e (the same eccentricity for all four propulsive units, for simplicity). Formally, we denote the ellipsoidal body of the microswimmer (the hull) as B h , the propulsive rotors as B p,i , i = 1, … , 4 , and we denote the swimmer body as B , which is the union of B h and of all the B p,i . The position vector with respect to the center gravity of B h of the i-th rotor is denoted by i and we will assume i = h 3 + e( i ) ⟂ , where ( i ) ⟂ is the unit vector obtained by normalizing the projection of i in the plane perpendicular to 3 , so that all the rotors have the same offset along the longitudinal axis h, and the same eccentricity e. We refer to Fig. 1 for a sketch of the geometry of the micro-swimmer.
For the rotors we will consider several different designs, including single helical tails, double helical tails of higher symmetry, and marine propellers, as discussed below. The fluid environment in which the robot moves is a Newtonian viscous fluid (air, water, glycerol, etc.) flowing at low Reynolds numbers, so that inertial forces can be neglected and its Here the volume of the fluid surrounding the body B extends to the whole space ℝ 3 (i.e., no confining walls are considered, for simplicity), ( ) ∈ ℝ 3 , p( ) ∈ ℝ are, respectively, the velocity and pressure of the fluid at a generic point ∈ ℝ 3 �B , while > 0 is the fluid dynamic viscosity and ∈ ℝ 3 is the velocity vector field of the solid swimmer. Equality between the velocity of the fluid and of the solid at points on the boundary of the swimmer, where the two are in contact, is the no-slip boundary condition. We consider a neutrally buoyant solid swimmer, and neglect inertial and gravity forces. The equations of motion reduce to the balance of generalized viscous forces (i.e., forces and moments, e.g., with respect to the center of mass of the swimmer's body) exerted by the fluid on the swimmer. The center of mass of the body of the swimmer is denoted by , and its orientation is described by a time-parametrized rotation that maps the orientation at time t = 0 , which we assume to coincide with the orientation of the lab frame, to the current orientation (the one of the body frame at time t) according to i (t) = ℚ(t) i (0) , with ℚ(0) = , the identity matrix. Translational and rotational velocities of the body of the swimmer are related to positional and orientational variables by the equations where we have denoted with ̌ the axial tensor associated with , i.e., the skew-symmetric tensor such that ̌ = × for every vector .
In view of the linearity of the Stokes system (1), the generalized viscous forces exerted by the fluid on the swimmer are linear functions of the translational and rotational velocities of the swimmer , and of the rotational velocities of the rotors i = i 3 , so that the equations of motion of the swimmer can be written as Here the 6 × 6 grand resistance matrix tot contains the hydrodynamic resistive coefficients for rigid motions of the swimmer, and the i-th column of the 6 × 4 propulsive matrix ℙ contains forces and torques generated by the i-th rotor turning at unit speed while all other rotors are still ( i = 1 , j = 0 for j ≠ i ) and the hull is not moving ( = = ). The procedure to construct these matrices is further described below. Following [10,11] we use additivity of resistance forces contributed by the individual components of the swimmer (this approximation arises from neglecting hydrodynamic interactions between the various parts of the composite swimmer, see the discussion in Sect. 7) and write where are the grand resistance matrices of the hull and the i-th rotor, respectively, each partitioned in 3 × 3 blocks containing translational ( ), coupling ( ℂ ) and rotational ( ) resistance coefficients, while where (ř i ) is the axial tensor associated with r i , the position vector of the i-th rotor with respect to . Matrices e,i contain the transport terms due to the eccentricity of the rotors, and are needed to account for the fact that the moments in the grand resistance matrices p,i are computed with respect to the center of mass of the i-th rotor, while the last three equations in (3) express the vanishing of the moment with respect to the center of mass of the swimmer body. Moreover, again due to the possible eccentricity of the rotors, a rotational motion around the center of the hull with angular velocity induces a translational velocity of the i-th rotor given by × i . Clearly, all these contributions vanish in the case i = 0 . Finally, recalling that i = i 3 , the i-th column of the 6 × 4 propulsion matrix ℙ is given by More in detail, the coefficients of the grand resistance matrices (5) are obtained from the following representation formulas for the hydrodynamic forces and moments on the hull and on the i-th rotor where h is the total force on the hull and h is the total moment with respect to the hull center of gravity of the forces acting on the hull, (i) represents the force acting on the i-th rotor, (i) is the moment computed with respect to the center of gravity of the i-th propeller p,i , and i (t) = ℚ(t) i (0) is the position vector of the i-th rotor with respect to the center of gravity of the hull. Computing the total force and the total moment with respect to of allthe forces in (8), and factoring out , , and i , we obtain the equations of motion (3), with coefficients given by (4)-(7).

An ideal swimmer with high symmetry
It is useful to consider an ideal swimmer (that is endowed with complete rotational symmetry, and for which the structure of the matrices containing the hydrodynamic coefficients is particularly simple, thanks to many cancellations) with the property that selected swimming problems can be solved analytically. These closed form solutions will be used as conceptual benchmarks for swimming problems formulated on more realistic swimmer geometries, that can be solved only numerically, as discussed later in Sect. 6.
The grand resistance matrix of the swimmer body (the hull) is assumed of the form Notice that the equation above gives the components of the grand resistance matrix with respect to the body frame, which we have chosen to be aligned with the principal axes of inertia of the swimmer body, whose longitudinal axis is parallel to 3 . These components are invariant with time. This is natural since viscous drag coefficients reflect the shape of the object, which is seen invariant in the body frame. We are assimilating here the hull to an ellipsoid of revolution, and we are neglecting the presence of rigid links anchoring the rotors to the swimmer body. In Sect. 6 we will model explicitly these links, and the structure of the matrix will be slightly perturbed. The grand resistance matrix of the i-th rotor is assumed of the form We will consider two cases. Either the rotors are all equal, or else the first two propellers have one chirality (e.g., they are both left-handed helical filaments), and the other two have opposite chirality (e.g., they are both right-handed helical filaments). We assume that the unit vectors 1 , 2 are parallel to the links between hull and propellers shown in Fig. 2. The first Fig. 2, case (a). All rotors have the same grand resistance matrix, given by The propulsive matrix for this case is then given by The second case is illustrated in Fig. 2, case (b). In this case we have in Eq. (10) (these coupling coefficients change sign because of the change of chirality) and the propulsive matrix is given by Again, we refer to Sect. 6 for examples of grand resistance and propulsion matrices arising for concrete examples of realistic rotors, shaped as either helical tails or as marine propellers, and computed by solving numerically outer Stokes problems of the type (1), with suitable assignments of the swimmer velocity field .

Swimming problems: solutions for the ideal swimmer
Substituting the hydrodynamic coefficients for the ideal swimmers given by Eq. (9) and either Eqs. (11)- (12) or Eqs. (13)-(14) into the equations of motion (3), we can formulate and solve analytically some physically relevant navigation problems that can be used as conceptual benchmarks for the discussion of more complex cases, for which only numerical simulation is feasible.

Direct swimming problem
The direct swimming problem consists in predicting the velocity of the swimmer in response to prescribed rotation rates of the rotors. In other words, given the components of the vector of rotor velocities as data, solve the equations of motion (3) in terms of the unknown translational and rotational velocities. Since the 6 × 6 matrix tot is invertible, there is a unique solution for any arbitrary assignment of the rotor velocities i , i = 1, … , 4 , namely, We are interested in particular to the following two concrete problems. For the ideal swimmer with identical rotors, case (a), consider equal rotor velocities Using Eqs. (9), (11), and (12), Eq. (15) becomes Its solution is given by U 1 = U 2 = Ω 1 = Ω 2 = 0 and or, more explicitly, We remark that, in this case, it is impossible to have a non-zero translational velocity accompanied by zero rotational velocity. Since it turns out that always C 33 ≪ K 33 , W 33 , zero rotational velocity con only occur if = 0 , and hence also the translational velocity must vanish. For the ideal swimmer with opposite rotors having opposite chirality, case (b), we are interested in the case and Eq. (15) becomes Its solution is given by U 1 = U 2 = Ω 1 = Ω 2 = Ω 3 = 0 and the only non-vanishing velocity component is We remark that, in this case, it is possible to have a non-zero translational velocity accompanied by zero rotational velocity. As discussed in more detail in Sect. 5, and as it is intuitively obvious, this leads to straight trajectories along the initial orientation of the longitudinal axis of the body swimmer 3 (0) , with no rotations along the body axis (purely translational motion along a straight line parallel to 3 (0)).

Inverse swimming problem
The inverse swimming problem consists in determining rotation rates i , i = 1, … , 4 of the rotors to best approximate a set of translational and rotational velocities and, ultimately, a target trajectory. Arranging the i in a vector = 1 , 2 , 3 , 4 T of unknowns, and the components of and in a vector = U 1 , U 2 , U 3 , Ω 1 , Ω 2 , Ω 3 T of data, we can rewrite the equations of motion (3) in the standard form = as follows The 6 × 4 matrix ℙ is clearly not invertible and there exist data for which Eq. (23) do not have exact solutions. One can then turn to approximate solutions, which are typically not unique. One possible notion of approximate solution of (23) is one that minimizes the error in the least squares sense The pseudo-inverse ℙ † of the coefficient matrix ℙ is useful to find the solution † of (24) with smallest norm ‖ ‖ = ( computed starting from the Singular Value Decomposition (SVD) of ℙ where , are two orthogonal matrices whose columns contain the left-singular and right-singular vectors of ℙ , respectively, and is a diagonal matrix whose entries are the singular values Σ ii of ℙ . The pseudo-inverse ℙ † of ℙ is then given by where the elements of the diagonal matrix † are given either by We then obtain Concretely, we are interested in the inverse problem of approximating a translational motion with velocity The interest of this problem is that it leads to a straight trajectory parallel to the initial orientation of the longitudinal axis of the swimmer, joining an initial and a target position along that axis with no energy losses that would accompany wobbling motions around the straight trajectory, and with no rotations around the body axis. Besides saving energy, the last feature would be useful for microswimmers equipped with a camera located inside the hull and pointing along the direction of the body axis. Similarly to what happens with drones hovering above a still target, having no rotations would imply that the still target would not appear rotating in the images collected by the camera, as is desirable. For the ideal swimmer with identical rotors, case (a), the approximate solution (27) corresponding to given by (28) (translational motion along longitudinal axis of the hull) is, in particular, Equation (29) above represents only an approximate solution for the inverse problem of pure translational swimming with velocity (28). As it is known from the discussion of the direct swimming problem for the ideal swimmer of type (a), equal rotation rates of the  (30) above represents also an exact solution for the inverse problem of pure translational swimming with velocity (28). As it is known from our previous discussion of the direct swimming problem for the ideal swimmer of type (b), if propellers with opposite chirality rotate with angular velocities of equal absolute value and opposite sign, then the motion of the swimmer is a pure translation along the longitudinal axis of the hull.

Trajectories
Once the translational and rotational velocities , of the swimmer are known, its trajectory and orientation can be recovered by integrating equations (2), see e.g. [17]. We consider here the case of constant rotor velocities. Since the components in the body frame { i } , i = 1, 2, 3 of the hydrodynamic matrices in the equations of motion (3) are time-invariant (either exactly or up to a good level of approximation), it is natural to write these equations in the body frame, and to solve for the components of , in this frame, which are denoted by U i and Ω i , i = 1, 2, 3 and are time-invariant. It follows from (2) that and using the skew-symmetric matrix we can rewrite the last equation as We assume that the body frame coincides with the laboratory frame at t = 0 , and the solution of Eq. (34) is easily expressed in terms of the exponential matrix (t) = exp( t) as we then integrate (31) to get the trajectory of the center of gravity. We set and, similarly, Moreover, we denote by = cos −1 (̂ 0 ⋅̂ 0 ) the (constant) angle between and , by ̂ 0 the unit vector perpendicular to both ̂ 0 and ̂ 0 such that ̂ 0 ×̂ 0 = sin( )̂ 0 and set ̂ 0 =̂ 0 ×̂ 0 , so that the triplet {̂ 0 ,̂ 0 ,̂ 0 } is an orthonormal basis. We obtain Equation (36) gives the trajectory in the lab frame of the center of gravity (t) starting from (0) at time t = 0 and it represents a circular helix with axis parallel to ̂ 0 . The straight trajectories of Sect. 4 are a special case corresponding to parallel to , leading to sin( ) = 0 . When and are orthogonal, cos = 0 and the helical trajectories (36) degenerate to circular ones (the swimmer swims in circles). We remark that helical trajectories are ubiquitous for biological microswimmers (in particular: unicellular ciliates or flagellates) that swim by beating periodically their cilia or flagella [18][19][20][21]. This is a direct consequence of the invariance (to translations and rotations) of the equations of motion of the swimmer, Eq. (3), written in the body frame. A more precise version of this statement can be found in [22,23], where it is called the Helix Theorem.

Navigation problems and trajectories for realistic geometries
In order to solve concretely the equations of motion (3), the hydrodynamic coefficients contained in the grand resistance matrices (of either the whole swimmer or of one of its parts; we contemplate both cases by referring to the grand resistance matrix of a generic body B ) and in the propulsion matrix ℙ need to be evaluated. The row index in these matrix identifies the type of generalized force (components of force along { i } , i = 1, 2, 3 for the first three rows, components of the moment with respect to center of gravity along { i } , i = 1, 2, 3 for the second three rows). The column index identifies the type of swimmer motion for which the corresponding generalized force is being evaluated: translation at unit speed along { i } , i = 1, 2, 3 for the first three columns of , rotation at unit speed around an axis through the center of gravity and parallel to { i } , i = 1, 2, 3 for the last three columns of , rotation at unit speed of the i − th rotor (while the rest of the swimmer is fixed) for the i − th column of ℙ . We follow here [14] and exploit the reciprocal theorem (see, e.g., [24]) to compute these hydrodynamic coefficients. More in detail, denoting by i the velocity field in the fluid region B ext ⧵B obtained by solving the outer Stokes problem (1) with prescribed velocity i on the solid B , the generic hydrodynamic coefficient ij in matrices and ℙ can be computed as where D is the symmetric part of the gradient of the velocity field . Here we are writing B ext rather than ℝ 3 since, in practice, we solve the outer Stokes problem in free space by replacing ℝ 3 with a large computational box B ext and setting ( ) = 0 on ∈ B ext . As the row index i and the column index j vary in (37), the Dirichlet data i in the outer Stokes problem (1) scan translations and rotations at unit speed, the j scan translations, rotations, and individual rotor rotation at unit speed, and ij scans all the coefficients of the matrices and ℙ . For example, the element ij of the 3 × 3 submatrix giving the components of viscous drag (i-th force component) opposing translations at unit speed in the direction of j is where i (respectively, j ) is the solution of the outer Stokes problem (1) corresponding to a prescribed velocity on the swimmer's boundary given by = i (respectively, = j ). In order to evaluate the hydrodynamic coefficients a CFD program was developed within an MSc thesis by the first author using Pyhton, based on the FEniCSx library, a widespread open-source computing platform that can be executed with good results on a personal computer (Dell XPS 9560 2.8GHz i7-7700HQ with 8 GB RAM in our case). The code has been validated against benchmark results from the literature, comparing both velocity fields with the case of a rotating helix [25], and drag coefficients [26,27] for a sphere and an ellipsoid in free space, showing that the computational box is large enough and mesh refinement is sufficient to guarantee relative errors of magnitude 10 −3 with the use of a personal computer. Use of an implementation on a High Performance Computing (HPC) platform allowing larger computational boxes and finer meshes showed smaller relative errors, not exceeding 10 −5 .
The HPC code was also written in Python, and based on the FEniCS [28] library, which uses an MPI implementation for parallelcomputations, and relied, for the solution of the linear systems resulting from the discretization, on the the PETSc library, which provides state-of-the art parallel iterative solvers. The code also allows for more complex configurations, such as the possibility of simulating the hydrodynamic flow around a swimmer in confined environment (e.g. a narrow channel), and of taking inertial effects into account, which becomes a necessity when the size of the swimmer is not microscopic. The Navier-Stokes equations are solved in place of the Stokes equations when needed, at moderate values of Reynolds number. This code was also validated against results from literature, from standard examples of Stokes [29] flow to more complex motions of falling discs and spheres in a fluid at moderate (but non-zero) Reynolds numbers [30,31]. We collect in the rest of this section the matrices containing the hydrodynamic coefficients for some concrete geometries, illustrated in the figures that follow, and use them to calculate the corresponding trajectories generated by specific choices of rotor velocities, as illustrated below. The body of the swimmer is the prolate spheroid shown in Fig. 3 with four links of equal length e protruding in directions perpendicualr to the axis of   . The corresponding grand resistance matrix is given in (39). Here, and in all the formulas that follow for the resistance matrices we have used the notation ex to denote ×10 x ; the value of For the rotors we will consider several alternative choices. For the case shaped as a single left-handed helix with axis length L, radius R, and pitch P, illustrated in Fig. 4, the grand resistance matrix p is given by Considering propellers shaped as double-helices, obtained from the previous ones by adding a second identical helix with the same axis but shifted by a phase-angle , see Fig. 5, we obtain a grand resistance matrix p given by  Experimentally tested rotors: double helical tail and marine propeller from [14], geometric data in Table 1. Surfaces used for mesh generation (left) and a detail of the trace of the computational mesh on the propeller blades (right) are shown The lack of symmetry in the propeller produces trajectories with more pronounced wobbling than in cases where the propellers have higher symmetry. This is true not only in the comparison with the ideal swimmers of Sect. 4, but also in going from the single-helix to the double-helix shape, as shown by the trajectories reported below.
The trajectories corresponding to the two choices of Finally, two additional rotor geometries were considered, one shaped as a double right-handed helical tail and another one shaped as a marine propeller, as shown in Fig. 9. These shapes included here for comparison, since they have been 3D-printed and their hydrodynamic coefficients have been computed on a HPC platform and validated experimentally, as reported in [14]. The grand resistance matrices for the two propellers shown in Fig. 9    The geometric data for the helical tail and the marine propeller of Fig. 9 are reported in Table 1. The resulting 3D trajectories are shown in Fig. 10 while their projection onto the plane XZ are shown in Fig. 11.
The swimming performance of the marine propeller is not good in terms of distance covered along the Z-axis, as expected. This is due to the fact that at low Reynolds numbers inertia is negligible, and thus the rotating blades are not able to generate significant lift. Hence, they produce small propulsive forces. This kind of propeller would fare better in a comparison at moderate Reynolds numbers, as is done in [14]. Such a study however is outside of the scope of the present work. For the single and double helix propellers, simulations are performed for each of the two configurations shown in Fig. 2, and the numerical values of the computed velocities are reported in Table 2 and Table 3. For the double helical tail and marine propellers, only case a) is considered since geometries of opposite chirality have not been studied yet with experiments. In case (b) the velocity of rotation Ω 3 is never exactly zero as in the ideal case of Sect. 3, but it is substantially reduced with respect to case (a), for a propeller of equivalent geometry. For propellers shaped as double-helices this rotational velocity is negligible, leading to good performance in applications requiring an on-board camera housed in the hull, and pointing along the direction of the longitudinal axis.

Discussion and perspectives
In conclusion, we have formulated and solved the problem of predicting velocities and trajectories of composite swimmers made of an ellipsoidal head and multiple eccentric rotors, as a consequence of prescribed rotation rates of the propellers (direct swimming problem). Conversely, we have discussed how to find the rotation rates of the propellers to best approximate desired translational and rotational velocities and, ultimately, target trajectories (inverse swimming problem). We have obtained analytical solutions for ideal swimmers with high symmetry, and numerical solutions for more realistic swimmers which do not exhibit full rotational symmetry. Future extensions of this work will examine the role of hydrodynamic interactions. To resolve them, the hydrodynamic coefficients in (3) must be evaluated directly for the whole composite swimmer, along the lines discussed in [13,14], and cannot be obtained by simply adding the individual contributions of the various parts of the composite swimmer, as in Eqs. (4)- (6). Another possible extension is to consider the presence of surrounding walls and the effect of inertial forces, similarly to what has been done in [14] for a robotic swimmer with a single non-eccentric rotor. The HPC code described in the previous section was developed with the aim of investigating these possible extensions of our work. A further direction worth exploring is to address optimal control problems for best matching of desired (target) trajectories, similarly to what is discussed in [32].