Over-Actuated Underwater Robots: Configuration Matrix Design and Perspectives

In general, for the configuration designs of underwater robots, the positions and directions of actuators (i.e., thrusters) are given and installed in conventional ways (known points, vertically, horizontally). This yields limitations for the capability of robots and does not optimize the robot’s resources such as energy, reactivity, and versatility, especially when the robots operate in confined environments. In order to optimize the configuration designs in the underwater robot field focusing on over-actuated systems, in the paper, performance indices (manipulability, energetic, reactive, and robustness indices) are introduced. The multi-objective optimization problem was formulated and analyzed. To deal with different objectives with different units, the goal-attainment method, which can avoid the difficulty of choosing a weighting vector to obtain a good balance among these objectives, was selected to solve the problem. A solution design procedure is proposed and discussed. The efficiency of the proposed method was proven by simulations and experimental results.


Introduction
The Actuation System (AS) is an important part of marine robots. The AS groups the different actuators carried by the system. Following the generic Navigation-Guidance-Control (NGC) structure, the AS is in charge of realizing the desired force (F d B ) provided by the control system (see Figure 1). Following Figure 1, the sensorial stage using sensors measurement and prior knowledge of the environment provides to the navigation system the necessary information to compute an estimation of the system state (η). Then, the guidance system uses this estimation and the reference system state (η d ) provided by the mission controller to compute the error function (ε). The control system is then in charge of computing the desired force (F d B ) in order to reduce the error function to zero. Note that classically, this desired force is expressed in the body-frame. Afterwards, the actuation system produces in the environment a resulting force (F B ), which should be as close as possible to F d B . Note that, in this paper, the desired force (F d B ) and resulting force (F B ) are (6 × 1) vectors and include force and torque elements. Inside the AS block, referring to Figure 2, the desired force (F d B ) is the output of the controller. Then, the dispatcher (D) considers the actuator allocation method (and eventually, redundancy management) to compute the desired actuator force (F d m ) that each actuator has to produce. The inverse actuator characteristics are then considered in order to compute the actuator inputs (c m ). Once applied, c m can produce actuator forces (F m ). The resulting force F B is produced with respect to the actuator configuration (A). The properties of the AS are indeed dependent on the actuator configuration (position and attitude of the actuators with respect to the body-frame), actuator dynamics (response characteristics), and dispatcher (control allocation, redundancy management) (see Figure 2) and afford the system different properties. Let us consider in the following that n is the number of Degrees of Freedom (DoFs) of the system and m is the number of actuators. If the system carries less actuators than DoFs, it is said to be underactuated (in that case, A will be an (n × m) matrix where n > m). Long-range Autonomous Underwater Vehicles (AUVs) and, for the terrestrial case, unicycle wheeled vehicles belong to this category [1]. In that case, specific nonlinear guidance strategies have to be used [2]. If the system carries more actuators than DoFs, it is said to be redundant (n < m). Then, there are different solutions (c m ) to produce an identical resulting force (F B ). Indeed, D is one of the multiple possible inverses of A, classically D = A + , where A + is the Moore-Penrose pseudo-inverse. The properties of the AS play a pivotal role in the system performances, in terms of achievable dynamics, maneuverability, robustness, and dependability. The properties of an overactuated system have been studied in aerospace control, where critical safety is required [3], and for marine vehicles [4], where the harsh oceanic condition may easily produce actuator failure. Redundancy was also used in [5] in order to compensate different and unknown actuator responses. The domain of robotic manipulators has also extensively studied this question of redundancy, especially with recent works on humanoid robotics, where the task function approach [6] has been used to concurrently achieve equilibriums [7], walking pattern following [8], and multicontact management [9]. For a global evaluation of an actuation system, we should of course consider many factors, including redundancy management, the control allocation method, the actuator characteristics (inverse and direct), and the actuator configuration. This paper focuses on the study of the actuator configuration; for other problems, the readers can refer to [5] and the references therein.
Different performance criteria related to the actuator configuration design have been proposed. For mobile manipulation, the manipulability index [10] measures the manipulation capability of the end-effector. Intuitively, this index regards the set of all end-effector velocities, which is realizable by joint velocities. This set is called the hyper-manipulability ellipsoid. This index is quantified by computing the hyper-manipulability ellipsoid's properties. Based on these properties, there are different ways to quantify the manipulability index, including the volume of the hyper-manipulability ellipsoid, the ratio of the minimum and maximum radii of the hyperellipsoid, and the minimum radius of the hyperellipsoid. The selection depends on the purpose of evaluation. When the uniformity of manipulating ability is important, the ratio of two radii of the hyperellipsoid is chosen (the optimal value will be close to one). Otherwise, the minimum radius of the hyperellipsoid is suited for the case where the minimum manipulating ability might be critical [11]. Another criterion, attainability [12][13][14], was studied using workspace volume estimation.
In the underwater robotics field, the manipulability index, energetic index, and force index were introduced in [15], and the manipulability index was applied in [16]. Specifically, the manipulability index is used to measure the system's ability to exert a desired force with a specific actuator configuration. Therefore, the closer to one this index is, the better the robot's isotropy is, i.e., the robot can exert the same forces/torques in any direction. The energetic index is a measurement of the variation of system energy when the direction of the desired force changes. This is evaluated by measuring the energy consumption when the direction of a normalized desired force changes over a 3D sphere. The basic idea of the energetic index is to keep the system's energy consumption constant and as low as possible when the direction of action changes. The force index is used to measure the ratio between the actual maximum value and the minimum value of realizing forces. However, these studies only considered a given and fixed actuator configuration. Regarding the design of the actuator configuration of an overactuated underwater robot, a general problem is how to achieve an optimal configuration considering different performance indices. This is a challenging issue that raises two specific questions:

