A Nonlinear Circuit Analysis Technique for Time-Variant Inductor Systems

Time-variant inductors exist in many industrial applications, including sensors and actuators. In some applications, this characteristic can be deleterious, for example, resulting in inductive loss through eddy currents in motors designed for high efficiency operation. Therefore, it is important to investigate the electrical dynamics of systems with time-variant inductors. However, circuit analysis with time-variant inductors is nonlinear, resulting in difficulties in obtaining a closed form solution. Typical numerical algorithms used to solve the nonlinear differential equations are time consuming and require powerful processors. This investigation proposes a nonlinear method to analyze a system model consisting of the time-variant inductor with a constraint that the circuit is powered by DC sources and the derivative of the inductor is known. In this method, the Norton equivalent circuit with the time-variant inductor is realized first. Then, an iterative solution using a small signal theorem is employed to obtain an approximate closed form solution. As a case study, a variable inductor, with a time-variant part stimulated by a sinusoidal mechanical excitation, is analyzed using this approach. Compared to conventional nonlinear differential equation solvers, this proposed solution shows both improved computation efficiency and numerical robustness. The results demonstrate that the proposed analysis method can achieve high accuracy.


Introduction
Time-variant inductors exist in many industrial applications such as electric motors [1], power electronics [2], magnetic bearings [3], Linear Variable Differential Transformers (LVDT) [4], piezoelectric actuators [5] and fluidic valves [6]. The changing inductance will make the flux inside a motor nonuniform, which causes eddy currents to occur. The eddy current loss (induction loss) is a major loss in addition to Ohmic losses in the copper, hysteresis loss and mechanical loss [7][8][9], especially if a motor is running at high speed [10,11]. Recent studies also indicate that the induction loss also introduces distortions in MRI imaging systems [12]. The industry has been expected to reduce induction loss in order to not only improve the power efficiency but to also cancel the steady state ripples resulting from the eddy currents. Previously, the researchers have focused on the this topic by topology or structural optimization approaches such as using an special spinning echo [13] and novel interface circuits [14].
On the other hand, many techniques utilize the eddy currents for sensing, such as metal detection sensors [15], motor fault diagnosis sensors [16] and non-destructive sensors [17]. Some automatic electrical braking systems use eddy currents for magnetic breaking [18]. In this circumstance, it tends to enlarge the eddy current effect to get a better signal to noise ratio or power distribution by using new materials or components [19].
Apart from the hardware modification approach, an interest has arisen in characterizing its principle from the analytical approach to either do compensation or signal enhancement [20]. However, solving a nonlinear circuit with a time-variant inductor is difficult. The typical way to solve this problem is considering the time-variant inductor as time-invariant, which will yield a linear circuit for analysis. However, when the eddy currents are not neglectable, the solutions are incorrect. With modern computer processing systems, numerical methods, such as the Runge-Kutta algorithm, are used to numerically solve these ordinary differential equations. Nevertheless, in many applications, such as model predictive control with high speed motors [21,22], it is still not fast enough or economical [23]. To overcome these challenges, an optimized, fast, sufficiently precise approximate solution for typical systems with variable inductance can be an alternative approach, while academia has still left it as blank.
This investigation proposes a nonlinear technique for analyzing systems containing time-variant inductors, such as motors and solenoids, that are driven by or can be approximated as being driven by a DC source. The method is given in the following steps: (a) Identifying the circuit model of a motor and transforming it into its Norton equivalent circuit model. (b) Extracting a time-varying formula for the inductance. (c) Calculating the current through the motor using a set of iterative equations. The rest of the paper provides the details.

