2.1. Split-step Fourier method

In this work, we model the optical network physical layer using the Manakov equation, given by (Millar et al., Reference Millar, Makovejs, Behrens, Hellerbrand, Killey, Bayvel and Savory2010)

(1) $$ {\displaystyle \begin{array}{r}\frac{\partial {E}_X}{\partial z}+\frac{\alpha }{2}{E}_X+\frac{j{\beta}_2}{2}\frac{\partial^2{E}_X}{\partial {t}^2}\\ {}-j\frac{8}{9}\gamma \left({\left|{E}_X\right|}^2+{\left|{E}_Y\right|}^2\right){E}_X=0,\\ {}\frac{\partial {E}_Y}{\partial z}+\frac{\alpha }{2}{E}_Y+\frac{j{\beta}_2}{2}\frac{\partial^2{E}_Y}{\partial {t}^2}\\ {}-j\frac{8}{9}\gamma \left({\left|{E}_X\right|}^2+{\left|{E}_Y\right|}^2\right){E}_Y=0,\end{array}} $$

where $ {\beta}_2 $ is the group velocity dispersion coefficient of the fiber, $ \alpha $ is the fiber attenuation coefficient, $ t $ and $ z $ denote the retarded time and the spatial coordinate along the direction of propagation respectively, $ {E}_X $ and $ {E}_Y $ are the electric field components along the X and Y axes, respectively, and $ \gamma $ is the fiber nonlinearity coefficient, given by

(2) $$ \gamma =\frac{2\pi {n}_2}{\lambda {A}_{\mathrm{eff}}}, $$

where $ {A}_{\mathrm{eff}} $ is the fiber effective mode area, $ {n}_2 $ is the nonlinear index coefficient, and $ \lambda $ is the center operating wavelength. $ {\beta}_2 $ relates to the dispersion coefficient $ D $ via $ D=-\frac{2\pi c}{\lambda^2}{\beta}_2 $ . Writing (1) in terms of its constituent linear and nonlinear parts can be shown to yield an exact solution given by (Millar et al., Reference Millar, Makovejs, Behrens, Hellerbrand, Killey, Bayvel and Savory2010)

(3) $$ \mathbf{E}\left(z+h,T\right)=\exp \left(h\left(\hat{D}+\hat{N}\right)\right)\mathbf{E}\left(z,T\right), $$

after propagation over distance $ h $ for time $ T $ , where

(4) $$ {\displaystyle \begin{array}{c}\hat{N}=-j\frac{8}{9}\gamma \hskip0.35em {\mathbf{E}}^{\mathrm{H}}\mathbf{E}-\frac{\alpha }{2},\\ {}\hat{D}=-\frac{j{\beta}_2}{2}\frac{\partial^2}{\partial {t}^2},\end{array}} $$

where $ {\mathbf{E}}^{\mathrm{H}} $ denotes the Hermitian transpose of $ \mathbf{E} $ . We can then assume (Millar et al., Reference Millar, Makovejs, Behrens, Hellerbrand, Killey, Bayvel and Savory2010)

(5) $$ \exp \left(h\left(\hat{D}+\hat{N}\right)\right)\mathbf{E}\left(z,T\right)\approx \exp \left(h\hat{D}\right)\exp \left(h\hat{N}\right)\mathbf{E}\left(z,T\right). $$

In this work, we utilize the symmetric SSFM, which involves taking a dispersive step of size $ \frac{h}{2} $ , followed by a nonlinear step of size $ h $ and an additional dispersive step of size $ \frac{h}{2} $ . This can be written as (Millar et al., Reference Millar, Makovejs, Behrens, Hellerbrand, Killey, Bayvel and Savory2010; Agrawal, Reference Agrawal2019)

(6) $$ \mathbf{E}\left(z+h,T\right)\approx \exp \left(\frac{h}{2}\hat{D}\right)\exp \left({L}_{\mathrm{eff}}\hat{N}\right)\exp \left(\frac{h}{2}\hat{D}\right)\mathbf{E}\left(z,T\right), $$

where $ {L}_{\mathrm{eff}}=\frac{1-\exp \left(-\alpha h\right)}{\alpha } $ is the fiber nonlinear effective length. It should be noted that here $ \alpha $ is measured in units of km $ {}^{-1} $ , rather than dB km $ {}^{-1} $ . Equation (6) can be solved iteratively using a fast Fourier transform to obtain an approximate solution for the transmission in the fiber (Ip and Kahn, Reference Ip and Kahn2008).

We also simulate polarization mode dispersion (PMD) using the waveplate model for birefringent fiber, such that each step of the SSFM includes a random polarization rotation and a frequency-dependent birefringence. We obtain the differential group delay between fast and slow axes $ \Delta {\tau}_p $ for a given fiber PMD $ \left\langle \Delta \tau \right\rangle $ as (Poole and Favin, Reference Poole and Favin1994)

(7) $$ \Delta {\tau}_p=\sqrt{\frac{3\pi }{8{n}_{\mathrm{sps}}}}\left\langle \Delta \tau \right\rangle, $$

where $ {n}_{\mathrm{sps}} $ is the number of steps per fiber span. To include PMD, a random polarization rotation is applied at the start of each symmetric SSFM step and the differential group delay is then included as part of each dispersive step of size $ \frac{h}{2} $ .

2.2. Digital backpropagation