1.
How do we define general and typical indices to evaluate an actuator configuration of an overactuated underwater robot? 2.
How do we solve the complex optimal problem, which is normally nonconvex and has some conflicting objectives?
This paper focuses on the design of the actuator configuration for an overactuated underwater robot with the contributions outlined below:

1.
We propose performance indices to evaluate the actuator configuration of underwater robots; 2.
We optimize an actuator configuration design of an overactuated underwater robot with respect to different performance indices simultaneously.
This paper focuses on the design of an actuator configuration of an overactuated underwater robot, which optimizes different performance indices. Mathematically, an actuator configuration is a mapping from an actuator force vector to a resulting force vector (note that these vectors include force and torque elements). Since we considered an underwater robot equipped with thrusters, the mapping is from a thruster force vector (F m space) to a body-frame vector (F B space) (see Figure 3). The mapping operator is a matrix, which has different names in the literature such as: control effectiveness matrix [4,17], static transformation matrix [18], geometrical distribution of thrusters [19], configuration matrix [16]. In this paper, the mapping of an actuator configuration is called a configuration matrix, denoted as A. The paper is organized as follows. The notations are given in Section 2. The problem formulation and performance indices are described in Section 3. The problem's solution is displayed in Section 4. Simulation results and analyses are depicted in Section 5. Real experiments are depicted in Section 6. Finally, conclusions and future works are discussed in Section 7.

Notation
This section provides most of notations used in the whole paper. However, further notations are introduced when needed. In order to illustrate the notations, a given robot configuration is shown in Figure 4, and detailed explanations are given in Table 1.  Volume of a space rank(.) Rank of a matrix

Problem Formulation
The relation between the desired force (F d B ) and resulting force (F B ) depends on different elements (see Figure 2). This paper only focuses on the actuator configuration. Therefore, three assumptions were considered:

1.
The inverse characteristics and direct characteristics of the actuators are perfectly known, i.e., F d m = F m ; 2.
The dispatcher is the Moore-Penrose pseudo-inverse of the actuators' configuration, i.e., if the actuators' configuration is the A matrix, the dispatcher is D = A + ; 3.
All actuators have the same characteristics.

Model of the Actuator Configuration
This part describes how to model an actuator configuration of an overactuated underwater robot equipped with thrusters. A thruster is modeled by the position and direction of the force produced with respect to the body-frame of the robot. The position of the ith thruster is described by a unit position vector r i and the distance d i to the system's Center of Mass (CM) in the body-frame. The direction of the ith thruster is represented by a unit direction vector u i with respect to the body-frame as in Figure 5, and the ith thruster induces a force with the magnitude denoted as F m,i . The relation of the thruster force vector and resulting force vector (note that this space includes force elements (F) and torque elements (Γ)) is described in Equation (1).
and m is the number of thrusters, m > 6. The configuration matrix A is described as: where A 1 and A 2 ∈ R 3×m are submatrices of A, which correspond to the force and torque elements, respectively. It is obvious to see that τ T i .u i = 0. This is one constraint of the configuration matrix.
In this paper, we assumed that all distances from the thruster positions to the center of the body-frame are the same, d i = d j = const, i, j = 1, . . . , m. Without loss of generality, we can assume that d i = 1, i = 1, . . . , m.

