A graph is used to represent the communication topology of an interacting system of agents. The communication graph is built upon a set of nodes V = { 1 , … , N } that represent agents. An edge denoted by the pair (k, j) exists between agent k ∈ V and j ∈ V if information flows from j to k. The set of all edges is denoted E ⊆ V 2 . The set of nodes V and the edges E define a graph G = ( V , E ) (Diestel, 2000). A sequence of edges {(k, k ₁), (k ₁, k ₂), …, (k _ξ , j)} with distinct nodes k l ∈ V , k _l ≠ k, k _l ≠ j, for l = 1, 2, …, ξ is called a path from node k to node j. A graph G is called undirected if ( k , j ) ∈ E implies ( j , k ) ∈ E . For a undirected graph, the set of neighbors to node k is denoted N k = { v ∈ V : ( k , v ) ∈ E } . If there exists a path between any pair of distinct nodes k , j ∈ V , then an undirected graph G is called connected. Edges are expressed using the adjacency matrix A ∈ R N × N , where the entry on the kth row and jth column is A k j = 1 if k , j ∈ E 0 otherwise .

The degree matrix D ∈ R N × N encodes how many unique edges are connected to each node and has nonzero elements on the diagonal, i.e., D k j = ∑ n = 1 N A k n if k = j 0 otherwise .

The symmetric and positive semi-definite Laplacian matrix L ∈ R N × N associated with the undirected graph G is L = D − A. The Laplacian is used to compute relative state information communicated between agents. The quadratic form z ^T L z ≥ 0, where z ∈ R N , is equal to zero if and only if z _k = z _j , for all k , j ∈ V .

The self-propelled particle model (Justh and Krishnaprasad, 2004) has often been used to describe the collective motion of N planar vehicles that move at a constant speed with steering controls inputs. The planar position of the kth particle with respect to the origin of the inertial frame is expressed using complex coordinates as r k = x k + i y k ∈ C , where k ∈ V . The dynamics of the kth particle are r ̇ k = v k e i θ k θ ̇ k = u k , (1) where, for the kth particle, v k ≜ x ̇ k 2 + y ̇ k 2 ∈ R is a constant speed, θ k ≜ a t a n ( y ̇ k / x ̇ k ) ∈ T is the orientation of the velocity (also called the phase of the particle), T ≜ S is the torus, and u k ∈ R is the steering control. The unit vector e i θ k is called the phasor of particle k and is aligned with its heading, whereas i e i θ k is perpendicular to the heading (Figure 2A). For a constant speed, v _k = v ₀, and a constant turn-rate, θ ̇ k = ω 0 , the particle moves on a circle with radius | v 0 ω 0 − 1 | and center c k = r k + i v 0 ω 0 − 1 e i θ k . This fixed-radius circle will later serve as a reference for stabilizing circular formations.

When referring to the positions, phase arrangement, reference circle centers, and control inputs of the collective of N particles, we use bold letters, i.e., r ≜ [ r 1 , … , r k ] T ∈ C N , θ ≜ [ θ 1 , … , θ k ] T ∈ T N , c ≜ [ c 1 , … , c k ] T ∈ C N , and u = [ u 1 , … , u k ] T ∈ R N , respectively. Similarly, e i θ ≜ [ e i θ 1 , … , e i θ k ] T ∈ C N . For complex numbers, z 1 , z 2 ∈ C , the inner product is defined as ⟨ z 1 , z 2 ⟩ = R e { z 1 ̄ z 2 } , where z 1 ̄ is the complex conjugate of z ₁. This inner product is equivalent to the standard inner product on R 2 . For complex vectors, z , y ∈ C N , the inner product is similarly defined as ⟨ z , y ⟩ = ∑ i = 1 N R e { z i ̄ y i } . The modulus of a complex number is denoted | ⋅ | = ⟨ ⋅ , ⋅ ⟩ .

