## 2 Materials and Methods
### Validation  
The first dataset consists of, as input data, the scattering parameters, \(|S_{11}|\) and \(|S_{21}|\), or the impedance parameters, \(\text{Re}(Z_{in})\) and \(\text{Im}(Z_{in})\). As output parameters, we have the distance \(d\) between the two coils and the configuration parameters \(\Delta x\), \(\Delta y\), and \(\theta\). With symbolic regression, the goal is to best fit these data with a multi-variable function:  
\[(d,\Delta x,\Delta y,\theta)=f_{1}(|S_{11},|S_{21}|)=f_{2}(Re(Z_{in}),\,Im(Z_{ in})) \tag{1}\]  
The second dataset consists of, as input data, the WPT frequency, the distance between the two coils, \(z_{0}\), and the load resistance, \(R\). As output parameters, we have the phase between the input current and voltage. With symbolic regression, the goal is to best fit these data with a multi-variable function:  
Figure 5: Data measured at fixed coils distance of \(z_{0}=4\) mm, for six different load conditions. Is it clear that data regularization is needed due to measurement error. Figure 6: Data measured at fixed load condition of \(R=100\)\(\Omega\), for three different distance coil conditions, with \(z_{0}=0,4\) and \(8\) mm. Is it clear that data regularization is needed due to the measurement error. \[phase=f_{3}(f,R,z_{0}) \tag{2}\]  
The SR algorithm pySR [46] is used. It is a Python 3.11.11 library that implements genetic algorithms to find the best formulas, with a cost function that, as previously mentioned, is inversely proportional to the MSE and the complexity of the formula. This tool allows you to define custom binary and unary operators, enabling the search for formulas with a given structure. Using pySR on a big dataset with more than 10,000 data points is not convenient since it significantly slows the research and increases the complexity of the problem. In the first simulation dataset, there are \(229\cdot 57\cdot 1001=13,066,053\) points, so the dataset is reduced, considering only the measured points at the resonant frequency of the WPT, 13.56 MHz. For pySR, however, that is still a large amount of data, so the analysis is performed in two reduced ways. First, a single configuration at a time is considered, so only 57 points in which only the distance variable \(d\) varies. The second analysis is performed at a given distance, with 229 points in which the configuration parameters vary. The number of populations and the number of iterations are both set to 10 in pySR. Increasing them to 50 can increase the accuracy by a few percentage points but at a significant computational cost. On the other hand, no particular operator is used, only exponential and logarithmic functions. But with some idea of the final formula to reach the research space, it can be drastically reduced. Due to how the dataset is built, it is convenient to use \(Z_{in}\) for the experimental measurements, and S in the simulated ones. Indeed, going from S to \(Z_{in}\) or vice versa, is possible but would increase the uncertainty due to numerical errors. The second dataset contains 3055 data points, with 18 configurations: 3 for the distance between the coils and 6 for the load conditions. The measured dataset has a lot of noise and also some invalid values as it can be seen from the figure at a fixed load value, Figure 6, and at a fixed distance value, Figure 5. For this reason, a data regularization is necessary. Standard regularization methods do not yield good results for high data irregularity, resulting in either overly strong regularization that also modifies good measurement points, or an inability to adequately regularize the data due to being too simple. Therefore, a custom adjustment function is defined for this dataset. The idea is to change only those points that differ significantly from their neighbors. Phase control of each point with the previous one is used, and for phase differences greater than 5 degrees, the phase of the point is adjusted to 0.9 times the phase of the previous point.