Background
Considering the inductor to be time-invariant, the voltage,V l , across it is: where L is the time-invariant inductor and I(t) is the current through the inductor. The prediction models do not have to perfectly match the actual system but a better model can improve performance. For a time-variant inductor in real practice, L(t), the voltage, V l , that across the inductor is: The time varying inductance makes the circuit a nonlinear system and linear circuit analysis techniques can therefore not be used to obtain a closed form solution. Linear circuit analysis can be applied by considering the inductor as a constant, where its value can be taken as the averaged value of L(t). It could be an issue if the application requires high precision when using the approximated linear analysis. Using linear circuit analysis has some disadvantages such as superpositioning cannot be used and the output frequency content differs from the input frequency content. Even though nonlinear differential equations can usually be solved using numerical techniques, they are time consuming and require high quality computing hardware and are difficult to implement, which can be key problems in predictive control. Furthermore, applications like predictive control usually require sampling and predicting the process as fast as it can, especially when dealing with a rapid response system, which conflicts with the limited calculation speed for this complicated nonlinear system. So, a simplified analysis can be used considering the case where the inductor circuit contains only resistors and DC sources, which is a typical configuration in inductive position sensors [4,24] and the steady state of a DC motor. An example circuit is shown in Figure 1. Figure 1a is the equivalent circuit of a DC motor, where V s is the steady state voltage, R 1 is the series resistor inside the motor, V l is the voltage generated by the rotor and load and L(t) is the variable inductor. It could be transformed to a Norton equivalent circuit for the sake of analysis as shown in Figure 1b. In this equivalent model, I sc and R n are the short circuit current source and the Norton equivalent resistor related to the original voltage input V s and resistor, I r is the current through the parallel resistor, R n and I l (t) is the current through the inductor L(t). The rest of the discussion is based on the standard Norton equivalent model.
The characterization of the inductors inside a motor can be done using either analytical methods [25,26] or finite element analysis (FEA) [27]. The circuit operation is described by: where V l is the voltage across the variable inductor. Though (3) fully characterizes the circuit's behavior, it is difficult to solve this nonlinear differential equation in practice and obtain a closed form solution.

Proposed Method
Because (3) is difficult to solve, an alternative iterative approach can be applied to obtain its solution [28,29]. The proposed method initiates with the following assumption that the variable inductor can be modelled as: where L 0 is the time-invariant base inductance and L 1 (t) is the time-variant part with a constraint that min(L 1 (t)) < L 0 to ensure that L(t) is always positive to obey the physical principle. It is also assumed that (4) is differentiable and can be obtained as: The first step is analysis of this circuit ignoring the effect of L 1 (t) and considering it as a linear circuit. The purpose of this step is to obtain the steady state of I l (t), denoted by I l0 (t). Then V l1 (t) can be calculated by (3) and (4). Correspondingly, the I l1 (t) term is: Additional terms can be calculated recursively: The I lk (t) shall be solved using as many terms as required to obtain sufficient precision. The overall equation of I l (t) can be solved as [28]:

Comparison with Numerical Methods
To address the advantage of the proposed method, solving the system by using the 4th order Runge-Kutta algorithm, a popular numerical method, is derived and analyzed. The first step would be reformulating (3) into a standard differential equation form: where f c (I l , L(t),L(t), V l ) represents the main body of the differential equation. The 4th order Runge-Kutta algorithm calculates the solution of I l (t) = I l (YT d ) iteratively from an initial condition is the Yth sample and T d is the computation interval. The iteration procedure can be described by: It can be observed that the recursive computation from K 1 to K 4 follows a pure serial flow, which means it cannot be accelerated by using parallel computation mechanisms. Moreover, advanced control systems for induction motors, such as model prediction control (MPC), require a behavior prediction in a finite horizon [30], which further exposes the flaw of this method. The other severe penalty of this explicit solver is that the numerical stability depends on the selection of T d , so T d has to be sufficiently small while It is a conflict to the computation efficiency. On the contrary, the proposed method can conduct an approximated but analytically accurate solution. The numerical stability is guaranteed and its computation can be optimized by a parallel form. A case study is given below to embody this advantage.