Manipulability Index
As mentioned before, the manipulability index was first introduced in [20] for manipulator mechanisms, and there are different ways to quantify the manipulability index. This paper focuses on the isotropic property of a marine robot. Then, the ratio between the maximum and minimum radii of the manipulability ellipsoid was chosen (see Figure 6). Because of the units' consistency, the matrices that relate to the force space, A 1 , and torque space, A 2 , were investigated separately. However, because of our assumption on d i , the manipulability index is defined as the condition number of the configuration matrix: where σ max and σ min are the maximum and minimum singular values of the configuration matrix, A, respectively. Following Figure 6, the manipulability index investigates the resulting force ellipsoid, which is realizable by the thruster forces (F m ) such that F m ≤ 1 (see Theorem in Appendix A). If I m = 1, the robot is isotropic, or if I m = ∞, the robot cannot act along at least one direction.

Energetic Index
Energy is very important for marine robots, and the energy consumption of robots depends on many factors such as the mechanical design, the environmental effects, and the specific mission. In order to evaluate the energy performance of an underwater robot, the energetic index was introduced in [15]. Being different from this, in our paper, the twonorm of the thruster force vector, p E = F m 2 , was used to quantify the energy that an underwater robot consumes to produce forces and torques, which can be calculated as follows in Equation (4): The energetic index is proposed to measure the variation of the energy consumption of an underwater robot when the direction of the desired force changes. It is quantified by computing the energy consumption when a unit desired force vector, (F d B ), changes over the unit hypersphere (see Figure 7 for a 3D sphere). Because of the units' consistency, however, the force and torque sphere were computed separately.
For the force sphere case, the unit desired force vector includes a unit vector of force elements and a zero vector of torque elements. For the torque sphere case, the unit desired force vector includes a zero vector of force elements and a unit vector of torque elements. Intuitively, this can be expressed as: where u s = [cos α cos β sin α cos β sin β] T is a unit vector in spherical coordinates with α ∈ [−π, π] and β ∈ [−π/2, π/2].
According to these two cases, the norm of the thruster force vector was also divided into two cases as follows: The energetic index is defined as: where S is the area of the three-dimensional sphere, p E f , p EΓ are the subvectors of p E corresponding to the force sphere and the torque sphere case, respectively, and w e f and w eΓ are the weighting coefficients. Note that the weighting coefficients were chosen to normalize the difference between the force and torque spheres. These choices depends on the robot's characteristics. However, because of the normalization of the spheres, it is normal to assign one to these coefficients.

Workspace Index
The term workspace volume was first introduced in [13] for manipulator mechanisms. In this paper, the work space index was used to measure the volume of the attainable regions of the resulting force space with respect to (w.r.t.) the body-frame. In general, the characteristics of thrusters always have limitations, namely saturations and dead-zones (in this index, the dead-zone is not considered). This yields the polytope of the thruster force space, the F m space, denoted as M. By properly choosing the configuration matrix, A = (A 1 A 2 ) T , the volume of the resulting force space for the force, the F F space, and the resulting force space for the torque, the F T space, can be maximized (see Figure 8). Note that the resulting spaces for the force and torque were studied separately because of the units' consistency.
In general, the set M of thruster forces is known (with the given saturations of thrusters), so M is a polytope and F F and F T are also polytopes (under the A 1 and A 2 linear transform actions). We define the workspace index as: where Vol is the volume of a polytope and ω w f and ω wτ are the weighting coefficients. In control perspectives, the larger the space's volumes are, the less control effort is need. The design objective was to maximize the workspace index, I w . Normally, the set M is convex and its vertices are known. It is easy to find the vertices of F F and F T . Of course, F F and F T are also convex sets (because of the linear transformation). This problem becomes a volume computation of convex polytopes.