DBP is a technique for the joint compensation of both linear and nonlinear impairments that accumulate in optical fibers. This is achieved by dividing the fiber link of interest up into small sections $ \Delta z $ and inversely solving the nonlinear Schrödinger equation (NLSE) which governs the propagation of optical signals through the optical fiber (Ip and Kahn, Reference Ip and Kahn2008). Specifically, we solve the Manakov equation using the SSFM, with the sign of the physical layer parameters reversed, that is, $ \alpha \to -\alpha $ , $ \gamma \to -\gamma $ , and $ {\beta}_2\to -{\beta}_2 $ . Thus, we simulate the propagation of the signal in reverse, in order to estimate the linear and nonlinear impairments accumulated in the span. We utilize a symmetric SSFM implementation, as outlined above. In this work, we utilized DBP with a fixed step size and a constant number of steps per span, although this can be adapted whilst performing DBP to improve the computational efficiency (Shao et al., Reference Shao, Liang and Kumar2014). All SSFM simulations, for both forward and reverse propagation, were performed with 100 steps per fiber span.

2.3. Gaussian process regression

GP regression is a probabilistic supervised ML method. GPs leverage Bayesian statistics to find the most likely function describing the relationship between a set of inputs and a set of outputs, given the data and a set of priors (Rasmussen and Williams, Reference Rasmussen and Williams2006). GP regression is a nonparametric ML method, meaning that rather than assuming a given parametric form and finding a set of parameters that describe the mapping, the space of functions is searched directly in a probabilistic way. The main benefit of using GPs over other supervised regression methods, such as neural networks (NNs), is that GPs provide a rigorous prediction error, the predictive variance. This can be used to judge the level of confidence associated with a given prediction, providing an extra layer of interpretability for the model user.

GPs are kernel methods, meaning that a kernel function is used to model the covariance between the data. Specifically, the kernel function specifies the covariance between the evaluation of the function being learned by the GP $ f\left(\mathbf{x}\right) $ at two input points (Rasmussen and Williams, Reference Rasmussen and Williams2006):

(8) $$ k\left({\mathbf{x}}_{\mathbf{i}},{\mathbf{x}}_{\mathbf{j}}\right)=\operatorname{cov}\left(f\left({\mathbf{x}}_{\mathbf{i}}\right),\hskip0.5em f\left({\mathbf{x}}_{\mathbf{j}}\right)\right). $$

In this work, we use the squared exponential kernel function, defined as (Rasmussen and Williams, Reference Rasmussen and Williams2006)

(9) $$ k\left({\mathbf{x}}_{\mathbf{i}},{\mathbf{x}}_{\mathbf{j}}\right)={h}_1^2\exp \left(\frac{-{\left\Vert {\mathbf{x}}_{\mathbf{i}}-{\mathbf{x}}_{\mathbf{j}}\right\Vert}^2}{2{{\mathbf{h}}_{\mathbf{2}}}^2}\right)+{h}_3^2{\delta}_{ij}, $$

where $ {x}_i $ and $ {x}_j $ are the input values of the points being compared, $ {h}_1 $ and $ {\mathbf{h}}_{\mathbf{2}} $ are the hyperparameters of the squared exponential kernel, $ {\delta}_{ij} $ is the Kronecker delta, that is, $ {\delta}_{ij}=1 $ if $ i=j $ and $ 0 $ otherwise, and $ \left\Vert \cdot \right\Vert $ represents the Euclidean distance. It should be noted that $ {\mathbf{h}}_{\mathbf{2}} $ is a vector with the same number of elements as the dimensionality of $ \mathbf{x} $ , whereas $ {h}_1 $ and $ {h}_3 $ are both scalars, as the targets have a single dimension. Commonly, the set of kernel hyperparameters are denoted by $ \theta $ , meaning for equation (9)), $ \theta \hskip0.5em \triangleq \hskip0.5em \left\{{h}_1,{\mathbf{h}}_{\mathbf{2}},{h}_3\right\} $ . $ {h}_1 $ acts as a global covariance scale for the GP, as it affects the magnitude of the covariances between all pairs of sample points $ {\mathbf{x}}_{\mathbf{i}},{\mathbf{x}}_{\mathbf{j}} $ . The role of $ {\mathbf{h}}_{\mathbf{2}} $ is to scale the distance between inputs $ {\mathbf{x}}_{\mathbf{i}},{\mathbf{x}}_{\mathbf{j}} $ that will yield a large covariance $ \operatorname{cov}\left(\hskip0.2em f\left({\mathbf{x}}_{\mathbf{i}}\right),f\left({\mathbf{x}}_{\mathbf{j}}\right)\right) $ , where each component of $ {\mathbf{h}}_{\mathbf{2}} $ scales a given input dimension. Therefore, larger values of the correlation length indicate that the function $ f\left(\mathbf{x}\right) $ varies less with respect to changes in a given component of $ \mathbf{x} $ , and vice versa. $ {h}_3 $ represents the noise of the observations—larger values of $ {h}_3 $ indicate that the GP is interpreting more of the variation in the data as noise, rather than the underlying signal. Using a kernel function to compute inner products in the original space is equivalent to using a large set of basis functions, facilitating efficient learning of complex functions (Rasmussen and Williams, Reference Rasmussen and Williams2006). When we fit the GP to data, we do so by maximizing the likelihood function with respect to the hyperparameters $ \theta ={h}_1 $ , $ {\mathbf{h}}_{\mathbf{2}} $ , and $ {h}_3 $ (Rasmussen and Williams, Reference Rasmussen and Williams2006). This corresponds to maximization of the probability of the targets given the data and hyperparameters—thus we find the underlying mapping between the inputs and targets in a probabilistic way. Specifically, this was achieved by computing the gradient of the log marginal likelihood with respect to the hyperparameters and performing a gradient-based maximization using the Limited memory Broyden–Fletcher–Goldfarb–Shanno Box constrained (L-BFGS) method (Liu and Nocedal, Reference Liu and Nocedal1989). This gradient is given by the expression (2.30) (Rasmussen and Williams, Reference Rasmussen and Williams2006, chap. 2):