Cooperative control laws for stabilizing the collective motion of identical, unit-speed, self-propelled particles in parallel or circular formations have been extended to include an external flow field (Paley, 2008), motion on spherical surfaces (Paley, 2009), and various communication topologies (Sepulchre et al., 2007; Sepulchre et al., 2008). For parallel formations, all particles are synchronized when they have equal and constant phase, θ = θ ₀ 1, where 1 = [1,…,1]^T is the N-by-1 vector of ones, for some constant θ 0 ∈ T . For synchronization, the relative positions of particles are arbitrary. For circular formations, all particles move in the same direction along the same circle, that is, θ ̇ = ω 0 1 for some constant ω ₀ and c = c ₀ 1 for some constant c 0 ∈ C . In a circular formation, the relative phases of the particles are arbitrary.

Parallel and circular formations may be achieved using Lyapunov-based control design to minimize a potential function for a desired formation. Consider the Laplacian parallel formation potential (Paley, 2007) U p θ ≜ 1 2 ⟨ e i θ , L e i θ ⟩ , (2) which is minimized when the agents are synchronized. Assume the Laplacian matrix L corresponds to a time-invariant, connected, and undirected graph G representing the communication topology of the agents. The time-derivative of U _p ( θ ) along trajectories of Eq. 1 is (Paley, 2007) U ̇ p θ = ∑ k = 1 N ∂ U p θ ∂ θ k ∂ θ k ∂ t = ∑ k = 1 N ⟨ i e i θ k , L k e i θ ⟩ u k , (3) where L _k is the kth row of the Laplacian matrix. The term L _k e ^{i θ} is the sum of the phasor of the kth agent relative to the phasors of all connected agents, i.e., L k e i θ = | N k | e i θ k − ∑ j ∈ N k e i θ j . Choosing the gradient control (Paley, 2007) u k = − K ⟨ i e i θ , L k e i θ ⟩ , (4) for K > 0, makes Eq. 3 negative semi-definite and drives U _p ( θ ) to zero so that agents converge to the set of synchronized parallel formations.

Similarly, to achieve a circular formation, the Laplacian circular formation potential (Paley, 2007) U c r , θ ≜ 1 2 ⟨ c , L c ⟩ (5) may be used. The potential U _c ( r , θ ) has a minimum value when the agents are in a circular formation. The time-derivative of U _c ( r , θ ) along trajectories of the self-propelled particle Eq. 1 is (Paley, 2007) U ̇ c r , θ = v 0 ω 0 − 1 ∑ k = 1 N v 0 − 1 ω 0 v k − u k ⟨ e i θ j , L j c ⟩ . (6)

Choosing the circular formation control (Paley, 2007) u k = v 0 − 1 ω 0 v k + K 0 ⟨ e i θ k , L k c ⟩ (7) makes Eq. 6 negative semi-definite and drives U _c ( r , θ ) towards zero so that the agents’ time-averaged circle centers coincide to a common point.

The Chaplygin sleigh is a canonical nonholonomic mechanical system consisting of a rigid body moving in the plane that is supported by two frictionless sliding points and a single knife edge that allows no motion perpendicular to its edge (Bloch, 2003). Previous studies have demonstrated that a fish robot driven by an internal reaction wheel can be modeled as a Chaplygin sleigh due to the nonholonomic constraint imposed by the Kutta condition (Kelly et al., 2012), (Lee et al., 2019), which constrains the fluid flow at the trailing edge. As the reaction wheel spins back and forth, it flaps the robot’s body, which interacts with the surrounding fluid to generate thrust.

Consider a system of N fish robots each modeled as a Chaplygin sleigh with the following dynamics in state-space form (Lee et al., 2019): r ̇ k = v k e i θ k θ ̇ k = ω k v ̇ k = l ω k 2 − d v k ω ̇ k = − m l v k b ω k − u k b , (8) where r k ∈ C is the position of the trailing edge of the fish robot (Figure 2B), v k ∈ R is the swimming speed, θ k ∈ T is the velocity orientation, ω k ∈ R is the angular rate of the kth fish, and u k ∈ R is the applied torque, where k = 1, …, N. Furthermore, d ≥ 0 is the drag coefficient, and m > 0, l > 0, and b > 0 are the mass, length, and moment of inertia, respectively. Unlike the self-propelled particle Eq. 1, the speed of the Chaplygin sleigh Eq. 8 is not constant and the control input is a torque rather than an angular rate.