Reactive Index
The reactive index quantifies how fast the actuation system is able to change the orientation of the resulting force F B (ideally, F d B ). Suppose that the robot is traveling in a direction with a set of thruster forces F m1 induced from desired force vector F d B1 . The robot wants to change to another direction (or the same direction with a different magnitude) with the desired force vector F d B2 , so the thrusters have to produce another set of thruster forces F m2 . The two-norm of the deviation of the thruster forces, F m = F m1 − F m2 = [ F m1 F m2 · · · F mm ] T , is considered as the reactive capability of the robot. Referring to the approximation of the characteristics of the thrusters, as Figure 9, the change from F m1 to F m2 is closer than that from F m1 to F m3 (in the linear section, the dead-zone of the thruster characteristics is neglected in this paper). Hence, we have: From Equation (11), the sensitivity of the thruster forces with respect to the desired forces, in other words the variation of the thruster forces w.r.t. the desired forces, is upperbounded by the norm of the pseudo-inverse of the configuration matrix, A + . We define the reactive index as: It is obvious to see that if this index is lower, the robot is more reactive. Then, the objective of the design process is to minimize the reactive index.

Robustness Index
This criterion measures the robustness level of the AS of an underwater robot. This means that if any thruster of the robot fails, the remaining ones can still perform the robot's mission. In particular, for any F d B vector, there always exists a F m vector to satisfy the equation F B = AF m , and F B is as close as possible to F d B . We have: where a i is the ith column of the matrix A and F m,i is the force magnitude of the ith thruster. When one or more thrusters completely fail, the value of F m,i = 0. Note that in the case that the ith thruster has partly failed, the value of F m,i remains small (not addressed in this paper). This is equivalent, as we considered that a corresponding column a i of the configuration matrix A equals the zero vector. Therefore, Equation (13) is equivalent to: where the A matrix is the A matrix with one or more corresponding columns equal to the zero vectors. We discuss hereafter two questions: What are the conditions of the matrix A to guarantee the robustness? What is the maximum number of failed thrusters?
To address these two questions, we supposed that k thrusters fail, and Equation (14) is a linear equation system with six equations (the dimensions of F B are 6 × 1) and (m − k) variables, because the matrix A is 6 × m, where k columns are zero vectors. It is obvious to see that if rank(A ) = 6, for a given F d B , there always exits F m such that F B = A F m and F B is as close as possible to F d B . This can be interpreted as m − k ≥ 6 or k ≤ m − 6. The condition on the configuration matrix and that on the maximum number of failed thruster that guarantee the robustness of an underwater robot are stated as:
Robustness condition: the rank of the configuration matrix always equals six, i.e., rank(A) = 6, if any columns, from one to a maximum of (m − 6), of the A matrix equal the zero vectors. If rank(A) < 6, the system becomes underactuated, and the guidance and control have to change to guarantee the robot's mission. This problem is not addressed in this paper.
We define the robustness index as: where A| ≤m−6 is the A matrix where the maximum number of columns being zero is (m − 6). This index is verified in the solving process of the problem.

Configuration Matrix Design Problem
With all the performance indices discussed above, we propose the following design problem: is the objective function vector. A is the feasible set of the configuration matrix (A) including the constraints of the configuration matrix (A) and the robustness index. The reciprocal of the workspace index, 1 I w , is in Equation (16), because we wanted to maximize the workspace index.
This is a multi-objective optimization problem, and the unique solution belongs to the convexity of each objective function in the objective vector and the feasible set, A. Note that this optimization problem is with respect to a matrix variable (matrix optimization), not a vector variable. However, the optimization techniques for vector variables (vector optimization) can be applied here because we do not lose the physical meaning when converting a matrix variable to a vector variable in this optimization problem (because of the independency of each column in the matrix derived from the independent positions and orientations of the thrusters).
Specifically, Equation (16) can be rewritten as: The problem (17) is to minimize an objective vector V(A), including the manipulability index, the energetic index, the reciprocal of the workspace index, and the reactive index, with respect to configuration matrix, A, and to satisfy the constraints of the matrix structure itself and the robustness index. It is clear that this is a nonconvex and multiobjective optimization problem, which normally has many solutions. In the following sections, we propose a mathematical analysis and a method for solving the multi-objective optimization problem.