(10) $$ \frac{\partial }{\partial {\theta}_j}\log p\left(\mathbf{y}|X,\theta \right)=\frac{1}{2}{\mathbf{y}}^T{K}^{-1}\frac{\partial K}{\partial {\theta}_j}{K}^{-1}\mathbf{y}-\frac{1}{2}\mathrm{tr}\left({K}^{-1}\frac{\partial K}{\partial {\theta}_j}\right), $$

where $ K $ is the matrix of the kernel evaluated at all pairs of input points and $ \mathbf{y} $ is the vector of targets. An in-depth discussion of this process is provided in Rasmussen and Williams (Reference Rasmussen and Williams2006, chap. 5, 5.3–5.7). Choosing a given kernel corresponds to making some assumptions about the data, as different kernels result in a different prior distribution over functions. We choose equation (9) because we expect the underlying functions to be well-described by a single-length scale with Gaussian noise, and we do not expect the data to exhibit properties such as periodicity and decay that may require more complex kernels (Rasmussen and Williams, Reference Rasmussen and Williams2006). We also tried using the more general Matérn kernel (Rasmussen and Williams, Reference Rasmussen and Williams2006, chap. 4), but saw no improvement relative to equation (9), and therefore we used the more simplistic kernel which was faster to fit. In this work, we utilize an implementation provided by the Multi-Output Gaussian Process Emulator library (Daub and Barlow, Reference Daub and Barlow2021).

The GP used in this work is constructed using a zero mean function and the covariance function given by equation (9). Therefore, we have a GP prior over our unknown function $ f\sim \mathrm{GP} $ , and we can write the prior over the function value at a specific input $ \mathbf{x} $ as:

(11) $$ P\left(f\left(\mathbf{x}\right)\right)=\mathcal{N}\left(f\left(\mathbf{x}\right);0,K\right), $$

where $ P $ denotes a probability, $ \mathcal{N} $ is a multivariate Gaussian distribution, and $ K $ is the kernel matrix, such that $ {K}_{i,j}=k\left({\mathbf{x}}_{\mathbf{i}},{\mathbf{x}}_{\mathbf{j}}\right) $ . It should be noted that the majority of applications of GPs utilize a zero mean function, as these models are able to learn a vast number of functions accurately within the range of the measured data (Rasmussen and Williams, Reference Rasmussen and Williams2006). As we have strong a priori bounds on the range over which we wish to make predictions, the incorporation of a specific mean function is not required. This GP model was observed to be of sufficient complexity to act as a high-accuracy surrogate model for this problem. For a detailed description of the theory of GPs, we refer the reader to Bishop (Reference Bishop2006, chap. 6) and Rasmussen and Williams (Reference Rasmussen and Williams2006). Moreover, we note that relevance vector machines (Candela and Hansen, Reference Candela and Hansen2004; Martino and Read, Reference Martino and Read2021) (RVMs) are very similar in formulation to GPs. There is an inherent trade-off between these two approaches—RVMs offer greater flexibility in the choice of basis functions than GPs, however, GPs have more interpretable predictive variance than RVMs. Exploring the use of RVMs in this technique forms part of the planned future work.

2.4. Simplicial homology global optimization

SHGO is a general-purpose global optimization algorithm based on algebraic topology and simplicial integral homology (Endres et al., Reference Endres, Sandrock and Focke2018). SHGO has been chosen to optimize the function learned by the GP in this work as it offers a number of theoretical and practical advantages. First, SHGO is theoretically proven to find the global optimum, even for the highly general case of noncontinuous, nonconvex, and nonsmooth functions (Endres et al., Reference Endres, Sandrock and Focke2018). In addition, SHGO has been shown to be highly performant compared to other state-of-the-art global optimization algorithms (Endres et al., Reference Endres, Sandrock and Focke2018). Also, SHGO accepts bounds on the variables being optimized, allowing us to input some a priori knowledge in the form of broad parameter ranges derived from fiber specification sheets and the literature. Finally, SHGO does not require the derivative of the objective function, eliminating an additional source of error given that the objective function is in this case learned by the GP. We utilized the SciPy implementation of SHGO in this work (Virtanen et al., Reference Virtanen, Gommers, Oliphant, Haberland, Reddy, Cournapeau, Burovski, Peterson, Weckesser, Bright, van der Walt, Brett, Wilson, Millman, Mayorov, Nelson, Jones, Kern, Larson, Carey, Polat, Feng, Moore, Van der Plas, Laxalde, Perktold, Cimrman, Henriksen, Quintero, Harris, Archibald, Ribeiro, Pedregosa and van Mulbregt2020).

5.1. Experimental setup

The experimental setup consists of transmission of a four-subcarrier 49.5 GBd polarization division multiplexed DP-256QAM superchannel across 1,010.75 km using a recirculating fiber loop. The experimental setup used is outlined in Figure 1.

Figure 1. Diagram of the experimental setup, consisting of a recirculating fiber loop that was used to transmit a 4 × 49.5 GBd superchannel over 13 spans of length 77.75 km. DAC, digital-to-analog converter; DP-IQ-MOD, dual-polarization IQ modulator; ECL, external cavity laser; MUX, multiplexer; PS, polarization scrambler.

Four external cavity lasers (ECLs) with <100 kHz linewidth and 50 GHz spacing were used as sources for two odd and two even subcarriers, which were modulated by two separate dual-polarization IQ modulators driven by two 33-GHz 92-GSa/s digital-to-analog converters (DACs). The 256-QAM driving signals were each spectrally shaped by matched root-raised cosine (RRC) filters with 1% roll-off. The frequency response of the transmitter components was compensated by applying digital pre-emphasis to the generated signals. The odd and even signals were modulated with independent data and combined with a 3-dB coupler to form a 4 × 49.5 GBd DP-256QAM superchannel with 200 GHz optical bandwidth. The superchannel was then amplified and launched into a recirculating fiber loop, which comprised of a polarization scrambler (PS), a span of 77.75 km of Corning^® SMF-28^® ULL fiber, two variable optical attenuators (VOAs), three erbium-doped fiber amplifiers (EDFAs) with 4.5 dB noise figure and a 2 nm optical bandpass filter to remove noise outside the signal bandwidth. The transmitted superchannel circulated through the loop 13 times, corresponding to a total transmission distance of 1,010.75 km. It was then mixed with a <100 kHz linewidth local oscillator (LO) laser using a 90° hybrid and detected by four balanced photo-detectors with 100 GHz electrical bandwidth. Finally, the received signals were digitized using a 256 GSa/s real-time digital sampling oscilloscope with 100 GHz analog electrical bandwidth.