Prior work has established that the Chaplygin-sleigh model exhibits limit-cycle dynamics under open-loop periodic control inputs (Pollard et al., 2019), as well as feedback control (Lee et al., 2019) (Free et al., 2020). Consider the feedback control (Lee et al., 2019) u k = b − K 1 ω k − K 2 ⁡ sin θ ̄ k − θ k , (9) where θ ̄ k is the desired heading angle of the kth fish, and K ₁, K ₂ > 0 are feedback gains. Substituting Eq. 9 into Eq. 8 yields the closed-loop system (Lee et al., 2019) r ̇ k = v k e i θ k θ ̇ k = ω k v ̇ k = l ω k 2 − d v k ω ̇ k = − m l b v k ω k + K 1 ω k + K 2 ⁡ sin θ ̄ k − θ k . (10)

The system Eq. 10 can be divided into a slow and fast subsystem (Lee et al., 2019), where the fast v _k -subsystem (Lee et al., 2019), v ̇ k = d l d ω k 2 − v k ,

converges to v k → l / d ω k 2 for a sufficiently large drag coefficient d. Let a = ml ²/bd > 0. The slow (θ _k , ω _k )-subsystem becomes (Lee et al., 2019) θ ̇ k = ω k ω ̇ k = − a ω k 3 + K 1 ω k + K 2 ⁡ sin θ ̄ k − θ k . (11)

Observe that Eq. 11 gives the equations of motion of a pendulum with nonlinear damping and natural frequency K 2 (Lee et al., 2019). The system Eq. 11 has two equilibrium points corresponding to a fish robot heading that is parallel ( θ k , ω k ) = ( θ ̄ k , 0 ) and anti-parallel ( θ k , ω k ) = ( θ ̄ k − π mod 2 π , 0 ) to the desired heading. Both equilibria are unstable and the system exhibits a stable limit cycle centered on ( θ ̄ k , 0 ) in the (θ _k , ω _k ) plane (Lee et al., 2019). The corresponding limit cycle of Eq. 10 is evident in the (v _k , ω _k ) plane as well. The limit cycle propels the robot in the desired direction by flapping; however, the limit cycle is achieved only for certain values of the control gains K ₁ and K ₂ (Lee et al., 2019). The average swimming velocity is proportional to K ₁, but if K ₁ is too large, then the angular rate in the resulting limit cycle does not switch signs and the robot spins in a circle (Lee et al., 2019). The control law Eq. 9 that enables each fish robot to swim in a desired direction can be modified, with interactions from neighboring fish, to achieve collective motion of the school, as described next.

Consider a collection of N identical fish robots modeled by the Chaplygin sleigh system Eq. 8. Assume a sufficiently large drag coefficient so that v k → ( l / d ) ω k 2 and the (θ _k , ω _k ) dynamics follow Eq. 11. For the purposes of control design, the simplified Chaplygin sleigh system Eq. 8 becomes r ̇ k = l / d ω k 2 e i θ k θ ̇ k = ω k ω ̇ k = − a ω k 3 − u k b . (12)

Inspired by the Laplacian parallel formation potential Eq. 2 for the self-propelled particle, consider the potential V p θ , ω = 1 2 ω T ω + 1 2 N K 2 ⟨ e i θ , L e i θ ⟩ . (13)

The time-derivative of V _p ( θ ) is V ̇ p θ , ω = ω ̇ T ω + 1 N K 2 ⟨ d d t e i θ , L e i θ ⟩ , (14) where, along trajectories of Eq. 12, ω ̇ T ω = ∑ k = 1 N − a ω k 3 − b − 1 u k ω k , (15) and 〈 d d t e i θ , L e i θ 〉 = ∑ k = 1 N 〈 i e i θ k , L k e i θ 〉 ω k . (16) By choosing the control u k = b − K 1 ω k + K 2 N ⟨ i e i θ k , L k e i θ ⟩ = b − K 1 ω k + K 2 N ∑ j ∈ N k sin θ j − θ k , (17) and substituting Eqs 15–17 into Eq. 14, V ̇ p ( θ , ω ) becomes V ̇ p θ , ω = ∑ k = 1 N − a ω k 2 + K 1 ω k 2 . (18)

