Nomenclature

$ a$

Sphere radius ( $ \mu m$ )

$ {\textrm{U}}_{i{\kern1pt}j}$

Distance between centers of sphere $ i\;$ and sphere $ j,$ $ i,j=1,\,2,\,3.$ ( $ {\mu}\textrm{m}$ )

$$\mu $$

Fluid viscosity ( $ \textrm{k}\textrm{g}.{\textrm{m}}^{-1}.{\textrm{s}}^{-1}$ )

$ \textbf{v}$

Fluid velocity ( $ m.\,{s}^{-1}$ )

p

Fluid pressure ( $ Pa$ )

$ {\textbf{f}}_{j}$ , $ {\textbf{v}}_{i}$

Force acting on sphere i ( $ N$ ), velocity of sphere i ( $ m.\,{s}^{-1}$ )

$ {\textbf{H}}_{i{\kern1pt}j}$

Hydrodynamic Oseen tensor ( $ m.\,{s}^{-1}.{N}^{-1}$ )

$ {\textbf{r}}_{i}$ , $ {\textbf{r}}_{i{\kern1pt}j}$

Position vector of sphere center i ( $ {\mu }\textrm{m}$ ), position vector of sphere center i relative to sphere center $ j$ ( $ {\mu }\textrm{m}$ )

$ \textbf{I}$

Identity matrix (-)

$ {x}_{i}$ , $ {y}_{i}$ , $ {z}_{i}$

$ x$ , $ y$ , and $ z$ positions of sphere $ i$ in three-dimensional Cartesian coordinate system ( $ {\mu }\textrm{m}$ )

$$\theta $$

Angle between arm connecting spheres 1 and 2 and $ x$ axis (radians)

$ {\left({x}_{1}\right)}_{d}$ , $ {\left({y}_{1}\right)}_{d}$

Desired $ x$ and $ y$ positions of sphere 1 ( $$\mu {\rm{m}}$$ )

$${\left( \theta \right)_d}$$

Desired $$\theta $$ (radians)

$ {k}_{p1}$ , $ {k}_{p2}$ , $ {k}_{p3}$

Proportional control coefficients ( $ {\textrm{s}}^{-1}$ )

$ \Delta \textrm{t}$

Simulation time step ( $ s$ )

$ J$

Optimization objective function

$ \textrm{L}$

Lagrange function

2. Governing Equations

The small scale of microswimmers results in a low-Reynolds flow regime. As a result, the Navier–Stokes equations reduce to the linear, time-independent Stokes equation. Therefore, fluid flow resulting from a swimming microrobot is governed by the Stokes and continuity equations [Reference Datt, Nasouri and Elfring23]. For an incompressible Newtonian fluid, the equations are as follows [Reference Happel, Brenner and Moreau28]:

(1) \begin{align} \boldsymbol{\nabla} \cdot {\bf{v}} = 0,\,\,\,\,\,\,\, - \boldsymbol{\nabla} {\rm{p}} + {\mu } {\nabla ^2}\,{\bf{v}} = 0, \end{align}

where $ \textbf{v}$ denotes the fluid velocity vector, p is the fluid pressure, and $ {\mu }$ is the fluid viscosity. In modeling the microrobot, dynamic interactions between the arms and the fluid are neglected, mainly because of their relative slenderness. Therefore, Eq. (1) need to be solved subject to no-slip conditions relative to the spheres and $ \textbf{v}=0$ at infinity. It is also assumed that forces acting on each sphere pass through its center.

The relation between each sphere’s velocity and the forces acting on each sphere can be presented as follows [Reference Ebrahimian, Yekehzare and Ejtehadi25]:

(2) \begin{align} {{\bf{v}}_i} = \sum _{j = 1}^3{{\bf{H}}_{i{\kern1pt}j}}{{\bf{f}}_j}, \end{align}

where $ {\textbf{v}}_{i}$ is the velocity of sphere i, $ {\textbf{f}}_{j}$ is the force acting on sphere j, and $ {\textbf{H}}_{i{\kern1pt}j}$ is the symmetric Oseen tensor, which depends on the geometry of the system. For a system of linked spherical swimmers, $ {\textbf{H}}_{i{\kern1pt}j}$ can be written as follows [Reference Alexander, Pooley and Yeomans29]:

(3) \begin{align}{{\bf{H}}_{i{\kern1pt}j}} = \left\{ \begin{array}{r@{\qquad}l} {{\;\;\;\;\;\;1 \over {6\pi \mu a}}\textbf{I},} & {{\rm{if}}\,\,i = j}\\[5pt] {{1 \over {8\pi \mu \left| {{{\bf{r}}_{ij}}} \right|}}(\textbf{I} + {{{{\bf{r}}_{i{\kern1pt}j}}{{\bf{r}}_{i{\kern1pt}j}}} \over {{{\left| {{{\bf{r}}_{i{\kern1pt}j}}} \right|}^2}}}),} & {{\rm{if}}\,\,i \ne j}\end{array} \right.\end{align}

where $ {\textbf{r}}_{i{\kern1pt}j}={\textbf{r}}_{i}-{\textbf{r}}_{j},$ and $ {\textbf{r}}_{i}$ is the position vector of sphere i. Here $ \textbf{I}$ is the identity matrix of size 3. Along with Eq. (3), the microswimmer’s motion is also subject to the conservation of linear and angular momentum [Reference Ebrahimian, Yekehzare and Ejtehadi25]:

(4) \begin{align} \sum\limits_{i = 1}^3 {{{\bf{f}}_i}} = 0,\,\,\,\,\,\,\,\,\,\,\,\sum\limits_{i = 1}^3 {{{\bf{r}}_i} \times {{\bf{f}}_i}} = 0. \end{align}

The set of equations in (2) and (4) forms the dynamic equations of the swimming predator.

2.1 Derivation of Oseen tensor

In order to use Eq. (2), we derive the explicit form of the hydrodynamic Oseen tensor presented in relation (3). The center of each sphere is presented in a three-dimensional Cartesian coordinate system, where $ {x}_{i},$ $ {y}_{i}$ , and $ {z}_{i}$ represent the coordinates of sphere i. The relative distances are shown in the following order:

(5) \begin{align} {z_{i{\kern1pt}j}} = {z_i} - {z_j},\,\,\,\,\,{y_{i{\kern1pt}j}} = {y_i} - {y_j},\,\,\,\,\,{x_{i{\kern1pt}j}} = {x_i} - {x_j}, \end{align}

(6) \begin{align} {{\bf{r}}_{i{\kern1pt}j}} = {{\bf{r}}_i} - {{\bf{r}}_j} = \left[\begin{matrix}{{x_{i{\kern1pt}j}}} \\ {{y_{i{\kern1pt}j}}} \\ {{z_{i{\kern1pt}j}}}\end{matrix} \right],\,\,\,\,\,\,\,\left| {{{\bf{r}}_{i{\kern1pt}j}}} \right| = \sqrt {{x_{i{\kern1pt}j}}^2 + {y_{i{\kern1pt}j}}^2 + {z_{i{\kern1pt}j}}^2}. \end{align}

Accordingly, for $ i\ne j,$ $ {\textbf{H}}_{i{\kern1pt}j}$ can be obtained in the following order:

(7) \begin{align} {{\bf{H}}_{i{\kern1pt}j}} = {1 \over {8\pi \mu \left| {{{\bf{r}}_{i{\kern1pt}j}}} \right|}}\left( {{\bf{I}} + {{{{\bf{r}}_{i{\kern1pt}j}}{{\bf{r}}_{i{\kern1pt}j}}} \over {{{\left| {{{\bf{r}}_{i{\kern1pt}j}}} \right|}^2}}}} \right) = {1 \over {8\pi \mu \left| {{{\bf{r}}_{i{\kern1pt}j}}} \right|}}\left( {\left[ \begin{matrix} 1\;\;\;\;\; & 0\;\;\;\;\; & 0 \\ 0\;\;\;\;\; & 1\;\;\;\;\; & 0 \\ 0\;\;\;\;\; & 0\;\;\;\;\; & 1 \\ \end{matrix} \right] + {{ \begin{bmatrix} {{x_{i{\kern1pt}j}}^2}\;\;\;\;\; & {{x_{i{\kern1pt}j}}{y_{i{\kern1pt}j}}}\;\;\;\;\; & {{x_{i{\kern1pt}j}}{z_{i{\kern1pt}j}}} \\ {{y_{i{\kern1pt}j}}{x_{i{\kern1pt}j}}}\;\;\;\;\; & {{y_{i{\kern1pt}j}}^2}\;\;\;\;\; & {{y_{i{\kern1pt}j}}{z_{i{\kern1pt}j}}} \\ {{z_{i{\kern1pt}j}}{x_{i{\kern1pt}j}}}\;\;\;\;\; & {{z_{i{\kern1pt}j}}{y_{i{\kern1pt}j}}}\;\;\;\;\; & {{z_{i{\kern1pt}j}}^2} \end{bmatrix}} \over {{x_{i{\kern1pt}j}}^2 + {y_{i{\kern1pt}j}}^2 + {z_{i{\kern1pt}j}}^2}}} \right) \end{align}

In the case where $ i=j:$

(8) \begin{align} {{\bf{H}}_{i{\kern1pt}j}} = {1 \over {6\pi \mu a}}{\bf{I}} = {1 \over {6\pi \mu a}}\left[ \begin{matrix} 1\;\;\;\;\; & 0\;\;\;\;\; & 0 \\ 0\;\;\;\;\; & 1\;\;\;\;\; & 0 \\ 0\;\;\;\;\; & 0\;\;\;\;\; & 1 \end{matrix} \right].\end{align}

In the case where a microswimmer’s motion is two-dimensional, Eqs. (7) and (8) can be simplified for two-dimensional motion:

(9) \begin{align} {{\bf{H}}_{ii}} = {1 \over {6\pi \mu a}}\left[ \begin{matrix} 1\;\;\; & 0 \\ 0\;\;\; & 1 \end{matrix} \right],\,\,\,\,\,\,\,\,{{\bf{H}}_{i{\kern1pt}j}} = {{\left[ \begin{matrix} {2{x_{i{\kern1pt}j}}^2 + {y_{i{\kern1pt}j}}^2}\;\;\;\;\; & {{x_{i{\kern1pt}j}}{y_{i{\kern1pt}j}}} \\ {{y_{i{\kern1pt}j}}{x_{i{\kern1pt}j}}}\;\;\;\;\; & {{x_{i{\kern1pt}j}}^2 + 2{y_{i{\kern1pt}j}}^2} \end{matrix} \right]} \over {8\pi \mu {{\left( {{x_{i{\kern1pt}j}}^2 + {y_{i{\kern1pt}j}}^2} \right)}^{{3 \over 2}}}}} \end{align}

Eq. (9) yields the Oseen tensor for a two-dimensional, three-sphere microsystem.

2.2 Derivation of the Governing Equations of Motion

The explicit form of the low-Reynolds predator’s equations of motion can be obtained by expanding Eqs. (2) and (4) in terms of the hydrodynamic Oseen tensor in (9). First, Eq. (2) needs to be expanded. To this end, distances between spheres i and j are represented by $ {\textrm{U}}_{i{\kern1pt}j},$ and forces exerted on sphere i are shown by $\left( \begin{matrix} {{{\rm{f}}_{xi}}} \\ {{{\rm{f}}_{yi}}} \end{matrix} \right).$ The position of each sphere center is represented according to Eqs. (5) and (6).

Merging Eqs. (2) and (9) results in the following matrix equation:

(10) \begin{align} \dot{x} = {\bf{Gf}}.\end{align}

Here, ${\bf\dot{x}} = {[\begin{matrix} {{{\dot x}_1}}\;\;\;\;\; & {{{\dot y}_1}}\;\;\;\;\; & {{{\dot x}_2}}\;\;\;\;\; & {{{\dot y}_2}}\;\;\;\;\; & {{{\dot x}_3}}\;\;\;\;\; & {{{\dot y}_3}} \end{matrix} ]^{\rm{T}}}$ and ${\bf{f}} = {[\begin{matrix} {{{\rm{f}}_{x1}}}\;\;\;\;\; & {{{\rm{f}}_{y1}}}\;\;\;\;\; & {{{\rm{f}}_{x2}}}\;\;\;\;\; & {{{\rm{f}}_{y2}}}\;\;\;\;\; & {{{\rm{f}}_{x3}}}\;\;\;\;\; & {{{\rm{f}}_{y3}}} \end{matrix} ]^{\rm{T}}}.$ $ \textbf{G}$ is a $ 6\times 6$ matrix which depends on the relative position of each sphere’s center:

(11) \begin{align} {\bf{G}} = \left[ \begin{matrix} \begin{matrix} {{1 \over {6\pi \mu a}}}\;\;\;\;\; & 0\;\;\;\;\; & {{{2{x_{12}}^2 + {y_{12}}^2} \over {8\pi \mu {{\rm{U}}_{12}}^3}}}\;\;\;\;\; & {{{{x_{12}}{y_{12}}} \over {8\pi \mu {{\rm{U}}_{12}}^3}}}\;\;\;\;\; & {{{2{x_{13}}^2 + {y_{13}}^2} \over {8\pi \mu {{\rm{U}}_{31}}^3}}}\;\;\;\;\; & {{{{x_{13}}{y_{13}}} \over {8\pi \mu {{\rm{U}}_{31}}^3}}} \\[5pt] 0\;\;\;\;\; & {{1 \over {6\pi \mu a}}}\;\;\;\;\; & {{{{y_{12}}{x_{12}}} \over {8\pi \mu {{\rm{U}}_{12}}^3}}}\;\;\;\;\; & {{{{x_{12}}^2 + 2{y_{12}}^2} \over {8\pi \mu {{\rm{U}}_{12}}^3}}}\;\;\;\;\; & {{{{y_{13}}{x_{13}}} \over {8\pi \mu {{\rm{U}}_{31}}^3}}}\;\;\;\;\; & {{{{x_{13}}^2 + 2{y_{13}}^2} \over {8\pi \mu {{\rm{U}}_{31}}^3}}} \\[5pt] {{{2{x_{12}}^2 + {y_{12}}^2} \over {8\pi \mu {{\rm{U}}_{12}}^3}}}\;\;\;\;\; & {{{{x_{12}}{y_{12}}} \over {8\pi \mu {{\rm{U}}_{12}}^3}}}\;\;\;\;\; & {{1 \over {6\pi \mu a}}}\;\;\;\;\; & 0\;\;\;\;\; & {{{2{x_{23}}^2 + {y_{23}}^2} \over {8\pi \mu {{\rm{U}}_{23}}^3}}}\;\;\;\;\; & {{{{x_{23}}{y_{23}}} \over {8\pi \mu {{\rm{U}}_{23}}^3}}} \\[5pt] {{{{y_{12}}{x_{12}}} \over {8\pi \mu {{\rm{U}}_{12}}^3}}}\;\;\;\;\; & {{{{x_{12}}^2 + 2{y_{12}}^2} \over {8\pi \mu {{\rm{U}}_{12}}^3}}}\;\;\;\;\; & 0\;\;\;\;\; & {{1 \over {6\pi \mu a}}}\;\;\;\;\; & {{{{y_{23}}{x_{23}}} \over {8\pi \mu {{\rm{U}}_{23}}^3}}}\;\;\;\;\; & {{{{x_{23}}^2 + 2{y_{23}}^2} \over {8\pi \mu {{\rm{U}}_{23}}^3}}} \\[5pt] {{{2{x_{13}}^2 + {y_{13}}^2} \over {8\pi \mu {{\rm{U}}_{31}}^3}}}\;\;\;\;\; & {{{{x_{13}}{y_{13}}} \over {8\pi \mu {{\rm{U}}_{31}}^3}}}\;\;\;\;\; & {{{2{x_{23}}^2 + {y_{23}}^2} \over {8\pi \mu {{\rm{U}}_{23}}^3}}}\;\;\;\;\; & {{{{x_{23}}{y_{23}}} \over {8\pi \mu {{\rm{U}}_{23}}^3}}}\;\;\;\;\; & {{1 \over {6\pi \mu a}}}\;\;\;\;\; & 0 \\[5pt] {{{{y_{13}}{x_{13}}} \over {8\pi \mu {{\rm{U}}_{31}}^3}}}\;\;\;\;\; & {{{{x_{13}}^2 + 2{y_{13}}^2} \over {8\pi \mu {{\rm{U}}_{31}}^3}}}\;\;\;\;\; & {{{{y_{23}}{x_{23}}} \over {8\pi \mu {{\rm{U}}_{23}}^3}}}\;\;\;\;\; & {{{{x_{23}}^2 + 2{y_{23}}^2} \over {8\pi \mu {{\rm{U}}_{23}}^3}}}\;\;\;\;\; & 0\;\;\;\;\; & {{1 \over {6\pi \mu a}}} \end{matrix} \hfill \end{matrix} \right].\end{align}

Eq. (10) can alternatively be written as

(12) \begin{align} {\bf{f}} = {{\bf{G}}^{\bf - 1}}\dot{\textbf{x}}.\end{align}

The motion of this microswimmer is controlled by changing the length of each arm using prismatic actuators. Therefore, we are interested in writing the force vector in terms of a combination of velocity and time derivatives of the length of each arm. This can be achieved by deriving the relation between the following two sets of variables:

(13) \begin{align} \dot{\textbf{x}} = \textbf{D}\dot{\textbf{z}}, \end{align}

where $\dot{\textbf{z}} = {\left[ \begin{matrix} {{{\dot x}_1}}\;\;\;\;\; & {{{\dot y}_1}}\;\;\;\;\; & {{{\dot{\theta } }_{12}}}\;\;\;\;\; & {{{{{\dot{\textrm{U}}}}}_{12}}}\;\;\;\;\; & {{{{{\dot{\textrm{U}}}}}_{23}}}\;\;\;\;\; & {{{{{\dot{\textrm{U}}}}}_{31}}} \end{matrix} \right]^{\rm{T}}}.$ Matrix $ \textbf{D}$ in relation (13) can be determined by obtaining time derivatives of geometric relations presented in Fig. 2. Using relationship between spheres coordinates and after some mathematical manipulation, we can obtain:

(14) \begin{align} \dot{\textbf{x}} = \left[ \begin{matrix} {{{\dot x}_1}} \\[5pt] {{{\dot y}_1}} \\[5pt] {{{\dot x}_2}} \\[5pt] {{{\dot y}_2}} \\[5pt] {{{\dot x}_3}} \\[5pt] {{{\dot y}_3}} \end{matrix} \right] = \left[ \begin{matrix} 1\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & 0 \\[5pt] 0\;\;\;\;\; & 1\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & 0 \\[5pt] {{{\rm{E}}_{11}}}\;\;\;\;\; & {{{\rm{E}}_{12}}}\;\;\;\;\; & {{{\rm{E}}_{13}}}\;\;\;\;\; & {{{\rm{E}}_{14}}}\;\;\;\;\; & {{{\rm{E}}_{15}}}\;\;\;\;\; & {{{\rm{E}}_{16}}} \\[5pt] {{{\rm{E}}_{21}}}\;\;\;\;\; & {{{\rm{E}}_{22}}}\;\;\;\;\; & {{{\rm{E}}_{23}}}\;\;\;\;\; & {{{\rm{E}}_{24}}}\;\;\;\;\; & {{{\rm{E}}_{25}}}\;\;\;\;\; & {{{\rm{E}}_{26}}} \\[5pt] 1\;\;\;\;\; & 0\;\;\;\;\; & { - {y_{31}}}\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & {{{{x_{31}}} \over {{{\rm{U}}_{31}}}}} \\[5pt] 0\;\;\;\;\; & 1\;\;\;\;\; & {{x_{31}}}\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & {{{{y_{31}}} \over {{{\rm{U}}_{31}}}}} \end{matrix} \right]\left[ \begin{matrix} \begin{matrix} {{{\dot x}_1}} \\[5pt] {{{\dot y}_1}} \\[5pt] {\dot \theta } \end{matrix} \\[5pt] \begin{matrix} {{{{{\dot{\textrm{U}}}}}_{12}}} \\[5pt] {{{{{\dot{\textrm{U}}}}}_{23}}} \\[5pt] {{{{{\dot{\textrm{U}}}}}_{31}}} \end{matrix} \end{matrix} \right] = {\bf D}\dot{\textbf{z}}, \end{align}

Figure 2. Position vector of each sphere center in x-y coordinate and angle $$\theta .$$

where $ \textbf{E}={\textbf{A}}^{-1}\textbf{B}\textbf{C}\textbf{},$ $$\matrix{ {{\bf{A}} = \left[ {\matrix{ {{x_{12}}\;\;\;\;\;{{\rm{y}}_{12}}} \cr {{x_{23}}\;\;\;\;\;{{\rm{y}}_{23}}} \cr } } \right],{\bf{B}} = \left[ {\matrix{ {{x_{12}}} & {{{\rm{y}}_{12}}} & 0 & 0 & { - {{\rm{U}}_{12}}} & 0 \cr 0 & 0 & {{x_{23}}} & {{{\rm{y}}_{23}}} & 0 & {{{\rm{U}}_{23}}} \cr } } \right],} \hfill \cr } $$ and ${\bf{C}} = \left[ \begin{matrix} 1\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; \\ 0\;\;\;\;\; & 1\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; \\ 1\;\;\;\;\; & 0\;\;\;\;\; & { - {y_{31}}}\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & {{{{x_{31}}} \over {{{\rm{U}}_{31}}}}} \\ 0\;\;\;\;\; & 1\;\;\;\;\; & {{x_{31}}}\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & {{{{y_{31}}} \over {{{\rm{U}}_{31}}}}} \\ 0\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & 1\;\;\;\;\; & 0\;\;\;\;\; & 0 \\ 0\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & 0 \;\;\;\;\; & 1\;\;\;\;\; & 0 \end{matrix} \right]$ . Here $ {\textrm{E}}_{i{\kern1pt}j}$ are the elements of matrix $ \textbf{E}.$ Hence, Eq. (12) can be written in the following form:

(15) \begin{align} {\bf{f}} = {{\bf{G}}^{ - 1}}\dot{\textbf{x}} = {\bf M}{\dot{\textbf{z}}},\end{align}

where $ \textbf{M}={\textbf{G}}^{-1}\textbf{D}.$ Now using Eqs. (4) and (15) we can obtain:

(16) \begin{align} {[\begin{matrix} {{{\bf{M}}_1}}\;\;\; & {{{\bf{M}}_2}}\;\;\; & {{{\bf{M}}_3}} \end{matrix} ]^{\rm{T}}}{[\begin{matrix} {{{\dot x}_1}}\;\;\; & {{{\dot y}_1}}\;\;\; & {\dot {\theta } }\;\;\; & {{{{{\dot{\textrm{U}}}}}_{12}}}\;\;\; & {{{{{\dot{\textrm{U}}}}}_{23}}}\;\;\; & {{{{{\dot{\textrm{U}}}}}_{31}}} \end{matrix} ]^{\rm{T}}} = 0,\end{align}

where $${{\bf{M}}_1} = \left[ {\matrix{ {1\;\;\;0\;\;\;1\;\;\;0\;\;\;1\;\;\;0} \cr } } \right]{\bf{M}},$$ $${{\bf{M}}_2} = \left[ {\matrix{ {0\;\;\;1\;\;\;0\;\;\;1\;\;\;0\;\;\;1} \cr } } \right]{\bf{M}}$$ , and $${{\bf{M}}_3} = \left[ { - {y_1}\;\;\;{x_1}\;\;\; - {y_2}\;\;\;{x_2}\;\;\; - {y_3}} \right.$$ $$\left. {{x_3}} \right]{\bf{M}}.$$ Relation (16), which consists of three equations and six variables, describes the dynamic properties of the microrobot. It can be used alongside Eq. (13) to determine the position of each sphere’s center.

2.3 Dynamic simulation and validation

If at each point in time, the length of each arm is known, Eq. (16) yields a system of three equations and three variables, which can be solved in each time step. Here, we have simulated an ordered motion consisting of six steps and two cycles. In the first cycle, each arm shrinks from $ \textrm{L}$ to $ {\textrm{L}}^{*}=\left(1-\epsilon\right)\textrm{L},$ with a constant rate. In the next cycle, each arm, in turn, returns to its initial length with the same rate. In the end, the swimmer has returned to its initial configuration. The order of steps is presented in Fig. 3.

Figure 3. The order by which each arm shrinks and then returns to its original configuration [Reference Ebrahimian, Yekehzare and Ejtehadi25].

After these six steps, the predator experiences a net center of mass displacement along with a net rotation. The displacement and rotation after a full cycle have been presented in Fig. 4, as a function of $ \epsilon$ .

Figure 4. Displacement and rotation for $ \textrm{L}=10a$ as a function of $ \epsilon$ . The marks on the plot represent data obtained by Ebrahimian et al. in [Reference Ebrahimian, Yekehzare and Ejtehadi25].

To validate our dynamic model, results from Ebrahimian et al. have also been illustrated in Fig. 4, which demonstrates that the results that have been obtained and presented in this paper are completely coincident with results from [Reference Ebrahimian, Yekehzare and Ejtehadi25]. In the proceeding sections, we will present control strategies for controlling the motion of the low-Reynolds predator.

3. Predator Position and Heading Control

In this section, our aim is to determine the control inputs in a way that the microswimmer’s position and orientation track the desired trajectories. Eq. (16) consists of three equations and six variables. The lengths of each arm are the variables used to control the predator’s motion. Consequently, one control strategy can involve controlling three variables such as $ {x}_{1},$ $ {y}_{1}$ , and $ \theta ,$ where these variables, respectively, are the first sphere x and y positions and the first arm orientation. In this case, the first sphere is essentially the predator and must be able to successfully hunt down the prey. These three variables are controlled to follow desired paths $ {x}_{1d},$ $ {y}_{1d}$ , and $ {\theta }_{d}$ using proportional control with feed forward terms. $ {x}_{1d}$ and $ {y}_{1d}$ are the prey’s x and y positions:

(17) \begin{align} {\dot x_1} = {\dot x_{1d}} - {k_{p1}}\left( {{x_1} - {x_{1d}}} \right) = {c_1}, \end{align}

(18) \begin{align} {\dot y_1} = {\dot y_{1d}} - {k_{p2}}\left( {{y_1} - {y_{1d}}} \right) = {c_2}, \end{align}

(19) \begin{align} \dot \theta = {\dot \theta _d} - {k_{p3}}\left( {\theta - {\theta _{\rm{d}}}} \right) = {c_3}. \end{align}

As a result, $ {\dot{x}}_{1}$ , $ {\dot{y}}_{1}$ , and $ \dot{\theta }$ are determined in each time step, and Eq. (16) can be merged with Eqs. (17) to (19) to produce the following system of equations:

(20a) $$\left[ {\matrix{ 1 & 0 & 0 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 & 0 & 0 \cr 0 & 0 & 1 & 0 & 0 & 0 \cr {{{\bf{M}}_1}\left( {{\rm{1}},{\rm{1}}} \right)} & {{{\bf{M}}_1}\left( {{\rm{1}},{\rm{2}}} \right)} & {{{\bf{M}}_1}\left( {{\rm{1}},{\rm{3}}} \right)} & {{{\bf{M}}_1}\left( {{\rm{1}},{\rm{4}}} \right)} & {{{\bf{M}}_1}\left( {{\rm{1}},{\rm{5}}} \right)} & {{{\bf{M}}_1}\left( {{\rm{1}},{\rm{6}}} \right)} \cr {{{\bf{M}}_2}\left( {{\rm{1}},{\rm{1}}} \right)} & {{{\bf{M}}_2}\left( {{\rm{1}},{\rm{2}}} \right)} & {{{\bf{M}}_2}\left( {{\rm{1}},{\rm{3}}} \right)} & {{{\bf{M}}_2}\left( {{\rm{1}},{\rm{4}}} \right)} & {{{\bf{M}}_2}\left( {{\rm{1}},{\rm{5}}} \right)} & {{{\bf{M}}_2}\left( {{\rm{1}},{\rm{6}}} \right)} \cr {{{\bf{M}}_3}\left( {{\rm{1}},{\rm{1}}} \right)} & {{{\bf{M}}_3}\left( {{\rm{1}},{\rm{2}}} \right)} & {{{\bf{M}}_3}\left( {{\rm{1}},{\rm{3}}} \right)} & {{{\bf{M}}_3}\left( {{\rm{1}},{\rm{4}}} \right)} & {{{\bf{M}}_3}\left( {{\rm{1}},{\rm{5}}} \right)} & {{{\bf{M}}_3}\left( {{\rm{1}},{\rm{6}}} \right)} \cr } } \right]\left[ {\matrix{ {{{\dot x}_1}} \hfill \cr {{{\dot y}_1}} \hfill \cr {\dot \theta } \hfill \cr {{{\dot U}_{12}}} \hfill \cr {{{\dot U}_{23}}} \hfill \cr {{{\dot U}_{31}}} \hfill \cr } } \right] = \left[ {\matrix{ {{c_1}} \cr {{c_2}} \cr {{c_3}} \cr 0 \cr 0 \cr 0 \cr } } \right]$$

.

Here, we have a system of six equations and six unknowns, which can be solved in each time step. Solving for $ {\dot{\textrm{U}}}_{12},$ $ {\dot{\textrm{U}}}_{23}$ , and $ {\dot{\textrm{U}}}_{31}$ yields:

(20b) \begin{align} \bf{W}\left[\begin{array}{c}{\dot{\textrm{U}}}_{12}\textrm{}\\ {\dot{\textrm{U}}}_{23}\\ {\dot{\textrm{U}}}_{31}\end{array}\right]=\textbf{Q}\left[\begin{array}{c}{\dot{x}}_{1d}-{k}_{p1}\left({x}_{1}-{x}_{1d}\right)\\ {\dot{y}}_{1d}-{k}_{p2}\left({y}_{1}-{y}_{1d}\right)\\ {\dot{\theta }}_{d}-{k}_{p3}\left(\theta -{\theta }_{\textrm{d}}\right)\end{array}\right]\ \end{align}

( ) $${{\bf{W}} = \left[ {\matrix{ {{{\bf{M}}_{\bf{1}}}\left( {{\rm{1}},{\rm{4}}} \right)\;\;\;\;\;{{\bf{M}}_{\bf{1}}}\left( {{\rm{1}},{\rm{5}}} \right)\;\;\;\;\;{{\bf{M}}_{\bf{1}}}\left( {{\rm{1}},{\rm{6}}} \right)} \cr {{{\bf{M}}_{\bf{2}}}\left( {{\rm{1}},{\rm{4}}} \right)\;\;\;\;\;{{\bf{M}}_{\bf{2}}}\left( {{\rm{1}},{\rm{5}}} \right)\;\;\;\;\;{{\bf{M}}_{\bf{2}}}\left( {{\rm{1}},{\rm{6}}} \right)} \cr {{{\bf{M}}_{\bf{3}}}\left( {{\rm{1}},{\rm{4}}} \right)\;\;\;\;\;{{\bf{M}}_{\bf{3}}}\left( {{\rm{1}},{\rm{5}}} \right)\;\;\;\;\;{{\bf{M}}_{\bf{3}}}\left( {{\rm{1}},{\rm{6}}} \right)} \cr } } \right],\quad {\bf{Q}} = - \left[ {\matrix{ {{{\bf{M}}_{\bf{1}}}\left( {{\rm{1}},{\rm{1}}} \right)\;\;\;\;\;{{\bf{M}}_{\bf{1}}}\left( {{\rm{1}},{\rm{2}}} \right)\;\;\;\;\;{{\bf{M}}_{\bf{1}}}\left( {{\rm{1}},{\rm{3}}} \right)} \cr {{{\bf{M}}_{\bf{2}}}\left( {{\rm{1}},{\rm{1}}} \right)\;\;\;\;\;{{\bf{M}}_{\bf{2}}}\left( {{\rm{1}},{\rm{2}}} \right)\;\;\;\;\;{{\bf{M}}_{\bf{2}}}\left( {{\rm{1}},{\rm{3}}} \right)} \cr {{{\bf{M}}_{\bf{3}}}\left( {{\rm{1}},{\rm{1}}} \right)\;\;\;\;\;{{\bf{M}}_{\bf{3}}}\left( {{\rm{1}},{\rm{2}}} \right)\;\;\;\;\;{{\bf{M}}_{\bf{3}}}\left( {{\rm{1}},{\rm{3}}} \right)} \cr } } \right]}$$

The schematic of the employed control method and the details of controller block are shown in Fig. 5.

Figure 5. Block diagram of the control system (a), controller block details (b).

After variables $ {\dot{x}}_{1},$ $ {\dot{y}}_{1}$ , $ \dot{\theta }$ , $ {\dot{\textrm{U}}}_{12}$ , $ {\dot{\textrm{U}}}_{23}$ , and $ {\dot{\textrm{U}}}_{31}$ are determined from, the velocity of each sphere (variables $ {\dot{x}}_{1},$ $ {\dot{y}}_{1}$ , $ {\dot{x}}_{2}$ , $ {\dot{y}}_{2}$ , $ {\dot{x}}_{3}$ , and $ {\dot{y}}_{3})$ can be determined using Eq. (13). Integrating the velocity of each sphere in each time step, results in the position vector of the next time step:

(21) $${{{\left[ {\matrix{ {\matrix{ {{x_1}\;\;\;{y_1}\;\;\;{x_2}\;\;\;} \cr } \matrix{ {{y_2}\;\;\;{x_3}\;\;\;{y_3}} \cr } } \cr } } \right]}_{k + 1}} = {{\left[ {\matrix{ {\matrix{ {{x_1}\;\;\;{y_1}\;\;\;{x_2}\;\;\;} \cr } \matrix{ {{y_2}\;\;\;{x_3}\;\;\;{y_3}} \cr } } \cr } } \right]}_k} + \Delta {\rm{t}}{{\left[ {\matrix{ {\matrix{ {{{\dot x}_1}\;\;\;{{\dot y}_1}\;\;\;{{\dot x}_2}\;\;\;} \cr } \matrix{ {{{\dot y}_2}\;\;\;{{\dot x}_3}\;\;\;{{\dot y}_3}} \cr } } \cr } } \right]}_k}}$$

For simulation, a time step of $ \Delta \textrm{t}=0.01\textrm{s}$ has been considered, which leads to convergent results over a simulation time span of 10 s. The following desired paths have been considered:

(22) \begin{align} {x_{1d}} = \left( { - 2 - 0.6t} \right)\left( {\mu m} \right), \end{align}

(23) \begin{align} {y_{1d}} = \left( { - 3.75 - 5.85{e^{ - o.1t}}{\rm{cos}}\,0.7t} \right)\left( {\mu m} \right), \end{align}

(24) \begin{align} {\theta _{\rm{d}}} = {\pi \over 6}\left( {rad} \right), \end{align}

Control coefficients $ {k}_{p1}$ and $ {k}_{p2}$ are set equal to 1, $ {k}_{p3}$ is set at 0.5, initial values for $ {x}_{1},$ $ {y}_{1}$ , and $ \theta $ are assumed to be zero, and the initial arm lengths are $ 10\;{\mu }\textrm{m}.$ The results are presented in the following figures.

Fig. 6 shows that the predator has been able to successfully hunt down its prey according to Eqs. (22), (23), and (24), and the position and angle errors have reduced to zero. The x and y positions have settled after almost 4 s, whereas the angle has taken almost 8 s to settle. This is due to the fact that the control coefficient for angle $ {\theta }$ has been set at half that of $ {x}_{1}$ and $ {y}_{1}.$ According to Eqs. (22), (23), and (24), the desired prey path for the predator lies on a curved line. This path, along with the predator’s path, has been presented in Fig. 7, which shows their convergence.

The necessary arm length inputs have also been presented in Fig. 8. These results demonstrate that the arm length connecting spheres 1 and 3 has increased to almost $ 30\;{\mu }\textrm{m}.$ While these results satisfy our control goals, it is desirable to decrease the amount and change of arm lengths in order to lower the swimmer’s energy consumption. This can be achieved using optimal control strategies.

4. Optimal Control

An alternative control strategy is optimal control of the microswimmer position. Accordingly, variables $ {x}_{1}$ and $ {y}_{1},$ which are the predator’s x and y positions, are controlled to follow reference paths $ {\left({x}_{1}\right)}_{d}$ and $ {\left({y}_{1}\right)}_{d},$ which are the prey’s x and y positions, according to Eqs. (17) and (18), while minimizing an objective function $ J$ in each time step. Therefore, in each time step, variables $ {\theta }_{k+1},$ $ {\textrm{U}}_{{12}_{k+1}}$ , $ {\textrm{U}}_{{23}_{k+1}}$ , and $ {\textrm{U}}_{{31}_{k+1}}$ must be determined such that they minimize function $ J$ while satisfying dynamical equations in relation (16). As a result, the Lagrange function to find minima of the objective function $ J,$ subject to the equality constraints (16), is the following:

(25) \begin{align} {\rm{L}} = {\rm{J}} + \left( {{\gamma _1}{{\bf{M}}_1} + {\gamma _2}{{\bf{M}}_2} + {\gamma _3}{{\bf{M}}_3}} \right)\dot{\textbf{z}}. \end{align}

Figure 6. (a) Position of $ {x}_{1}$ and error as a function of time for control strategy involving control of position and angle. (b) Position of $ {y}_{1}$ and error as a function of time for control strategy involving control of position and angle. (c) Angle $ \theta $ and error as a function of time for control strategy involving control of position and angle.

Here parameters $ {\gamma }_{i}$ are the Lagrange multipliers. The equations for optimal control in each time step are the following:

(26) \begin{align} \left[\begin{array}{c}{\textbf{M}}_{1}\\ {\textbf{M}}_{2}\\ {\textbf{M}}_{3}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{\dot{x}}_{1}\\ {\dot{y}}_{1}\\ \dot{\theta }\end{array}\\ \begin{array}{c}{\dot{\textrm{U}}}_{12}\textrm{}\\ {\dot{\textrm{U}}}_{23}\\ {\dot{\textrm{U}}}_{31}\end{array}\end{array}\right]=\frac{\left[\begin{array}{c}{\textbf{M}}_{1}\\ {\textbf{M}}_{2}\\ {\textbf{M}}_{3}\end{array}\right]\left({{\left[\begin{array}{c}{x}_{1}\\ {y}_{1}\\ {\theta }\\ {\textrm{U}}_{12}\\ {\textrm{U}}_{23}\\ {\textrm{U}}_{31}\end{array}\right]}_{k+1}-\left[\begin{array}{c}{x}_{1}\\ {y}_{1}\\ {\theta }\\ {\textrm{U}}_{12}\\ {\textrm{U}}_{23}\\ {\textrm{U}}_{31}\end{array}\right]}_{k}\right)}{\Delta \textrm{t}}=0, \end{align}

(27) \begin{align} {x_1}_{k + 1} = {x_1}_k + \Delta {\rm{t}}{c_1},\quad {y_1}_{k + 1} = {y_1}_k + \Delta {\rm{t}}{c_2}, \end{align}

(28) \begin{align} {{\partial {\rm{L}}} \over {\partial {\theta _{k + 1}}}} = {{\partial {\rm{L}}} \over {\partial {{\rm{U}}_{{{12}_{k + 1}}}}}} = {{\partial {\rm{L}}} \over {\partial {{\rm{U}}_{{{23}_{k + 1}}}}}} = {{\partial {\rm{L}}} \over {\partial {{\rm{U}}_{{{31}_{k + 1}}}}}} = 0. \end{align}

Relations (26) to (28) consist of 9 equations that are used to find $ {x}_{{1}_{k+1}},$ $ {y}_{{1}_{k+1}}$ , $ {\theta }_{k+1}$ , $ {\textrm{U}}_{{12}_{k+1}}$ , $ {\textrm{U}}_{{23}_{k+1}}$ , $ {\textrm{U}}_{{31}_{k+1}}$ , $ {\gamma }_{1}$ , $ {\gamma }_{2}$ , and $ {\gamma }_{3}.$ After that, positions of other spheres of the microswimmer can be found according to Eqs. (13) and (21).

Figure 7. The predator and prey paths for control strategy involving control of position and angle. The position of the first sphere (predator) and the prey at t=0, 2, 5, and 10 s has also been illustrated.

Figure 8. Arm length between spheres as a function of time for control strategy involving control of position and angle.

If the objective in each time step is to minimize fluctuations in the control efforts which are the swimmer arm lengths, the following objective function should be used:

(29) \begin{align} {{\rm{J}}_1} = {1 \over 2}\left[ {{{\left( {{{\rm{U}}_{{{12}_{k + 1}}}} - {{\rm{U}}_{{{12}_k}}}} \right)}^2} + {{\left( {{{\rm{U}}_{{{23}_{k + 1}}}} - {{\rm{U}}_{{{23}_k}}}} \right)}^2} + {{\left( {{{\rm{U}}_{{{31}_{k + 1}}}} - {{\rm{U}}_{{{31}_k}}}} \right)}^2}} \right]. \end{align}

If it is desired to minimize arm lengths in the next time step, the following objective function is suitable:

(30) \begin{align} {{\rm{J}}_2} = {1 \over 2}\left[ {{{\left( {{{\rm{U}}_{{{12}_{k + 1}}}} - {\rm{b}}} \right)}^2} + {{\left( {{{\rm{U}}_{{{23}_{k + 1}}}} - {\rm{b}}} \right)}^2} + {{\left( {{{\rm{U}}_{{{31}_{k + 1}}}} - {\rm{b}}} \right)}^2}} \right]. \end{align}

Reference arm length b should not be set too low; otherwise, the swimmer’s arm lengths might decrease below $ 2\;{\mu }\textrm{m},$ resulting in collision between its spheres. Implementing $ {\textrm{J}}_{2}$ results in undesirable shocks and impulses in the beginning of motion. In order to reduce the intensity of these shocks, a linear combination of $ {\textrm{J}}_{1}$ and $ {\textrm{J}}_{2}$ is used instead:

(31) \begin{align} {\textrm{J}}_{3}={\textrm{J}}_{1}+\alpha {\textrm{J}}_{2} \end{align}

The constant $ \alpha $ shoud be small in order to avoid large impulses.

The swimmer’s motion has been simulated for objective function $ {\textrm{J}}_{1}$ and the following prey paths (22) and (23). Control coefficients $ {k}_{p1}$ and $ {k}_{p2}$ are set equal to 1, and initial values are similar to the previous simulation. The results are presented in the following figures:

Figure 9. Arm length between spheres as a function of time for optimal control strategy with objective function $ {J}_{1}.$

The x and y positions of the first sphere are the same as the previous simulation in Fig. 6(a) and (b). Therefore, the predator has again been able to achieve the goal of tracking its prey, and both x and y positions of the first sphere (predator) have settled by $ t=4s.$ In this case, however, the changes in arm lengths have decreased and the maximum arm length has been lowered by $ 5\;{\mu }\textrm{m},$ which can be seen in Fig. 9. The predator and prey paths for this simulation are the same as the previous simulation presented in Fig. 7.

The swimmer’s motion has also been simulated for objective function $ {\textrm{J}}_{3}$ and desired paths (22) and (23). Control coefficients $ {k}_{p1}$ and $ {k}_{p2}$ are set equal to 1, and initial conditions are the same as the previous simulations. In Eq. (31), the constant $ \alpha $ is set at 0.1, and reference arm length b is set at $ 7\;{\mu }\textrm{m}.$ The results are presented in the following figure:

Figure 10. Arm length between spheres as a function of time for optimal control strategy with objective function $ {J}_{3}.$

As before, the predator, which is the swimmer’s first sphere, has been able to successfully hunt down the prey. The x and y positions of the first sphere are similar to previous simulations. Figure 10 illustrates the arm length inputs for the swimmer’s motion. In this case, the arm length between spheres 2 and 3 has decreased quickly near to $ 5\;{\mu }\textrm{m},$ while peak arm length has only slightly increased compared to the previous simulation. The desired and actual paths for the predator and prey are similar to previous simulations.

5. Optimal Control with Constraints on the Arm Lengths

The optimal control strategies presented thus far do not place any constraints on the variations in arm lengths. In order to avoid spheres colliding with each other, or to limit the maximum increase of arm lengths, constraints must be considered for the minimum and maximum admissible values of the arm lengths. Controlling variables $ {x}_{1}$ and $ {y}_{1}$ to follow reference prey paths $ {x}_{1d}$ and $ {y}_{1d},$ while minimizing objective function J with arm length constraints in each time step, is possible with quadratic programming. In order to use the quadratic programming technique, objective function J, system equations, and constraint inequalities can be written in matrix formulation [Reference Tonon, Aronna and Kalise30]:

(32) \begin{align} {\rm{min}}{1 \over 2}{\bf{x}}^T \bf{Hx} + {\bf{f}}^T \textbf{x},\quad {\rm{S}}.{\rm{t}}.{\textbf{A}_{eq}}\textbf{x} = {\textbf{b}_{eq}},\quad {\rm{S}}.{\rm{t}}.\textbf{lb} \le \textbf{x} \le \textbf{ub}. \end{align}

In Eq. (32), $ \,\textbf{x}$ represents the vector containing variables that need to be determined. The term $ \frac{1}{2}{\textbf{x}}^{T}\textbf{Hx}+{\textbf{f}}^{T}\textbf{x}$ is the quadratic objective or cost function which has been written in matrix formulation. The goal of quadratic programming is to minimize the objective function with constraints $ {\textbf{A}}_{eq}\textbf{x}={\textbf{b}}_{eq}$ and $ \textbf{lb}\le \textbf{x}\le \textbf{ub}.$ The term $ {\textbf{A}}_{eq}\textbf{x}={\textbf{b}}_{eq}$ is the equality constraints of the system written in matrix formulation, and $ \textbf{lb}$ and $ \textbf{ub}$ are the lower and upper bounds of the variables.

Here, we have presented the objective function for minimizing arm length change, $ {\textrm{J}}_{1}$ , in matrix form. $ {\textrm{J}}_{1}$ can be expanded in the following form:

(33) \begin{align} {{\rm{J}}_1} = {1 \over 2}\left[ {{{\rm{U}}_{{{12}_{k + 1}}}}^2 + {{\rm{U}}_{{{23}_{k + 1}}}}^2 + {{\rm{U}}_{{{31}_{k + 1}}}}^2} \right] - \left[ {{{\rm{U}}_{{{12}_{k + 1}}}}{{\rm{U}}_{{{12}_k}}} + {{\rm{U}}_{{{23}_{k + 1}}}}{{\rm{U}}_{{{23}_k}}} + {{\rm{U}}_{{{31}_{k + 1}}}}{{\rm{U}}_{{{31}_k}}}} \right]. \end{align}

In Eq. (33), constant terms have been eliminated in the expansion of $ {\textrm{J}}_{1}.$ Therefore, $ \textbf{}\textbf{x}$ , $ \textbf{H}$ , and $ \textbf{f}$ are in the following order:

(34) $$\matrix{ {{\bf{x}} = {{\left[ {\matrix{ {{x_1}} \cr {{y_1}} \cr \theta \cr {{{\rm{U}}_{12}}} \cr {{{\rm{U}}_{23}}} \cr {{{\rm{U}}_{31}}} \cr } } \right]}_{k + 1}},\quad {\bf{H}} = \left[ {\matrix{ 0 & 0 & 0 & 0 & 0 & 0 \cr 0 & 0 & 0 & 0 & 0 & 0 \cr 0 & 0 & 0 & 0 & 0 & 0 \cr 0 & 0 & 0 & 1 & 0 & 0 \cr 0 & 0 & 0 & 0 & 1 & 0 \cr 0 & 0 & 0 & 0 & 0 & 1 \cr } } \right],\quad {\bf{f}} = \left[ {\matrix{ 0 \cr 0 \cr 0 \cr { - {{\rm{U}}_{12}}} \cr { - {{\rm{U}}_{23}}} \cr { - {{\rm{U}}_{31}}} \cr } } \right].} \hfill \cr } $$

The dynamic equations for the system are Eqs. (16), (17), and (18), which if merged together, can be written in the following form:

(35) $$\left[ {\matrix{ 1 & 0 & 0 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 & 0 & 0 \cr {{{\bf{M}}_1}\left( {{\rm{1}},{\rm{1}}} \right)} & {{{\bf{M}}_1}\left( {{\rm{1}},{\rm{2}}} \right)} & {{{\bf{M}}_1}\left( {{\rm{1}},{\rm{3}}} \right)} & {{{\bf{M}}_1}\left( {{\rm{1}},{\rm{4}}} \right)} & {{{\bf{M}}_1}\left( {{\rm{1}},{\rm{5}}} \right)} & {{{\bf{M}}_1}\left( {{\rm{1}},{\rm{6}}} \right)} \cr {{{\bf{M}}_2}\left( {{\rm{1}},{\rm{1}}} \right)} & {{{\bf{M}}_2}\left( {{\rm{1}},{\rm{2}}} \right)} & {{{\bf{M}}_2}\left( {{\rm{1}},{\rm{3}}} \right)} & {{{\bf{M}}_2}\left( {{\rm{1}},{\rm{4}}} \right)} & {{{\bf{M}}_2}\left( {{\rm{1}},{\rm{5}}} \right)} & {{{\bf{M}}_2}\left( {{\rm{1}},{\rm{6}}} \right)} \cr {{{\bf{M}}_3}\left( {{\rm{1}},{\rm{1}}} \right)} & {{{\bf{M}}_3}\left( {{\rm{1}},{\rm{2}}} \right)} & {{{\bf{M}}_3}\left( {{\rm{1}},{\rm{3}}} \right)} & {{{\bf{M}}_3}\left( {{\rm{1}},{\rm{4}}} \right)} & {{{\bf{M}}_3}\left( {{\rm{1}},{\rm{5}}} \right)} & {{{\bf{M}}_3}\left( {{\rm{1}},{\rm{6}}} \right)} \cr } } \right]\left[ {\matrix{ {{{\dot x}_1}} \hfill \cr {{{\dot y}_1}} \hfill \cr {\dot \theta } \hfill \cr {{{\dot U}_{12}}} \hfill \cr {{{\dot U}_{23}}} \hfill \cr {{{\dot U}_{31}}} \hfill \cr } } \right] = \left[ {\matrix{ {{c_1}} \cr {{c_2}} \cr 0 \cr 0 \cr 0 \cr } } \right]$$

Eq. (35) can be rewritten as the following expression:

(36) $$\matrix{ {{\bf{R}}{{\left[ {\matrix{ {\matrix{ {{x_1}\;\;\;{y_1}\;\;\;\theta \;\;\;} \cr } \matrix{ {{{\rm{U}}_{12}}\;\;\;{{\rm{U}}_{23}}\;\;\;{{\rm{U}}_{31}}} \cr } } \cr } } \right]}_{k + 1}}^{\rm{T}} = {\bf{R}}{{\left[ {\matrix{ {\matrix{ {{x_1}\;\;\;{y_1}\;\;\;\theta \;\;\;} \cr } \matrix{ {{{\rm{U}}_{12}}\;\;\;{{\rm{U}}_{23}}\;\;\;{{\rm{U}}_{31}}} \cr } } \cr } } \right]}_k}^{\rm{T}} + \Delta {\rm{t}}{{\left[ {\matrix{ {{c_1}\;\;\;{c_2}\;\;\;\matrix{ {0\;\;\;0\;\;\;0} \cr } } \cr } } \right]}^{\rm{T}}},} \hfill \cr } $$

where $${\bf{R}} = \left[ {\matrix{ 1 & 0 & 0 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 & 0 & 0 \cr {{{\bf{M}}_1}\left( 1 \right)} & {{{\bf{M}}_1}\left( 2 \right)} & {{{\bf{M}}_1}\left( 3 \right)} & {{{\bf{M}}_1}\left( 4 \right)} & {{{\bf{M}}_1}\left( 5 \right)} & {{{\bf{M}}_1}\left( 6 \right)} \cr {{{\bf{M}}_2}\left( 1 \right)} & {{{\bf{M}}_2}\left( 2 \right)} & {{{\bf{M}}_2}\left( 3 \right)} & {{{\bf{M}}_2}\left( 4 \right)} & {{{\bf{M}}_2}\left( 5 \right)} & {{{\bf{M}}_2}\left( 6 \right)} \cr {{{\bf{M}}_3}\left( 1 \right)} & {{{\bf{M}}_3}\left( 2 \right)} & {{{\bf{M}}_3}\left( 3 \right)} & {{{\bf{M}}_3}\left( 4 \right)} & {{{\bf{M}}_3}\left( 5 \right)} & {{{\bf{M}}_3}\left( 6 \right)} \cr } } \right].$$ As a result, $ \textbf{A}_{eq}$ and $ \textbf{b}_{eq}$ are the following:

(37) $$\matrix{ {{{\bf{A}}_{eq}} = {\bf{R}},\quad {{\bf{b}}_{eq}} = {\bf{R}}{{\left[ {\matrix{ {\matrix{ {{x_1}\;\;\;{y_1}\;\;\;\theta \;\;\;} \cr } \matrix{ {{{\rm{U}}_{12}}\;\;\;{{\rm{U}}_{23}}\;\;\;{{\rm{U}}_{31}}} \cr } } \cr } } \right]}_k}^{\rm{T}} + \Delta {\rm{t}}{{\left[ {\matrix{ {{c_1}\;\;\;{c_2}\;\;\;\matrix{ {0\;\;\;0\;\;\;0} \cr } } \cr } } \right]}^{\rm{T}}}.} \hfill \cr } $$

For simulation, the constraints for arm lengths are set between $ 8$ and $ 25\;{\mu }\textrm{m}.$ Therefore, the upper and lower bounds in (32) are the following:

(38) $${{\bf{lb}} = {{\left[ {\matrix{ {\matrix{ { - \infty \;\;\;\;\; - \infty \;\;\;\;\; - \infty \;\;\;\;\;} \cr } \matrix{ {8\;\;\;\;\;8\;\;\;\;\;8} \cr } } \cr } } \right]}^{\rm{T}}},\quad {\bf{ub}} = {{\left[ {\matrix{ {\matrix{ { + \infty \;\;\;\;\; + \infty \;\;\;\;\; + \infty \;\;\;\;\;} \cr } \matrix{ {25\;\;\;\;\;25\;\;\;\;\;25} \cr } } \cr } } \right]}^{\rm{T}}}.}$$

After $ {x}_{{1}_{k+1}},$ $ {y}_{{1}_{k+1}}$ , $ {\theta }_{k+1}$ , $ {\textrm{U}}_{{12}_{k+1}}$ , $ {\textrm{U}}_{{23}_{k+1}}$ , and $ {\textrm{U}}_{{31}_{k+1}}$ are determined in each time step from quadratic programming, the position vector of sphere centers, in the next time step, is determined using Eq. (13).

The swimmer’s motion has been simulated for objective function $ {\textrm{J}}_{1}$ and desired paths (22) and (23), with constraints according to (38). Control coefficients $ {K}_{p1}$ and $ {K}_{p2}$ are set equal to 1, and initial conditions are the same as the previous simulations. The results are presented in the Fig. 11. This figure presents arm length variations for this simulation. The results from Fig. 11 show that arm lengths have remained within their defined boundaries. However, as a result, overall arm length change has slightly increased, which, in turn, means more energy consumption. The predator and prey paths are similar to previous simulations.

The results from various proposed methods have been compared in Table I. According to this table, the last control strategy, optimal control with constraints on the arm lengths, not only minimizes the microswimmer length changes but also maintains the arms lengths in the desired interval.

Table I. The average of objective function $ {J}_{1},$ average of objective function $ {J}_{3}$ for the first second of motion, maximum arm length, minimum arm length, angle control ability, constraints on arm lengths, and settling time, for various simulations

Figure 11. Arm length between spheres as a function of time for optimal control strategy with objective function $ {J}_{1}$ along with arm length constraints.