# Over-Actuated Underwater Robots: Configuration Matrix Design and Perspectives

Huu-Tho Dang, Lionel Lapierre, René Zapata, Benoît Ropars, Pascal Lepinay

► To cite this version:

Huu-Tho Dang, Lionel Lapierre, René Zapata, Benoît Ropars, Pascal Lepinay. Over-Actuated Underwater Robots: Configuration Matrix Design and Perspectives. Sensors, 2021, 21 (22), pp.7729-7760. 10.3390/s21227729. lirmm-03587052

HAL Id: lirmm-03587052

https://hal-lirmm.ccsd.cnrs.fr/lirmm-03587052

Submitted on 24 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Over-Actuated Underwater Robots: Configuration Matrix Design and Perspectives †

Tho Dang <fnref n="1" />*, Lionel Lapierre <fnref n="1" />*, Rene Zapata <fnref n="1" />, Benoit Ropars <fnref n="2" /> and Pascal Lepinay <fnref n="1" />

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Dang, T.; Lapierre, L.; Zapata, R.; Ropars, B.; Lepinay, P. Over-Actuated Underwater Robots: Configuration Matrix Design and Perspectives. *Sensors* **2021**, *21*, 7729. https://doi.org/10.3390/s21227729

Academic Editors: Reza Ghabcheloo
and Antonio M. Pascoal

Received: 15 October 2021

Accepted: 18 November 2021

Published: 20 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** 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.

**Keywords:** over-actuated underwater robots; multi-objective optimization; underwater robots; performance indices

## 1. 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_B^d$) provided by the control system (see Figure 1).

<figure><img src="image_4.png" /><figcaption>**Figure 1.** NGC structure augmented with the actuation system and sensorial stage.</figcaption></figure>

<footnote n="1">¹ Laboratory of Informatics, Robotics and MicroElectronics (LIRMM) (UMR 5506 CNRS—UM), Université Montpellier, 161 rue Ada, CEDEX 5, 34392 Montpellier, France; zapata@lirmm.fr (R.Z.); lepinay@lirmm.fr (P.L.)</footnote>

<footnote n="2">² Reeds Company, 199 rue Hélène Boucher, 34170 Castelnau-Le-Lez, France; ropars@lirmm.fr</footnote>

<footnote>* Correspondence: danghuu@lirmm.fr (T.D.); lapierre@lirmm.fr (L.L.)</footnote>