The feedback control law Eq. 17 relies only on relative-state measurements between agents and does not include feedback linearization of the agents’ dynamics. Recall K ₁, K ₂ > 0 are control gains. Since Eq. 18 is a summation of quartic functions with roots at ω _k = 0 and | ω k | = K 1 / a , then V ̇ < 0 outside Ω p = { ( θ , ω ) ∈ T N × R N : | ω k | ≤ K 1 / a ∀ k ∈ V } . Therefore, all trajectories are trapped in Ω _p . The gain K ₂ in Eq. 17 is chosen to ensure forward flapping motion for Eq. 8, as discussed in Section 2.3.

The control law Eq. 17 is illustrated by numerical simulation using control Eq. 17 and the full dynamics Eq. 8. The simulation was conducted for 150 s with N = 8 robots using the parameters listed in Table 1. The robots were initialized with random headings and zero linear and angular velocities. A communication range of 3 m determined the communication topology, which remained invariant during the simulation based on the agent’s random initial positions. Figures 3A,B show all N robots converging to the same limit cycle in the (θ _k , ω _k ) and (v _k , ω _k ) planes. As a result, all robots move in the same direction (on average), as shown in Figure 3C. The parallel potential, V _p (t), initially decreases (Figure 3D), but instead of converging to zero, it oscillates around a fixed value as the robots converge to different phases on the same limit cycle.

The parallel formation control Eq. 17 is based on the forward swimming control Eq. 9; however, a desired heading was not prescribed, but rather the average heading emerged through interactions among agents depending on their initial conditions. Similarly, a circular formation control is proposed here that drives the fish robots to continuously adjust their heading at a known average rate, while aligning the center position of the nominal circles to an (average) consensus value. Virtual fish can be introduced to achieve a reference heading or position (Sepulchre et al., 2007, 2008). Consider the following circular formation potential inspired by Eq. 5: V c t , r , θ , ω = 1 2 ω T ω − K 2 γ T 1 + 1 2 K 3 ⟨ c , L c ⟩ , (19) where γ = [ γ 1 , … , γ k ] T with γ _k (t) = cos(θ _k − ω ₀ t) and K ₂, K ₃ > 0. The time-derivative of V _c along trajectories of Eq. 12 is V ̇ c t , r , θ , ω = ω ̇ T ω − K 2 γ ̇ T 1 + K 3 ⟨ c ̇ , L c ⟩ , (20) where γ ̇ T 1 = − ∑ k = 1 N ⁡ sin ( θ k − ω 0 t ) ( ω k − ω 0 ) , and ω ̇ T ω is given in Eq. 15.

When averaged over time, the motion of a fish robot resembles that of a self-propelled particle. Recall that the reference circle center for a self-propelled particle is c k = r k + v 0 ω 0 − 1 i e i θ k . Since the average swimming speed of the robots is bK ₁/ml (Lee et al., 2019), then by setting v ₀ = bK ₁/ml, the parameter ω ₀ may be chosen to yield an average turn rate and reference circle with radius | v 0 ω 0 − 1 | . However, ω ₀ should be sufficiently small to ensure that ω _k switches signs along the limit-cycle so that the robots flap. The third term in Eq. 20, evaluated along trajectories of Eq. 12, is 〈 c ̇ , L c 〉 = ∑ k = 1 N l d ω k − v 0 ω 0 ω k 〈 e i θ k , L k c 〉 . (21) Therefore Eq. 20 may be rewritten as V ̇ c t , r , θ , ω = ∑ k = 1 N − a ω k 3 − b − 1 u k ω k + K 2 ⁡ sin θ k − ω 0 t ω k − ω 0 + K 3 l d ω k − v 0 ω 0 ω k ⟨ e i θ k , L k c ⟩ . (22) Choosing the control u k = b − K 1 ω k + K 2 ⁡ sin θ k − ω 0 t + K 3 l d ω k − v 0 ω 0 ⟨ e i θ k , L k c ⟩ ) , (23) the derivative Eq. 22 becomes V ̇ c = ∑ k = 1 N − a ω k 2 + K 1 ω k 2 − K 2 ω 0 ⁡ sin θ k − ω 0 t ≤ ∑ k = 1 N − a ω k 2 + K 1 ω k 2 + K 2 | ω 0 | . (24)