The offline digital signal processing (DSP) chain was conducted as outlined in Wakayama et al. (Reference Wakayama, Sillekens, Galdino, Lavery, Killey and Bayvel2019) using QPSK pilot symbols. Specifically, the pilot symbols are used for equalization and frame synchronization, as well as for carrier phase estimation (CPE). A block diagram of the DSP is given in Figure 2. After the signal is received, it is converted to the electrical domain using an analog-to-digital converter (ADC). It is then deskewed and normalized. Following this, if desired, DBP is performed on the entire four-subcarrier superchannel. If we do not perform DBP, that is, to quantify the DBP gain, we perform chromatic dispersion compensation (CDC) instead. Following either of these blocks, we then shift the signal by 25 GHz to move the third subchannel to baseband. This is then followed by resampling to two samples per symbol, RRC-matched filtering to extract the third subchannel, and then frame synchronization. The filtered subchannel is equalized by an adaptive FIR filter, which utilizes two algorithms for training the equalizer filter weights; recursive least squares for coarse-tuning and least mean squares for fine-tuning, which are both trained on the received pilot symbols (Wakayama et al., Reference Wakayama, Sillekens, Galdino, Lavery, Killey and Bayvel2019). Having trained the equalizer, the optimized filter weights are then used to equalize the payload symbols. The signal is then demultiplexed into its two constituent orthogonal polarizations before we perform frequency offset estimation and CPE. The CPE algorithm is implemented as a two-step process, consisting of initial linear interpolation followed by more detailed residual phase estimation (RPE), as outlined in detail in Wakayama et al. (Reference Wakayama, Sillekens, Galdino, Lavery, Killey and Bayvel2019). Following this, we apply the Gram–Schmidt organization procedure (GSOP) and then estimate the SNR. This estimation consists of taking the noisy symbols outputted by the DSP chain and calculating the SNR as

(12) $$ \mathrm{SNR}=\frac{\unicode{x1D53C}\left[{\left|X\right|}^2\right]}{\unicode{x1D53C}\left[{\left|X-Y\right|}^2\right]}, $$

where $ \unicode{x1D53C}\left[\cdot \right] $ denotes the expectation, $ X $ are the transmitted symbols, and $ Y $ is the received symbols after DSP, and a k-means-based algorithm was used to extract the constellation centroid.

Figure 2. Schematic of the DSP chain, adapted from Wakayama et al. (Reference Wakayama, Sillekens, Galdino, Lavery, Killey and Bayvel2019). CDC, chromatic dispersion compensation; CPE, carrier-phase estimation; DBP, digital backpropagation; FOE, frequency offset estimation; GSOP, Gram–Schmidt orthogonalization procedure; LO, local oscillator; RPE, residual phase estimation; RRC, root-raised cosine.

5.2. Experimental a priori knowledge

The data sheet of the Corning® SMF-28® ULL fiber used in this experiment provides us with some a priori knowledge. First, $ \alpha $ is stated to have a maximum value of 0.17 dB km $ {}^{-1} $ . Thus, we consider a range of $ \alpha $ values between 0.14 and 0.17 dB km $ {}^{-1} $ . As GPs are a kernel-based method, they estimate the similarity between neighboring data to infer the underlying function. As a result, the model accuracy is limited near the edge of the range of training data, as there are fewer neighboring data. Thus, we want to avoid the ground truth value lying too close to the edge of the range, where the model accuracy will suffer. If any of the estimated parameters were observed to fall close to any of the corresponding a priori bounds, the bounds should be increased. The lower bound in this case is chosen to be lower than the world record of 0.142 dB km $ {}^{-1} $ (Hasegawa et al., Reference Hasegawa, Tamura, Sakuma, Kawaguchi, Yamamoto and Koyan2018), to avoid the estimated value falling too close to the lower bound. It was observed that the upper bound was sufficiently conservative, such that the estimated values were not close to this value. Additionally, this fiber type has a specified expected nonlinear refractive index component value $ {n}_2=2.1\times {10}^{-20} $ , however, no range is given. Similarly, the expected fiber effective area is given as 83 $ {\unicode{x03BC} \mathrm{m}}^2 $ , with no range or tolerance. However, the value of $ \gamma $ has significant uncertainty for deployed fibers, which may suffer from aging effects, or in some cases, the fiber type may even be unknown (Seve et al., Reference Seve, Pesic, Delezoide, Giorgetti, Sgambelluri, Sambo, Bigo and Pointurier2019). Moreover, variations during the manufacturing process can modify the value of $ \gamma . $ Thus, we define a broad range for $ \gamma $ , from 0.9 to 1.2 W $ {}^{-1} $ km $ {}^{-1} $ . As for the dispersion, we have stronger a priori knowledge as it is possible to obtain a measured estimate of the dispersion from the adaptive equalizer in the DSP chain (Corsini et al., Reference Corsini, Peracchi, Matarazzo, Foggi, Nijhof, Meloni, Potì, Magri and Ciaramella2013, eq. 4). Thus, we use a tight range around this estimated value, from 15.9 to 16.1 ps nm $ {}^{-1} $ km $ {}^{-1} $ , to allow the SHGO optimizer to tweak the dispersion to the value that yields optimal gain. It should be noted that, for deployed systems, the uncertainty surrounding these parameters is significantly higher than in laboratory experiments (Pointurier, Reference Pointurier2017, Reference Pointurier2021).