System Modeling and Analysis
Solenoids are electromagnetic actuators, which have a similar formula to inductance motors. A solenoid has a variable inductor and can be modelled as: where x is the displacement of the armature, µ 0 is the permeability of free space, µ 1 is the relative permeability of the dielectric material between the coil and armature, A is the cross-sectional area of the core, N is the number turns of the coil, x o is the initial air gap between the armature and the backside of the frame and d is the additional initial air gap related to the solenoid's geometry [31,32].
In the applications of magnetic field measurement [33] and motion control [34], the solenoid is stimulated by a external mechanical sinusoidal displacement input: where y 1 is a small magnitude and ω is the stimulating frequency. If x x 0 , the inductor's expression can be approximated as: where y 0 is the time-invariant part when x = 0: The derivative of L(t) is:L (t) = y 1 ωcos(ωt).
Considering the Norton equivalent circuit model for the solenoid circuit with a constant current source I sc , then Based on (15) and (16), apply (3) to obtain v l1 : v l1 = I sc y 1 ωcos(ωt).

Computational Complexity Analysis
The first clue that the derived formula from (18) to (21) might be complex, however, is that its computation can be significantly accelerated by using a parallel type architecture. Taking (19) as an example, the terms of −I sc y 0 y 1 ω 2 sin(ωt)/R 2 n and I sc y 2 1 ω 2 cos(2ωt)/R 2 n can be computed simultaneously. Similarly, the components in I l3 and I l4 can be divided into parallel branches to decrease the total time consumption. Though advanced sinusoidal solvers like Coordinate Rotation Digital Computer (CORDIC) have improved the computation efficiency, the time cost of each sinusoidal operation is still 1 to 2 orders of magnitude higher than addition/multiplication [35]. Thus, the number of serial stages of sinusoidal operations is dominating the total time cost in this case.
Recalling the system as the standard differential Equation (9) and reorganizing it into a standard equation, then replacing the variables with (13) and (15) in this case: f c (I l (t), L(t),L(t), V l (t)) = V l − I l y 1 ωcos(ωt) which can be used for the 4th order Runge-Kutta method to recursively solve the differential equation as shown in (10). Comparing the optimized time consumption of the proposed method and 4th order Runge-Kutta, it can be observed that the proposed method only takes 1 stage of sinusoidal operation, while 4th order Runge-Kutta 4 stages because of the recursive procedures. Another significance of this proposed method would be the in-dependency of the solutions at different times such that I l (t) for ∀t are always solvable without involving the previous states. By contrast, the Runge-Kutta method always requires an initial condition I l (0) at t = 0 and spends t/T d cycles to reach the solution of I l (t), which could potentially hurt the performance of MPC with induction motors type of applications.
It is notable that the proposed approximated methods do require more hardware resources (adders and multipliers) to get the optimized computation efficiency than the Runge-Kutta method, especially if more higher order terms are desired to improve the accuracy. However, as the embedded computation platforms like Field Programmable Gate Arrays (FPGAs) and portable General Purpose Graphic Processor Units (GPGPUs) are mushrooming, these additional costs become insignificant [36,37].

Simulation Study
To verify the feasibility of the proposed technique, a series of simulations was performed. A commercial solenoid was chosen as the target device. Its resistance was 18.7 Ω and its inductance varied between 18.6 mH and 64.8 mH depending on the position of its stroke. A DC power supply was connected with a tuned voltage of 6 V, which indicates that its Norton equivalent current source I sc was set to 0.32 A. The time-invariant y 0 was 43.2 mH while y 1 was 1.05 mH and the stimulating frequency was 20 Hz.
At the same time, the differential equation describing the system (13) with the above parameters was built in MATLAB Simulink and the iterative solution (8) with four different high order terms was obtained for comparison. In Figure 2, the left side is the system (13) and the right side is the iterative solutions with different orders. In Figure 3-5 the waveform with the caption "Simulink" presents the output of the Simulink model based on (3) and (13) by using 4th order Runge-Kutta Method with a step size T d = 10 µs, the waveforms with caption l 1 , l 1 -l 2 , l 1 -l 3 and l 1 -l 4 indicate the iterative solutions with different higher order terms. Figure 3 shows the time response of the system and the analytical solutions with up to 4th order terms. The results generated by the proposed method match the theoretical result well in the steady state when the 4th order terms are induced. Figure 4 presents the errors between the different orders configuration's analytical solutions and the simulated system. In the worst case, the steady states errors are less than ±5 mA with one term and the error decreases rapidly with more terms. Figure 5 shows the zoomed in errors in percentage scale. It is shown that the final steady errors are less than ±3% with one term and less than ±0.1% with four terms. All three figures demonstrate the increase in accuracy with additional higher order terms.
The numerical robustness is essential to the proposed method as analyzed in the analytic sections. The comparison between the proposed method with 4 terms and the 4th order Runge-Kutta method with different T d is illustrated in Figure 6. The reference signal is generated from the 4th order Runge-Kutta solver with T d = 10 µs, which is as same as the previous simulation results. It is obvious that the accuracy of the Runge-Kutta method is degraded even with a 0.2 ms step size variation and becomes unstable when it is larger than a threshold. As claimed before that the convergence of the proposed method is guaranteed, the results are performed as expected even with a large 3 ms step size.