Thus, V ̇ c < 0 outside Ω c = { ( θ , ω ) ∈ T N × R N : ω k 2 ≤ ( K 1 + K 1 2 + 4 a K 2 | ω 0 | ) / 2 a ∀ k ∈ V } , which implies the system is driven to a bounded set containing the desired limit cycle. The K ₃ term in Eq. 23 biases the torque to align the reference circles, whereas the remaining terms produce a flapping motion with a given average turn rate. Although Eq. 24 remains negative for any positive K ₃, this gain must be chosen sufficiently small to ensure flapping. As with Eq. 17, the feedback control Eq. 23 relies only on relative-state measurements between agents and does not include feedback linearization of the agents’ dynamics.

The circular formation feedback control Eq. 17 is numerically illustrated by simulating Eq. 8 with parameters from Table 1 and using ω ₀ = 0.05 rad/s. The simulation was conducted for ten minutes to demonstrate circular motion. Figure 4A shows all N robots converge to the same limit cycle in the (v _k , ω _k ) plane. The orbits in Figure 4B resemble the limit cycle in Figure 3B; however, due to the time-varying term in Eq. 19, they translate along the perimeter of the phase cylinder. The net result is motion along a circle whose center position is determined by the initial conditions of the agents (Figure 4C). Since the oscillating centers are aligned when all robots have identical position and phase, the controller drives all the robots to one side of the circle. In ongoing work, we seek to stabilize symmetric circular formations (Sepulchre et al., 2007; Sepulchre et al., 2008), which evenly distribute the agents around the formation. The circular potential exhibits low and high frequency oscillations (Figure 4D), which correspond to motion around the reference circle and flapping, respectively.

The fish-inspired soft robots used in the experiments (Figure 1) are each driven by a Pololu 12 V DC motor (with a 4.4 to 1 gear ratio) that oscillates a reaction wheel located at each robot’s center of mass. Each fish robot measures its orientation with an onboard BNO055 inertial measurement unit (IMU) sensor that features a built-in extended Kalman filter. A micro-SD card is used to store sensor data and a 900 MHz xBee radio enables each fish robot to communicate with a ground station and other fish robots on the water’s surface. Each fish robot performs onboard sensor processing and control with a Teensy 3.2 microcontroller. For a more detailed discussion of the fish robots’ design, refer to (Lee et al., 2019).

The experiments were conducted in a 367,000 gallon water tank at the Neutral Buoyancy Research Facility at the University of Maryland, College Park. An overhead camera is mounted above the experimental area to record the true position of each robotic fish; however, the position data is not used in real time by the fish for formation control. Instead, during formation control experiments, each robot exchanges orientation data from their onboard IMU with other fish in the school using the xBee radios. Since the position of each robotic fish is not computed in real-time, an invariant complete communication graph was used in the experiments; whereas, a proximity-based communication graph was used in simulation.^AT The overhead camera images are post-processed after each experiment to visualize the trajectory of each robotic fish. The image processing uses MATLAB’s built-in corner/object detection based on a minimum eigenvalue algorithm (The MathWorks, 2019) and a built-in constant velocity Kalman filter for object tracking. Since all of the sensing and control law computations occur onboard, the school of robotic fish are a self-contained system when performing formation control experiments. The block diagram in Figure 5 gives an overview of the experimental testbed.