Thus, as well as being a useful research tool that avoids time-consuming and nonoptimal grid search, the proposed method provides a systematic way to obtain a high accuracy characterization of deployed systems from broad a priori ranges. To demonstrate the feasibility of using this methodology in deployed networks, we assume relatively broad a priori bounds, defined from our approximate a priori knowledge.

6.1. Parameter estimates

We first test the efficacy of this technique using simulated traces, generated for a simulated system using the SSFM. This system was designed to emulate the experimental system outlined in Figure 1, with a superchannel consisting of four subcarriers which are used to transmit a 49.5 GBd DP-256QAM signal over 13 spans of length 77.75 km, for a total transmitted distance of 1,010.75 km. As we know exactly the parameters that we used to generate the traces, this provides us with a ground truth that we can use to verify the accuracy of the physical layer parameter estimation. The ground truth values of $ D $ , $ \gamma $ , and $ \alpha $ are shown in Table 1, along with the a priori search ranges used. The same ranges were used for the experimental traces. It should be noted that we model the effects of polarization mode dispersion (PMD) in our simulation, as outlined above, using a value of 0.04 ps km $ {}^{-\frac{1}{2}} $ . This value is the link design value taken from the manufacturer specification sheet of the Corning^® SMF-28^® ULL fiber used in this experiment, a statistical upper limit on the total link PMD.

Table 1. A priori physical layer parameter ranges.

During training of the GP, we take 100 additional samples from the parameter space and record the SNR after application of DBP and DSP, for verification of the trained models. We compare the predictions of the GP models at launch powers of 1, 2, 3, and 4 dBm with the corresponding DBP chain SNRs. The mean absolute error (MAE) values were 0.009, 0.016, 0.010, and 0.012 and the MSE values were 0.013, 0.022, 0.013, and 0.017, respectively for the traces at 1, 2, 3, and 4 dBm. As the GP is trained on SNR in linear units, these errors are given in terms of linear SNR. This indicates that 500 samples are sufficient for training a high-accuracy GP model, which learns the functional mapping from the physical layer parameters to the SNR gain with a small prediction error. This is crucial, as the GP is used as a computationally cheap surrogate model for the DBP chain, which is used by the SHGO optimizer to find the parameters that maximize the DBP gain. Exploring the precise minimum number of samples required for a good estimation forms part of the planned future work. It is important to note the motivation for using the GP to learn a surrogate model for the SNR after DBP in conjunction with SHGO, rather than simply using SHGO to optimize the SNR after DBP directly. Specifically, the bottleneck in terms of performance is running the DBP chain. It was observed that training the GP surrogate model required running the DBP chain 500 times, whereas SHGO required approximately 1,500 function evaluations to converge. Thus, using SHGO directly on the DBP chain would be more computationally expensive, as the DBP chain would need to be run a larger number of times. In addition, a trained GP surrogate model is itself a useful tool, providing an interpretable model of the SNR as a function of $ D $ , $ \gamma $ , and $ \alpha $ that can be used to perform a sensitivity analysis. This is explored further in Section 7.2.

We use the simulated traces at launch powers per channel of 1, 2, 3, and 4 dBm to perform the estimation, and the optimal parameter values obtained are shown in Table 2. For the traces at 1, 2, and 3 dBm, one global optimum was found, whereas for the trace at 4 dBm two optima with equivalent performance were found. These traces were chosen as they correspond to the four values around the peak in the experimental system after DBP is applied. We can see that we are able to correctly estimate the ground truth value of $ D $ to a precision of 4 significant figures. Moreover, we are able to estimate $ \gamma $ with a maximum error of 0.02 W $ {}^{-1} $ km $ {}^{-1} $ . As for $ \alpha $ , the estimation is accurate to within an error of 0.001 dB km $ {}^{-1} $ . Thus, this technique is able to obtain highly accurate estimates of $ D $ , $ \alpha $ , and $ \gamma $ from broad a priori parameter ranges for a simulated system. Sources of error in the simulated system include PMD, numerical errors introduced by the DSP and interpolation errors introduced by the GP surrogate model. As a representative example, Figure 3 shows the comparison of the SNR after decoding the simulated traces both with and without DBP for using the parameters $ D=16.00 $ ps nm $ {}^{-1} $ km $ {}^{-1} $ , $ \gamma =1.02 $ W km $ {}^{-1} $ , and $ \alpha =0.156 $ dB km $ {}^{-1} $ , estimated using the trace at 4 dBm. We show the achieved SNR in Figure 3a and the SNR gain at each launch power in Figure 3b. The gain in maximum SNR was 1.14 dB with all five sets of parameters found. Moreover, the sum of the gain for each trace, that is, the difference between the performance with and without DBP for each trace, was 10.10 dB for the parameters found using all traces. Thus, as expected, these parameter sets give practically equivalent performance.

Table 2. Parameter estimates for simulated traces.

Figure 3. Simulated traces comparison with and without DBP using parameters $ D=16.00 $ ps nm $ {}^{-1} $ km $ {}^{-1} $ , $ \gamma =1.02 $ W $ {}^{-1} $ km $ {}^{-1} $ , and $ \alpha =0.156 $ dB km $ {}^{-1} $ , estimated using Algorithm 1 using the simulated trace at a launch power of 4 dBm. The SNR versus launch power curves with and without DBP are shown, as well as the DBP gain at each launch power.

6.2. GP model parameter sweeps

