Compressive synthetic aperture sonar imaging with distributed optimization.

Synthetic aperture sonar (SAS) provides high-resolution acoustic imaging by processing coherently the backscattered acoustic signal recorded over consecutive pings. Traditionally, object detection and classification tasks rely on high-resolution seafloor mapping achieved with widebeam, broadband SAS systems. However, aspect- or frequency-specific information is crucial for improving the performance of automatic target recognition algorithms. For example, low frequencies can be partly transmitted through objects or penetrate the seafloor providing information about internal structure and buried objects, while multiple views provide information about the object's shape and dimensions. Sub-band and limited-view processing, though, degrades the SAS resolution. In this paper, SAS imaging is formulated as an ℓ1-norm regularized least-squares optimization problem which improves the resolution by promoting a parsimonious representation of the data. The optimization problem is solved in a distributed and computationally efficient way with an algorithm based on the alternating direction method of multipliers. The resulting SAS image is the consensus outcome of collaborative filtering of the data from each ping. The potential of the proposed method for high-resolution, narrowband, and limited-aspect SAS imaging is demonstrated with simulated and experimental data.


I. INTRODUCTION
Synthetic aperture sonar (SAS) combines coherently the backscattered echoes from successive acoustic pulses (pings) for high-resolution seafloor imaging with application in mine countermeasures, underwater archaeology, or inspection of underwater installations. 1,2 Specifically, an active sonar transmits a short pulse to insonify the seafloor and records the backscattered waves repeatedly while it moves along a predefined (usually linear) path to form a synthetic aperture. SAS imaging refers to the inverse problem of reconstructing the seafloor reflectivity by coherently processing the recorded signals from the synthetic aperture. The range resolution in SAS imaging is determined by matched filtering, i.e., cross-correlating the recorded and the transmitted signal, hence it is inversely proportional to the bandwidth of the transmitted pulse. The crossrange resolution achieved with SAS systems is independent of frequency and range, i.e., it is not limited by the physical dimensions of a real aperture which is characterized by a constant angular resolution determined by its beampattern resulting in a range and frequency dependent cross-range resolution. 3,4 Utilizing widebeam systems and broadband pulses at a high ping repetition rate, synthetic aperture processing can achieve a centimetric range and cross-range resolution resulting in optical-like images of the seafloor reflectivity. 5,6 Conventional target recognition tasks rely on highresolution seafloor mapping. 7 However, in applications such as unexploded ordnance remediation or mine countermeasures, aspect-and frequency-specific features are useful for reducing the false alarm rate of automatic target recognition (ATR) algorithms, [8][9][10] as well as for increasing the detection rate in challenging environments. For example, low frequencies can excite acoustic resonances in objects with structural symmetries or detect buried objects due to improved seafloor penetration, [11][12][13][14] while multi-aspect imaging increases the information on the object's shape and dimensions. 15 However, sub-band or sub-view processing reduces the resolution and the signal-to-noise ratio (SNR) in conventional SAS imaging. Sub-band and sub-view feature extraction for classification is based on image denoising methods such as spatial filtering followed by deconvolution, 16 wavenumber domain filtering, 17 or sparse feature selection through wavelet shrinkage. 18 In this paper, we formulate the SAS imaging problem within the compressive sensing (CS) framework, which asserts that underlying sparse signals can be reconstructed from very few measurements with convex optimization. 19 In ocean acoustics, sparse signal reconstruction is shown to improve radically the resolution in direction-of-arrival estimation 20,21 and matched field processing. 22 In synthetic aperture radar compressive imaging has been applied for feature enhancement, 23 autofocusing, 24 and deception jamming suppression. 25 In SAS imaging for ultrasound applications in air, CS has been applied to overcome the ambiguities arising from spatial under-sampling. 26 We develop a compressive SAS imaging method for underwater applications which, to the best of the authors' knowledge, has not been reported in the literature yet. In particular, the main contributions of this study are: • SAS imaging is formulated in the frequency domain, in Sec. IV, as a least-squares optimization problem regularized with an ' 1 -norm penalty which promotes sparse solutions. 27 • Since the backscattered data are naturally distributed over multiple pings, we use the alternating direction method of multipliers 28 (ADMM) to solve the optimization problem in a distributed way. The resulting consensus optimization problem involves minimizing a data fitting term per ping with the constraint that the local solutions agree to a global model for the SAS image, while it allows for decentralized storage and processing of the large amount of data resulting from aperture synthesis. • Guidelines on regularization parameter selection for tuning the optimization problem are provided in Sec. V B. • The high-resolution imaging and denoising capabilities of the proposed method are demonstrated in Sec. V with simulations and validated with experimental data for narrowband and multi-view imaging.
The paper makes heavy use of convex optimization theory. The reader is advised to consult the Appendix for a summary on the basic notions on convex optimization and the mathematical background of the ADMM algorithm.