Problem Solution
Our final objective was to find a distribution (position and orientation) of the thrusters of an underwater robot. This means obtaining the u i and r i vectors for i = 1, 2, . . . , m. These vectors can be extracted from configuration matrix A, which is the solution of the problem (17). Recall that our problem (17) is the multi-objective optimization problem with nonconvexity, and theoretically, this problem has infinitely many Pareto-optimal solutions. Our objective was to find one Pareto-optimal solution to build the robot. Analyzing the underlying mathematical properties of the problem helped us simplify the solving process. Thus, the mathematical analysis of the problem is shown in the next section.

Mathematical Analysis
The configuration matrix A has the form: We have: B is an m × m symmetric matrix where each element is denoted as b ij . We have: where λ i is the ith eigenvalue of matrix B. From Equations (19) and (20), we have: In the case of manipulability index optimization, the condition of configuration matrix A is one, cond(A) = 1. This means that the maximum singular value equals the minimum singular value, σ max = σ min . Note that the matrix A is the n × m matrix with n < m. The matrix A has n nonzero singular values (we have to guarantee that rank(A) = n), then the matrix B has n nonzero eigenvalues and m − n zero eigenvalues.
In the optimization case of the manipulability index, cond(A) = 1 ⇒ σ max = σ min , we have λ i = λ max = λ min = λ (σ = √ λ). Equation (21) is rewritten as: Considering the fact that τ i 2 ≤ 1, we have: Therefore, we have λ max = 2 m n when τ i 2 = 1. In the singular-value decomposition of a matrix, when cond(A) = 1, the matrix A can be written as: where U ∈ R n×n , V ∈ R m×m are orthogonal matrices, The pseudo-inverse of matrix A is A + and can be written as: where Our objective for the reactive index was to minimize A + . From Equation (25), the reactive index I re = A + = 1 σ , and the minimum value of the reactive index is equivalent to the maximum value of σ. This leads to the equality of Equation (23).
In order to minimize the reactive index and the manipulability index, the configuration matrix A has the following structure: where S(n × m) is like-diagonal and σ = √ λ = 2 m n ; U(n × n) and V(m × m) are orthogonal matrices (UU T = I, VV T = I). This result can be used as the initial value of the numerical optimization process and is useful for solving the problem.
We continue to discuss the energetic index. First, we introduce a proposition as follows: Proof. We have: With M = PΣQ T : We have: where I is a q × q identity matrix. Replacing Equation (29) to (28), we have: The energetic index is stated as: Choosing w e f = w eΓ = 1 (because the desired force vectors, F d B ( f ), F d B (τ), are units), we have: In the case in which the reactive index and the manipulability index are minimum, the configuration matrix A(n × m) has the form of Equation (26); therefore, the pseudoinverse matrix A + (m × n, m > n) has the following structure: where V, U are orthogonal matrices. It is clear that matrix A + satisfies the condition of Proposition 1. Applying this proposition, we have: . Therefore, Equation (32) becomes: From the aforementioned mathematical analysis of the energetic index, we can see that the energetic index belongs to the norm of the pseudo-inverse of the configuration matrix, I re = 2 A + , when the configuration matrix A has the form of (26).
We then discuss the upper-bound of the workspace index. For the units' consistency, the workspace index for the force space and that for the torque space were investigated separately, denoted as I w f and I wτ , respectively. Recall that the objective of the workspace index is to maximize the volume of the resulting force space (F B space), including the resulting space for the force and the resulting space for the torque with given thrusters' force (the F m space).
The fact that for all vectors F m ∈ R m , AF m ≤ A F m . The volume of the resulting force space is maximum when the equality holds.
Following Figure 10, the volumes of the resulting force spaces (F B ) (the force and torque spaces) are always less than the volumes of the exterior hypersphere of F B spaces of the force and torque (this may be the circumscribed spheres or not). This means that:  The volume of a Euclidean ball of radius R in n-dimensional Euclidean space is [21]: where (2k + 1)!! = 1.3.5, . . . , (2k − 1).(2k + 1).
Proof. We have: On the other hand: From (37) and (38), we have: where I 1 and I 2 are partitioned matrices of matrix I. From (39) and the uniqueness of singular-value decomposition [22], it is obvious to see that the structures of A 1 and A 2 are the same as (26) with different dimensions. Therefore, cond(A 1 ) = cond(A 1 ) = 1 and A 1 = A 2 = σ.
From (35) and (36) and Proposition 2, it is obvious to obtain the upper-bound of the resulting spaces of the force and torque of the system and then the upper-bound of the workspace index. Normally, the weighting coefficients in the workspace index are chosen as one because of our assumption of d i .