Experimental Verification
The solenoid was fixed on a LDS800-440 large shaker, which could generate the sinusoidal input stimulation, as shown in Figure 7. The current actually through the solenoid was measured using a LEM LTS 6-NP current transducer and sampled using National Instrument's A/D board numbered 9223. A low-pass filter was applied to decrease the noise level, because the switching power supply was noisy and the current transducer was imperfect. Figure 8 shows the time response of the system under test and the analytical solutions with up to 4th order terms. The experimental results match the theory well. Figure 9 presents the errors between the different orders configuration's analytical solutions and the experimental system. In the worst case, the steady states errors are less than ±4.5 mA with one term. Figure 10 shows the zoomed in errors in percentage scale. It is shown that the final steady errors are less than ±1.5% with one term and less than ±0.5% with four terms. Compared to simulation studies, due to measurement uncertainty and noise, the experimental data is less precises than the results predicted by computer modeling but they still match the theory well and prove again that more terms can enhance the accuracy effectively.
These results demonstrate that the proposed method can be used to predict the states of a circuit with a time-variant inductor and has the potential to improve the prediction model in predictive controllers.  Though the measurement results validate the feasibility with satisfied accuracy (error ±0.5%) as expected, it can be observed that the accuracy has reached a limit and cannot be further improved by adding more higher order terms. This barrier is generated from multiple sources: (1) Impedance measurement error. The prerequisite of this method requires the accurate values of L(t) and R n but the inevitable parasitic impedance and instrumental noise will either bias or re-scale it.
(2) Motional error in the shaker system. In this experimental setup, the sinusoidal motion was provided system. The vibration system do have a closed-loop energy control system which means the overall amplitude maintains as a constant but the perfection of the sinusoidal trajectory is not guaranteed.
(3) Phase delay effects in the low-pass filters. Both the data acquisition module and post-process software have low-pass filters to reject the band noise that out-of-interest. The side effect would be introducing a phase delay in the measurement results.
Thus, the imperfections above led to residual amplitude and phase errors that can not be eliminated by adding more higher order terms. These can be decreased by optimizing the interface circuits and compensating the phase delays by adding phase-lead type of algorithms.
In summary, the proposed method can effectively solve this nonlinear circuit system with time-variant equation with correct system identification and accuracy can be improved by adding more higher order higher order terms. Its numerical robustness and computation efficiency perform better than the Runge-Kutta method. Nevertheless, this method requires more hardware resources to be implemented as a parallel architecture to fulfill its efficiency such as FPGAs or high capacity digital signal processors (DSPs) as a major limitation.

Conclusions
A nonlinear analysis technique was proposed for obtaining an approximate closed form model for time-variant inductor based systems consisting of the inductor, DC sources and resistors. First, the Norton equivalent circuit model is obtained. Then the approximate closed form solution for the current through the inductor is obtained using an iterative analysis approach. A case study with simulation and experimental verification demonstrated the technique's effectiveness in modelling the system. The theoretical study and the detailed derivation show the enhanced numerical robustness and computational efficiency. The results using the proposed technique compared to the numerical differential equation solution resulted in an error of less than 0.1% in simulation and 0.5% in experiments with a four term solution. Including additional higher order terms will further enhance the accuracy. This method is practical and straight forward to implement and can improve the applications which contain time-variant inductors.