A. Mathematical notation
Herein, vectors are represented by bold lowercase letters and matrices by bold uppercase letters. The symbols T , H denote the transpose and the Hermitian (i.e., conjugate transpose) operator, respectively, on vectors and matrices. The ' p -norm of a vector x 2 C n is defined as kxk p ¼ ð P n i¼1 jx i j p Þ 1=p . The infinity norm of x is defined as the maximum vector element in absolute value kxk 1 ¼ max i2n jx i j. The Frobenius norm of a matrix X 2 C mÂn is defined as where Trð:Þ denotes the matrix trace.

II. STRIP-MAP SAS GEOMETRY
The SAS geometry in strip-map mode 2 is depicted in Fig.  1. Let the seafloor lie on the xy-plane and be characterized by the complex reflectivity function sðr s Þ; r s ¼ ðx; yÞ. The platform carrying a SAS system, comprising a configuration of transmitters and receivers, moves along a linear path, parallel to the y axis at height h from the seafloor plane. The reconstruction of the seafloor reflectivity is depicted on the slant-plane, i.e., in slant-range-cross-range coordinates. In strip-map mode, the antenna is focused toward broadside, i.e., the central axis of the real-aperture beampattern is perpendicular to the platform path.
The active sonar transmits a short pulse of duration s q , onwards referred to as ping, and records the backscattered echoes repeatedly as the platform moves along the track. The stop-and-hop approximation postulates that the platform is stationary during each ping transmission and reception, before it jumps instantaneously to the next position. 2 Hence, the position of the platform is discretized according to the ping number p as (1) where v p is the constant speed of the platform and s rec is the duration of the recording which defines the ping repetition period. The insonified area per ping is determined by the radiation pattern of the transmitting antenna. The total imaging area is determined by the synthetic aperture length and the recording duration s rec allowing a maximum swath of where c is the speed of sound in water.
Multi-element receiver arrays are employed to allow longer imaging ranges without violating the spatial sampling condition for a moderate platform speed (a few knots). 1 With the phase center approximation 29 (PCA), the resulting single-transmitter/multiple-receivers multistatic configuration is replaced with a virtual array of monostatic elements located at the middle of the distance between each transmitter-receiver pair on the actual array. The resulting PCA virtual array has half the length of the receiver array L.
The resolution of a SAS system is defined in the range and cross-range directions, indicated in Fig. 1 as d x and d y , respectively. The range resolution, obtained with matched filtering, is determined by the bandwidth Df of the transmitted ping, : The cross-range (angular) resolution depends on the apparent synthetic aperture length L SA , For a given transducer size D, the corresponding synthetic aperture length is proportional to the wavelength and the range of the target, L SA % kr 0 =D, resulting in a cross-range resolution which is independent of frequency and range, d y % D=2. 1