Problem Solution
Based on the above mathematical analysis, the goal-attainment method was chosen to solve the problem with the given desired values. The idea of this method is to minimize the deviation of the desired values and the obtained values. One advantage of the goalattainment method is that the problem does not need to be normalized to a dimensionless problem. The solution of this method has been proven to be Pareto-optimal. This method is also suitable when the feasible objective set is nonconvex [23]. All Pareto-optimal solutions may be found by changing the attainment vector; however, this depends on the properties of the problem.
Our problem using the goal-attainment method becomes: is the desired objective vector, w is an attainment vector, which can be chosen. The goal-attainment method with two objective functions is illustrated in Figure 11. By altering w vector, we searched for the Pareto-optimal solutions. Therefore, our solving process included two phases: 1. Phase 1: Find one Pareto-optimal solution of the configuration matrix with the goalattainment method; 2.
Phase 2: Check the robustness index of the chosen solution in Phase 1.
The optimization toolbox in the MATLAB environment was used to solve our problem. Note that our problem was formulated as a multi-objective optimization problem. One objective has one desired value excluding the robustness index, which is as a constraint, and therefore, the desired vector is set up. The goal-attainment method was used to solve the problem. An attainment vector was chosen as a trade-off between the underachievement and overachievement of the objective functions. In multi-objective optimization, an optimal solution depends on a decision maker. Theoretically, there is no method for this choice. In our work, this vector was selected by trial and error. In particular, for the manipulability, reactivity, and energetic indices, we know exactly the desired values, so the corresponding values in the attainment vector were chosen as zero, which means that these are hard constraints. For the workspace index, we only know the upper-bound of the desired value; therefore, a positive value was chosen for underachievement. Figure 11. Goal-attainment method with two objective functions.

Simulation Results
We designed an overactuated underwater robot with m = 8 thrusters and n = 6 degrees of freedom. Two cases were simulated: the general case and the given position case. In the general case, we have to identify both the positions and orientations of eight thrusters optimizing the performance indices. In the given position case, the thrusters are installed at the corners of a cube, and we only have to determine the directions of the thrusters. In this simulation, the thruster characteristics were chosen as in [5], then the maximum and minimum values of the thrusters forces were F imax = 1.1N and F imin = −0.4N, respectively. The desired values of the performance indices were subsequently I d m = 1, I d e = 1.2248, I d wF = 597.7, I d wT = 597.7, I d re = 0.6124 (σ max = 2 m n = 1.6330; see Table 2 for more details).

General Case
In this case, the robot is called a ball robot, and the positions and orientations of the thrusters are not known. The problem (40) is solved as follows. to choose this attainment vector. It was chosen by trial and error. By our approach, we found that some values of this vector can be assigned zero (this imposes hard constraints), except the workspace index (because of the upper-bound value).
The simulation results are shown in Figures 12 and 13a,b. The configuration matrix A and optimal values are shown in Table 3. Specifically, in Figure 12, the positions of the thrusters are at the top of the blue line, and the orientations of the thrusters are shown as the red arrow. Furthermore, we can see that the isotropic property of the robot is guaranteed (see Figure 13a,b) with the sphere shapes of the attainable spaces of the forces and torques. From Table 3, the obtained values of the manipulability index, the energetic index, and the reactive index were almost the same as the desired values. However, the obtained workspace index was smaller than the desired one.   In this phase, the robustness index was checked. The optimal configuration matrix A in Table 3 satisfies the robustness constraint. Specifically, the maximum number of thrusters that are acceptable (for failures) is two.

Given Position Case
In this case, the robot is called the cube robot, and the positions of thrusters are at the corners of the cube. We only had to find their orientations. The number of variables in the problem (40) was reduced. The desired objective vector and attainment vector were the same as in general case. The results are presented in the sequel.

Phase 1
The optimization toolbox was used to solve our problem, and the simulation results are shown in Figures 14 and 15a,b and Table 4. The directions of the thrusters are depicted as red arrows in Figure 14. Being similar to the general case, the isotropic property is also guaranteed in this case (see Figure 15a,b). One Pareto-optimal configuration matrix is shown in Table 4. We can see that the obtained objective values in Table 4 are the same as the general case.  The optimal configuration matrix A in Table 4 satisfies the conditions of the robustness index. Similarly, the maximum number of thrusters that can fail is two.