In order to further interpret the trained GP model, we perform parameter sweeps and plot the predictive mean and predictive variance. Specifically, this is done by sweeping each parameter in turn with the others held fixed at the optimum found, for each of the GP models trained. The results of these sweeps are shown for $ D $ , $ \gamma $ , and $ \alpha $ in Figures 4–6, respectively. The ranges for each sweep are the same as those defined in Table 1, meaning that we are sweeping across the full a priori search range learned by each GP. A confidence region corresponding to one predictive standard deviation is also shown, where the predictive standard deviation is simply given by the square root of the predictive variance of the GP model. From the confidence bands, we can see that the GP model has the lowest uncertainty with respect to $ D $ , followed by $ \gamma $ and $ \alpha $ . Additionally, from the variation of the predictive mean of the GP across the range, we can see that $ D $ has the largest effect on the SNR after DBP, followed by $ \gamma $ and then $ \alpha $ . Specifically, the ranges of variation are 0.58, 0.17, and 0.03 dB for $ D $ , $ \gamma $ , and $ \alpha $ , respectively. We also consider these values as a fraction of the DBP gain achieved for the trace used in the estimation. For the trace at 4 dBm shown, the ranges for $ D $ , $ \gamma $ , and $ \alpha $ correspond to 16.7, 4.8, and 0.8% of the achieved gain respectively. However, it is important to contextualize these results—the objectives of this methodology are to estimate the parameters $ D $ , $ \alpha $ , and $ \gamma $ at the receiver, as well as to provide a methodology for systematically finding the optimal parameters in terms of DBP gain. Even if the effect on the DBP is small for $ \alpha $ and $ \gamma $ , estimating these parameters with high accuracy is crucial for building high-accuracy network models, which in turn can be used to unlock extra capacity by reducing network margins.

Figure 4. $ D $ GP model sweep, in which we sweep $ D $ across the a priori search range with $ \gamma =1.03 $ W $ {}^{-1} $ km $ {}^{-1} $ and $ \alpha =0.156 $ dB km $ {}^{-1} $ . The hyperparameters for this model are $ {h}_1=21.5, $ $ {\mathbf{h}}_{\mathbf{2}}=\left(0.01\;\mathrm{ps}\;{\mathrm{nm}}^{-1}\;{\mathrm{km}}^{-1},0.84\;{\mathrm{W}}^{-1}\;{\mathrm{km}}^{-1},0.11\;\mathrm{dB}\;{\mathrm{km}}^{-1}\right) $ , and $ {h}_3=0.01 $ .

Figure 5. $ \gamma $ GP model sweep, in which we sweep $ \gamma $ across the a priori search range with $ D=16.00 $ ps nm $ {}^{-1} $ km $ {}^{-1} $ and $ \alpha =0.156 $ dB km $ {}^{-1} $ . The hyperparameters for this model are $ {h}_1=21.5 $ , $ {\mathbf{h}}_{\mathbf{2}}=\left(0.01\;\mathrm{ps}\;{\mathrm{nm}}^{-1}\;{\mathrm{km}}^{-1},0.84\;{\mathrm{W}}^{-1}\;{\mathrm{km}}^{-1},0.11\;\mathrm{dB}\;{\mathrm{km}}^{-1}\right) $ , and $ {h}_3=0.01 $ .

Figure 6. $ \alpha $ GP model sweep, in which we sweep $ \alpha $ across the a priori search range with $ D=16.00 $ ps nm $ {}^{-1} $ km $ {}^{-1} $ and $ \gamma =1.03 $ W $ {}^{-1} $ km $ {}^{-1} $ . The hyperparameters for this model are $ {h}_1=21.5 $ , $ {\mathbf{h}}_{\mathbf{2}}=\left(0.01\;\mathrm{ps}\;{\mathrm{nm}}^{-1}\;{\mathrm{km}}^{-1},0.84\;{\mathrm{W}}^{-1}\;{\mathrm{km}}^{-1},0.11\;\mathrm{dB}\;{\mathrm{km}}^{-1}\right) $ , and $ {h}_3=0.01 $ .

7.1. Parameter estimates

Having confirmed the efficacy of Algorithm 1 for simulated traces, we demonstrate the DBP parameter estimation process for experimentally measured traces, using the experimental system outlined in Figure 1. We perform the parameter estimation process outlined in Algorithm 1 for the experimental traces measured at launch powers 1.22, 2.27, 3.24, and 4.25 dBm, denoted $ {P}_1 $ , $ {P}_2 $ , $ {P}_3 $ , and $ {P}_4 $ , respectively. These particular traces were chosen due to their proximity to the optimum launch power, once DBP is applied.

In order to validate the GP model learned, we again draw an additional 100 validation samples from the DBP chain and compare these samples to the predictions of the trained GP model at these sample inputs. This yielded a linear SNR mean absolute error (MAE) of 0.015, 0.032, 0.012, and 0.010 for the GPs trained on traces at $ {P}_1 $ , $ {P}_2 $ , $ {P}_3 $ , and $ {P}_4 $ , respectively. The corresponding linear SNR RMSE values were 0.020, 0.049, 0.019, and 0.013, respectively. Thus, the GP trained on the trace at $ {P}_2 $ has the highest validation error, whereas the model trained at $ {P}_4 $ has the lowest validation error. The magnitude of the validation errors are very small in practical terms, indicating that 500 samples is sufficient for the GP to learn an accurate probabilistic mapping from $ D $ , $ \gamma $ , and $ \alpha $ to the SNR after DBP.