III. SAS MODEL
The transmitted ping q(t) is expressed in a general form as where w(t) is a window function producing a slowly-varying envelop, and /ðtÞ is the phase as a function of time t. For instance, a linear frequency modulated (LFM) chirp of duration s q sweeping linearly from a low frequency f l to a high frequency where the signal's instantaneous frequency f is a linear function of time t 2 ½0; s q with slope a.
With the stop-and-hop assumption and the PCA, the backscattered wave from an isotropic point scatterer located at r s within the insonified area I is a replica of the transmitted pulse q(t) delayed by the travel time for the two-way distance between the scatterer and the virtual transceiver, multiplied by the scatterer's complex reflectivity sðr s Þ. Then, with the Born approximation, 30 the signal dðr r ; tÞ recorded at a receiver located at r r is the superposition of the backscattered echoes from all scatterers within the insonified area I, dðr r ; tÞ ¼ ð I q t À 2kr s À r r k 2 c sðr s Þdr s : The frequency response of the recorded signal [Eq. (7)] is obtained through the Fourier transform as where we have included, for completeness, an amplitude factor Gðr r ; r s ; f Þ ¼ G t ðr s ; f ÞG r ðr s ; f Þ=kr s À r r k 2 2 accounting for the shading function of the transmitter array G t ðr s ; f Þ and the receiver array G r ðr s ; f Þ and the attenuation due to the spherical spreading over the two-way wave propagation distance.
The matched filtered frequency response is computed by multiplying the recorded signal [Eq. (8)] by the complex conjugate of the transmitted pulse, IV. SAS IMAGING Discretizing the reflectivity field into N pixels, the data model [Eq. (9)] can be written in a matrix-vector formulation, where d 2 C M is the vector of the matched filtered measurements at frequency f for all M receivers comprising the real aperture at ping p, s 2 C N is the unknown vector of the complex reflectivity values over a two-dimensional (2D) grid of N pixels, and n 2 C M is the additive noise vector. The matrix A 2 C MÂN maps the unknown reflectivity s to the observations d and has as columns the steering vectors, which describe the propagation delay from the sth scatterer to all the M sensors on the real aperture at ping p. In this study, the matrix A is considered known and fixed. Note that we have incorporated the gain factors, kqðf Þk 2 2 and Gðr r ; r s ; f Þ, into the reflectivity vector s as they can be easily accounted for in a calibrated system. SAS imaging refers to the inverse problem of reconstructing the reflectivity field s, given the sensing matrix A and a set of measurements d over a range of frequencies and pings.

A. CBF
An estimateŝ of the reflectivity field can be obtained by spatial filtering the array data d. Conventional (delay-andsum) beamforming (CBF) uses the steering vectors [Eq. (11)] as spatial weights to combine the sensor outputs coherently, enhancing the signal at the look-direction from the ubiquitous noise. In SAS imaging, CBF provides the reflectivity estimate, by combining coherently the sensor outputs over P pings and F frequencies.
To improve readability, we use the indices i and j for the ping p and frequency f dependency, respectively, of the quantities in Eq. (10), i.e., dðp; f Þ d ij and Aðp; f Þ A ij .