<footnote>† The manuscript is extended by our publication Dang, HT.; Lapierre, L.; Zapata, R.; Lepinay, P.; Ropars, B. Configuration Matrix Design of Over-Actuated Marine Systems. In Proceedings of the OCEANS 2019, Marseille, France, 17–20 June 2019; pp. 1–10, doi:10.1109/OCEANSE.2019.8867445.</footnote>
Following Figure 1, the sensorial stage using sensors measurement and prior knowl-
edge of the environment provides to the navigation system the necessary information to
compute an estimation of the system state ($\hat{\eta}$). Then, the guidance system uses this esti-
mation and the reference system state ($\eta^d$) provided by the mission controller to compute
the error function ($\epsilon$). The control system is then in charge of computing the desired force
($F_B^d$) 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_B^d$. Note that, in this paper,
the desired force ($F_B^d$) 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_B^d$) is the
output of the controller. Then, the dispatcher ($\mathcal{D}$) considers the actuator allocation method
(and eventually, redundancy management) to compute the desired actuator force ($F_m^d$) 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, $\mathcal{D}$ is one of the multiple possible
inverses of **A**, classically $\mathcal{D} = \mathbf{A}^+$, where $\mathbf{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 *overactu-
ated* 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].

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Actuation system scheme.</figcaption></figure>

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 manipula-
tion 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–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<sub>m</sub> space) to a body-frame vector (F<sub>B</sub> 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.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Actuator configuration mapping.</figcaption></figure>
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.

## 2. 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.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** A given robot configuration.</figcaption></figure>

<table><caption>**Table 1.** Notations.</caption><tr><td><b>A</b></td><td>Configuration matrix</td></tr><tr><td><b>A</b><sup>+</sup></td><td>Moore-Penrose pseudo-inverse of <b>A</b> matrix</td></tr><tr><td><b>u</b><sub>i</sub></td><td>(3 × 1) unit vector of the direction of the <i>i</i>th thruster</td></tr><tr><td><b>r</b><sub>i</sub></td><td>(3 × 1) unit vector of the position of the <i>i</i>th thruster</td></tr><tr><td><b>F</b><sub>m</sub></td><td>(<i>m</i> × 1) force vector of <i>m</i> thrusters</td></tr><tr><td><b>F</b><sub>m</sub><sup>d</sup></td><td>(<i>m</i> × 1) desired force vector of <i>m</i> thrusters</td></tr><tr><td><i>F</i><sub>m,i</sub></td><td>Force magnitude of the <i>i</i>th thruster</td></tr><tr><td><b>F</b><sub>B</sub><sup>d</sup></td><td>(6 × 1) desired force (force and torque elements) w.r.t. the body-frame</td></tr><tr><td><b>F</b><sub>B</sub> = (<b>F</b><sub>B</sub>/<b>Γ</b>)</td><td>(6 × 1) resulting force (force and torque elements) w.r.t. the body-frame</td></tr><tr><td><b>c</b><sub>m</sub></td><td>(<i>m</i> × 1) input vector of thrusters</td></tr><tr><td>⊗</td><td>Cross product</td></tr><tr><td>||·||</td><td>Euclidean norm</td></tr><tr><td>||·||<sub>p</sub></td><td>p-norm</td></tr><tr><td><i>m</i></td><td>Number of thrusters</td></tr><tr><td><i>n</i></td><td>Number of Degrees of Freedom (DoFs)</td></tr><tr><td><b>F</b></td><td>(3 × 1) vector of force elements in the resulting force <b>F</b><sub>B</sub></td></tr><tr><td><b>Γ</b></td><td>(3 × 1) vector of torque elements in the resulting force <b>F</b><sub>B</sub></td></tr><tr><td><b>D</b></td><td>Dispatcher</td></tr><tr><td><i>d</i><sub>i</sub></td><td>Distance of the <i>i</i>th thruster to the center of the body-frame</td></tr><tr><td>cond(<b>A</b>)</td><td>Condition number of the matrix <b>A</b></td></tr><tr><td>Vol(.)</td><td>Volume of a space</td></tr><tr><td>rank(.)</td><td>Rank of a matrix</td></tr></table>
### 3. Problem Formulation

The relation between the desired force ($F_B^d$) 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_m^d = F_m$;

2. The dispatcher is the Moore-Penrose pseudo-inverse of the actuators' configuration, i.e., if the actuators' configuration is the $\mathbf{A}$ matrix, the dispatcher is $\mathbf{D} = \mathbf{A}^+$;

3. All actuators have the same characteristics.

#### 3.1. Model of the Actuator Configuration

This part describes how to model an actuator configuration of an overactuated under-water 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 $i$th thruster is described by a unit position vector $\mathbf{r}_i$ and the distance $d_i$ to the system's Center of Mass (CM) in the body-frame. The direction of the $i$th thruster is represented by a unit direction vector $\mathbf{u}_i$ with respect to the body-frame as in Figure 5, and the $i$th 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 ($\mathbf{F}$) and torque elements ($\mathbf{\Gamma}$)) is described in Equation (1).

<figure><img src="image_8.png" /><figcaption>Figure 5. Actuator configuration model.</figcaption></figure>

$$ \mathbf{F}_B = \mathbf{A}\mathbf{F}_m = \begin{pmatrix} \mathbf{F} \\ \mathbf{\Gamma} \end{pmatrix} \qquad (1) $$

where $\mathbf{F}_B = [F_u \ F_v \ F_w \ F_p \ F_q \ F_r]^T \in \mathbb{R}^6$, $\mathbf{A} \in \mathbb{R}^{6 \times m}$, $\mathbf{F}_m = [F_{m,1} \ F_{m,2} \ \dots \ F_{m,m}]^T \in \mathbb{R}^m$, and $m$ is the number of thrusters, $m > 6$. The configuration matrix $\mathbf{A}$ is described as:

$$ \mathbf{A} = \begin{pmatrix} \mathbf{u}_1 & \mathbf{u}_2 & \cdots & \mathbf{u}_m \\ d_1 \mathbf{r}_1 \otimes \mathbf{u}_1 & d_2 \mathbf{r}_2 \otimes \mathbf{u}_2 & \cdots & d_m \mathbf{r}_m \otimes \mathbf{u}_m \end{pmatrix} \\ = \begin{pmatrix} \mathbf{u}_1 & \mathbf{u}_2 & \cdots & \mathbf{u}_m \\ \tau_1 & \tau_2 & \cdots & \tau_m \end{pmatrix} = \begin{pmatrix} \mathbf{A}_1 \\ \mathbf{A}_2 \end{pmatrix} \qquad (2) $$

where $\mathbf{A}_1$ and $\mathbf{A}_2 \in \mathbb{R}^{3 \times m}$ are submatrices of $\mathbf{A}$, which correspond to the force and torque elements, respectively. It is obvious to see that $\tau_i^T \cdot \mathbf{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 = \text{const}, i, j = 1, \dots, m$. Without loss of generality, we can assume that $d_i = 1, i = 1, \dots, m$.

## 3.2. 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, $\mathbf{A}_1$, and torque space, $\mathbf{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:

$$I_m = \text{Cond}(\mathbf{A}) = \frac{\sigma_{max}}{\sigma_{min}} \qquad (3)$$

where $\sigma_{max}$ and $\sigma_{min}$ are the maximum and minimum singular values of the configuration matrix, $\mathbf{A}$, respectively.

<figure><img src="image_9.png" /><figcaption>Figure 6. Manipulability ellipsoid with mapping.</figcaption></figure>

Following Figure 6, the manipulability index investigates the resulting force ellipsoid, which is realizable by the thruster forces ($\mathbf{F}_m$) such that $\|\mathbf{F}_m\| \le 1$ (see Theorem in Appendix A). If $I_m = 1$, the robot is isotropic, or if $I_m = \infty$, the robot cannot act along at least one direction.

## 3.3. 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 two-norm of the thruster force vector, $p_E = \|\mathbf{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):

$$p_E = \|\mathbf{F}_m\|_2 = \sqrt{\sum_{i=1}^{m} F_{mi}^2} = \|\mathbf{A}^+ \cdot \mathbf{F}_B^d\|_2 \qquad (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, ($\mathbf{F}_B^d$), 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:

$$ \mathbf{F}_B^d = \begin{pmatrix} \mathbf{F} \\ \mathbf{\Gamma} \end{pmatrix} = \begin{cases} (\mathbf{u}_s), & \text{for the force sphere} \\ (\mathbf{0}), & \text{for the torque sphere.} \end{cases} \qquad (5) $$

where $\mathbf{u}_s = [\cos \alpha \cos \beta \quad \sin \alpha \cos \beta \quad \sin \beta]^T$ is a unit vector in spherical coordinates with $\alpha \in [-\pi, \pi]$ and $\beta \in [-\pi/2, \pi/2]$.

According to these two cases, the norm of the thruster force vector was also divided into two cases as follows:

$$ p_E = \begin{cases} p_{Ef} = \|\mathbf{A}^+(\mathbf{u}_s)\|, & \text{for the force sphere case} \\ p_{EG} = \|\mathbf{A}^+(\mathbf{0})\|, & \text{for the torque sphere case.} \end{cases} \qquad (6) $$

<figure><img src="image_10.png" /><figcaption>Figure 7. The rotation of the unit desired vector in the 3D sphere.</figcaption></figure>

The energetic index is defined as:

$$ I_e = \frac{1}{S} \int_S (w_{ef} p_{Ef} + w_{e\Gamma} p_{EG}) dS \qquad (7) $$

where S is the area of the three-dimensional sphere, $p_{Ef}$, $p_{EG}$ are the subvectors of $p_E$ corresponding to the force sphere and the torque sphere case, respectively, and $w_{ef}$ and $w_{e\Gamma}$ 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.

### 3.4. Workspace Index

The term workspace volume was first introduced in [13] for manipulator mechanisms. In this paper, the workspace 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 $\mathbf{F}_m$ space, denoted as $\mathbf{M}$. By properly choosing the configuration matrix, $\mathbf{A} = (\mathbf{A}_1 \mathbf{A}_2)^T$, the volume of the resulting force space for the force, the $\mathbf{F}_F$ space, and the resulting force space for the torque, the $\mathbf{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 $\mathbf{M}$ of thruster forces is known (with the given saturations of thrusters), so $\mathbf{M}$ is a polytope and $\mathbf{F}_F$ and $\mathbf{F}_T$ are also polytopes (under the $\mathbf{A}_1$ and $\mathbf{A}_2$ linear transform actions). We define the workspace index as:

$$ I_w = \omega_{wf} \mathrm{Vol}(\mathbf{F}_F) + \omega_{w\Gamma} \mathrm{Vol}(\mathbf{F}_T) \qquad (8) $$
where Vol is the volume of a polytope and $\omega_{wf}$ and $\omega_{wT}$ are the weighting coefficients.

In control perspectives, the larger the space's volumes are, the less control effort is needed. The design objective was to maximize the workspace index, $I_w$. Normally, the set $\mathcal{M}$ is convex and its vertices are known. It is easy to find the vertices of $\mathbb{F}_F$ and $\mathbb{F}_T$. Of course, $\mathbb{F}_F$ and $\mathbb{F}_T$ are also convex sets (because of the linear transformation). This problem becomes a volume computation of convex polytopes.

<figure><img src="image_11.png" /><figcaption>Figure 8. Space mapping ($v_i$ is denoted as the vertex).</figcaption></figure>

### 3.5. Reactive Index

The reactive index quantifies how fast the actuation system is able to change the orientation of the resulting force $\mathbb{F}_B$ (ideally, $\mathbb{F}_B^d$). Suppose that the robot is traveling in a direction with a set of thruster forces $\mathbb{F}_{m1}$ induced from desired force vector $\mathbb{F}_{B1}^d$. The robot wants to change to another direction (or the same direction with a different magnitude) with the desired force vector $\mathbb{F}_{B2}^d$, so the thrusters have to produce another set of thruster forces $\mathbb{F}_{m2}$. The two-norm of the deviation of the thruster forces, $\Delta \mathbb{F}_m = \mathbb{F}_{m1} - \mathbb{F}_{m2} = [\Delta \mathbb{F}_{m1} \Delta \mathbb{F}_{m2} \cdots \Delta \mathbb{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 $\mathbb{F}_{m1}$ to $\mathbb{F}_{m2}$ is closer than that from $\mathbb{F}_{m1}$ to $\mathbb{F}_{m3}$ (in the linear section, the dead-zone of the thruster characteristics is neglected in this paper). Hence, we have:

$$ \Delta \mathbb{F}_m = \mathbf{A}^+ (\mathbb{F}_{B1}^d - \mathbb{F}_{B2}^d) = \mathbf{A}^+ \Delta \mathbb{F}_B^d \quad (9) $$

$$ \| \Delta \mathbb{F}_m \| = \| \mathbf{A}^+ \Delta \mathbb{F}_B^d \| \le \| \mathbf{A}^+ \| \| \Delta \mathbb{F}_B^d \| \quad (10) $$

$$ \frac{\| \Delta \mathbb{F}_m \|}{\| \Delta \mathbb{F}_B^d \|} \le \| \mathbf{A}^+ \| \quad (11) $$

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 upper-bounded by the norm of the pseudo-inverse of the configuration matrix, $\|\mathbf{A}^+\|$. We define the reactive index as:

$$ I_{re} = \| \mathbf{A}^+ \| \quad (12) $$

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.
<figure><img src="image_12.png" /><figcaption>Figure 9. Thruster characteristics’ approximation.</figcaption></figure>

### 3.6. 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_B^d$ 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_B^d$.

We have:

$$F_B = AF_m = \sum_{i=1}^{m} a_i F_{m,i} \qquad (13)$$

where $a_i$ is the $i$th column of the matrix $A$ and $F_{m,i}$ is the force magnitude of the $i$th thruster.

When one or more thrusters completely fail, the value of $F_{m,i} = 0$. Note that in the case that the $i$th 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:

$$F_B = A'F_m \qquad (14)$$

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 \times 1$) and $(m-k)$ variables, because the matrix $A'$ is $6 \times m$, where $k$ columns are zero vectors. It is obvious to see that if $rank(A') = 6$, for a given $F_B^d$, there always exits $F_m$ such that $F_B = A'F_m$ and $F_B$ is as close as possible to $F_B^d$. This can be interpreted as $m-k \ge 6$ or $k \le 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:

1. The maximum failed thrusters: $m-6$;

2. 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:

$$I_{ro} = rank(A|_{\le m-6}) = 6 \qquad (15)$$
where $A|_{\le 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.

### 3.7. Configuration Matrix Design Problem

With all the performance indices discussed above, we propose the following design problem:

$$ \begin{array}{ll} \min_{\mathbf{A}} & \mathbf{V}(\mathbf{A}) = \min_{\mathbf{A}} [I_m I_e \frac{1}{I_w} I_{re}]^T \\ \text{s.t} & \mathbf{A} \in \mathcal{A} \end{array} \qquad (16) $$

where $\mathbf{V}(\mathbf{A})$ is the objective function vector. $\mathcal{A}$ is the feasible set of the configuration matrix ($\mathbf{A}$) including the constraints of the configuration matrix ($\mathbf{A}$) and the robustness index. The reciprocal of the workspace index, $\frac{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, $\mathcal{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:

$$ \begin{array}{ll} \min_{\mathbf{A}} & \mathbf{V}(\mathbf{A}) = \min_{\mathbf{A}} [I_m I_e \frac{1}{I_w} I_{re}]^T \\ \text{s.t} & \|\mathbf{u}_i\| = 1, i = 1, 2, \dots, m \\ & \|\boldsymbol{\tau}_i\| \le 1, i = 1, 2, \dots, m \\ & \boldsymbol{\tau}_i^T \mathbf{u}_i = 0, i = 1, 2, \dots, m \\ & I_{re} = \mathrm{rank}(\mathbf{A}|_{\le m-6}) = 6 \end{array} \qquad (17) $$

The problem (17) is to minimize an objective vector $\mathbf{V}(\mathbf{A})$, including the manipulability index, the energetic index, the reciprocal of the workspace index, and the reactive index, with respect to configuration matrix, $\mathbf{A}$, and to satisfy the constraints of the matrix structure itself and the robustness index. It is clear that this is a nonconvex and multi-objective 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.

## 4. Problem Solution

Our final objective was to find a distribution (position and orientation) of the thrusters of an underwater robot. This means obtaining the $\mathbf{u}_i$ and $\mathbf{r}_i$ vectors for $i = 1, 2, \dots, m$. These vectors can be extracted from configuration matrix $\mathbf{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.

### 4.1. Mathematical Analysis

The configuration matrix $\mathbf{A}$ has the form:

$$ \mathbf{A} = \begin{pmatrix} \mathbf{u}_1 & \mathbf{u}_2 & \cdots & \mathbf{u}_m \\ \tau_1 & \tau_2 & \cdots & \tau_m \end{pmatrix} \qquad (18) $$
We have:

$$ \mathbf{B} = \mathbf{A}^T \mathbf{A} = \begin{pmatrix} \mathbf{u}_1 & \mathbf{u}_2 & \cdots & \mathbf{u}_m \\ \tau_1 & \tau_2 & \cdots & \tau_m \end{pmatrix}^T \begin{pmatrix} \mathbf{u}_1 & \mathbf{u}_2 & \cdots & \mathbf{u}_m \\ \tau_1 & \tau_2 & \cdots & \tau_m \end{pmatrix} \quad (19) $$

B is an $m \times m$ symmetric matrix where each element is denoted as $b_{ij}$. We have:

$$ \begin{aligned} \mathrm{Tr}(\mathbf{B}) &= \sum_{i=1}^{m} b_{ii} \\ &= \sum_{i=1}^{m} \lambda_i \end{aligned} \quad (20) $$

where $\lambda_i$ is the ith eigenvalue of matrix B.

From Equations (19) and (20), we have:

$$ \begin{aligned} \sum_{i=1}^{m} \lambda_i &= \sum_{i=1}^{m} \mathbf{u}_i^T \mathbf{u}_i + \tau_i^T \tau_i \\ &= \sum_{i=1}^{m} \|\mathbf{u}_i\|^2 + \|\tau_i\|^2 \end{aligned} \quad (21) $$

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, $\sigma_{max} = \sigma_{min}$. Note that the matrix A is the $n \times 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 \Rightarrow \sigma_{max} = \sigma_{min}$, we have $\lambda_i = \lambda_{max} = \lambda_{min} = \lambda (\sigma = \sqrt{\lambda})$. Equation (21) is rewritten as:

$$ \begin{aligned} n\lambda &= m + \sum_{i=1}^{m} \|\tau_i\|^2 \\ \lambda &= \frac{m}{n} + \frac{1}{n} \sum_{i=1}^{m} \|\tau_i\|^2 \end{aligned} \quad (22) $$

Considering the fact that $\|\tau_i\|^2 \le 1$, we have:

$$ \lambda \le 2 \frac{m}{n} \quad (23) $$

Therefore, we have $\lambda_{max} = 2 \frac{m}{n}$ when $\|\tau_i\|^2 = 1$.

In the singular-value decomposition of a matrix, when $cond(A) = 1$, the matrix A can be written as:

$$ \mathbf{A} = \mathbf{U}\mathbf{S}\mathbf{V}^T = \mathbf{U}[\sigma]_{n \times m}\mathbf{V}^T \quad (24) $$

where $\mathbf{U} \in \mathbb{R}^{n \times n}$, $\mathbf{V} \in \mathbb{R}^{m \times m}$ are orthogonal matrices, $\mathbf{S} = [\sigma]_{n \times m} = \begin{pmatrix} \sigma & 0 & \cdots & 0 \\ \vdots & \sigma & \cdots & 0 \\ 0 & \cdots & \sigma & 0 \end{pmatrix} \in \mathbb{R}^{n \times m}$

The pseudo-inverse of matrix A is $A^+$ and can be written as:

$$ \mathbf{A}^+ = \mathbf{V}\mathbf{S}^+\mathbf{U}^T = \mathbf{V}[\frac{1}{\sigma}]_{m \times n}\mathbf{U}^T \quad (25) $$
$$ \text{where } \mathbf{S}^{+} = [\frac{1}{\sigma}]_{m \times n} = \begin{pmatrix} \frac{1}{\sigma} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \frac{1}{\sigma} & 0 \\ 0 & 0 & \frac{1}{\sigma} \\ 0 & \cdots & 0 \end{pmatrix} \in \mathbb{R}^{m \times n} $$

Our objective for the reactive index was to minimize $||\mathbf{A}^+||$. From Equation (25), the reactive index $I_{re} = ||\mathbf{A}^+|| = \frac{1}{\sigma}$, and the minimum value of the reactive index is equivalent to the maximum value of $\sigma$. 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:

$$ \mathbf{A} = \mathbf{U}\mathbf{S}\mathbf{V}^T \\ = \mathbf{U} \begin{pmatrix} \sigma & 0 & \cdots & 0 & 0 & 0 \\ 0 & \sigma & 0 & \cdots & 0 & 0 \\ 0 & 0 & \sigma & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \sigma & 0 & 0 \end{pmatrix} \mathbf{V}^T \quad (26) $$

where $\mathbf{S}(n \times m)$ is like-diagonal and $\sigma = \sqrt{\lambda} = \sqrt{2\frac{m}{n}}$; $\mathbf{U}(n \times n)$ and $\mathbf{V}(m \times m)$ are orthogonal matrices ($\mathbf{U}\mathbf{U}^T = \mathbf{I}, \mathbf{V}\mathbf{V}^T = \mathbf{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:

**Proposition 1.** Let $\mathbf{M}$ be a $p \times q$ matrix ($p \ge q$), $\mathbf{M} \in \mathbb{R}^{p \times q}$. For all $\mathbf{x} \in \mathbb{R}^q$, if $\mathbf{M} = \mathbf{P}\mathbf{\Sigma}\mathbf{Q}^T$,

where $\mathbf{P} \in \mathbb{R}^{p \times p}$, $\mathbf{Q} \in \mathbb{R}^{q \times q}$ are orthogonal matrices, $\mathbf{\Sigma} = \begin{pmatrix} \mu & 0 & \cdots & 0 \\ 0 & \mu & \cdots & 0 \\ 0 & \cdots & \mu & 0 \\ 0 & \cdots & 0 & \mu \\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 \end{pmatrix} \in \mathbb{R}^{p \times q}$, then
$$ ||\mathbf{M}\mathbf{x}|| = ||\mathbf{M}|| ||\mathbf{x}||. $$

**Proof.** We have:

$$ ||\mathbf{M}\mathbf{x}||^2 = (\mathbf{M}\mathbf{x})^T (\mathbf{M}\mathbf{x}) = \mathbf{x}^T \mathbf{M}^T \mathbf{M}\mathbf{x} \quad (27) $$

With $\mathbf{M} = \mathbf{P}\mathbf{\Sigma}\mathbf{Q}^T$:

$$ \begin{aligned} ||\mathbf{M}\mathbf{x}||^2 &= \mathbf{x}^T (\mathbf{P}\mathbf{\Sigma}\mathbf{Q}^T)^T (\mathbf{P}\mathbf{\Sigma}\mathbf{Q}^T) \mathbf{x} \\ &= \mathbf{x}^T \mathbf{Q}\mathbf{\Sigma}^T \mathbf{P}^T \mathbf{P}\mathbf{\Sigma}\mathbf{Q}^T \mathbf{x} \\ &= \mathbf{x}^T \mathbf{Q}\mathbf{\Sigma}^T \mathbf{\Sigma}\mathbf{Q}^T \mathbf{x} \end{aligned} \quad (28) $$
We have:

$$
\begin{aligned}
\Sigma^T \Sigma &= \begin{pmatrix}
\mu & 0 & \cdots & 0 \\
0 & \mu & \cdots & 0 \\
0 & \cdots & \mu & 0 \\
0 & \cdots & 0 & \mu \\
\vdots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & 0
\end{pmatrix}^T \begin{pmatrix}
\mu & 0 & \cdots & 0 \\
0 & \mu & \cdots & 0 \\
0 & \cdots & \mu & 0 \\
0 & \cdots & 0 & \mu \\
\vdots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & 0
\end{pmatrix} \\
&= \begin{pmatrix}
\mu^2 & 0 & \cdots & 0 \\
0 & \mu^2 & \cdots & 0 \\
\vdots & \vdots & \vdots & \vdots \\
0 & \cdots & 0 & \mu^2
\end{pmatrix} = \mu^2 \mathbf{I}
\end{aligned}
\qquad (29)
$$

where **I** is a q × q identity matrix.

Replacing Equation (29) to (28), we have:

$$
\begin{aligned}
\|\mathbf{M}\mathbf{x}\|^2 &= \mathbf{x}^T \mathbf{V} \mu^2 \mathbf{I} \mathbf{V}^T \mathbf{x} \\
&= \mu^2 \mathbf{x}^T \mathbf{x} = \|\mathbf{M}\|^2 \|\mathbf{x}\|^2
\end{aligned}
\qquad (30)
$$

Therefore, $||\mathbf{M}\mathbf{x}|| = ||\mathbf{M}||||\mathbf{x}||$. $\square$

The energetic index is stated as:

$$
I_e = \frac{1}{S} \int_S (w_{ef} \|\mathbf{A}^+(\mathbf{F}_B^d(f))\| + w_{e\Gamma} \|\mathbf{A}^+ \mathbf{F}_B^d(\Gamma)\|) dS \qquad (31)
$$

Choosing $w_{ef} = w_{e\Gamma} = 1$ (because the desired force vectors, $\mathbf{F}_B^d(f)$, $\mathbf{F}_B^d(\tau)$, are units), we have:

$$
I_e = \frac{1}{S} \int_S (\|\mathbf{A}^+ \mathbf{F}_B^d(f)\| + \|\mathbf{A}^+ \mathbf{F}_B^d(\Gamma)\|) dS \quad (32)
$$

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 pseudo-
inverse matrix **A**+(m × n, m > n) has the following structure:

$$
\mathbf{A}^{+} = \mathbf{V}\mathbf{S}^{+}\mathbf{U}^{T} = \mathbf{V} \begin{pmatrix} \frac{1}{\sigma} & 0 & \cdots & 0 \\ 0 & \frac{1}{\sigma} & \cdots & 0 \\ 0 & \cdots & \frac{1}{\sigma} & 0 \\ 0 & \cdots & 0 & \frac{1}{\sigma} \\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 \end{pmatrix} \mathbf{U}^{T} \quad (33)
$$

where **V**, **U** are orthogonal matrices.

It is clear that matrix **A**<sup>+</sup> satisfies the condition of Proposition **1**. Applying this proposition, we have: ||**A**<sup>+</sup>**F**<sub>B</sub><sup>d</sup>(f)|| = ||**A**<sup>+</sup>||**F**<sub>B</sub><sup>d</sup>(f)|| and ||**A**<sup>+</sup>**F**<sub>B</sub><sup>d</sup>(Γ)|| = ||**A**<sup>+</sup>||**F**<sub>B</sub><sup>d</sup>(Γ)||. Therefore, Equation (32) becomes:

$$
\begin{align*}
I_e &= \frac{1}{S} \int_S (\|\mathbf{A}^+\| \|\mathbf{F}_B^d(f)\| + \|\mathbf{A}^+\| \|\mathbf{F}_B^d(\Gamma)\|) dS \\
&= \frac{1}{S} \|\mathbf{A}^+\| \int_S (\|\mathbf{F}_B^d(f)\| + \|\mathbf{F}_B^d(\Gamma)\|) dS \\
&= 2 \|\mathbf{A}^+\| \tag{34}
\end{align*}
$$
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\|\mathbf{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_{wf}$ 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 \in \mathbb{R}^m$, $\|AF_m\| \le \|\mathbf{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:

<figure><img src="image_13.png" /><figcaption>Figure 10. Upper-bound of the resulting force space.</figcaption></figure>

$$ 
\begin{aligned}
I_{wf} &\le Vol(\mathbf{B}(R1)) \\
I_{wT} &\le Vol(\mathbf{B}(R2))
\end{aligned}
\qquad (35) $$

where $\mathbf{B}(R1)$ and $\mathbf{B}(R2)$ are Euclidean balls of radius $R1 = \|\mathbf{A}(\mathbf{1}: \mathbf{3}, : )\|\|F_m\| = \|\mathbf{A}_1\|\|F_m\|$ and of radius $R2 = \|\mathbf{A}(\mathbf{4}: \mathbf{6}, : )\|\|F_m\| = \|\mathbf{A}_2\|\|F_m\|$, respectively; $\mathbf{A}(\mathbf{1}: \mathbf{3}, : )$ is composed of the first three rows of $\mathbf{A}$, and $\mathbf{A}(\mathbf{4}: \mathbf{6}, : )$ is composed of the last three rows of $\mathbf{A}$.

The volume of a Euclidean ball of radius $R$ in $n$-dimensional Euclidean space is [21]:

$$ V_n(R) = \begin{cases} \frac{\pi^k}{k!} R^{2k}, & \text{if } n = 2k \\ \frac{2^{k+1} \pi^k}{(2k+1)!!} R^{2k+1}, & \text{if } n = 2k+1. \end{cases} \qquad (36) $$

where $(2k+1)!! = 1.3.5, \dots, (2k-1)(2k+1)$.

**Proposition 2.** If the configuration matrix **A** has the form of (26), then $\mathrm{cond}(\mathbf{A}_1) = \mathrm{cond}(\mathbf{A}_2) = 1$ and $\|\mathbf{A}_1\| = \|\mathbf{A}_2\| = \sigma$.

**Proof.** We have:

$$ 
\begin{aligned}
\mathbf{A}\mathbf{A}^T &= (\mathbf{U}\mathbf{S}\mathbf{V}^T)(\mathbf{U}\mathbf{S}\mathbf{V}^T)^T = \mathbf{U}\mathbf{S}\mathbf{V}^T\mathbf{V}\mathbf{S}^T\mathbf{U}^T \\
&= \mathbf{U}\mathbf{S}\mathbf{S}^T\mathbf{U}^T = \sigma^2\mathbf{I}
\end{aligned}
\qquad (37) $$

On the other hand:
$$ 
\begin{aligned}
\mathbf{A}\mathbf{A}^T &= \begin{pmatrix} \mathbf{A}_1 \\ \mathbf{A}_2 \end{pmatrix} \begin{pmatrix} \mathbf{A}_1 \\ \mathbf{A}_2 \end{pmatrix}^T = \begin{pmatrix} \mathbf{A}_1 \\ \mathbf{A}_2 \end{pmatrix} \begin{pmatrix} \mathbf{A}_1^T & \mathbf{A}_2^T \end{pmatrix} \\
&= \begin{pmatrix} \mathbf{A}_1 \mathbf{A}_1^T & \mathbf{0} \\ \mathbf{0} & \mathbf{A}_2 \mathbf{A}_2^T \end{pmatrix}
\end{aligned} 
\quad (38) $$

From (37) and (38), we have:

$$ 
\begin{aligned}
\mathbf{A}_1 \mathbf{A}_1^T &= \sigma^2 \mathbf{I}_1 \\
\mathbf{A}_2 \mathbf{A}_2^T &= \sigma^2 \mathbf{I}_2
\end{aligned} 
\quad (39) $$

where $\mathbf{I}_1$ and $\mathbf{I}_2$ are partitioned matrices of matrix $\mathbf{I}$.

From (39) and the uniqueness of singular-value decomposition [22], it is obvious to see that the structures of $\mathbf{A}_1$ and $\mathbf{A}_2$ are the same as (26) with different dimensions. Therefore, $\mathrm{cond}(\mathbf{A}_1) = \mathrm{cond}(\mathbf{A}_1) = 1$ and $\|\mathbf{A}_1\| = \|\mathbf{A}_2\| = \sigma$. $\square$

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$.

## 4.2. 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 goal-attainment 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:

$$ 
\begin{array}{ll}
\min_{\mathbf{A},\gamma} & \gamma \\
\text{s.t} & \mathbf{A} \in \bar{\mathbf{A}} \\
& \mathbf{V}(\mathbf{A}) - \mathbf{w}\gamma \leq \mathbf{V}_{\text{goal}}
\end{array} 
\quad (40) $$

where $\bar{\mathbf{A}} = \mathbf{A} \setminus I_{ro}$, i.e., the feasible set of configuration matrices, $\mathbf{A}$, without robustness index $I_{ro}$, $\gamma$ is a slack vector variable, and $\mathbf{V}_{\text{goal}} = [I_m^d \quad I_e^d \quad \frac{1}{I_{re}^d} \quad I_{re}^d]$ is the desired objective vector, $\mathbf{w}$ is an attainment vector, which can be chosen. The goal-attainment method with two objective functions is illustrated in Figure 11. By altering $\mathbf{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 goal-attainment 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><img src="image_14.png" /><figcaption>**Figure 11.** Goal-attainment method with two objective functions.</figcaption></figure>

## 5. 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_m^d = 1$, $I_e^d = 1.2248$, $I_{wF}^d = 597.7$, $I_{wT}^d = 597.7$, $I_{re}^d = 0.6124$ ($\sigma^{max} = \sqrt{2\frac{m}{n}} = 1.6330$; see Table 2 for more details).

<table><caption>**Table 2.** Desired values of the indices.</caption><thead><tr><th>Index</th><th>Optimal Formula and Condition</th><th>Desired Value</th></tr></thead><tbody><tr><td>$I_m^d$</td><td>$\sigma_{max} = \sigma_{min}$</td><td>1</td></tr><tr><td>$I_e^d$</td><td>$2 ||\mathbf{A}^+||$</td><td>1.2248</td></tr><tr><td>$\frac{1}{I_w^d}$</td><td>see Equations (35) and (36) and $\frac{1}{I_w^d} = \frac{1}{I_{wF}^d} + \frac{1}{I_{wT}^d}$</td><td>0.0033</td></tr><tr><td>$I_{re}^d$</td><td>$\frac{1}{\sigma^{max}}$</td><td>0.6124</td></tr></tbody></table>

### 5.1. 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.

#### 5.1.1. Phase 1

The optimization toolbox was used to solve the problem (40) with the desired goal vector, and the constraints were $\mathbf{V}_{goal} = [I_m^d \ I_e^d \ \frac{1}{I_w^d} \ I_{re}^d] = [1 \ 1.2248 \ 0.0033 \ 0.6124]^T$, the constraint set $\bar{\mathbf{A}} = \{\mathbf{A} \in \mathbb{R}^{6 \times 8} / ||\mathbf{u}_i|| = 1, ||\boldsymbol{\tau}_i|| \le 1, \boldsymbol{\tau}_i^T \mathbf{u}_i = 0\}$, and the attainment vector $\mathbf{w} = [0 \ 0 \ 0 \ 0.0036]^T$. The attainment vector allows setting the overachievement or underachievement of the individual goals. At the moment, there is no general method
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.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Positions and directions of the thrusters (general case) (XI-axis = **u**-axis; YI-axis = **v**-axis; ZI-axis = **w**-axis).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Attainable spaces in the general case (Fx-axis = **u**-axis; Fy-axis = **v**-axis; Fz-axis = **w**-axis). (a) Force space. (b) Torque space.</figcaption></figure>
<table><caption>**Table 3.** Configuration matrix in the general case.</caption><thead><tr><th></th><th colspan="8">Configuration Matrix</th><th>Optimal Value</th><th>Attainment Factor</th></tr></thead><tbody><tr><td>A =</td><td>(-0.8891<br>-0.0985<br>0.4471<br>-0.4308<br>0.5107<br>-0.7441)</td><td>(-0.3645<br>-0.3036<br>0.8803<br>0.4701<br>0.7561<br>0.4554)</td><td>(0.5438<br>-0.5911<br>0.5957<br>-0.8386<br>-0.4103<br>0.3583)</td><td>(0.9879<br>-0.0608<br>0.1429<br>0.0379<br>0.9868<br>0.1577)</td><td>(0.3134<br>-0.9493<br>0.0260<br>-0.1336<br>-0.0712<br>-0.9885)</td><td>(0.0148<br>0.0515<br>0.9986<br>0.5628<br>-0.8259<br>0.0342)</td><td>(0.0495<br>0.8919<br>0.4495<br>-0.9972<br>0.0690<br>-0.0272)</td><td>(0.6090<br>0.7158<br>0.3417<br>0.4758<br>0.0149<br>-0.8794)</td><td>Fval =</td><td>(1.0000<br>1.2200<br>0.0050<br>0.6124)</td><td>0.3896</td></tr></tbody></table>

### 5.1.2. Phase 2

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.

### 5.2. 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.

#### 5.2.1. 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.

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Robot design with the directions of the thrusters (given position case) (XI-axis = **u**-axis; YI-axis = **v**-axis; ZI-axis = **w**-axis).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Attainable spaces in the given position case (Fx-axis = **u**-axis; Fy-axis = **v**-axis; Fz-axis = **w**-axis). (a) Force space. (b) Torque space.</figcaption></figure>

<table><caption>**Table 4.** Configuration matrix in the given position case.</caption><thead><tr><th></th><th colspan="8">Configuration Matrix</th><th>Optimal Value</th><th>Attainment Factor</th></tr></thead><tbody><tr><td>A =</td><td>0.0836</td><td>0.6616</td><td>-0.8122</td><td>0.4785</td><td>-0.6616</td><td>-0.0836</td><td>-0.4785</td><td>-0.8122</td><td rowspan="5">Fval = <br> (1.0000) <br> 1.2200 <br> 0.0050 <br> (0.6124)</td><td rowspan="5">0.3868</td></tr><tr><td></td><td>0.7452</td><td>0.7452</td><td>-0.3337</td><td>0.3337</td><td>0.7452</td><td>0.7452</td><td>0.3337</td><td>-0.3337</td></tr><tr><td></td><td>0.6616</td><td>-0.0836</td><td>-0.4785</td><td>-0.8122</td><td>0.0836</td><td>-0.6616</td><td>0.8122</td><td>-0.4785</td></tr><tr><td></td><td>-0.8122</td><td>0.4785</td><td>-0.0836</td><td>-0.6616</td><td>-0.4785</td><td>0.8122</td><td>0.6616</td><td>-0.0836</td></tr><tr><td></td><td>-0.3337</td><td>-0.3337</td><td>0.7452</td><td>0.7452</td><td>-0.3337</td><td>-0.3337</td><td>0.7452</td><td>-0.7452</td></tr><tr><td></td><td>0.4785</td><td>0.8122</td><td>0.6616</td><td>-0.0836</td><td>-0.8122</td><td>-0.4785</td><td>0.0836</td><td>0.6616</td><td></td><td></td></tr></tbody></table>

### 5.2.2. Phase 2

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.

### 5.3. 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).

$$A_1 = \begin{pmatrix} 0 & 1 & 0 & 0 & 0 & -1 & 0 \\ 1 & 0 & 0 & -1 & 1 & 0 & 0 \\ 0 & 0 & -1 & 0 & 0 & 1 & 0 \\ 0.27 & 0 & -0.27 & 0.27 & 0.27 & 0 & 0.27 \\ 0 & -0.27 & 0.27 & 0 & 0.27 & -0.27 & -0.27 \\ 0.27 & -0.27 & 0 & 0.27 & 0.27 & 0 & 0 \end{pmatrix} \quad (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).

$$A_2 = \begin{pmatrix} 0.6616 & -0.8122 & 0.4785 & 0.0836 & -0.0836 & -0.4785 & -0.8122 & -0.6616 \\ 0.7452 & 0.3337 & 0.3337 & 0.7452 & 0.7452 & 0.3337 & -0.3337 & 0.7452 \\ -0.0836 & -0.4785 & -0.8122 & 0.6616 & -0.6616 & 0.8122 & -0.4785 & 0.0836 \\ 0.1608 & 0.0111 & -0.2459 & -0.3708 & 0.3642 & 0.2015 & 0.0011 & -0.1658 \\ -0.0989 & 0.3556 & 0.3633 & -0.0989 & -0.1056 & 0.3508 & -0.3456 & -0.1056 \\ 0.3906 & 0.2292 & 0.0044 & 0.1583 & -0.1649 & -0.0254 & 0.2392 & -0.3708 \end{pmatrix} \quad (42)$$
T200 Thruster: Thrust vs. PWM Input to ESC

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Thruster characteristics (BlueRobotics) [24].</figcaption></figure>

Note that the configuration matrices $\mathbf{A}_1$ and $\mathbf{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 $\mathbf{C}^1$ and $\mathbf{C}^2$ are illustrated in Figure 17a,b. It is obvious that the $\mathbf{C}^2$ configuration is more isotropic than the $\mathbf{C}^1$ configuration. However, for some specific points of the attainable force and torque spaces, the $\mathbf{C}^1$ configuration is better than the $\mathbf{C}^2$ configuration.

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Attainable spaces for different configurations (X-axis = **u**-axis; Y-axis = **v**-axis; Z-axis = **w**-axis). (a) $\mathbf{C}^1$(blue), $\mathbf{C}^2$(red). (b) $\mathbf{C}^1$(blue), $\mathbf{C}^2$(red).</figcaption></figure>

Thanks to the properties of matrices $\mathbf{A}_1$ and $\mathbf{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 $\mathbf{C}^2$ configuration are better than $\mathbf{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><caption>**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)).</caption><thead><tr><th>No.</th><th>Indices</th><th>C<sup>1</sup></th><th>C<sup>2</sup></th></tr></thead><tbody><tr><td>1</td><td>$I_m$</td><td>7.12</td><td>2.559</td></tr><tr><td>2</td><td>$I_e$</td><td>3.32</td><td>2.09</td></tr><tr><td>3</td><td>$I_w$</td><td>6,511,536.45</td><td>10,919,428.13</td></tr><tr><td>4</td><td>$I_{re}$</td><td>4.05</td><td>1.56</td></tr><tr><td>5</td><td>$I_{ro}$</td><td>0</td><td>2</td></tr></tbody></table>

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.

<figure><img src="image_21.png" /><figcaption>**Figure 18.** The simulation of the cube rotation about the **u-axis** for $C^1$ and $C^2$ (X-axis = **u-axis**). (a) Applied torque about the **u-axis**. (b) PWM inputs of $C^1$. (c) PWM inputs of $C^2$.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 19.** The simulation of the cube rotation about the **v-axis** for **C**¹ and **C**² (**Y-axis** = **v-axis**). (a) Applied torque about the **v-axis**. (b) PWM inputs of **C**¹. (c) PWM inputs of **C**².</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** The simulation of the cube rotation about the **Z-axis** for **C**¹ and **C**² (**Z-axis** = **w-axis**). (a) Applied torque about the **w-axis**. (b) PWM inputs of **C**¹. (c) PWM inputs of **C**².</figcaption></figure>
**6. 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.

<figure><img src="image_24.png" /><figcaption>Figure 21. $C^1$ of the cube robot.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 22. $C^2$ of the cube robot.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure 23.** Cube robot in the water https://www.youtube.com/watch?v=RKiWUOxDKdw (accessed on 18 October 2019)</figcaption></figure>

## 6.1. 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).

<figure><img src="image_27.png" /><figcaption>**Figure 24.** The cube rotates about the **u-axis** for $C^1$ and $C^2$ (X-axis = **u-axis**).</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Figure 25. The cube rotates about the **v**-axis for **C**¹ and **C**² (Y-axis = **v**-axis).</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure 26. The cube rotates about the Z-axis for **C**¹ and **C**² (Z-axis = **w**-axis).</figcaption></figure>

## 6.2. Energetic Validation

In this section, we verify the energy consumption during these experiments for the two configurations. An energy-like criterion is proposed:
$$ \mathbf{E} = \sum_{i=1}^{m} \int_{t=0}^{T} |PWM^{i}(t) - 1500| dt \qquad (43) $$

where *m* is the number of thrusters, *T* is the time of the experiment, and PWM<sup>*i*</sup>(*t*) is the PWM inputs of the *i*th 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**<sup>2</sup> was lower than that of **C**<sup>1</sup>. For **v**-axis and **w**-axis rotations, the duration of the experiments of **C**<sup>2</sup> was longer than that of **C**<sup>1</sup>, and the energy consumption, therefore, was higher.

<table><caption>**Table 6.** Energy consumption of the two configurations.</caption><thead><tr><th>No.</th><th>Rotation</th><th>E<sub>C1</sub></th><th>E<sub>C2</sub></th></tr></thead><tbody><tr><td>1</td><td><b>u</b></td><td>7.2303 × 10<sup>4</sup></td><td>6.9603 × 10<sup>4</sup></td></tr><tr><td>2</td><td><b>v</b></td><td>7.5480 × 10<sup>4</sup></td><td>1.0590 × 10<sup>5</sup></td></tr><tr><td>3</td><td><b>w</b></td><td>3.1637 × 10<sup>4</sup></td><td>7.4350 × 10<sup>4</sup></td></tr></tbody></table>

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**<sup>2</sup> was lower than that of **C**<sup>1</sup>. However, for the **w**-axis, the energy value of **C**<sup>2</sup> was higher. This happened because the robot dived deeper in **C**<sup>2</sup> in the experiment of the **w**-axis rotation, and the robot had to deliver more power to maintain a greater constant depth.

<table><caption>**Table 7.** Energy consumption of the two configurations with the same time duration.</caption><thead><tr><th>No.</th><th>Rotation</th><th>E<sub>C1</sub></th><th>E<sub>C2</sub></th></tr></thead><tbody><tr><td>1</td><td><b>v</b></td><td>7.5480 × 10<sup>4</sup></td><td>7.2715 × 10<sup>4</sup></td></tr><tr><td>2</td><td><b>w</b></td><td>3.1637 × 10<sup>4</sup></td><td>3.3312 × 10<sup>4</sup></td></tr></tbody></table>

### 6.3. Robustness and Reactive Validation

This section validates the robustness and reactivity of the optimal configuration (C<sup>2</sup>) in comparison to the normal one (C<sup>1</sup>). For robustness, the robot performed a mission, and one or two thrusters were turned off. For the normal configuration C<sup>1</sup>, the mission would fail, and for the optimal configuration C<sup>2</sup>, the mission would be guaranteed. Specifically, for the robustness index, we carried out the following experiments:

1. 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).

<figure><img src="image_30.png" /><figcaption>**Figure 27.** Depth control for $C^1$ and $C^2$ with one and two motors stopped. (a) Depth control of two configurations with one motor stopped. (b) Depth control of two configurations with two motors stopped.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Figure 28.** Depth control for $C^1$ and $C^2$ with three motors stopped.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure 29. PWM evaluation for $C^1$ and $C^2$ with 3 motors stopped. (a) PWM of $C^1$. (b) PWM of $C^2$.</figcaption></figure>

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).

<figure><img src="image_33.png" /><figcaption>Figure 30. Angular velocity evaluation for $C^1$ and $C^2$: diving, rotating about the u-axis, and rotating about the diagonal-axis ($Wx = p$; $Wy = q$; $Wz = r$). (a) Angular velocities of $C^1$. (b) Angular velocities of $C^2$.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 31. Angular velocity evaluation for $C^1$ and $C^2$: diving, rotating about the u-axis, and rotating about the v-axis ($Wx = p$; $Wy = q$). (a) Angular velocities of $C^1$. (b) Angular velocities of $C^2$.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>**Figure 32.** Angular velocity evaluation for **C**¹ and **C**²: diving, rotating about the **u**-axis, and rotating about the **v**-axis (**W**x = **p**; **W**y = **q**).</figcaption></figure>

## 7. 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.

**Author Contributions:** Conceptualization, T.D., L.L., and R.Z.; methodology, T.D., L.L., and R.Z.; software, T.D., B.R., and P.L.; validation, T.D., L.L., R.Z., B.R., and P.L.; writing—original draft preparation, T.D.; writing—review and editing, L.L.; supervision, L.L. and R.Z.; project administration, L.L.; funding acquisition, L.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project was supported by the LabEx NUMEV (ANR-10-LABX-0020) within the I-SITE MUSE (ANR-16-IDEX-0006) and the Region Occitanie (french FEDER funds).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**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.
## Abbreviations

The following abbreviations are used in this manuscript:

AS Actuation System

NGC Navigation-Guidance-Control

DoFs Degrees of Freedom

AUVs Autonomous Underwater Vehicles

CM Center of Mass

PWM Pulse-Width Modulation

## Appendix A

**Theorem A1.** The image of the unit hypersphere under any $n \times m$ matrix is a hyperellipsoid.

**Proof.** Let $\mathbf{A}$ be an $n \times m$ matrix with rank $r$. Let $\mathbf{A} = \mathbf{USV}^T$ be a singular-value decomposition of $\mathbf{A}$. The left and right singular vectors of $\mathbf{A}$ are denoted as $\mathbf{u}_1, \mathbf{u}_2, \dots, \mathbf{u}_n$ and $\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_m$, respectively. Since $\text{rank}(\mathbf{A}) = r$, the singular values of $\mathbf{A}$ have the properties: $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_r > 0$ and $\sigma_{r+1} = \sigma_{r+2} = \dots = \sigma_m = 0$.

Let $\mathbf{x} = \begin{pmatrix} x_1 \\ \vdots \\ x_m \end{pmatrix}$ be an unit vector in $\mathbb{R}^m$. Because $\mathbf{V}$ is an orthogonal matrix and $\mathbf{V}^T$ is also, we have that $\mathbf{V}^T \mathbf{x}$ is an unit vector (it is easy to see that $\|\mathbf{V}^T \mathbf{x}\| = \|\mathbf{x}\|$). Therefore, $(\mathbf{v}_1^T \mathbf{x})^2 + (\mathbf{v}_2^T \mathbf{x})^2 + \dots + (\mathbf{v}_m^T \mathbf{x})^2 = 1$.

On the other hand, we have $\mathbf{A} = \sigma_1 \mathbf{u}_1 \mathbf{v}_1^T + \sigma_2 \mathbf{u}_2 \mathbf{v}_2^T + \dots + \sigma_r \mathbf{u}_r \mathbf{v}_r^T$. Therefore:

$$ \mathbf{A} \mathbf{x} = \sigma_1 \mathbf{u}_1 \mathbf{v}_1^T \mathbf{x} + \sigma_2 \mathbf{u}_2 \mathbf{v}_2^T \mathbf{x} + \dots + \sigma_r \mathbf{u}_r \mathbf{v}_r^T \mathbf{x} $$

$$ \begin{aligned} \mathbf{A} \mathbf{x} &= \sigma_1 \mathbf{u}_1 \mathbf{v}_1^T \mathbf{x} + \sigma_2 \mathbf{u}_2 \mathbf{v}_2^T \mathbf{x} + \dots + \sigma_r \mathbf{u}_r \mathbf{v}_r^T \mathbf{x} \\ &= (\sigma_1 \mathbf{v}_1^T \mathbf{x}) \mathbf{u}_1 + (\sigma_2 \mathbf{v}_2^T \mathbf{x}) \mathbf{u}_2 + \dots + (\sigma_r \mathbf{v}_r^T \mathbf{x}) \mathbf{u}_r \\ &= \mathbf{y}_1 \mathbf{u}_1 + \mathbf{y}_2 \mathbf{u}_2 + \dots + \mathbf{y}_r \mathbf{u}_r \\ &= \mathbf{U} \mathbf{y} \end{aligned} \quad (A1) $$

where $\mathbf{y}_i$ denotes the $\sigma_i \mathbf{v}_i^T \mathbf{x}$ and $\mathbf{y} = \begin{pmatrix} y_1 \\ \vdots \\ y_r \end{pmatrix}$.

From (A1), we have: $\|\mathbf{A}\mathbf{x}\| = \|\mathbf{U}\mathbf{y}\| = \|\mathbf{y}\|$ (since $\mathbf{U}$ is an orthogonal matrix). Moreover, $\mathbf{y}$ has the following property:

$$ \begin{aligned} & \left(\frac{y_1}{\sigma_1}\right)^2 + \left(\frac{y_2}{\sigma_2}\right)^2 + \dots + \left(\frac{y_r}{\sigma_r}\right)^2 = \\ &= (\mathbf{v}_1^T \mathbf{x})^2 + (\mathbf{v}_2^T \mathbf{x})^2 + \dots + (\mathbf{v}_r^T \mathbf{x})^2 \le 1 \end{aligned} \quad (A2) $$

Specifically:

1. If $r = m$ (of course, we must have $m \le n$), the equality in Equation (A2) holds, and the image of the unit hypersphere forms the surface of a hyperellipsoid;

2. If $r < m$, the image of the unit hypersphere corresponds to a solid hyperellipsoid.

This completes the proof. □

## References

1. Lapierre, L. Robust diving control of an AUV. *Ocean Eng.* **2009**, *36*, 92–104. [CrossRef]

2. Lapierre, L.; Jouvencel, B. Robust nonlinear path-following control of an AUV. *IEEE J. Ocean. Eng.* **2008**, *33*, 89–102. [CrossRef]

3. Levine, W.S. *The Control Systems Handbook: Control System Applications*; CRC Press: Boca Raton, FL, USA, 2010.

4. Johansen, T.A.; Fossen, T.I. Control allocation—A survey. *Automatica* **2013**, *49*, 1087–1103. [CrossRef]

5. Ropars, B.; Lapierre, L.; Lasbourygues, A.; Andreu, D.; Zapata, R. Redundant actuation system of an underwater vehicle. *Ocean Eng.* **2018**, *151*, 276–289. [CrossRef]
6. Nakamura, Y.; Hanafusa, H.; Yoshikawa, T. Task-priority based redundancy control of robot manipulators. *Int. J. Robot. Res.* **1987**, *6*, 3–15. [CrossRef]

7. Adorno, B.V.; Fraisse, P.; Druon, S. Dual position control strategies using the cooperative dual task-space framework. In *Proceedings of the IROS’10: International Conference on Intelligent Robots and Systems, Taipei, Taiwan, 18–22 October 2010*; pp. 3955–3960.

8. Agravante, D.J.; Sherikov, A.; Wieber, P.B.; Cherubini, A.; Kheddar, A. Walking pattern generators designed for physical collaboration. In *Proceedings of the ICRA 2016: 2016 IEEE International Conference on Robotics and Automation (ICRA)*, Stockholm, Sweden, 16–21 May 2016; pp. 1573–1578.

9. Pham, T.H.; Caron, S.; Kheddar, A. Multicontact Interaction Force Sensing From Whole-Body Motion Capture. *IEEE Trans. Ind. Inform.* **2018**, *14*, 2343–2352. [CrossRef]

10. Yoshikawa, T. Dynamic manipulability of robot manipulators. In *Proceedings of the 1985 IEEE International Conference on Robotics and Automation*, St. Louis, MO, USA, 25–28 March 1985; Volume 2, pp. 1033–1038. [CrossRef]

11. Yoshikawa, T. *Foundations of Robotics: Analysis and Control*; MIT Press: Cambridge, MA, USA, 1990.

12. Kumar, A.; Waldron, K. The workspaces of a mechanical manipulator. *J. Mech. Des.* **1981**, *103*, 665–672. [CrossRef]

13. Paden, B.; Sastry, S. Optimal kinematic design of 6R manipulators. *Int. J. Robot. Res.* **1988**, *7*, 43–61. [CrossRef]

14. Park, F.C.; Brockett, R.W. Kinematic dexterity of robotic mechanisms. *Int. J. Robot. Res.* **1994**, *13*, 1–15. [CrossRef]

15. Pierrot, F.; Benoit, M.; Dauchez, P. Optimal thruster configuration for omni-directional underwater vehicles. SamoS: A Pythagorean solution. In *Proceedings of the OCEANS ’98 Conference Proceedings*, Nice, France, 28 September–1 October 1998; Volume 2, pp. 655–659. [CrossRef]

16. Kharrat, H. Optimization of Thruster Configuration for Swimming Robots. Master’s Thesis, Rice University, Houston, TX, USA, 2015.

17. Stephan, J.; Fichter, W. Fast Exact Redistributed Pseudoinverse Method for Linear Actuation Systems. *IEEE Trans. Control Syst. Technol.* **2017**, *27*, 451–458. [CrossRef]

18. Grechi, S.; Caiti, A. Comparison between Optimal Control Allocation with Mixed Quadratic & Linear Programming Techniques. In *Proceedings of the 10th IFAC Conference on Control Applications in Marine Systems* CAMS, Trondheim, Norway, 13–16 September 2016; Volume 49, pp. 147–152. doi: 10.1016/j.ifacol.2016.10.335. [CrossRef]

19. Ropars, B.; Lasbouygues, A.; Lapierre, L.; Andreu, D. Thruster’s dead-zones compensation for the actuation system of an underwater vehicle. In *Proceedings of the Control Conference (ECC)*, Linz, Austria, 15–17 July 2015; pp. 741–746.

20. Yoshikawa, T. Manipulability of robotic mechanisms. *Int. J. Robot. Res.* **1985**, *4*, 3–9. [CrossRef]

21. Olver, F.W.; Lozier, D.W.; Boisvert, R.F.; Clark, C.W. *NIST Handbook of Mathematical Functions* Hardback and CD-ROM; Cambridge University Press: Cambridge, UK, 2010.

22. Trefethen, L.N.; Bau, D., III. *Numerical Linear Algebra*; Siam: Philadelphia, PA, USA, 1997; Volume 50.

23. Gembicki, F.; Haimes, Y. Approach to performance and sensitivity multiobjective optimization: The goal-attainment method. *IEEE Trans. Autom. Control* **1975**, *20*, 769–771. [CrossRef]

24. BlueRobotics. Available online: https://bluerobotics.com/ (accessed on 20 December 2018).