A Comparison of Two Configurations
In this section, a comparison of two configurations is illustrated. The choice of the configurations corresponds to a real robot (cube robot), which is used in the experiments in the next section. The first one is a normal configuration (denoted as C 1 ) in which the thrusters are distributed vertically or horizontally (in practice, this configuration is easier to install, as Figure 21 shows). The configuration matrix of the C 1 configuration, denoted as A 1 , is shown in Equation (41).
The second one (denoted as C 2 ) is an optimal configuration, denoted as A 2 , which is a solution of the optimization problem (given position case) thanks to the thruster characteristics of BlueRobotics (Figure 16), and the optimal configuration matrix is shown in Equation (42). Note that the configuration matrices A 1 and A 2 were calibrated with the corresponding geometrical properties of the real cube robot at the LIRMM Institute, Montpellier University. The attainable force space and torque space corresponding to the two configurations C 1 and C 2 are illustrated in Figure 17a,b. It is obvious that the C 2 configuration is more isotropic than the C 1 configuration. However, for some specific points of the attainable force and torque spaces, the C 1 configuration is better than the C 2 configuration. Thanks to the properties of matrices A 1 and A 2 (Equations (41) and (42)) and the thruster characteristics (Figure 16), Table 5 shows the values of the performance indices for both configurations. The performances of the C 2 configuration are better than C 1 . Because of the calibration (the distance d i is different between the motors), the manipulability index (I m ) is larger than one (Note that, theoretically, the distances of all thrusters with respect to the center of mass were assumed the same (without loss of generality, they were assigned one). However, in practice, for our cube robot, these distances were not completely the same, and we had to calibrate the configuration matrix.). Table 5. Comparison between the two configurations (I ro shows the maximum number of thrusters that can fail to make sure that rank(A = 6)).

No.
Indices In order to verify the attainability of the two configurations (workspace index), incremental torques were applied about the u-, v-, and w-axis, respectively (Figures 18a, 19a and 20a), and the corresponding Pulse-Width Modulation (PWM) inputs (c m ) of the eight thrusters were computed. The results are shown in Figures 18b,c, 19b,c and 20b,c, in which the two PWM saturation values of the thrusters (upper saturation value: 1900, lower saturation value: 1100) are plotted with two bold lines. We can see that the performances of the robot with the two configurations are almost the same for the rotation about the u-and v-axis. However, the C 2 configuration showed better performance for the rotation about the w-axis. In fact, the thrusters with the C 1 configuration reached saturations very earlier in comparison with the thrusters with the C 2 configuration (Figure 20b,c).
In order to validate the robustness of the optimal configuration (C 2 ) in comparison with the normal configuration (C 1 ), the rank of matrices A 1 and A 2 was checked when one or two arbitrary columns have been nullified. When the resulting matrices are rank deficient, this means that the robustness is not guaranteed because one direction is not actuated. Therefore, we cannot control all 6 DoFs independently. The robustness index in Table 5 shows the checking results. In particular, when the fifth thruster of the C 1 configuration fails, the robustness is not guaranteed.

Experimental Results
Experiments were carried out on the cube robot to compare between the two configurations, C 1 (see Figure 21) and C 2 (see Figure 22), in the swimming pool at Montpellier University. The cube in the water and a video link for the cube's operations can be seen in Figure 23.

Attainability Validation
The incremental torques about the u-axis, v-axis, and w-axis were applied to cube robot respectively, and angular velocities and PWM input values were stored to evaluate these two configurations. For safety, the experiments were stopped when one thruster reached the saturation values. The experimental results are shown in Figures 24-26. For rotating about the u-axis (Figure 24), the attainability of configurations C 1 and C 2 was almost the same: all thrusters operated in a feasible region. Otherwise, for rotating about the v-axis or w-axis, the attainability of configuration C 2 was better than that of C 1 . In particular, with the v-axis experiment (Figure 25), the cube robot with C 1 stopped the mission earlier than with C 2 (at Time Step 771) because one thruster reached saturation. The same thing happened with the w-axis experiment (at Time Step 451) (see Figure 26).