B. Sparse reconstruction with the ADMM
In the case that there are only a few strong scatterers in the reflectivity field (K ( N), SAS imaging can be solved as a sparse model fitting problem, where l > 0 is a regularization parameter which controls the relative importance between the quadratic data-fitting term and the sparsity promoting ' 1 -norm regularization term. The ' 1 -norm regularized least-squares problem (13) is known as the least absolute shrinkage and selection operator (lasso), 27 and can be solved with convex optimization methods. 31 In a statistical framework, the lasso is equivalent to maximum a posteriori estimation of a linear model with Gaussian noise model and Laplacian prior on the model parameters. 32 Interior point solvers, implemented for instance in the Python-embedded CVXPY package, 33,34 provide a highly accurate solution to the lasso problem for applications with relatively few model parameters such as in directionof-arrival estimation. 20,35 However, in SAS imaging the single-ping, single-frequency lasso problem (13) is highly underdetermined, i.e., the number of measurements from the M sensors on the real array is much smaller than the number of pixels with unknown reflectivity N. Hence, the mapping A ij is highly coherent, i.e., it is characterized by linearly dependent columns (see Sec. V B). The coherence of A ij makes the solution obtained with interior point solvers, highly inaccurate for l > 0. 19,20 Besides, the formulation (14) does not account for the common sparsity profile among pings and frequencies identified by the location of strong scatterers in the reflectivity scene. To reduce the coherence of the mapping A ij , problem (13) could be reformulated by concatenating the real aperture data vectors for all pings to form the synthetic aperture data vector d SA j 2 C MP with the corresponding sensing matrix A SA j . The resulting solution, is more accurate than Eq. (14), yet it does not account for the common sparsity profile among frequencies while it is computationally impractical due to the increased problem dimensions and requires a centralized collection of the data before processing.
Here, we solve problem (13) with the ADMM which is an efficient convex optimization solver for problems with a very large number of model parameters N; see Appendix 2 for a full derivation of the iterative ADMM solution. Accounting for separable measurements over pings and frequencies, the lasso problem (13) is reformulated as Problem (16) is a consensus problem, since the constraint imposes that the local reflectivity variables s ij related to the data set from each ping and each frequency should be equal to the global variable z. A common reflectivity profile z over pings and frequencies implies isotropic scattering for the considered view angle and frequency band. Isotropic scattering is a typical assumption in SAS, where the backscattered waves from consecutive pings are coherently integrated to achieve high-resolution imaging. Consensus optimization promotes collaborative filtering as the P data sets over the frequency band of interest are collaborating to develop a global model for z. Note that Eq. (16) is a separable, equality constrained, convex optimization problem which can be derived from the generic formulation (A8) with f ðsÞ ¼ where q > 0 is the augmented Lagrangian regularization parameter. A closed-form solution for the s ij and z update steps is obtained with partial differential calculus by solving, respectively, and The overline operator in Eq. (19) denotes averaging over the P pings and F frequencies, v ¼ 1=ðPFÞ P P i¼1 P F j¼1 v ij and the subgradient @ z is a generalization of the partial differential operator for functions that are not differentiable everywhere (Ref. 31, p. 338).
The resulting ADMM iterative algorithm for solving problem (16) is where S : ð:Þ denotes a soft-thresholding operator defined as S j ðaÞ ¼ ða=jajÞmaxðjaj À j; 0Þ for a 2 C; j 2 R. 28 Note that the s ij -update is an ' 2 -norm regularized least-squares solution (ridge regression), so the ADMM is essentially solving the lasso problem by an iterative ridge regression procedure. The sand u-variable update steps in Eq. (20) are independent for each ping i 2 ½1; P and for each frequency j 2 ½1; F. Hence, the ADMM provides the solution to a large global problem by coordinating the solutions of smaller decentralized local problems, i.e., through a broadcast and gather procedure. In the following, the ADMM algorithm (20) is initialized with s 0 ij ¼ A H ij d ij ; z 0 ¼ s 0 and u 0 ij ¼ 0 for i 2 ½1; P and j 2 ½1; F. We defer the discussion on regularization parameter selection until Sec. V. The stopping criterion is based on the ' 2 -norm of the primal and dual residual (see Appendix 3), where kþ1 r ¼ abs ffiffiffi ffi N p þ rel maxðk s kþ1 k 2 ; kz kþ1 k 2 Þ and kþ1 y ¼ abs ffiffiffi ffi N p þ rel kq u kþ1 k 2 are the corresponding convergence tolerances. 28 Herein, the absolute and relative tolerance bounds are abs ¼ rel ¼ 10 À3 . Making use of the matrix inversion lemma identity ðA þ BCDÞ À1 ¼ A À1 À A À1 BðC þDA À1 BÞ À1 DA À1 which holds for matrices A, B, C, D with conformable dimensions when all the inverses exist, 28 the matrix inverse in the s ij -minimization step can be replaced by The identity [Eq. (22)] allows caching the factorization of a smaller matrix which is computationally advantageous in SAS imaging where M ( N.

V. RESULTS
The potential of the proposed sparse reconstruction method over CBF is demonstrated first on sub-band SAS imaging with synthetic data and then on narrowband and limited-aspect SAS imaging with experimental measurements.

A. Simulations
For the simulations, we employ a SAS system, as in Fig.  1, with a uniform linear receiver array with 12 sensors at spacing d ¼ 0.03 m. The transmitted signal is a LFM pulse [Eq. (6)] of duration s q ¼ 10 ms with bandwidth 3-20 kHz and a Tukey (tapered cosine) window w(t). The PCA is used to replace the single-transmitter/multiple-receivers multistatic configuration with a monostatic one. The array displacement between pings is Dy ¼ 0:09 m and the cross-range distance from À2 to 2 m is used for aperture synthesis so that the whole imaging area is within the À3 dB beam width of the transmitter which is 40 . The sound speed is c ¼ 1500 m/s and the imaging grid comprises N ¼ 31 Â 41 ¼ 1271 pixels. With the selected system parameters and assuming that the SAS platform's velocity is 3 knots, the stop-and-hop approximation is valid; see Eq. (51) in Ref. 36.
To describe the SAS system, both in terms of resolution and artifacts from sidelobes due to matched filtering and beamforming, we first calculate the point spread function (psf) which is the response of an imaging system to an isotropic point scatterer. 36 Figure 2 shows the psf of CBF and sparse reconstruction implemented with the ADMM for a point scatterer at r 0 ¼ ð8:5; 0Þ m with the considered SAS system for full-band (3-20 kHz) and sub-band (7-10 kHz) processing. The range resolution of CBF deteriorates for sub-band processing as expected from Eq. (3). The crossrange resolution also deteriorates since the synthetic aperture length remains unaltered for full-band and low frequency sub-band processing; see Eq. (4). Sparse reconstruction with consensus optimization offers high-resolution imaging, free of sidelobes in both cases. Figure 3 shows the SAS reconstruction for sub-band processing when a second isotropic scatterer is added to the reflectivity field with power À20 dB relative to the power of In the latter case, the artifacts in sparse reconstruction which are adjacent to the strong scatterer (at the order of À12 dB), where the spatial coherence of the psf is high, are due to the fact that the ADMM converges to a solution with moderate accuracy; see Appendix 3. Nevertheless, Fig. 4 demonstrates that consensus optimization with the ADMM algorithm provides more accurate sparse reconstruction in SAS imaging than interior point solvers. Specifically, since the lasso problem in Eq. (15) is solved independently for each frequency, the coherence in slant-range is still high (low slant-range resolution) resulting in artifacts along the slant-range. Besides, solving Eq. (15) with interior point solvers requires 8.8 min due to the size of the problem while the ADMM provides the solution in Fig. 4

(d) in less than 40 s (8 iterations).
Consider now a more realistic scenario where the seafloor is modeled as a random field of isotropic scatterers 37 with characteristic length 0.2 m and exponential covariance with zero-mean and standard deviation À15 dB relative to the reflectivity of an embedded strongly scattering disk as depicted in Fig. 5(a). The SNR due to additive noise is 20 dB. Sub-band, low frequency (7-10 kHz) processing degrades the range resolution of conventional SAS imaging [ Fig. 5(b)] while sparse reconstruction provides a more accurate estimate of the strongly scattering disk [ Fig. 5(c)].

B. Regularization parameter selection
It is important to note that the solution of the ADMM consensus optimization [Eq. (20)] is tuned with two regularization parameters: the regularization parameter for the augmented Lagrangian term q and the regularization parameter for the ' 1 -norm constraint on the model features l. Even though the best choice for the values of the regularization parameters is problem-dependent, we provide general guidelines on regularization parameter selection.  To make the analysis generally applicable to SAS imaging, i.e., independent of a particular system or model, we normalize the sensing matrix by its Frobenius norm, A ij =kA ij k F , such that the N eigenvalues of the associated Gram matrix sum up to one P eigðA H ij A ij Þ ¼ 1. Besides, we normalize the data d ij =kd ij k 2 to focus on the phase differences between the recorded data on the physical array.
The parameter q regularizes the least-squares solution for the s-variable update in Eq. (20) in terms of diagonal loading of the Gram matrices A H ij A ij to stabilize its inverse. Strip-map SAS requires receiver arrays with broad beam widths to maximize the effective length of aperture synthesis providing fine cross-range resolution which in turn requires fine imaging grids. The wide beampattern of the (physical) receiver array in combination with the fine spatial sampling of the SAS image results in rank-deficient Gram matrices A H ij A ij . Figure 6 shows that less than 0.1% of the eigenvalues have a significant value. In practice, there is one non-zero eigenvalue approximately equal to 1 and setting q ¼ 1 regularizes the least-squares s-variable update sufficiently. Bear in mind that, even though the Gram matrix for a single ping and frequency has highly coherent columns, the resulting coherence pattern of the total SAS processing is better described by the psf; see Fig. 2.
The parameter l controls the sparsity of the solution, hence it depends on the specific setup. For large values of l, the solution is very sparse but the data fit is poor. In fact for l > l max ¼ k P P i¼1 P F j¼1 A H ij d ij k 1 the solution z of the distributed lasso problem (20) is zero. 28 As l decreases toward zero the solution becomes gradually less sparse improving the data fit. For l ¼ 0 the ADMM solution becomes equal to the CBF solution z ¼ŝ CBF . The best value for l for a specific setup is determined through the lasso path, 21,38 i.e., by solving problem (20) for different values of l in the interval 0 < l l max . Figure 7 shows the ADMM solution for the setup in

C. Convergence
Convergence is another important aspect of the iterative ADMM solution as it determines both the computational efficiency of the algorithm and the accuracy of the optimal solution. Figure 8 shows the evolution of the primal and dual residual norms of the ADMM algorithm for the sparse reconstruction in Fig. 5 (q ¼ 1, l ¼ 0:5l max ). The stopping criterion, as determined by the tolerance bounds [Eq. (21)], is satisfied after eight iterations, while each iteration takes less than 5 s to compute. The fast convergence of the ADMM iterative solution is very useful for practical applications.   Figure 9 depicts the primal and dual variable updates as the iterative algorithm (20) approaches convergence. The consensus problem (16) gradually constrains the local variables s ij , which are initialized with the very low resolution CBF output for a single ping and frequency bin, to match the global variable z. Besides, the global variable z, which is initialized with the total SAS beamforming output, is optimized toward a sparser representation by the ' 1 -norm, while the dual variable is updated with the running sum of the residuals The convergence of the ADMM algorithm is highly determined by the choice of the regularization parameter q; see Sec. 3.4.1 in Ref. 28. The ADMM update Eqs. (17) indicate that large values of q favor primal feasibility, i.e., the fulfillment of the constraints in the optimization problem (16). Hence, a large q reduces the number of iterations needed for the primal residual norm [Eq. (21)] to reach its tolerance bound. The fast convergence of the primal residual norm comes at the expense of a less sparse optimal solution for z which tends to increase the dual residual norm. Contrarily, small values of q induce  opposite trends in the convergence of the primal and dual residual norms as demonstrated in Fig. 10. The convergence rate achieved for compressive SAS imaging with the suggested choice of q ¼ 1 (less than 10 iterations irrespective of SNR for all examples in this study) provides a good balance between computational efficiency and optimization accuracy.

D. Experimental results
The measurements are from the target-oriented railbased experiment (TORHEX'18) 39 collected in July 2018 at the port of La Spezia, Italy, with a SAS system similar to the one used for the simulations. A steel air-filled sphere of 1 m diameter and a concrete cylinder of 2 m length are used to exemplify the potential frequency-and aspect-specific backscattering features, respectively, of man-made objects. The objects are lying on a muddy seafloor at shallow water between 6 and 10 m depth. The broadband 3-20 kHz, fullview [À10 , 10 ] SAS imaging of the two objects is depicted in Fig. 11 as a reference. Figure 12 shows the frequency-dependent, full-view [À10 , 10 ] SAS reconstruction with CBF and compressive SAS imaging for the spherical object. With narrowband, low frequency processing CBF exhibits limited resolution, while compressive SAS imaging indicates clearly the frequencydependent scattering features of the shell, suppressing completely the background reverberation. Notably, at frequencies around f c ¼ 4 kHz the sound waves are mainly scattered from the convex part of the shell, while at frequencies around f c ¼ 7 kHz it is the concave part which scatters the sound.
Similarly, Fig. 13 shows the aspect-dependent SAS reconstruction with CBF and compressive SAS imaging for the cylindrical object. In this case, the processing is broadband (3-20 kHz) while sub-views of the cylinder are achieved by limiting the length of the synthetic aperture. The cylinder is slightly tilted towards negative angles which results in strong specular backscattering by its full length at sub-views [À10 , À3 ] and [À3 , 3 ], while at [3 , 10 ] the main scatterer is the cylinder's base.

VI. CONCLUSION
Frequency-and aspect-specific SAS imaging can provide important information for the classification of man-made objects. However, sub-band or limited-view processing comes at the expense of reduced resolution and increased sidelobe levels in CBF. In this paper, SAS imaging is formulated as a sparse signal reconstruction problem. The resulting optimization problem is solved efficiently in a distributed way with the ADMM. Sparse reconstruction with distributed optimization is shown to improve significantly the resolution and the SNR in SAS processing with limited data indicating the great potential of the method in SAS imaging. Specifically, we demonstrate both with synthetic and measured data that compressive SAS imaging improves the resolution compared to conventional delay-and-sum beamforming and offers a sidelobe-free reconstruction. Future efforts will focus on combining the developed compressive SAS imaging with ATR algorithms.

ACKNOWLEDGMENTS
This work was performed under the Project SAC000905-High Resolution Low Frequency Synthetic Aperture Sonar (HRLFSAS) of the STO-CMRE Programme of Work, funded by the NATO Allied Command Transformation. The authors would like to thank Dr. S. Fioravanti, A. Sapienza, A. Carta, and F. Aglietti for the data collection during the TORHEX'18 trials.

APPENDIX: DISTRIBUTED CONVEX OPTIMIZATION
This appendix summarizes the basic formulations in distributed convex optimization as analytically presented in Refs. 28 and 31. First, the primal and dual formulation of a generic equality-constrained convex optimization problem is presented, followed by a concise description of the ADMM, which is an algorithm for distributed convex optimization particularly well-suited for solving large-scale problems. Basically, this section provides the mathematical framework for the derivation of the consensus optimization problem (16), the iterative ADMM solution [Eq. (17)], and the stopping criterion [Eq. (21)].

Equality-constrained convex optimization
Consider the equality-constrained convex optimization problem, min s f ðsÞ subject to Hs ¼ b; (A1) where s 2 C N is the optimization variable, f : C N ! R is the convex objective (or cost) function, and the constraints are an affine transformation with H 2 C qÂN and b 2 C q . The optimal value p Ã of the optimization problem (A1), achieved at the optimal variable s Ã , is The Lagrangian for problem (A1) is obtained by augmenting the objective function with a weighed sum of the constraints, where m 2 C q is the dual variable of the problem (A1). The corresponding dual function is the minimum value of the Lagrangian (A3) over s, gðmÞ ¼ inf s Lðs; mÞ: (A4) The dual function, which is concave even when problem (A1) is not convex, 31 yields lower bounds on the optimal value p Ã (A2), since gðmÞ ¼ inf s Lðs; mÞ Lðs;mÞ f 0 ðsÞ for every feasible points.  3 ], and (c) and (f) [3 , 10 ] relative to the point of closest approach.