## General Near-Field Synthesis Based on the Optimization of Relevant Figures of Merit
### Direct Layout Optimization Algorithm
#### 2.2.2 Multi-frequency Direct Layout Optimization  
The proposed multi-frequency DLO consists of an iterative procedure. Starting with the initial layout, given by the values contained in vector \(\vec{\xi}\), the NF generated by the reflectarray is calculated at a number of frequencies. From this NF, the initial magnitude and phase FoM are calculated for each orthogonal polarization and frequency. The next step is to calculate the cost function. To this end, three elements are needed: the value of the current FoM (\(\Delta_{k,\Omega_{m}}^{m/p}\)), the specifications that are to be met (\(\bar{\Delta}_{k,\Omega_{m}}^{m/p}\)), and the weight functions to leverage the error of each FoM (\(W\)). Then, the cost function is defined as:  
\[r_{f,k,m}=W_{f,k,m}\Big{(}\bar{\Delta}_{k,\Omega_{m}}^{m/p}-\Delta_{k,\Omega_{ m}}^{m/p}\Big{)}, \tag{7}\]  
where \(f=1,2,\ldots,N_{f}\) is the index for the dependence on the frequency, with \(N_{f}\) the total number of frequencies at which the optimization is carried out; \(k=1,2,\ldots,K\) is the index for the number of FoM (\(K\) in total) defined in region \(\Omega_{m}\); and \(m=1,2,\ldots,M\) is the index for the different regions where FoM may be defined. Please note that it is not necessary to optimize FoM in all regions. Since the LM algorithm is employed, the cost function defined in (7) is also known as the residual [38].  
The next step is to compute the Jacobian matrix. The size of this matrix depends on the total number of DoF taken into account in the current iteration (which may vary from one to another) and the total number of FoM that are optimized. The number of DoF define the number of columns, while the number of FoM define the number of rows. Thus, the size of the matrix is \(4N_{f}\times N_{P}\), where \(N_{P}\leq P=LN\). Any element of the Jacobian matrix may be calculated as:  
\[J(q,s)=\frac{\partial r_{q}}{\partial\vec{\xi}_{s}}, \tag{8}\]  
where \(q=1,2,\ldots,4N_{f}\), and \(s=1,2,\ldots,N_{P}\). In (8), the dependence on indices \(f\), \(k\) and \(m\) has been dropped to alleviate the notation. In addition, although not directly specified, the residual depends on the optimizing variables since \(\vec{\xi}\) defines the electromagnetic response of the reflectarray unit cell that will shape the NF from which the FoM are calculated to obtain the residual. In addition, the partial derivative in (8) may be numerically evaluated by means of finite differences [38; 42], that, for the particular case of gradient-based algorithms for array antenna synthesis, can be efficiently computed using the differential contributions technique [43].  
Once the Jacobian matrix has been computed, the normal equations are formed:  
\[\Big{[}J_{i}^{T}I_{i}+\mu_{i}\operatorname{diag}\Big{(}J_{i}^{T}I_{i}\Big{)} \Big{]}\delta_{i}=-J_{i}^{T}r_{i}, \tag{9}\]  
which can be compactly written as:  
\[A_{i}\delta_{i}=b_{i}, \tag{10}\]  
where subindex \(i\) denotes the current iteration (see Figure 2), \(J^{T}\) is the transpose of the Jacobian matrix, \(\mu\) is a real positive number, \(\operatorname{diag}(\cdot)\) is the diagonal matrix, and \(\delta_{i}\) is the updating vector that satisfies the equality. Since matrix \(A_{i}\) is at least semi-positive definite, an efficient Cholesky factorization can be used to solve (10) for \(\delta_{i}\).  
After the normal equations are solved, the solution is updated as:  
\[\vec{\xi}_{i+1}=\vec{\xi}_{i}+\delta_{i}. \tag{11}\]  
With the updated layout \(\vec{\xi}_{i+1}\), the NF and FoM are computed again and checked against the specifications. If the requirements are met, the algorithm stops and \(\vec{\xi}_{i+1}\) is the final layout. Otherwise, the number of DoF in the optimization may be modified, as well as the specifications, and the process is repeated again until convergence or a given maximum number of iterations is reached.