In Figure 7, we show the estimation results for each trace. SHGO returned 10, 3, 2, and 2 local optima respectively for the traces at $ {P}_1 $ , $ {P}_2 $ , $ {P}_3 $ , and $ {P}_4 $ . These local optima were observed to give equivalent performance in terms of DBP gain, to within a maximum tolerance of 0.006 dB across all the experimentally measured traces. Specifically, this tolerance was calculated by performing DBP for all the experimental traces with each set of parameters and calculating the difference between the mean SNR after DBP is applied for the best and worst-performing parameter sets. This highlights the limitations of considering approaches such as gradient descent, which are prone to getting stuck in a given local optima. It should be noted that $ D $ was estimated to be $ 16.05 $ ps nm $ {}^{-1} $ km $ {}^{-1} $ for all optima using traces $ {P}_2 $ , $ {P}_3 $ , and $ {P}_4 $ , whereas $ D=16.06 $ $ \mathrm{ps}\;{\mathrm{nm}}^{-1}\;{\mathrm{km}}^{-1} $ for all optima estimated using the trace measured at $ {P}_1 $ . It is clear that there is a correlation between $ \alpha $ and $ \gamma $ : we can achieve equivalent performance with higher $ \gamma $ if $ \alpha $ is increased, and vice versa. The reason for using a search range that extends down to 0.140 dB km $ {}^{-1} $ is that GPs are a kernel-based method. Thus, the accuracy of the model will be much lower at the boundaries of the training data range, and so we want to avoid the optimum falling too close to the edge of the range, where the model accuracy will be reduced. If the optimum fell near the edge of the range, it would indicate to the user that the a priori ranges should be increased, such that the optimum is not too close to the edge of the range. In addition, we indicate that lower values of $ \gamma $ are more likely to reflect the true system values. This is because the fiber specifications list an effective area of 83 $ \mu \mathrm{m} $ and $ {n}_2=2.1\times {10}^{-20} $ , which corresponds to $ \gamma =1.03 $ W $ {}^{-1} $ km $ {}^{-1} $ . Thus, we would expect parameter sets with a smaller deviation from this value to be more likely. For clarity, we show the effective area values for each $ \gamma $ value on Figure 7. Thus, we consider parameter sets with the lowest value of $ \gamma $ in the feasible $ \alpha $ range to be the most likely parameter estimates. For the traces measured at launch powers $ {P}_2 $ , $ {P}_3 $ , and $ {P}_4 $ , there is one clear choice for the most likely set of parameters. For the trace measured at $ {P}_1 $ , there are three traces which we consider to be equally likely. These parameter best estimates are shown in Table 3 for each of the traces, including the three estimates from the trace at $ {P}_1 $ .

Figure 7. Plot of local optima obtained from the experimental traces at launch powers $ {P}_1=1.22 $ , $ {P}_2=2.27 $ , $ {P}_3=3.24 $ , and $ {P}_4=4.25 $ dBm. All optima have equivalent performance in terms of DBP gain, to a maximum tolerance of 0.006 dB. For traces at $ {P}_2 $ , $ {P}_3 $ , and $ {P}_4 $ , $ D=16.05 $ ps nm $ {}^{-1} $ km $ {}^{-1} $ , whereas for $ {P}_1 $ , $ D=16.06 $ ps nm $ {}^{-1} $ km $ {}^{-1} $ . The independently measured fiber loss of $ \alpha >0.156 $ dB km $ {}^{-1} $ is indicated on the plot.

Table 3. Most likely parameter estimates for the experimental system.

Additionally, estimated parameter values that are more consistent across multiple traces are more likely to be correct—we can see that a grouping of similar values emerges for all the traces in the range for $ \gamma $ values from 1.07 to 1.10 W $ {}^{-1} $ km $ {}^{-1} $ and $ \alpha $ values from 0.151 to 0.157 $ \mathrm{dB}\;{\mathrm{km}}^{-1} $ . In deployment, parameter estimates from measurements of multiple lightpaths occupying the same fiber segment could be used to cross-reference and identify more likely parameter estimates. We also notice that the number of local optima is observed to increase as we decrease the launch power. This is due to the fact that, at higher launch power, we have more nonlinear interference noise. Thus, the sensitivity of the SSFM to $ \gamma $ increases, making it easier to estimate, as fewer values of $ \gamma $ have equivalent performance. As a demonstrative example, we performed the parameter estimation process using the trace measured at −9.75 dBm, which returned 129 local optima. In this linear regime, $ \gamma $ has negligible effect, so SHGO finds a large number of local optima with equivalent performance.

To consider how these parameters perform across all the experimentally measured traces, we record the DBP gain, both in terms of the maximum SNR achieved and the SNR gain for each trace. All sets of estimated parameters yield a gain in maximum SNR of 0.83 dB. As for the total gain, taken as the sum of individual gain at each trace, gain values of 10.92, 10.95, 10.94, and 10.95 dB were achieved using the parameters estimated using the traces at $ {P}_1 $ , $ {P}_2 $ , $ {P}_3 $ , and $ {P}_4 $ , respectively. Here, 10.92 dB was achieved for each of the local optima estimated using the trace at $ {P}_1 $ . Thus, there is a small difference in the achieved gain, with the parameters estimated using the $ {P}_4 $ trace giving the highest gain, and the parameters estimated using the trace at $ {P}_1 $ giving the lowest gain. As an representative example, Figure 8 shows the DBP gain given by the parameters estimated for the trace at $ {P}_4 $ , for each of the experimental traces. The SNR achieved is compared with and without DBP in Figure 8a and the SNR gain is presented in Figure 8b. As in simulation, sources of error include PMD, DSP errors, and inaccuracies in the GP surrogate model. For the experimental case, there are additional sources of error arising from the nonideal nature of the experimental components. These sources of error lead to increased variation in the estimation of the physical layer parameters $ D $ , $ \gamma $ , and $ \alpha $ , as compared to the simulated case. Moreover, due to these error sources, in the experimental case we obtain multiple local optima with equivalent gain for all traces, whereas in simulation, only one trace returned more than one optimum.

Figure 8. Experimental traces comparison with and without DBP using parameters $ D=16.05 $ ps nm $ {}^{-1} $ km $ {}^{-1} $ , $ \gamma =1.07 $ W $ {}^{-1} $ km $ {}^{-1} $ , and $ \alpha =0.153 $ dB km $ {}^{-1} $ , estimated using Algorithm 1 using the experimental trace at a launch power of 4.25 dBm. The SNR versus launch power curves with and without DBP are shown, as well as the DBP gain at each launch power.