The nonlinear control laws Eq. 17 or Eq. 23 generate a reference torque command that steers the fish robots to their respective formations. However, this torque cannot be commanded directly since the control input into the reaction wheel system is the voltage applied to the DC motor generated by the motor driver and Teensy micro-controller. Furthermore, the only available measurement is the angular velocity of the motor’s shaft measured by an encoder. Thus, to track the reference torque a linear-quadratic-Gaussian (LQG) controller and estimator was implemented. The LQG controller assumes the following motor dynamics (Kim, 2017): Λ = μ d i d t + R i + e (25) e = K e ψ ̇ (26) τ = K τ i (27) J ψ ̈ = τ − ζ m ψ ̇ , (28) where Λ is the voltage input, R is the motor’s electrical resistance, μ is the inductance, i is the current, e is the motor’s back electromotive force (EMF), ψ ̇ is the angular rate of the output shaft, J is the sum of inertias between the reaction wheel and motor’s output shaft, τ is the torque applied by the motor, ζ _m is the internal damping friction applied to output shaft, and K _e and K _τ are the motor’s generator and torque constants, respectively. To identify the values of the DC motor and reaction wheel parameters, a series of system identification experiments were conducted, as described next.

The resistance of the motor R was measured directly with an ohm meter, and the sum of the motor and reaction wheel inertias J was approximated using analytical expressions for the moment of inertia of a cylinder about its axis of symmetry with a known mass and diameter. To determine the motor generator constant K _e , substitute Eq. 26 into Eq. 25 and examine the system at steady current state, i.e., di/dt = 0. Solving for ψ ̇ yields ψ ̇ = 1 K e Λ − R i . (29)

Using a Rigol DP711 power supply and quadrature encoder attached to the motor’s output shaft, a series of constant voltages were applied to the motor. For each voltage, the current through the motor and the angular velocity of the output shaft at steady state were recorded. A linear regression was used to determine the quantity 1/K _e , which is the slope of the line in Eq. 29 and Figure 6A.

The internal frinction coefficient ζ _m was found by setting the torque to zero, i.e., τ = 0 in Eq. 28, and solving the resulting differential equation to obtain ln ψ ̇ t = − ζ m J t + ln ψ ̇ 0 , (30) where ψ ̇ ( 0 ) is the initial shaft angular velocity. Notice that Eq. 30 is a linear equation that describes ln ( ψ ̇ ) as a function of time. A similar technique as described previously can be used to find − ζ _m /J from a time-history of torque-free motor data. To identify the damping parameter, the motor was initialized with a constant nonzero angular velocity and the input voltage was removed. The motor’s angular velocity was recorded as it decays under internal friction. As before, a linear regression was used to determine the slope − ζ _m /J in Eq. 30 from Figure 6B and, hence, ζ _m can be inferred since J is known.

To determine the torque constant K _τ , substitute Eq. 27 into Eq. 28 at steady state ( ψ ̇ constant) and solve for ψ ̇ to obtain ψ ̇ = K τ ζ m i . (31)

By repeating the procedure used to determine K _e , using the value of ζ _m , the value of K _τ is found from the slope K _τ /ζ _m in Eq. 31 and Figure 6C.

Since the inductance μ only plays a role in the transient response of the motor, which are sufficiently fast, we estimate this parameter heuristically. The simulated response of the motor to a sinusoidal voltage input is visually compared to the actual response of the motor under the same input. The process is repeated while adjusting the value of μ to obtain a similar response, as shown in Figure 6D.

Lastly, the motor was found to exhibit a range of deadband voltages near zero that resulted in the motor being unresponsive. To determine the range of this deadband, a series of incrementally increasing voltages were applied to the motor, giving an approximate deadband range of ± 4V. The parameter values determined through this system identification process are summarized in Table 2.

To implement the LQG torque controller and state estimator for the DC motor, convert Eqs 25–28 into state space form where q = τ , ψ ̇ T denotes the states of the motor and the output is ϒ = ψ ̇ . The continuous-time state space equations take the form q ̇ = A q + B Λ (32) ϒ = C q (33) where A = − R μ − K e K τ μ 1 J − ζ m J , B = K τ μ 0 , and C = 0 1 . (34)