Energetic Validation
In this section, we verify the energy consumption during these experiments for the two configurations. An energy-like criterion is proposed: where m is the number of thrusters, T is the time of the experiment, and PW M i (t) is the PWM inputs of the ith thruster. Table 6 shows the energy consumption of the robot during the three rotation experiments. For u-axis rotation, the attainability of the two configurations was the same, but the energy consumption of C 2 was lower than that of C 1 . For v-axis and w-axis rotations, the duration of the experiments of C 2 was longer than that of C 1 , and the energy consumption, therefore, was higher.  Table 7 shows the comparison of the energy consumption of the two configurations with the same time duration. For the v-axis rotation, the energy value of C 2 was lower than that of C 1 . However, for the w-axis, the energy value of C 2 was higher. This happened because the robot dived deeper in C 2 in the experiment of the w-axis rotation, and the robot had to deliver more power to maintain a greater constant depth.

Robustness and Reactive Validation
This section validates the robustness and reactivity of the optimal configuration (C 2 ) in comparison to the normal one (C 1 ). For robustness, the robot performed a mission, and one or two thrusters were turned off. For the normal configuration C 1 , the mission would fail, and for the optimal configuration C 2 , the mission would be guaranteed. Specifically, for the robustness index, we carried out the following experiments: The cube robot dives to a predefined depth with all motors being in the normal operating conditions; 2.
The cube robot dives to the same predefined depth with one vertical motor being stopped; 3.
The cube robot dives to the same predefined depth with two vertical motors being stopped; 4.
The cube robot dives to the same predefined depth with three motors being stopped (two vertical motors and one arbitrary motor); 5.
The cube robot simultaneously dives to the same predefined depth and rotates about the w-axis with three motors being stopped (two vertical motors and one horizontal motor) For the reactive index, we measured how fast the robot changed missions. The following experiments were carried out:

1.
The cube robot goes down to the predefined depth and goes up to another predefined depth, then goes down again to the former predefined depth; 2.
In the sequel, the cube robot goes down to the predefined depth, rotates about the u-axis, and after that, rotates about the v-axis. The rotation time of each axis should be 60 s or longer; 3.
Next, the cube robot goes down to the predefined depth, rotates about the u-axis, and after that, rotates about the diagonal-axis (diagonal of the cube robot). The rotation time of each axis should be 60 s or longer.
The experimental results for the robustness validation of C 1 and C 2 are shown in Figures 27-29. In the case of one or two motors stopped, the depth control performances of C 1 and C 2 were almost the same (see Figure 27). The differences are clear in the case of three thrusters stopped (Figure 29): the performance of C 1 was not guaranteed ( Figure 28) and violations of the PWM values occurred (see Figure 29a).  The results for the reactive validation are shown in Figures 30-32. We measured the reactive time of the angular velocities when the directions of the cube's actions changed. It is clear that the reactive time of C 2 was faster than that of C 1 . Specifically, the reactive time is the region formed by the vertical dashed lines in Figures 30-32. It is obvious that the reactive time of C 2 was smaller than that of C 2 (see Figures 31 and 32).

Conclusions and Future Work
In this paper, an approach for designing an optimal configuration matrix (which depends on the positions and directions of the thrusters) of overactuated underwater robots was presented. The performance indices (related to manipulability, energy, workspace, reactivity, and robustness) were proposed and analyzed. Specifically, the manipulability index shows the isotropic properties of a robot; the energetic index minimizes the energy consumption under some assumptions; the workspace index is related to the attainable spaces (i.e., the force and torque spaces) of the robot; the reactive index presents how fast the robot changes the direction of the resulting actuation force; finally, the robustness index is related to the capacity of the robot to maintain its performance in the case of failures (i.e., some thrusters are completely stopped). It was formulated as a multi-objective optimization problem. Because the different indices exhibit different magnitudes and physical meanings, the goal-attainment method was chosen to find one Pareto-optimal solution. Simulation and experimental results showed that the performances of the optimal configuration were better than a "normal" configuration, which is often used (thrusters are installed vertically or horizontally). Because of the nonconvexity of the problem, finding all Pareto-optimal solutions, the Pareto front, remains a challenging problem and will be future work. Moreover, a design problem relaxing the assumptions (i.e., perfectly known characteristics of the actuators, pseudo-inverse dispatcher) is also an interesting direction for future research.  Acknowledgments: The authors would like to thank Numev Labex, MUSE, Montpellier University; Region Occitanie; and FEDER for supporting this research.

Conflicts of Interest:
The authors declare no conflict of interest.