7.2. GP model parameter sweeps

We again perform parameter sweeps of the trained GP model for the experimental case, sweeping each parameter in turn, with the other two parameters fixed at the optimal values. These sweeps are shown in Figures 9–11 for $ D $ , $ \gamma $ , and $ \alpha $ , respectively. The hyperparameters found to maximize the marginal likelihood during the GP fitting process are also reported. As with the simulated traces, $ D $ has the largest variation in SNR across the a priori search range, followed by $ \gamma $ and $ \alpha $ . Thus, as in simulation, $ D $ is the most important parameter in terms of maximizing the DBP gain, followed by $ \gamma $ and then $ \alpha $ . Specifically, the ranges found are 0.76, 0.13, and 0.02 dB for $ D $ , $ \gamma $ , and $ \alpha $ , respectively. We also consider these values as a fraction of the DBP gain achieved for the trace used in the estimation. For the trace at 4.25 dBm shown, the ranges for $ D $ , $ \gamma $ , and $ \alpha $ correspond to 22.4, 3.9, and 0.7% of the achieved gain respectively. Also, we observe that the predictive variance is smallest for $ D $ , followed by $ \gamma $ , and then $ \alpha $ . An intuitive explanation for this is that the larger the sensitivity of the SNR to a given parameter, the more confident the GP model is in the learned underlying function.

Figure 9. $ D $ GP model sweep, in which we sweep $ D $ across the a priori search range with $ \gamma =1.07 $ W $ {}^{-1} $ km $ {}^{-1} $ and $ \alpha =0.153 $ dB km $ {}^{-1} $ . The hyperparameters for this model are $ {h}_1=13.6, $ $ {\mathbf{h}}_{\mathbf{2}}=\left(0.03\;\mathrm{ps}\;{\mathrm{nm}}^{-1}{\mathrm{km}}^{-1},0.26\;{\mathrm{W}}^{-1}{\mathrm{km}}^{-1},0.08\;\mathrm{dB}\;{\mathrm{km}}^{-1}\right) $ , and $ {h}_3=0.01 $ .

Figure 10. $ \gamma $ GP model sweep, in which we sweep $ \gamma $ across the a priori search range with $ D=16.05 $ ps nm $ {}^{-1} $ km $ {}^{-1} $ and $ \alpha =0.153 $ dB km $ {}^{-1} $ . The hyperparameters for this model are $ {h}_1=13.6 $ , $ {\mathbf{h}}_{\mathbf{2}}=\left(0.03\;\mathrm{ps}\;{\mathrm{nm}}^{-1}{\mathrm{km}}^{-1},0.26\;{\mathrm{W}}^{-1}{\mathrm{km}}^{-1},0.08\;\mathrm{dB}\;{\mathrm{km}}^{-1}\right) $ , and $ {h}_3=0.01 $ .

Figure 11. $ \alpha $ GP model sweep, in which we sweep $ \alpha $ across the a priori search range with $ D=16.05 $ ps nm $ {}^{-1} $ km $ {}^{-1} $ and $ \gamma =1.07 $ W $ {}^{-1} $ km $ {}^{-1} $ . The hyperparameters for this model are $ {h}_1=13.6 $ , $ {\mathbf{h}}_{\mathbf{2}}=\left(0.03\;\mathrm{ps}\;{\mathrm{nm}}^{-1}{\mathrm{km}}^{-1},0.26\;{\mathrm{W}}^{-1}{\mathrm{km}}^{-1},0.08\;\mathrm{dB}\;{\mathrm{km}}^{-1}\right) $ , and $ {h}_3=0.01 $ .

GP sweeps can also be used to illuminate the issues associated with manually optimizing $ D $ , $ \alpha $ , and $ \gamma $ in terms of DBP. Crucially, if one or more of the other parameters is set to a nonoptimal value, the sweeps can change significantly, resulting in an incorrect optimum value. This is because $ \alpha $ and $ \gamma $ jointly contribute to the nonlinear noise, meaning that there will not be uniquely optimal values of $ \alpha $ and $ \gamma $ with respect to the SNR after DBP. This can be seen in Chen and Shieh (Reference Chen and Shieh2010, eqs. 31–Reference Wakayama, Sillekens, Galdino, Lavery, Killey and Bayvel33). Thus, performing a one-at-a-time grid search is likely to lead to an incorrect optimum. To demonstrate the benefits of the multi-parameter optimization framework proposed as compared to single-parameter-at-a-time grid search, Figure 12 shows parameter sweeps over $ \alpha $ for the 4.25 dBm experimental trace with $ \gamma $ set to the nonoptimal value of $ \gamma =1.00 $ . We can see that for $ \gamma =1.00 $ W $ {}^{-1} $ km $ {}^{-1} $ , the optimum value of $ \alpha $ occurs at 0.143 dB km $ {}^{-1} $ , with a lower maximum SNR than the global optimum reported in Table 3. This highlights the limitations of doing a one-at-a-time grid search—the optimum found in one dimension is not in general independent of the values of the other parameters.

Figure 12. Demonstrative example $ \alpha $ GP model sweep, in which we sweep $ \alpha $ across the a priori search range with $ D=16.05 $ ps nm $ {}^{-1} $ km $ {}^{-1} $ and $ \gamma =1.00 $ W $ {}^{-1} $ km $ {}^{-1} $ —a deliberately nonoptimal value. This highlights the limitations of one-at-a-time grid search optimization: the optimal value for a single parameter sweep depends on the values of the other parameters.

It is clear that having an incorrect value of $ \gamma $ results in an incorrect optimum at 0.143 dB km $ {}^{-1} $ , at a lower SNR than that obtained using SHGO.