Novel Method of Coupling Coefﬁcient Estimation Based on the Bifurcation Phenomena in Inductive Power Transfer

: Inductive power transfer (IPT) applications, such as stationary charging of electric vehicles (EVs), at least moderate coupling between the coils to achieve high efﬁciency, but the coefﬁcient k typically varies between of 0.1 to 0.4, depending on the displacement of the coils according to SAE J2954. Thus, accurate and reliable methods for estimation of k are required for positioning of the EV to achieve optimal alignment with the charging pad. Additionally, in IPT, numerous control strategies are available for regulating output power and optimizing system efﬁciency that require an accurate estimate of the mutual inductance or k . However, existing estimation methods tend to require detailed a-priori information of a large number of circuit parameters, or they need measurement of currents or voltages in both primary and secondary sides. This paper presents a preliminary evaluation of a novel, primary-side method to estimate k , which is based solely on the frequency response of the input phase while operating the system in bifurcation. The method does not require any additional measurements of the system parameters. The theoretical background of the method is presented together with the description of the measurement procedure. The method is experimentally veriﬁed and compared with two currently used estimation methods. According to the presented experimental evaluation, the proposed method estimates k with an error of 3.62% with respect to the reference over the evaluated range of 0.08 to 0.36. In addition, we demonstrate that the presented method is resilient to detuning.