For implementation onboard the micro-controller, the continuous system Eqs 32, 33 is converted into a discrete-time system (with an addition of torque process noise and heading measurement noise): q k + 1 = Φ q k + Γ Λ k + ξ w k (35) ϒ k = C q k + η k , (36) where (Crassidis and Junkins, 2011) Φ = e A Δ t , Γ = ∫ 0 Δ t e A τ d τ B , ξ = 1 0 , (37) k is an integer indexing the discrete state, Δt = 100Hz is the time-step of the microcontroller, w _k is zero-mean, Gaussian, additive process noise with variance σ w 2 , and η _k is zero-mean, Gaussian, additive measurement noise with variance σ η 2 . We adopt the standard approach outlined in (Crassidis and Junkins, 2011) and (Brown and Hwang, 1997) to implement a discrete-time Kalman filter for the system Eqs 35, 36. The estimated state q ̂ k is used in a linear-quadratic-regulator (LQR) feedback control law of the form Λ k = κ T ( q d − q ̂ k ) , where the κ is a gain matrix found by solving the discrete-time algebraic Riccati equation and q _d are the desired state values. Numerical values used in the LQG controller and estimator design are given in Table 3.

To evaluate the LQG torque controller/estimator, the control law Eq. 9 was used to generate a reference torque for a step input change in desired heading. The performance of the torque tracking controller is shown in Figure 7A, and the heading trajectory of the fish robot during this experiment is examined in the next section.

We first examine the accuracy of the closed-loop directional controller Eq. 9, where θ _d is a desired heading. After testing Eq. 9 on a single robotic fish, we found that the DC motor’s deadband does not allow the average heading to completely converge to θ _d . When the fish robot’s heading is close to the desired heading, the voltage required to close the gap is too small and falls within the dead zone making the motor unresponsive. The accuracy of Eq. 9 is limited on our testbed due to this voltage deadband, but can also be mitigated by choice of the values of the nonlinear control gains K ₁ and K ₂. In both Eq. 9 and Eq. 17, increasing the value of K ₂ makes the control laws more sensitive to relative heading errors and increases the torque required to minimize the error. Similarly K ₁ amplifies the torque and, therefore the voltage input required to keep the fish swimming. Tuning these control gains increases the accuracy of the aforementioned control laws. The experimental results of this test with a step input for θ _d is shown in Figure 7B. The fish robot’s heading oscillates about a desired heading with a small persistent error between the average and desired headings. The control gains used in this experiment are K ₁ = 0.5 and K ₂ = 7.

To validate the theoretical results of the parallel formation control law Eq. 17, the results from six experiments are presented. The fish robots’ micro-controller uses IMU measurements for the heading and xBee radios for communication within the school. An overhead camera observes the positions of each fish robot; these were not used by the school during the experiments.

Four fish robots were initialized with random initial positions and orientations at the beginning of each parallel formation control experiment (implemented with gains K ₁ = 3 and K ₂ = 5). Consensus was achieved by the four fish with a small phase shift (up to 20 degrees) of the mean heading for each fish. This offset may be attributed to excess noise in angular velocity measurements and the aforementioned voltage dead zone. An example heading time-history from one of the experiments is shown in Figure 8A.

The excess noise in the angular velocity measurements caused large spikes to appear when computing the Lyapunov potential Eq. 13 from experimental data. Thus, to illustrate the convergence of this potential, the 1 2 ω T ω term was removed, leaving only the heading alignment term K 2 2 N < e i θ , L e i θ > plotted in Figure 8B for each of the six experiments. As mentioned in Section 3, the potential function does not decrease completely to zero when consensus is achieved because each robot converges to a different phase on the limit cycle.

An example trajectory of the fish robots in Figure 9A and the Supplementary Material S1 shows the fish robots achieving consensus and swimming in a parallel formation. However, since the ground truth orientations of the fish robots are not used, errors in the IMU measurements sometimes cause a subset of the school to swim in a different direction than the rest of the school, as shown in Figure 9B. The error in IMU measurements may be caused by an erroneous calibration or magnetic anomalies in the water tank interfering with the onboard magnetometer.