Introduction
Inductive power transfer (IPT) is a developing technology for simplifying and streamlining the electric vehicle (EV) charging experience [1][2][3]. Based on the movement of the vehicle, this application can be divided in dynamic, semi-dynamic or stationary charging [4,5]. In dynamic charging, an EV is charged while it travels over an electrified road [6]. Multiple scenarios of the electrification track are possible: it can be formed by a single elongated primary coil (similar to materials handling systems), or by segmented track, where each segment powered or switched individually, or by a sequence of pads similar to those used for EV charging [6][7][8].
Semi-dynamic charging is derived from the dynamic charging, but instead of powering vehicles in-move, they are powered at the intersections, traffic lights or bus stops (electrified public transport vehicles), where the vehicles are temporary stationary [4].
Stationary IPT chargers replace conductive chargers and can be fully integrated in the parking space, so that no part of the charger is visible [9]. Several standards were developed for stationary EV charging, e.g., SAE J2954 [10] or IEC 61980 [11]. The modern IPT systems for EV charging achieve similar levels of efficiency (>90%) and transferred power (tens and hundreds of kW) as the conductive chargers [2,[12][13][14].
Unfortunately, the power transfer capability and efficiency of IPT systems are dependent on the coil alignment [15]. Thus, to ensure sufficient efficiency, the coils need to be correctly positioned [10], which requires feedback of their position or, preferably, the coupling coefficient k. Additionally, several control algorithms have been proposed, such as [16,17], which aim to provide power regulation and optimize efficiency under coil displacement. However, they require the value of the coupling coefficient k or mutual inductance M for operation.
There are several methods for estimating the coupling coefficient k [18][19][20][21], but they usually require knowledge of a large number circuit parameters [18,20,21] or require voltage or current information from both the transmitter and receiver [19]. If these parameters are measured inaccurately or they change (e.g., due to the displacement or temperature), the estimate quality will decrease. Also, some of the methods are based on complicated procedures requiring increased computational power [19,20].
In this paper, we present a preliminary evaluation of a novel primary-sided method for estimating the coupling coefficient, which is based solely on the frequency response of the input phase while the IPT system is intentionally operated in bifurcation. The proposed method also requires only a minimal number of calculations.
Bifurcation is a phenomenon during which the single frequency at which the input current is in phase with voltage, splits into three such frequencies. These frequencies are called the zero-phase angle (ZPA) frequencies. The original ZPA frequency is the resonance frequency f 0 , and the additional ZPA frequencies are denoted f ϕ1 and f ϕ2 in this work.
The ZPA frequency f ϕ1 is lower than f 0 , and f ϕ2 is higher than f 0 . The onset of bifurcation can be initiated for example by appropriately adjusting the load or coupling coefficient k. As bifurcation increases, the ZPA frequencies f ϕ1 and f ϕ2 approach their respective asymptotes f 1a and f 2a , which are functions solely of the coupling coefficient k and f 0 .
The proposed method is developed and applied for the series-series (SS) compensated IPT system. The most practical way to achieve deep bifurcation for the SS compensation is to short-circuit the output. Under such a condition, the frequency response of the input phase ϕ in is measured to find the ZPA frequencies f ϕ1 ≈ f 1a and f ϕ2 ≈ f 2a , which are then used to calculate the coupling coefficient k.
Since the method requires shorting the load, the method would not be used while power is being transferred to the EV. Instead, this method is applicable for initial estimation of k during start-up (i.e., before the charging begins) or intermittent estimation of k (i.e., in between periods of charging). In practice, the load could be short circuited using an active rectifier (in which both of the bottom transistors are turned on), through the DC/DC converter between the rectifier and EV battery, or a dedicated switch at the DC side after the rectifier.
In comparison with the current estimation methods, the proposed method requires (coarse) measurement of the frequency response of the input phase. However, its main advantage is that it does not require a-priori information about the system (except for an approximate frequency range for the frequency response measurement). It also requires only the primary side information (i.e., no wireless feedback between the secondary and primary is needed). The preliminary measurement results presented in this paper also show that the method is resilient to the detuning due to change of inductances in displacement. The method's estimate of k only differed from the reference value by an average absolute error of 3.62% over the range k = 0.08 to k = 0.36.
The second chapter of the paper describes the initial AC equivalent circuit model, of which the method is based. The third chapter follows by deriving of the asymptotes of the ZPA frequencies f ϕ1 and f ϕ2 . The fourth chapter outlines the presented method itself, which is then experimentally verified in the fifth chapter and compared to the selected currently used methods for coupling coefficient estimation. The practical aspects of the selected methods are further discussed in the sixth chapter.

AC Equivalent Circuit
The proposed k-parameter estimation method is based on the AC equivalent circuit model of the IPT system depicted in Figure 1. The arrangement represents a single transmitter and receiver with the series-series compensation topology. The primary and secondary coil inductances L p and L s , respectively, are compensated by capacitances C p and C s to form series resonant circuits with resonant frequencies f p and f s . The frequencies f p and f s are set to the same value of f 0 , which is the common resonant frequency f 0 . The capacitances C p and C s are then calculated as: where ω 0p = ω 0s = ω 0 = 2π f 0 for the primary tuned to the secondary.

AC Equivalent Circuit
The proposed -parameter estimation method is based on the AC equivalent circuit model of the IPT system depicted in Figure 1. The arrangement represents a single transmitter and receiver with the series-series compensation topology. The primary and secondary coil inductances and , respectively, are compensated by capacitances and to form series resonant circuits with resonant frequencies and . The frequencies and are set to the same value of , which is the common resonant frequency . The capacitances and are then calculated as: where 2 for the primary tuned to the secondary. The coils are coupled by the mutual inductance . The circuit is loaded by a purely resistive load, which is represented in Figure 1 by the equivalent load resistance . The equivalent series resistances (ESRs) representing losses of each side are summarized in and . The circuit is powered from a source, which is considered sinusoidal (first harmonic approximation). Therefore, all the voltages and currents are represented by phasors, which are marked as . The symbol without this marking refers to phasor amplitude. The input voltage is described by its amplitude , its frequency (the operating frequency of the circuit) and its phase with respect to the input current , which is taken as the reference.
is considered constant. The IPT system output is characterized by the output voltage and current . The input impedance of the circuit in Figure 1 is described by:

Theoretical Background of the Estimation Method
As discussed in the introduction, the phase bifurcation is the splitting of a single ZPA frequency into three. Specifically, the ZPA frequencies are described by condition: The formulas describing the phase bifurcation frequencies can be obtained from (4). To simplify the calculation, the normalization process presented by Wang et al. in [22] can be advantageously used. In the normalization, the ESRs and are neglected. The IPT model is normalized as follows. First, the normalized operating frequency is defined as the operating frequency divided by : The coils are coupled by the mutual inductance M. The circuit is loaded by a purely resistive load, which is represented in Figure 1 by the equivalent load resistance R L . The equivalent series resistances (ESRs) representing losses of each side are summarized in R p and R s . The circuit is powered from a source, which is considered sinusoidal (first harmonic approximation). Therefore, all the voltages and currents are represented by phasors, which are marked asX. The symbol X without this marking refers to phasor amplitude. The input voltageV in is described by its amplitude V in , its frequency f (the operating frequency of the circuit) and its phase ϕ in with respect to the input currentÎ in , which is taken as the reference. V in is considered constant. The IPT system output is characterized by the output voltageV out and currentÎ out .
The input impedance of the circuit in Figure 1 is described by:

Theoretical Background of the Estimation Method
As discussed in the introduction, the phase bifurcation is the splitting of a single ZPA frequency into three. Specifically, the ZPA frequencies are described by condition: The formulas describing the phase bifurcation frequencies can be obtained from (4). To simplify the calculation, the normalization process presented by Wang et al. in [22] can be advantageously used. In the normalization, the ESRs R p and R s are neglected. The IPT model is normalized as follows. First, the normalized operating frequency u is defined as the operating frequency f divided by f 0 : Second, the imaginary part of the input impedance Im Ẑ is divided by the real part of the reflected impedance at resonance Re Ẑ r0 to obtain normalized Im Ẑ n : A series of substitutions is then applied after which Im Ẑ is expressed solely in terms of u and the loaded quality factors Q Lp and Q Ls , which are defined as: In the last step, Equation (4) is transformed into the phase bifurcation equation by substituting Im Ẑ n for Im Ẑ : The detailed derivation of the phase bifurcation equation is described in [22]. When solved for u, Equation (9) has six roots, three positive ones and three negative ones. The negative roots are neglected as u cannot be negative. The first root u ϕ0 = 1 describes the resonance frequency f 0 , and second and third roots describe the ZPA frequencies u ϕ1 and u ϕ2 , respectively: and when the roots u ϕ1 and u ϕ2 become real, the phase bifurcation occurs. Analyzing Equations (10) and (11) for the ZPA frequencies u ϕ1 and u ϕ2 shows that these frequencies are asymptotic and that the asymptotes depend solely on the coupling coefficient k. The asymptotic behavior can be demonstrated by decreasing the equivalent load resistance R L towards zero. Figure 2 for IPT system described in Section 5.1. shows the course of the ZPA frequencies u ϕ0 , u ϕ1 and u ϕ2 with decreasing R L while the remaining equivalent circuit parameters remain constant. As it can be seen, the ZPA frequencies u ϕ1 and u ϕ2 begin close to the resonance frequency u 0 . As R L decreases further, u ϕ1 and u ϕ2 recede until they stabilize at the asymptotes. The deep bifurcation, i.e., very low R L , is required for the roots of the bifurcation equation approach the asymptotes closely. The easiest way to achieve this is to short-circuit the output of the IPT system. To obtain the equation for the coupling coefficient calculation, it is necessary to obtain the formulas for the asymptotes. In the first step, the loaded quality factors and given by (8) are substituted in the roots and of the bifurcation equation, which are given by (10) and (11), respectively. Then their asymptotic values are calculated in the limit as goes to zero. In the next step the normalized frequency given by (5) is substituted together with the coupling coefficient : The final results are the equations for the asymptotes and : This result is interesting, because the asymptotes of the frequencies describing the frequency splitting of the output power maxima can be approximated by a similar set of formulas [23,24]. Our examination of the frequencies of the input current amplitude maxima showed that they also have the same asymptotes. The coupling coefficient can be calculated from the asymptotes and given by (13) and (14), respectively, as: (15) Note that the issue of coupling coefficient extraction is also important in distributedelement filter design; and, as discussed in [25], an equation identical to (16) is used to characterize coupled RF/microwave resonators based on the splitting of the scattering parameter .

Proposed Estimation Method (Bifurcation Asymptotes Method)
Equation (15) can be advantageously used for estimation of the coupling coefficient . Thus, we propose the following procedure based on the bifurcation asymptotes, which estimates using the frequency response of the input phase . The bifurcation asymptotes method has the following steps: (1) Short-circuit the output of the IPT system.
(2) Find the approximate values of the frequencies and by searching for either zero crossings of or input current amplitude maxima with a large step of . To obtain the equation for the coupling coefficient k calculation, it is necessary to obtain the formulas for the asymptotes. In the first step, the loaded quality factors Q Lp and Q Ls given by (8) are substituted in the roots u ϕ1 and u ϕ2 of the bifurcation equation, which are given by (10) and (11), respectively. Then their asymptotic values are calculated in the limit as R L goes to zero. In the next step the normalized frequency u given by (5) is substituted together with the coupling coefficient k: The final results are the equations for the asymptotes f 1a and f 2a : This result is interesting, because the asymptotes of the frequencies describing the frequency splitting of the output power P out maxima can be approximated by a similar set of formulas [23,24]. Our examination of the frequencies of the input current amplitude I in maxima showed that they also have the same asymptotes.
The coupling coefficient k can be calculated from the asymptotes f 1a and f 2a given by (13) and (14), respectively, as: Note that the issue of coupling coefficient extraction is also important in distributedelement filter design; and, as discussed in [25], an equation identical to (16) is used to characterize coupled RF/microwave resonators based on the splitting of the scattering parameter S 21 .

Proposed Estimation Method (Bifurcation Asymptotes Method)
Equation (15) can be advantageously used for estimation of the coupling coefficient k. Thus, we propose the following procedure based on the bifurcation asymptotes, which estimates k using the frequency response of the input phase ϕ in .
The bifurcation asymptotes method has the following steps: (1) Short-circuit the output of the IPT system.  The individual steps are examined in further detail in the following discussion. In order to use (15), it is necessary to achieve such values of Q Lp and Q Ls that f ϕ1 and f ϕ2 would approach f 1a and f 2a , respectively. While there are multiple ways to achieve such a change of Q Lp and Q Ls , e.g., by decreasing f 0 , minimizing of R L by short-circuiting the output is the most practical.
Because the method is based on measuring the frequency response, a knowledge of a frequency range to examine is necessary. However, this range depends on k based on (13) and (14), which is unknown. We find that it is most convenient is to start at f 0 , and then extend the search for f 1a and f 2a both above and below f 0 . Therefore, an approximate value of f 0 is the only information about the IPT system required for the method. An IPT system with f 0 of 100 kHz was used for the method verification.
Measurement of the full frequency response with a small step of frequency is time consuming. This could render the method unusable for more dynamic applications. However, because we are only looking for two frequency values (i.e., f 1a and f 2a ), it is also unnecessary. To reduce the number of steps, two measures were taken.
At first, a rough approximation of frequencies f 1a and f 2a is generated by employing a large measurement step-in testing of the method a step of 1 kHz was used. In this search, the positions of f 1a and f 2a in the frequency range can be evaluated based on multiple criteria. For example, if the evaluation is based on the input phase ϕ in , the positions of f 1a and f 2a can be localized based on the change of sign of ϕ in . Frequencies f 1a and f 2a can be also approximated by searching for the frequencies corresponding to maxima of the input current amplitude I in , as these have the same asymptotes as the ZPA frequencies f ϕ1 and f ϕ2 . Alternatively, if the primary side DC link current is measured, the frequencies of its maxima can be used as well, as they match I in maxima. In the verification of the presented method, the frequencies of I in maxima were used to generate a first approximation of f 1a and f 2a .
The second measure is to use the linear approximation to find f 1a and f 2a with increased precision. Figure 3 shows an example of the frequency response of the input phase ϕ in for deep bifurcation (i.e., f ϕ1 ≈ f 1a and f ϕ2 ≈ f 1a ) with marked f 0 , f 1a and f 2a . As the figure shows, the shape of the frequency response ϕ in is close to linear near the ZPA frequencies f ϕ1 and f ϕ2 .
Thus, the values of ϕ in are measured near the approximated values of f 1a and f 2a obtained in the first step with a smaller frequency step. In the verification of the method, we measured the ±1 kHz surroundings of the approximated frequencies f 1a and f 2a with a step of 0.5 kHz. From the resulting data series, the frequencies and phases of the points nearest the zero crossing of ϕ in were selected for each asymptote: the point f ia−1 , ϕ ia−1 with phase below zero and the point f ia+1 , ϕ ia+1 with phase above zero. Both points are marked by hollow circles in Figure 3. From these two points, the values of asymptote frequencies are approximated by: The obtained values of f 1a and f 2a (full circles in Figure 3) are then used to calculate the coupling coefficient according to (15). Electronics 2021, 10, 2548 7 of 20 Figure 3. Frequency response of the input phase for the IPT system in deep bifurcation, with marked measured points, linear approximation and interpolated zero crossings. Based on the circuit parameters from Table 1.
Thus, the values of are measured near the approximated values of and obtained in the first step with a smaller frequency step. In the verification of the method, we measured the ±1 kHz surroundings of the approximated frequencies and with a step of 0.5 kHz. From the resulting data series, the frequencies and phases of the points nearest the zero crossing of were selected for each asymptote: the point , with phase below zero and the point , with phase above zero. Both points are marked by hollow circles in Figure 3. From these two points, the values of asymptote frequencies are approximated by: (16) The obtained values of and (full circles in Figure 3) are then used to calculate the coupling coefficient according to (15).
The number of measured frequency points depends on . For example, if 0.3 is measured and 100 kHz, then 87.70 kHz and 119.52 kHz according to (13) and (14), respectively. The first approximation of and with a step of 1 kHz then requires measurement of 34 points and the second refining for the linear approximation requires 10 steps if the ±1 kHz surroundings of frequencies and are measured with a step of 0.5 kHz. In total, 44 points are measured. However, if the method is used to estimate increasing from a small value (e.g., parking of an EV over a charging station), the number of points can be further reduced by using and estimated from the previous step in the first approximation instead of .
Because the asymptotes (13) and (14) are the same for phase bifurcation, frequency splitting of the output power maxima and splitting of input current amplitude maxima [23,24], alternative approaches to finding the zero crossings of can be used for asymptote approximation, e.g., finding the frequencies of the maxima of or . However, the linear approximation cannot be used to estimate the positions of the maxima; thus, the step of the frequency response would need to be much smaller, resulting in increased measurement time. The number of measured frequency points depends on k. For example, if k = 0.3 is measured and f 0 = 100 kHz, then f 1a = 87.70 kHz and f 2a = 119.52 kHz according to (13) and (14), respectively. The first approximation of f 1a and f 2a with a step of 1 kHz then requires measurement of 34 points and the second refining for the linear approximation requires 10 steps if the ±1 kHz surroundings of frequencies f 1a and f 2a are measured with a step of 0.5 kHz. In total, 44 points are measured. However, if the method is used to estimate k increasing from a small value (e.g., parking of an EV over a charging station), the number of points can be further reduced by using f 1a and f 2a estimated from the previous step in the first approximation instead of f 0 .
Because the asymptotes f 1a (13) and f 2a (14) are the same for phase bifurcation, frequency splitting of the output power P out maxima and splitting of input current amplitude I in maxima [23,24], alternative approaches to finding the zero crossings of ϕ in can be used for asymptote approximation, e.g., finding the frequencies of the maxima of I in or I out . However, the linear approximation cannot be used to estimate the positions of the maxima; thus, the step of the frequency response would need to be much smaller, resulting in increased measurement time.

Experimental Verification
The estimates of k provided by the bifurcation asymptotes method are evaluated with respect to a reference value measured with a calibrated LRC bridge and compared with the results of other active (based on the measurement of currents and voltages) estimation methods. A direct method based on Faraday's law (i.e., measuring the induced electromotive force (EMF)) and the method introduced by Jiwariyavej et al. in [18] are used for comparison.
In the first part of the chapter, the measurement setup is described together with the measured sequence of coupling coefficients, which is used for evaluation of all three estimation methods. Then, the procedures of the Faraday's law method and Jiwariyavej's method are outlined, and the estimation results are compared with the reference for all three selected methods. In the last part of the chapter, the match to the reference of all three methods is evaluated.

Measurement Setup
The measured setup used for verification is depicted in Figure 4 (photographs) and in Figure 5 (schematics). The setup was designed to closely resemble the AC equivalent circuit in Figure 1.
The circuit is powered by a linear amplifier driven by the function generator BK Precision 4063 which generates a sinusoidal voltage. The linear amplifier is supplied by power sources Diametral P230R51D. Both the primary and the secondary use DD pads [15,26]. Each coil inductance is compensated by a capacitive decade to provide control over f 0p and f 0s . Nominal resonance frequency was selected as 100 kHz.
The output is loaded by a custom-built resistive decade, designed for high frequency operation. However, for lower values of R L it has a slight inductive character, but this is considered in L s and the calculation of C s , which compensates it. The currents are calculated from the measured voltages at current shunts. One of two available oscilloscopes Rohde & Schwarz RTH1004 measures the voltages with RT-ZI10 probes and the shunt voltages corresponding with currents with P2220 probes set to 1:1 attenuation. The voltages are measured at designated points at the switchboards, which connect together all components of each side. The function generator and oscilloscope are controlled by an acquisition MATLAB script running on the laptop. The setup was designed with a focus on the IPT transformer operation and the bifurcation phenomena; thus, the linear amplifier and the resistor decade are used in place of an inverter and output rectifier, respectively. The component parameter values were measured with the LRC bridge Hameg HM8118 with Kelvin probe. The system parameters for the aligned coils (position x = 0 mm) are listed Table 1, except R L , which is specific for each method.

Experimental Verification
The estimates of provided by the bifurcation asymptotes method are evaluated with respect to a reference value measured with a calibrated LRC bridge and compared with the results of other active (based on the measurement of currents and voltages) estimation methods. A direct method based on Faraday's law (i.e., measuring the induced electromotive force (EMF)) and the method introduced by Jiwariyavej et al. in [18] are used for comparison.
In the first part of the chapter, the measurement setup is described together with the measured sequence of coupling coefficients, which is used for evaluation of all three estimation methods. Then, the procedures of the Faraday's law method and Jiwariyavej's method are outlined, and the estimation results are compared with the reference for all three selected methods. In the last part of the chapter, the match to the reference of all three methods is evaluated.

Measurement Setup
The measured setup used for verification is depicted in Figure 4 (photographs) and in Figure 5 (schematics). The setup was designed to closely resemble the AC equivalent circuit in Figure 1.   The circuit is powered by a linear amplifier driven by the function generator BK Precision 4063 which generates a sinusoidal voltage. The linear amplifier is supplied by power sources Diametral P230R51D. Both the primary and the secondary use DD pads [15,26]. Each coil inductance is compensated by a capacitive decade to provide control over and . Nominal resonance frequency was selected as 100 kHz. The output is loaded by a custom-built resistive decade, designed for high frequency operation. However, for lower values of it has a slight inductive character, but this is considered in and the calculation of , which compensates it. The currents are calculated from the measured voltages at current shunts. One of two available oscilloscopes Rohde & Schwarz RTH1004 measures the voltages with RT-ZI10 probes and the shunt voltages corresponding with currents with P2220 probes set to 1:1 attenuation. The voltages are measured at designated points at the switchboards, which connect together all components of each side. The function generator and oscilloscope are controlled by an acquisition MATLAB script running on the laptop. The setup was designed with a focus   The circuit is powered by a linear amplifier driven by the function generator BK Precision 4063 which generates a sinusoidal voltage. The linear amplifier is supplied by power sources Diametral P230R51D. Both the primary and the secondary use DD pads [15,26]. Each coil inductance is compensated by a capacitive decade to provide control over and . Nominal resonance frequency was selected as 100 kHz. The output is loaded by a custom-built resistive decade, designed for high frequency operation. However, for lower values of it has a slight inductive character, but this is considered in and the calculation of , which compensates it. The currents are calculated from the measured voltages at current shunts. One of two available oscilloscopes Rohde & Schwarz RTH1004 measures the voltages with RT-ZI10 probes and the shunt voltages corresponding with currents with P2220 probes set to 1:1 attenuation. The voltages are measured at designated points at the switchboards, which connect together all components of each side. The function generator and oscilloscope are controlled by an acquisition MATLAB script running on the laptop. The setup was designed with a focus

Measured Sequence Used for Method Evaluation
The presented method is evaluated together with selected existing methods by estimating the coupling coefficient k during the secondary pad displacement. The secondary pad is horizontally moved in the x-axis (see Figure 6) in 10 mm steps from the initial position of 0 to 200 mm. The step is decreased to 5 mm near the DD pads null position (see Figure 9c in [15]).

Measured Sequence Used for Method Evaluation
The presented method is evaluated together with selected existing methods by estimating the coupling coefficient during the secondary pad displacement. The secondary pad is horizontally moved in the x-axis (see Figure 6) in 10 mm steps from the initial position of 0 to 200 mm. The step is decreased to 5 mm near the DD pads null position (see Figure 9c in [15]). The coupling coefficient values used as a reference for method evaluation were measured by the following procedure using the calibrated LRC bridge.
(1) Both pads are disconnected from the rest of the circuit.  [10]. As Figure 7b shows, the movement of the secondary coil changes the coil inductances, which in turn detunes the primary and secondary circuits. The impacts of detuning will be discussed for each method individually. The coupling coefficient values used as a reference for method evaluation were measured by the following procedure using the calibrated LRC bridge.
(1) Both pads are disconnected from the rest of the circuit.
(2) The primary coil inductance L 1a is measured while the secondary coil is open.
(3) The secondary coil is short circuited. The primary inductance is measured again and labeled as L 1b . L 1b is related to L 1a via: (4) The coupling coefficient k is then calculated from (17) as: Figure 7 depicts the resulting k values (a) together with the primary and secondary inductances L 1 and L 2 (b). The figure clearly shows the null point of the DD pads near x = 75 mm. The resulting sequence of k covers the interval of 0.08 to 0.35, which is an aggregated range of k, which the IPT systems presented in SAE J2954 as examples for EV stationary charging gain under displacement [10]. As Figure 7b shows, the movement of the secondary coil changes the coil inductances, which in turn detunes the primary and secondary circuits. The impacts of detuning will be discussed for each method individually.
The methods for estimating k are compared with the reference measurement presented in Figure 7a. For each displacement position, the relative error with respect to the reference value is calculated according to: (19) where k M is the value estimated by the method and k re f is the reference value. The methods for estimating are compared with the reference measurement presented in Figure 7a. For each displacement position, the relative error with respect to the reference value is calculated according to: where is the value estimated by the method and is the reference value.

Bifurcation Asymptotes Method
The proposed method for estimating was implemented and tested according to the procedure outlined in Chapter 4. The method employed the setup shown in Figure 5, only the load was short-circuited at its connection point at the secondary switchboard. For each displacement position, the input current amplitude was measured with a step size of 1 kHz to find the approximate values of and . These were then refined by measurement of the input phase measured with a step size of 0.5 kHz in interval ±1 kHz around the previously approximated frequencies and , and consequent linear interpolation. The obtained estimates of "BA 0 mm" and "BA AT" are compared in Figure 8a with the reference "ref". Label "BA 0 mm" marks the measurement with the system tuned to 100 kHz in the initial position x = 0 mm and "BA AT" refers to the measurement when the system was manually tuned by adjusting the primary and secondary capacitances for each misalignment position to 100 kHz (AT stands for always tuned).

Bifurcation Asymptotes Method
The proposed method for estimating k was implemented and tested according to the procedure outlined in Section 4. The method employed the setup shown in Figure 5, only the load was short-circuited at its connection point at the secondary switchboard. For each displacement position, the input current amplitude I in was measured with a step size of 1 kHz to find the approximate values of f 1a and f 2a . These were then refined by measurement of the input phase ϕ in measured with a step size of 0.5 kHz in interval ±1 kHz around the previously approximated frequencies f 1a and f 2a , and consequent linear interpolation. The obtained estimates of k "BA 0 mm" and "BA AT" are compared in Figure 8a with the reference "ref". Label "BA 0 mm" marks the measurement with the system tuned to f 0 = 100 kHz in the initial position x = 0 mm and "BA AT" refers to the measurement when the system was manually tuned by adjusting the primary and secondary capacitances for each misalignment position to f 0 = 100 kHz (AT stands for always tuned). The methods for estimating are compared with the reference measurement presented in Figure 7a. For each displacement position, the relative error with respect to the reference value is calculated according to: where is the value estimated by the method and is the reference value.

Bifurcation Asymptotes Method
The proposed method for estimating was implemented and tested according to the procedure outlined in Chapter 4. The method employed the setup shown in Figure 5, only the load was short-circuited at its connection point at the secondary switchboard. For each displacement position, the input current amplitude was measured with a step size of 1 kHz to find the approximate values of and . These were then refined by measurement of the input phase measured with a step size of 0.5 kHz in interval ±1 kHz around the previously approximated frequencies and , and consequent linear interpolation. The obtained estimates of "BA 0 mm" and "BA AT" are compared in Figure 8a with the reference "ref". Label "BA 0 mm" marks the measurement with the system tuned to 100 kHz in the initial position x = 0 mm and "BA AT" refers to the measurement when the system was manually tuned by adjusting the primary and secondary capacitances for each misalignment position to 100 kHz (AT stands for always tuned).  Relative error with respect to the reference over k. "BA 0 mm" refers to system tuned for x = 0 mm. "BA AT" refers to system tuned for each misalignment position.
For the system tuned at x = 0 mm, the method works well until approximately k = 0.17, as it has a relative error of approximately −2% (see Figure 8b), which is calculated according to (19). The negative sign means that the method underestimates k in comparison with the reference.
As k decreases less than k = 0.17, the relative error becomes positive, and its value increases up to 10%. For k smaller than approximately 0.08, the method fails altogether. However, for most practical systems [10], k is larger than 0.08 (typically, k > 0.2 is expected for high efficiency) so the method's breakdown at low values of k is not expected to pose a limitation in practice.
The increase of the relative error for k below 0.17 and the method's failure for k lower than 0.08 is primarily caused by detuning. Figure 9 shows that with increasing displacement the primary resonance frequency f 0p increases above the secondary resonance frequency f 0s . This detuning alters the frequency response of ϕ in . Because f 0p > f 0s , the frequency response moves below zero in comparison with the case f 0p = f 0s -see Figure 10a for the tuned system and Figure 10b for the detuned system. according to (19). The negative sign means that the method underestimates in comparison with the reference.
As decreases less than 0.17, the relative error becomes positive, and its value increases up to 10%. For smaller than approximately 0.08, the method fails altogether. However, for most practical systems [10], k is larger than 0.08 (typically, k > 0.2 is expected for high efficiency) so the method's breakdown at low values of k is not expected to pose a limitation in practice.
The increase of the relative error for below 0.17 and the method's failure for lower than 0.08 is primarily caused by detuning. Figure 9 shows that with increasing displacement the primary resonance frequency increases above the secondary resonance frequency . This detuning alters the frequency response of . Because , the frequency response moves below zero in comparison with the case -see Figure  10a for the tuned system and Figure 10b for the detuned system.  Under the displacement, the frequency responses change both with detuning and with changing , as illustrated in Figure 11. In deep bifurcation, the frequency response of is relatively sharp near the zero crossings and it approaches its limiting values of −90° and +90° between the crossings if is sufficiently high (see Figure 11a,b). Thus, the impacts of detuning are smaller for high , as even with the shift the frequency response according to (19). The negative sign means that the method underestimates in comparison with the reference.
As decreases less than 0.17, the relative error becomes positive, and its value increases up to 10%. For smaller than approximately 0.08, the method fails altogether. However, for most practical systems [10], k is larger than 0.08 (typically, k > 0.2 is expected for high efficiency) so the method's breakdown at low values of k is not expected to pose a limitation in practice.
The increase of the relative error for below 0.17 and the method's failure for lower than 0.08 is primarily caused by detuning. Figure 9 shows that with increasing displacement the primary resonance frequency increases above the secondary resonance frequency . This detuning alters the frequency response of . Because , the frequency response moves below zero in comparison with the case -see Figure  10a for the tuned system and Figure 10b for the detuned system.  Under the displacement, the frequency responses change both with detuning and with changing , as illustrated in Figure 11. In deep bifurcation, the frequency response of is relatively sharp near the zero crossings and it approaches its limiting values of −90° and +90° between the crossings if is sufficiently high (see Figure 11a,b). Thus, the impacts of detuning are smaller for high , as even with the shift the frequency response Under the displacement, the frequency responses change both with detuning and with changing k, as illustrated in Figure 11. In deep bifurcation, the frequency response of ϕ in is relatively sharp near the zero crossings and it approaches its limiting values of −90 • and +90 • between the crossings if k is sufficiently high (see Figure 11a,b). Thus, the impacts of detuning are smaller for high k, as even with the shift the frequency response resembles a straight line at the zero crossing and can be approximated as a linear function there (see Figure 11b). However, as k decreases due to displacement, the frequency response curve becomes less pronounced, and it begins to cross zero near its peak, where it does not resemble a straight line (see Figure 11c). This increases the relative error. With a further decrease of k and increase of detuning, former zero crossings are shifted below zero and only a single zero crossing occurs- Figure 11d. Thus, for low k and increased detuning, the method fails. resembles a straight line at the zero crossing and can be approximated as a linear function there (see Figure 11b). However, as decreases due to displacement, the frequency response curve becomes less pronounced, and it begins to cross zero near its peak, where it does not resemble a straight line (see Figure 11c). This increases the relative error. With a further decrease of and increase of detuning, former zero crossings are shifted below zero and only a single zero crossing occurs- Figure 11d. Thus, for low and increased detuning, the method fails. The method's performance for low can be improved by maintaining the tuned system for specific inductance values under displacement by adjusting the primary and secondary compensation capacitances. The compensation adjustment can be done by e.g., switched capacity, or manually, as in our case. By removing the detuning, it was possible to obtain an approximate value of in positions 65, 70, 180, and 190 mm-see "BA AT" in Figure 8a. However, even though as decreases the measurement error grows up to 30%, as shown in Figure 8b. The high error is caused by the fact that the bifurcation is not deep enough such that frequencies and are not adequately meeting their asymptotes and , respectively. After decreases under 0.04, not even tuning for a specific displacement position could help to achieve any meaningful results.
In the remaining points (75, 80, and 200 mm), only one zero crossing of the input phase was measured, even when the measurement step was decreased to 0.25 kHz. This may be caused either because the lack of bifurcation or the measurement step is still too coarse. The method's performance for low k can be improved by maintaining the tuned system for specific inductance values under displacement by adjusting the primary and secondary compensation capacitances. The compensation adjustment can be done by e.g., switched capacity, or manually, as in our case. By removing the detuning, it was possible to obtain an approximate value of k in positions 65, 70, 180, and 190 mm-see "BA AT" in Figure 8a. However, even though as k decreases the measurement error grows up to 30%, as shown in Figure 8b. The high error is caused by the fact that the bifurcation is not deep enough such that frequencies f ϕ1 and f ϕ2 are not adequately meeting their asymptotes f 1a and f 2a , respectively. After k decreases under 0.04, not even tuning for a specific displacement position could help to achieve any meaningful results.
In the remaining points (75, 80, and 200 mm), only one zero crossing of the input phase was measured, even when the measurement step was decreased to 0.25 kHz. This may be caused either because the lack of bifurcation or the measurement step is still too coarse.
In applying the bifurcation asymptotes method, it is necessary to consider the impacts of the bifurcation phenomena on the input impedanceẐ. When short-circuiting the output, the input impedance is relatively high at f 0 (e.g., Z = 99.3 Ω for x = 0 mm). However, near the ZPA frequencies f ϕ1 and f ϕ2 the impedance is very low (e.g., Z = 0.73 Ω), which is determined by the ESRs R p and R s . Thus, to limit the currents to 5 to 7 A in the primary side corresponding with the regular operation, the input voltage V in was decreased to 5 V when applying the method.

Faraday's Law Method
For comparison, we consider another common approach for extracting k by directly measuring the induced EMF in the secondary coil and solving for M using Faraday's law of electromagnetic induction. The adjusted measurement setup is depicted in Figure 12. Both pads are placed in their positions. The secondary coil is open, and the primary coil is in series with a resistor and a capacitor, which provide control of the coil current I 1 . The primary coil current was limited between 5.8 and 6.6 A during the measurement process. output, the input impedance is relatively high at (e.g., 99. 3 Ω for x = 0 mm).
However, near the ZPA frequencies and the impedance is very low (e.g., 0.73 Ω), which is determined by the ESRs and . Thus, to limit the currents to 5 to 7 A in the primary side corresponding with the regular operation, the input voltage was decreased to 5 V when applying the method.

Faraday's Law Method
For comparison, we consider another common approach for extracting by directly measuring the induced EMF in the secondary coil and solving for using Faraday's law of electromagnetic induction. The adjusted measurement setup is depicted in Figure 12. Both pads are placed in their positions. The secondary coil is open, and the primary coil is in series with a resistor and a capacitor, which provide control of the coil current . The primary coil current was limited between 5.8 and 6.6 A during the measurement process. The primary coil is excited with current , which is measured with a precision current shunt. The value is controlled with a resistor and a capacitor, which reduces the inductive character. The value of is adjusted to equal the value of during regular operation, to consider ferrite saturation and similar effects. The voltage is measured across the secondary coil terminals, and the mutual inductance is calculated as: The coupling coefficient is calculated from and the measured coil inductances , according to (12). The estimated values of obtained by Faraday's law method are compared with the reference measurements in Figure 13a. The relative error depicted in Figure 13b is very low even for low values. The distribution of error shows that the method generally estimates approximately by 1.55% above the reference. Faraday's law method is unaffected by detuning. The primary coil is excited with current I 1 , which is measured with a precision current shunt. The I 1 value is controlled with a resistor and a capacitor, which reduces the inductive character. The value of I 1 is adjusted to equal the value of I in during regular operation, to consider ferrite saturation and similar effects. The voltage V 2 is measured across the secondary coil terminals, and the mutual inductance M is calculated as: The coupling coefficient k is calculated from M and the measured coil inductances L 1 , L 2 according to (12).
The estimated values of k obtained by Faraday's law method are compared with the reference measurements in Figure 13a. The relative error depicted in Figure 13b is very low even for low k values. The distribution of error shows that the method generally estimates k approximately by 1.55% above the reference. Faraday's law method is unaffected by detuning.

Method of Jiwariyavej
The second method used for comparison was presented by Jiwariyavej et al. in [18]. The variant for a single receiver was used. This method does not directly measure but the mutual inductance . The calculation of is based on the equation for the real part of input impedance: (21) where is the secondary side reactance: 1 (22)

Method of Jiwariyavej
The second method used for comparison was presented by Jiwariyavej et al. in [18]. The variant for a single receiver was used. This method does not directly measure k but the where X s is the secondary side reactance: The mutual inductance M is then calculated as: The real part of input impedance Re{Z in } is calculated from the input current and voltage amplitudes I in and V in , respectively, and the input phase: The coupling coefficient k is then calculated according to (12). Obtaining the value of k according to this method requires preliminary measurement of R L , R p , R s , L p , L s , and C s .
The method was tested using the measurement setup shown in Figure 5, with the resistive load set to 31.8 Ω. Due to the higher number of preliminary parameter measurements, the method was measured in reduced number of points. The primary and secondary capacitances were adjusted to keep the system tuned at resonance frequency f 0 = 100 kHz against variations in coil inductance; see the illustration in Figure 7b. The specific values of the circuit parameters measured for each displacement position were used to calculate M according to (23) and consequently k according to (12).
Values of k obtained using the method of Jiwariyavej are compared to the reference values in Figure 14a. In addition to the plotted values, the method was also applied at displacement of 75, 190, and 200 mm, but did not provide any useful results. Thus, despite tuning for each displacement position, the method of Jiwariyavej suffers a similar limitation for small values of k as the bifurcation asymptote method, which was tuned only for x = 0 mm. The relative error of the method of Jiwariyavej is depicted in Figure 14b. For high values (above approximately 0.23), this method exhibits good agreement with the reference. On the other hand, as is decreased, the error grows significantly and reaches 19.3% at 0.08.

Measurement Comparison
A comparative evaluation of the methods for estimating k is summarized in Figure  15. Figure 15a shows how the methods compare directly to the reference and Figure 15b shows their relative errors. The relative error of the method of Jiwariyavej is depicted in Figure 14b. For high k values (above approximately 0.23), this method exhibits good agreement with the reference. On the other hand, as k is decreased, the error grows significantly and reaches 19.3% at k = 0.08.

Measurement Comparison
A comparative evaluation of the methods for estimating k is summarized in Figure 15. Figure 15a shows how the methods compare directly to the reference and Figure 15b shows their relative errors. The relative error of the method of Jiwariyavej is depicted in Figure 14b. For high values (above approximately 0.23), this method exhibits good agreement with the reference. On the other hand, as is decreased, the error grows significantly and reaches 19.3% at 0.08.

Measurement Comparison
A comparative evaluation of the methods for estimating k is summarized in Figure  15. Figure 15a shows how the methods compare directly to the reference and Figure 15b shows their relative errors. The evaluation of the methods is divided in two intervals ∈ ⟨0.08, 0.23 and ∈ 〈 0.23,0.36〉 based on the relative error results shown in Figure 15b. The bottom limit of 0.08 was selected as the final value, at which point all three methods yielded a result.
Thus, only the bifurcation asymptote method with tuning at x = 0 mm was evaluated. For each interval, the method is evaluated based on the average error ∆ and average absolute error ∆ , which are calculated as: The evaluation of the methods is divided in two intervals k ∈ 0.08, 0.23 and k ∈ 0.23, 0.36 based on the relative error results shown in Figure 15b. The bottom limit of k = 0.08 was selected as the final value, at which point all three methods yielded a result. Thus, only the bifurcation asymptote method with tuning at x = 0 mm was evaluated. For each interval, the method is evaluated based on the average error ∆k r and average absolute error ∆k r , which are calculated as: while ∆k r quantifies the disagreement between the method's results from the reference in the selected interval. ∆k r describes the shift of the estimate with respect to the reference, i.e., if ∆k r → 0 then the values are distributed evenly around the reference. The results are summarized in Table 2. All three estimation methods show a good match to the reference for k above 0.23. From these values it is evident that the Jiwariyavej's method exhibits the best match to the reference ( ∆k r = 0.20%), with the estimates which are distributed evenly around the reference (∆k r = −0.04%). The estimates of the method based on Faraday's law are shifted above the reference (∆k r = 1.91%), while the estimates of the bifurcation asymptotes method are shifted below the reference (∆k r = −2.01%).
However, with decreasing k all methods except Faraday's law method decrease in precision (interval k ∈ 0.08, 0.23 ). Faraday's law method is the most accurate ( ∆k r = 1.39%), but its estimates are again shifted above the reference (∆k r = 1.39%).
The second-best method is the bifurcation asymptotes method ( ∆k r = 4.36%) which is now shifted above the reference (∆k r = 3.87%). The Jiwariyavej's method performs much worse ( ∆k r = 10.38%) and its value begins to move further below the reference with decreasing k-see Figure 15b, ∆k r = −9.30%.
Across the whole measured interval of k ∈ 0.08, 0.36 , Faraday's law method is the most accurate ( ∆k r = 1.55%), followed by the bifurcation asymptote method ( ∆k r = 3.62%) and then Jiwariyavej's method ( ∆k r = 6.31%). The estimates of the two first methods are on average shifted above the reference ( ∆k r = 1.55% and ∆k r = 3.62%, respectively), while the third is shifted below the reference (∆k r = −5.6%).
The shifts of the estimates with respect to the reference could be partially explained by the different measurement currents used in each method (see Table 3) and the reference, which was measured with the LRC meter and thus with very small currents. The high currents could result in detuning caused by the dependance of the coil inductances on the currents through the coils and ferrite non-linearity [27], as the ferrites are used as the magnetic cores of the DD pads. This effect could especially impact the Jiwariyavej's method, as the primary currents were significantly higher in the second interval. Table 3. Currents used for testing of the estimation methods.

I in I out
Bifurcation asymptote method 5 to 7 A 5 to 10 A Faraday's law method 5 to 7 A Open Jiwariyavej's method 4.5 to 6 A (k ∈ 0.23, 0.36 ) 10.5 A (k ∈ 0.08, 0.23 ) 0.5 to 1 A (k ∈ 0.23, 0.36 ) 0.5 A (k ∈ 0.08, 0.23 ) Some of the differences could also be explained as a result of measurement errors due to limited instrument accuracy, loading effects, and environmental factors such as temperature, etc.

Discussion
As evaluated in the previous chapter, the bifurcation asymptote method provides relatively accurate estimates ( ∆k r = 3.62%) across the interval of k ∈ 0.08, 0.36 . However, the estimate precision is not the only consideration for evaluating the usability of the method. Table 4 summarizes some of the practical aspects of all three evaluated methods, which should also be considered. None to obtain M, L p , L s to obtain k R L , R p , R s , L s , and C s to obtain M, plus L p to obtain k

Other requirements
Ability to short-circuit the load. Source with variable frequency Decreased source voltage Open secondary coil Air-core coils, or L p and L s must be remeasured after displacement While the method based on Faraday's law provides very accurate estimates, it requires the ability to disconnect the secondary coil from the rest of the circuit and measure the coil voltage across its open terminals. Thus, this method is more suitable for characterizing the coils in the laboratory setting than performing real-time estimations of k for practical applications such as EV charging.
Jiwariyavej's method for a single receiver also provides accurate estimates, but only for higher k values: k ∈ 0.08, 0.36 . Another disadvantage is that it requires additional preliminary measurements including L s to calculate the value of M. Due to this, use of Jiwariyavej's method during operation is limited to air-core coils in which the inductance does not vary with displacement. The method is not very useful for the coils with magnetic cores, as their inductances change with displacement.
The main downside of the bifurcation asymptotes method is that it requires shortcircuiting of the load; therefore, it cannot provide estimation of k during the operation of an IPT system. However, the method can be used to estimate k intermittently during a charging cycle and at critical times such as at start-up. In addition, the method involves measuring the frequency response; thus, a variable frequency source is required, and a certain amount of time given by transients must be allotted for the extraction process to occur.
In return, however, the bifurcation asymptotes method provides significant advantages. Namely, the method requires no a-priori information of the system component values, it is based purely on primary-side measurements, it estimates k with good accuracy ( ∆k r = 3.62% for k ∈ 0.08, 0.36 ), and it is resilient to the detuning caused by the displacement of the coils with magnetic cores. Thanks to this, the bifurcation asymptotes method presents an interesting option for estimating k that can be readily employed by control strategies that require an accurate estimation of this variable. The method is especially suitable for implementation at the startup of stationary charging applications. For example, the method can be used to estimate k immediately after an EV is parked over a stationary wireless charger to determine if the EV is suitably positioned to allow for efficient wireless charging before high power is applied to the primary.

Conclusions
This paper presented a novel bifurcation asymptotes method for accurately estimating the coupling coefficient k for SS compensated IPT systems. The method involves temporarily operating the system in bifurcation by short-circuiting the load and measuring the zero phase angle frequencies at the primary. The method can be used by control algorithms that require the value of k, and is especially suitable for implementation at the startup of certain applications such as EV stationary charging. This paper presents a preliminary evaluation of the method, which contains the theoretical background and description of the method, experimental verification, and a comparison with existing k-estimation methods.
The proposed method provides good estimates (average absolute error of 3.62%) for a wide range of coupling coefficients (k ∈ 0.08, 0.36 ).
For future work, we intend to evaluate the method performance in a non-idealized setup by implementing it in an inverter. We also intend to explore how this method can be extended to other compensation topologies such as series-parallel, LCC-series, etc. Finally, we intend to perform more comparisons with current k-estimation methods (e.g., [18,19] cited in the Introduction).