Online Monitoring Method for Transformer Winding Deformation Based on Three-Dimensional Lissajous Curves

1  Introduction

Power transformers are fundamental components of electrical grids, primarily responsible for power transmission and voltage conversion [1]. As one of the most critical and costly pieces of equipment in power systems, their failure can severely disrupt grid reliability and operational stability [2,3]. Among various failure mechanisms, winding deformation is recognized as one of the most destructive. Studies indicate that faults attributable to winding deformation account for approximately one-third of all transformer failures [4]. In its early stages, winding deformation often exhibits no apparent symptoms; however, under continuous operational stress, the damage tends to accumulate progressively, potentially leading to irreversible deterioration of the transformer’s integrity over time [5,6]. Therefore, the development and application of effective technical methods for accurately detecting winding deformation in power transformers are essential to ensuring the secure and reliable operation of power networks.

Conventional techniques for detecting winding deformation primarily include the short-circuit impedance method [7], frequency response analysis (FRA) [8], and the low-voltage impulse (LVI) method [9]. While widely adopted in maintenance practices, these approaches share a fundamental limitation: they are offline diagnostic methods that require the transformer to be de-energized and taken out of service. Consequently, their application is restricted to periodic inspections, rendering them unsuitable for continuous condition monitoring or early fault warning. Furthermore, the diagnostic accuracy of these techniques can be compromised by factors such as measurement repeatability, inconsistent test setup, and varying environmental conditions. These factors often lead to ambiguous or inconclusive interpretations, particularly in cases of incipient or marginal faults.

To overcome the limitations of offline techniques, several online monitoring methods have been proposed in recent years, including online frequency response analysis [10], vibration-based monitoring [11], and ultrasonic detection [12]. Online FRA has demonstrated potential for detecting winding deformation; however, its diagnostic performance is often constrained by the insufficient physical interpretability of frequency-domain features and high sensitivity to operating conditions [13,14]. Vibration-based methods suffer from low signal-to-noise ratios in sub-stations and are highly susceptible to environmental noise and harmonic interference, which can severely degrade diagnostic accuracy [15,16]. Ultrasonic techniques, although repeatable under controlled conditions, are strongly influenced by background acoustic noise, complex propagation paths, and sensor installation positions, leading to unstable echo characteristics and limited reliability in field applications [17,18].

To avoid deploying additional equipment while achieving real-time, accurate online monitoring of transformer winding deformation, reference [19] proposed a method that constructs a two-dimensional (2D) Lissajous curve from the terminal voltage and current waveforms of the transformer. Winding deformation is diagnosed by analyzing morphological changes in the resulting curve. Building on this approach, reference [20] extracted four characteristic parameters from the curve and established mathematical models that link these parameters to specific fault types and their severity levels. Reference [21] investigated the influence of signal frequency on the curve’s parameters and incorporated cloud computing to develop a remote monitoring system. To advance diagnostic automation, reference [22] integrated the 2D Lissajous curve with a support vector machine (SVM), enabling early-stage fault diagnosis with a classification accuracy of 82%. Reference [23] introduced a fault diagnosis framework based on 2D Lissajous trajectory features and artificial neural networks (ANN), demonstrating accurate identification and localization of five types of internal faults.

Nevertheless, these methods are limited to two-dimensional graphical analysis and fail to exploit the multidimensional information inherent in terminal electrical signals fully. Moreover, they lack a systematic consideration of the impact of complex operating conditions on affect monitoring performance during actual service. Consequently, their diagnostic accuracy and operational robustness remain limited in practical engineering applications. To further enhance monitoring performance, reference [24] proposed an online data-driven method based on a 3D Lissajous curve, introducing the current differential value as the third dimension. This extension constructs a richer state-feature space, enabling accurate short-circuit impedance identification through a data-driven algorithm. However, this research remains confined to determining fault occurrence from the curve and lacks a practical approach for distinguishing fault types. Furthermore, the rich spatial and geometric information embedded in the 3D Lissajous curve remains underexploited, which limits the further application of this method for precise diagnosis.

To address the aforementioned challenges, this paper proposes an online monitoring method for transformer winding deformation based on 3D Lissajous Curves. First, a 3D Lissajous curve is constructed using the primary current i1(t), the derivative of the primary current di1(t)/dt, and the voltage difference between the primary and secondary sides Δu(t) as coordinate axes. Corresponding spatial curves for typical winding faults are simulated under both sinusoidal and non-sinusoidal excitation. To thoroughly explore their intrinsic geometric characteristics, feature parameters are extracted using the covariance matrix and an adaptive particle swarm optimization (APSO) algorithm with adaptively adjusted inertia weight and learning factors. These parameters include the lengths of the major and minor axes, eccentricity, the inclination angle of the major axis relative to the XY-plane, and planar parameters, thereby quantifying morphological variations in the curves. Subsequently, the extracted feature parameters are input into an SVM for fault classification, and the model’s adaptability is evaluated under various operating conditions. Comparative experiments against conventional fault diagnosis methods based on 2D Lissajous curves demonstrate that the proposed approach achieves superior classification accuracy and reliability. Furthermore, this method requires no additional detection equipment, offering the advantages of low cost and ease of implementation. It thus provides a novel technical pathway for the online monitoring and classification of winding deformation faults.

2  Principle and Analysis of the Lissajous Curve

2.1 Basic Principle

When two sinusoidal signals are superimposed along perpendicular axes, they form a closed trajectory known as a Lissajous curve, provided that their frequencies are integer multiples of one another. If the frequency ratio of two sine signals, x(t), y(t) is 1:1, their functional expressions can be defined as follows:

x(t)=A1cos⁡(ωt+θ1)y(t)=A2cos⁡(ωt+θ2)(1)

where A1 and A2 are the amplitude of the signal x(t), y(t), respectively; ω is the signal frequency; θ1 and θ2 are the phase angle of the signal x(t), y(t), respectively.

After eliminating the parameter ωt in the above formula, the following expression is obtained:

x2A12+y2A22−2xyA1A2cos⁡(θ1−θ2)=sin2⁡(θ1−θ2)(2)

In the Cartesian coordinate system, the Lissajous curve assumes different forms depending on the value of the phase difference θ1−θ2. When θ1−θ2 equals 0 or π, the curve forms a straight line; when θ1−θ2 lies between 0 and π, it takes the shape of an ellipse. During normal operation of a transformer, the terminal voltage and current both vary sinusoidally at the same frequency. However, due to the influence of magnetization characteristics and load impedance, a certain phase deviation exists between the electrical quantities on the two sides. As a result, the curve constructed from the terminal voltage and current forms an ellipse.

2.2 Construction of the 3D Lissajous Curve for Transformers

The T equivalent circuit of the secondary side of the single-phase double-winding transformer reduced to the primary side is shown in Fig. 1 [19]. In the diagram, u1(t), u2(t), i1(t), i2(t) are the voltage and current of the primary and secondary sides of the transformer; R1, R2, L1, L2 are the equivalent resistance and inductance of the primary and secondary windings, Δu(t) is the voltage difference between the primary and secondary sides, and ZL is the secondary side load.

Figure 1: T equivalent circuit of a transformer.

From the equivalent circuit diagram in Fig. 1, we obtain:

Δu(t)=u1(t)−u2(t)=R1i1(t)+L1di1(t)/dt+R2i2(t)+L2di2(t)/dt(3)

In practical engineering applications, the excitation current is sufficiently small to be neglected; consequently, i1(t)=i2(t), and Eq. (3) can be simplified as:

Δu(t)=u1(t)−u2(t)=Ri1(t)+Ldi1(t)/dt(4)

where R=R1+R2 is the equivalent total resistance of the windings, and L=L1+L2 is the equivalent total leakage inductance of the windings.

Given the difficulty of directly obtaining the differential of the measured value, it can be solved by an integral operation on both sides of Eq. (4):

{∫t1t2Δu(t)dt=R∫t1t2i1(t)dt+L[i1(t2)−i1(t1)]⋮∫tn−1tnΔu(t)dt=R∫tn−1tni1(t)dt+L[i1(tn)−i1(tn−1)]}(5)

When the time step is T, the trapezoidal integration yields:

(T2)∑k=1n−1(△uk+1+△uk)=R(T2)∑k=1n−1(i1(k+1)+i1(k))+L∑k=1n−1(i1(k+1)−i1(k))(6)

According to Section 2.1, this paper combines two 2D Lissajous curves Δu(t)−i1(t) and Δu(t)−di1(t)/dt to construct 3D Lissajous curves with i1(t), di1(t)/dt, Δu(t) as x, y, z axes respectively, as shown in Fig. 2. It is observed that the resulting 3D Lissajous curve lies entirely within a specific plane in space and is mathematically represented by Eq. (6) [24]. Since the parameters of transformer winding resistance and leakage inductance are closely related to the winding structure and relative positions, any winding deformation that alters the geometric dimensions will cause corresponding changes in the winding resistance and leakage inductance, thereby leading to detectable variations in the spatial curve.

Figure 2: 3D Lissajous curve.

3  Transformer Winding Deformation Simulation

3.1 Transformer Model Construction

Reference [19] indicates that transformer windings can be represented by distributed RLC equivalent circuits, as shown in Fig. 3. In Fig. 3, Rs and Ls represent the winding resistance and inductance in the distributed equivalent circuit, respectively. Cs and Cg represent the capacitance between the winding pancakes and the winding to ground capacitance, respectively. When the transformer is subjected to mechanical damage such as winding deformation, core displacement, or external short-circuit faults, the parameters of resistance, capacitance, and inductance in its equivalent circuit will change accordingly [25]. Axial displacement alters the relative positions of the high-voltage and low-voltage windings, reducing the winding inductance. Radial deformation causes the high-voltage winding to expand outward, increasing its inductance and decreasing inter-disk capacitance. Inter-disk short circuits alter the transformer ratio, with faults of varying severity represented by different numbers of short-circuited RLC loops in the winding equivalent model.

Figure 3: Distributed parameter RLC equivalent circuit of a transformer winding.

To investigate winding deformation faults, a distributed parameter transformer model was established using PSIM circuit simulation software, as shown in Fig. 3. The modeled transformer has a rated voltage of 35/10.5 kV and consists of 60 disks. Simulations were performed for inter-disk short-circuits, axial displacement, and radial deformation. The fault severity was varied by setting the proportion of faulty disks to 30% and 60% of the total winding. To facilitate online monitoring, a 50 Hz sinusoidal voltage source was employed as the excitation signal.

3.2 Fault Simulation under Sinusoidal Input

The simulated data were processed in MATLAB using trapezoidal integration to generate 3D Lissajous curves for both standard windings and various fault conditions, followed by spatial curve analysis. Fig. 4a presents the 3D Lissajous curves for 30% and 60% inter-disk short-circuit faults in the transformer winding. As illustrated, with increasing fault severity, the fault curves exhibit significant deviation from the standard winding curve. In three-dimensional space, the fault curve progressively moves away from the plane associated with the regular winding and shifts notably toward the XY plane. Concurrently, both the major and minor axes of the fault curve expand substantially, resulting in a clear increasing trend in the enclosed area.

Figure 4: 3D Lissajous curves for different fault types at 30% and 60% severity under sinusoidal input conditions: (a) Inter-disk short-circuit fault; (b) Axial displacement fault; (c) Radial deformation fault.

Fig. 4b displays the 3D Lissajous curves for transformer windings with 30% and 60% axial displacement. As illustrated, as the fault severity increases, the fault curve shows slight deviations relative to the standard winding curve. It gradually departs from the reference plane of the regular winding in three-dimensional space and shifts moderately toward the XY plane. The major and minor axes of the fault curve increase modestly, leading to a gentle expansion of the curve.

The 3D Lissajous curves for transformer windings with 30% and 60% radial deformation are shown in Fig. 4c. As illustrated, with increasing fault severity, the fault curves demonstrate slight deviations from the standard winding curves. The fault curve gradually leaves the plane corresponding to the regular winding in three-dimensional space and moves slightly away from the XY plane. Notably, the major axis of the fault curve decreases, a trend that differs from those observed in short-circuit faults and axial displacement faults.

3.3 Fault Simulation under Non-Sinusoidal Input

With the rapid development of modern power systems toward intelligent and distributed paradigms, the large-scale integration of renewable energy, widespread deployment of distributed generation, and proliferation of power electronic devices across industrial, transportation, and residential applications have progressively established non-linear loads as dominant components in power networks. These loads achieve energy conversion through switching operations, exhibiting distinctly non-linear voltage-current characteristics that cause severe distortion in grid voltage waveforms, thereby posing multifaceted challenges to system stability and equipment security. To investigate the impact of harmonic-distorted input voltages on the proposed 3D Lissajous curves, the primary-side input voltage was constructed as a superposition of a fundamental-frequency component and a third-order harmonic component, expressed as:

u1(t)=U1sin⁡(ωt)+U3sin⁡(3ωt+φ3)(7)

where U1 is the amplitude of the fundamental component, U3 is the amplitude of the third-order harmonic component, and φ3 is the harmonic phase angle. In this study, the third-order harmonic amplitude was set to U3=0.1U1, and the harmonic phase angle was set to φ3=30∘, representing a severely distorted operating condition. Under this configuration, the total harmonic distortion (THD) of the input voltage is 10%. Fig. 5 illustrates the corresponding 3D Lissajous curves under various fault conditions in the presence of harmonic interference.

Figure 5: 3D Lissajous curves for different fault types at 30% and 60% severity under non-sinusoidal input conditions: (a) Inter-disk short-circuit fault; (b) Axial displacement fault; (c) Radial deformation fault.

As shown in Fig. 5, the presence of harmonic components significantly alters the geometric shape of the 3D Lissajous curve. Specifically, Fig. 5a shows that, under inter-disk short-circuit faults with harmonic distortion, the fault curve progressively approaches the XY plane as the fault severity increases, accompanied by significant expansion of both the major and minor axes and a marked enlargement of the enclosed area. Fig. 5b demonstrates that, under the same harmonic conditions, as axial displacement increases, the fault curve gradually moves toward the XY plane, with both axes lengthening and the curve exhibiting an outward expansion. As seen in Fig. 5c, for radial displacement faults in the presence of harmonics, the curve gradually shifts away from the XY plane as the fault severity increases, its major axis shortens, and the overall curve contracts inward.

4  Spatial Geometric Feature Extraction

The preceding sections conducted simulation studies on typical transformer faults under both sinusoidal and non-sinusoidal excitation conditions, providing a preliminary qualitative analysis of the impacts of different fault types based on the spatial morphology of the 3D Lissajous curves. The simulation results indicate that different fault types produce distinguishable variations in the trajectory’s shape, deflection angle, and spatial distribution. However, qualitative analysis alone has limitations in accurately identifying specific fault types. To quantitatively characterize these curve variations, this paper calculates a set of feature parameters derived from the 3D Lissajous curve. In the 2D Lissajous curve analysis method, geometric parameters—particularly the lengths of the major and minor axes, eccentricity, and the angle—have been extensively validated as effective indicators for diagnosing winding deformation [19–21]. These features sensitively capture geometric distortions in 2D Lissajous curves resulting from winding deformation, while also demonstrating strong stability and repeatability. This paper extends the methodology to three-dimensional space, aiming to retain the advantages of the two-dimensional approach while further enhancing the ability to identify winding deformation faults by leveraging the richer information provided by the additional dimension.

To fully exploit the unique spatial structural characteristics of 3D Lissajous curves, this study introduces planar parameters as supplementary features. The selection of these features over other geometric features is primarily motivated by the following considerations. First, the major axis, minor axis, eccentricity, and angle exhibit clear physical correlations with the magnitude and phase relationships of voltage and current, directly reflecting the overall scale and orientation of the curve. As a result, they are susceptible to variations in electrical parameters caused by winding deformation faults. Second, planar parameters effectively characterize the curve’s spatial position in three-dimensional space and provide insight into changes in the winding leakage reactance. In contrast, localized geometric features—such as curvature—are more susceptible to harmonic interference and demonstrate comparatively lower stability. Consequently, they are not adopted as primary identification features in this work.

4.1 Analysis of Covariance Matrix

Due to digital sampling, the 3D Lissajous curve can be regarded as a spatial point set composed of i1(t), di1(t)/dt, Δu(t). To quantify the shape characteristics of the 3D Lissajous curve, this paper employs covariance matrix analysis for feature extraction. The core idea of this method is to regard a series of discrete point sets as a multivariate data distribution. By calculating the covariance matrix of the distribution and analyzing the eigenvalues and eigenvectors, the main distribution direction and stretching degree of the data point set in space are revealed, to describe the geometric characteristics of the 3D Lissajous curve [26].

Firstly, the data of N sampling points are averaged and centralized:

Cx=∑j=1NFxjN,Cy=∑j=1NFyjN,Cz=∑j=1NFzjN(8)

Fxj′=Fxj−CxFyj′=Fyj−CyFzj′=Fzj−Cz(9)

where, Cx, Cy, Cz are the center of the ellipse, Fxj, Fyj, Fzj, are the primary current, derivative of the primary current and the voltage difference between the primary and secondary sides obtained by the trapezoidal approximation method in Eq. (6). Fxj′, Fyj′, Fzj′ are the values after centralization.

The calculation of the spatial curve covariance matrix A is as follows:

A=[∑j=1NFxj′2∑j=1NFxj′Fyj′∑j=1NFxj′Fzj′∑j=1NFxj′Fyj′∑j=1NFyj′2∑j=1NFyj′Fzj′∑j=1NFxj′Fzj′∑j=1NFyj′Fzj′∑j=1NFzj′2]/N(10)

The characteristic decomposition of the covariance matrix A is performed to solve the characteristic equation:

Avi=λivi,i=1,2,3(11)

where λ1≥λ2≥λ3 are eigenvalues, v1, v2, v3 are the corresponding eigenvectors. The eigenvalue represents the degree of extension of the ellipse in each principal axis direction. The maximum eigenvalue λ1 corresponds to the square of the major axis of the ellipse, and the intermediate eigenvalue λ2 corresponds to the square of the minor axis of the ellipse. The eigenvectors of the two eigenvalues point to the major axis and the minor axis of the ellipse, respectively. The eigenvector corresponding to the minimum eigenvalue λ3 is the normal direction of the ellipse. Then the major axis length la and the minor axis length lb of the ellipse are as follows:

la=λ1,lb=λ2(12)

Based on the above geometric significance, this paper defines the ellipse’s major and minor axes la and lb, the eccentricity e, and the angle ψ between the ellipse’s major axis and the XY plane as spatial geometric features. The calculation formulas of e and ψ are as follows:

e=la−lbla+lbψ=arcsin⁡(|v1z|‖v1‖)(13)

where v1z is the component of the eigenvector v1 in the z-axis.

4.2 Extraction of Plane Parameters Based on APSO Algorithm

From the above, it can be seen that the 3D Lissajous curve lies in a space plane, which transforms the plane recognition problem into an optimization problem, and that the APSO algorithm is used to calculate the planar parameters.

As a classical heuristic optimization algorithm, particle swarm optimization (PSO) is inspired by the foraging behavior of birds. Its core idea is to find the optimal solution the problem by using the mutual cooperation and information sharing among individuals in the group [27]. Suppose that there are m particles in the search space, xtn=(x1,x2,⋯,xm) denotes the position of m particles in the search space at the tn iteration. vtn=(v1,v2,⋯,vm) denotes the velocity of m particles at the tn iteration. pbest represents the optimal position searched by the particle itself; gbest represents the optimal position searched by the group. The velocity and position update formula of the particle is:

vtn=vtn−1⋅w+c1⋅rand1⋅(pbest−xtn−1)+c2⋅rand2⋅(gbest−xtn−1)(14)

xtn=xtn−1+vtn(15)

where w is the inertia weight, c1 is self-learning factor and c2 is social learning factor; rand1 and rand2 are random numbers that change between [0, 1].

However, the fixed parameters w, c1, and c2 in the PSO algorithm limit the optimization ability of the particles, making it difficult to fly to a better position in time due to excessive inertia, and eventually falling into a local optimum [28]. To achieve accurate identification across different working conditions, an adaptive inertia weight and an adaptive learning factor control strategy are introduced into the traditional PSO algorithm, constructing a dynamic mapping between the iterative process and the search behavior to adapt to different problem types and search environments. APSO has dual advantages: it can maintain the particles’ wide-area search ability in the early stages of iteration and avoid premature convergence of the algorithm. With each iteration, the algorithm places greater emphasis on local search for the optimal solution, thereby improving its quality. Therefore, this paper updates the inertia weight w based on the number of iterations according to Eq. (16), causing w to gradually decrease.

w=wmin+(wmax−wmin)⋅e−10×(tnTmax)2(16)

where wmax is the maximum inertia weight and wmin is the minimum inertia weight, respectively, and Tmax is the maximum number of iterations.

To achieve the optimal balance between global search and local development, update the learning factor according to Eq. (17) as the number of iterations increases.

c1=cmax−(cmax−cmin)⋅tnTmaxc2=cmin+(cmax−cmin)⋅tnTmax(17)

where cmax is the maximum learning factor and cmin is the minimum learning factor.

In the APSO algorithm, each particle represents a potential solution to a planar parameter problem. The motion state of the particle is dynamically updated by the velocity vector (determining the moving step) and the position vector (determining the moving direction). This process depends on the particle’s current optimal position and its optimal position in the group, and the fitness function evaluates the particle’s position. Based on the above Eq. (6), this paper establishes a planar parameter identification model with the minimum sum of squared residuals as the fitness function J, as shown in Eq. (18):

J=min{(T2)∑k=1n−1(△uk+1+△uk)−R(T2)∑k=1n−1(i1(k+1)+i1(k))−L∑k=1n−1(i1(k+1)−i1(k))}2(18)

According to the above analysis, the steps of planar parameter identification based on the APSO algorithm are as follows:

(1)  Define the parameters to be identified, particle dimension, and fitness function. The parameters to be identified are R and L, that is, the particle dimension is 2-dimensional, and the search space is 2-dimensional. Set the appropriate population size and number of iterations.

(2)  Initialize the position and velocity of the particle swarm, and calculate the fitness value of the particle swarm according to Eq. (18), and obtain the optimal position of the individual and the population.

(3)  Iterative update. The inertia weight and learning factor are updated by Eqs. (16) and (17). The velocity and position of each particle are updated by Eqs. (14) and (15). The fitness value of the particle is calculated and compared with the optimal position of the historical individual and the population. Update yourself on the optimal position of the population.

(4)  Determine whether the maximum number of iterations has been reached. If the requirements are met, the identification is completed and the optimal positions of the population, namely the R and L parameters, are output; otherwise, the iteration continues.

The process of parameter identification of the APSO algorithm is shown in Fig. 6. A summary of all APSO parameters is provided in Table 1.

Figure 6: APSO flowchart.

5  Feature Extraction Results Analysis

5.1 Comparative Analysis of Algorithms

To verify the effectiveness and superiority of the APSO algorithm used in this paper for the planar parameter identification problem, this paper compares it with the traditional PSO algorithm and the Harris- hawks optimization algorithm (HHO), which has strong optimization ability and few controllable parameters [29]. Due to the low particle dimension of the method proposed in this paper, given the different update strategies across algorithms, the population size is set to 50, the maximum number of iterations to 100, and the sum of residual squares as the fitness function. According to the previous analysis, the inter-disk short-circuit fault under non-sinusoidal input has the most significant effect on the curve, so this working condition is selected to verify the performance of the three optimization algorithms. The convergence curves of fitness values for different optimization algorithms are shown in Fig. 7.

Figure 7: Convergence curves of fitness values for different optimization algorithms.

As shown in Fig. 7, the HHO algorithm has the fastest convergence speed, the APSO algorithm is second, and the PSO algorithm is the slowest. The iterative convergence times of the three algorithms are 10, 27, and 38, respectively. Although the HHO algorithm converges rapidly, its fitness value is 5.0238E−4, the PSO algorithm is 2.5339E−5, and the APSO algorithm is 1.7058E−6, which is the smallest of the three. The results show that the APSO algorithm achieves a good balance between convergence speed and optimization accuracy, delivers better parameter identification within a reasonable number of iterations, and exhibits higher comprehensive performance.

5.2 Analysis of Spatial Geometric Features

To quantitatively analyze the influence of different fault types on the shape of 3D Lissajous curve, based on the analysis method established in Section 4, the geometric characteristic parameters of spatial curve under various working conditions are extracted in MATLAB, include the lengths of the major and minor axes, eccentricity, the inclination angle of the central axis relative to the XY-plane and parameters (R, L) describing the plane where the curve is located. Under sinusoidal input, the calculation results of characteristic parameters corresponding to different fault types are shown in Table 2.

As shown in Table 2, during an inter-disk short-circuit, the major and minor axes of the curve increase significantly with fault severity. In contrast, eccentricity increases, and the angle decreases. The curve exhibits noticeable spatial deflection, while the planar parameters decrease substantially. In the case of axial displacement faults, both major and minor axis lengths increase slightly, though the overall change remains limited; the eccentricity tends to decrease, the angle decreases, and while the planar parameter R shows no clear trend, parameter L decreases. For radial deformation faults, the curve contracts along the major axis, with no significant change in minor axis length, and the eccentricity decreases slightly; the angle increases, and the planar parameters rise accordingly. The above analysis demonstrates that different fault types exhibit distinct patterns of variation in the geometric parameters of the 3D Lissajous curve. Consequently, by extracting and analyzing features such as the major and minor axes, eccentricity, angle, and planar parameters of the curve, winding fault types can be effectively identified and distinguished.

Under non-sinusoidal input conditions, the characteristic parameters for each fault type are summarized in Table 3. As shown in the table, although harmonic components are present in the input voltage, the overall variation trends of the characteristic parameters for each fault type remain largely consistent with those under sinusoidal excitation. However, the introduction of harmonics disturbs some parameters. Specifically, under non-sinusoidal excitation, the characteristic parameters related to the major axis show an increase compared to those under sinusoidal conditions. In contrast, the variation in the minor axis parameters is relatively small and shows no significant differences. This asymmetric change between the major and minor axes further increases the eccentricity. To quantify the influence of harmonics on the characteristic parameters, the 60% fault-severity case is used as an example for comparative analysis. For an inter-disk short-circuit fault, the major axis under non-sinusoidal input is 0.82% higher than under sinusoidal input. For axial displacement faults, the major axis increases by 0.56%, and for radial deformation faults, it increases by 0.52%. Although the increase in the major axis is not substantial, in practical engineering applications, such minor changes could interfere with the identification of early-stage or slight winding deformation faults. It is worth noting that the characteristics of the minor axis, angle, and planar parameters remain essentially consistent under both input conditions and are less affected by harmonics, demonstrating strong robustness. Therefore, in harmonic-rich environments, parameters with minimal harmonic interference can be preferentially selected as the primary basis for identifying winding fault types, thereby enhancing the adaptability and reliability of the diagnostic method.

6  Fault Diagnosis and Quantitative Assessment

6.1 Construction of an SVM-Based Diagnostic Model

To precisely identify the types of winding deformation faults, a fault diagnosis model employing SVM is established in this paper. SVM is a machine learning algorithm based on statistical learning theory, which maps input vectors into a high-dimensional feature space via kernel functions to construct an optimal hyperplane. This approach demonstrates significant advantages in addressing pattern recognition problems characterized by small samples and nonlinearity [30]. The radial basis function (RBF) is selected as the kernel function, with the penalty coefficient set to c=8 and the kernel parameter to γ=10. As described in Section 4, the six feature parameters extracted from the 3D Lissajous curve constitute the input feature vector X=(la,lb,e,ψ,R,L). The model output corresponds to three winding states: Inter-disk short-circuit fault, axial displacement, and radial deformation. The specific classification process is illustrated in Fig. 8.

Figure 8: Fault diagnosis flowchart.

To validate the effectiveness of the proposed model, this study utilized a total of 177 sample datasets. The collected samples encompass fault severity levels ranging from 2% to 98%, with an equal distribution of 59 samples per fault category. The datasets were randomly partitioned into a training set (124 samples) and a test set (53 samples) at a 7:3 ratio. Model performance was quantified using accuracy, precision, and confusion matrices. The calculation formulas for each metric are detailed in [31]. All experimental procedures were implemented in MATLAB R2023b.

6.2 Diagnostic Findings Analysis

6.2.1 Classification Results under Sinusoidal Input

Under standard simulation conditions with a rated load, no measurement errors, and a sampling frequency of 100 kHz, the classification results obtained using the aforementioned feature parameters and the SVM model are shown in Fig. 9.

Figure 9: Confusion matrix under sinusoidal input.

As shown in Fig. 9, the proposed method achieves a classification accuracy of 98.15%. Simulation results demonstrate that the diagnostic model based on 3D Lissajous features and an SVM can accurately distinguish interturn short circuits, radial deformation, and axial displacement faults with fault severity exceeding 2%.

6.2.2 Classification Results under Non-Sinusoidal Input

To evaluate the model’s robustness under non-ideal grid conditions, harmonic interference with a Total Harmonic Distortion (THD) of 10% was introduced into the simulation power supply, as outlined in Section 3.3. The fault data generated under these non-sinusoidal conditions served as the test set, which was then used to train the SVM model established in Section 6.1 for classification. The corresponding confusion matrix is shown in Fig. 10.

Figure 10: Confusion matrix under non-sinusoidal input.

The analysis indicates that, even under significant harmonic interference, the model maintains an overall classification accuracy of 98.15%. This result further validates the proposed method’s excellent adaptability and robustness to non-sinusoidal inputs, underscoring its strong applicability in practical environments with harmonic disturbances.

6.2.3 Classification Results under Load Variations

To verify the effectiveness of the proposed method under varying load conditions, tests were conducted over a fault severity range of 2% to 98% at 20% intervals. Corresponding load ratios of 40%, 70%, 90%, 120%, and 150% were configured for each fault level. The simulation data generated under all the aforementioned operating conditions were compiled into a test dataset, which was then fed into the pre-trained SVM model from Section 6.1 for classification. The confusion matrix for the classification results under varying load conditions is shown in Fig. 11.

Figure 11: Confusion matrix under load fluctuation.

As shown in Fig. 11, under test conditions encompassing a wide range of load variations (from 40% to 150%), the overall classification accuracy of the proposed method remains as high as 96.30%. This result demonstrates that the method maintains excellent fault-discriminative capability and robustness across diverse load conditions.

6.2.4 Classification Results at Different Sampling Frequencies

Variations in sampling frequency affect the numerical integration results of winding voltage and current signals, leading to proportional scaling of the geometric dimensions of the corresponding 3D Lissajous curves. This scaling effect can compromise classification accuracy. Notably, some parameters—such as the inclination angle, eccentricity, and planar attributes of these curves—are scale-invariant. To eliminate the scale-related interference introduced by varying sampling frequencies and to ensure feature consistency across different sampling conditions, a frequency-based normalization method is proposed. Specifically, the major and minor axis lengths derived at any sampling frequency fs are proportionally scaled relative to a reference sampling frequency (fs0=100 kHz). The normalization formula is defined as follows:

la,norm=la⋅(fsfs0)lb,norm=lb⋅(fsfs0)(19)

where la,norm and lb,norm represent the equivalent axial lengths normalized to the reference frequency. This process effectively eliminates scale inconsistencies caused by differences in sampling frequency.

To validate the impact of sampling frequency on classification performance, fault datasets were generated at simulation sampling frequencies of 100 kHz, 50 kHz and 10 kHz. Feature normalization was then applied using the formula described above. The normalized data were used as the test set and input into an SVM model trained on reference data sampled at 100 kHz for classification verification. The comparison of classification accuracy across different sampling frequencies is presented in Table 4.

As shown in Table 4, after data normalization, the classification accuracy experienced a slight decrease when the sampling frequency was reduced to 10 kHz, yet it remained at a high level of 96.30%. This indicates that the proposed method maintains robust fault-discrimination accuracy across varying sampling frequencies.

6.2.5 Classification Results Considering Measurement Errors

In the practical construction of the 3D ellipse, measurement errors from instruments and rounding errors in computation are inherently present. To evaluate the impact of measurement errors on fault classification, ±0.5% amplitude errors were introduced into both the voltage and current signals in the dataset. The data incorporating these errors was then used as a test set and input into the SVM model established in Section 6.1 for classification. The classification accuracy under different error scenarios is presented in Table 5.

As shown in Table 5, even with a ±0.5% amplitude measurement error, the diagnostic model’s overall accuracy remains consistently high at 98.15%. This indicates that the influence of such measurement inaccuracies on classification performance is negligible.

6.2.6 Comparative Analysis with the 2D Lissajous Feature Method

To verify the superior performance of the proposed 3D Lissajous curve features in diagnosing winding deformation faults, a comparative analysis was conducted against a diagnostic method based on 2D Lissajous features and SVM. For the comparative method, key parameters of the 2D Lissajous curve—including its major axis (la,2D), minor axis (lb,2D), eccentricity (e2D), and angle (ψ2D)—were extracted using a covariance matrix to form the feature vector X2D=(la,2D,lb,2D,e2D,ψ2D). This vector was then input into an SVM classifier identical to the one described in Section 6.1. The comparison utilized the same training set (124 samples) and test set (53 samples), and was performed under standard sinusoidal operating conditions. The classification confusion matrix for the 2D Lissajous method is shown in Fig. 12. The proposed 3D Lissajous-based method achieves an overall accuracy of 98.15%, which is markedly higher than the 90.74% achieved by the conventional 2D Lissajous method. This represents a significant improvement of 7.41 percentage points.

Figure 12: Confusion matrix for 2D Lissajous curves.

The classification performance of the 2D Lissajous-based method under various non-ideal operating conditions is summarized in Table 6. As shown, the overall classification accuracy of the method using 2D Lissajous features ranges from 87.04% to 90.74% across different test scenarios. In comparison, the accuracy achieved by the proposed 3D Lissajous feature-based method ranges from 96.30% to 98.15%. This comparative result demonstrates that, under both ideal testing conditions and a variety of complex operational environments, the proposed 3D Lissajous feature-based approach exhibits significantly superior diagnostic accuracy and enhanced model robustness relative to the conventional 2D Lissajous method.

6.3 Computational Efficiency and Real-Time Analysis

To ensure that the proposed method satisfies the real-time requirements for online monitoring of transformer winding faults, an evaluation of the temporal consumption across each processing stage was conducted. The tests were performed on a platform running MATLAB R2023b, with a hardware configuration comprising an Intel(R) Core i5-8265U processor and 8 GB of RAM. A dataset of 59 simulation records was used, each sampled at 100 kHz and corresponding to one power-frequency cycle (20 ms). The time distribution for the feature extraction stage is illustrated in the accompanying Fig. 13. As shown, the predominant computational cost in this stage lies in the APSO module, whose complexity is mainly governed by the population size (Ns) and the maximum number of iterations (Tmax). To investigate the optimal trade-off between computational efficiency and diagnostic accuracy, comparative tests were performed on different parameter combinations, with the results detailed in Table 7.

Figure 13: Distribution of processing time per cycle in the feature extraction stage.

As shown in Table 7, the computational time of the proposed method is overwhelmingly spent in the APSO-optimized feature extraction stage, accounting for 98% of the total processing time. In contrast, the covariance analysis and SVM classification processes demonstrate exceptional efficiency, each requiring less than 0.2 ms; thus, their impact on the overall system’s real-time performance is negligible. When the parameters were optimized from (Ns = 50, Tmax = 100) to (Ns = 20, Tmax = 70), the average processing time per diagnostic instance decreased significantly from 33.03 to 10.51 ms, representing a 68.18% improvement in computational efficiency. Concurrently, diagnostic accuracy remained stable, consistently exceeding 96.30%. Compared to the maximum accuracy of 98.15% achieved under the previous high-parameter configuration, the maximum observed precision degradation was limited to 1.85 percentage points. This indicates that the method retains its core fault-identification capability even after parameter simplification.

These results demonstrate that, in engineering practice, judicious reductions in APSO optimization parameters can substantially enhance algorithmic execution efficiency while incurring only a minimal compromise in diagnostic accuracy. This strategy effectively ensures that the real-time operational requirements of online monitoring systems are met.

7  Conclusion

This paper presents a novel online monitoring method for transformer winding deformation based on 3D Lissajous curves. Spatial geometric features are extracted from the constructed curves using covariance matrix analysis, while the corresponding spatial plane parameters are determined via an APSO algorithm. These components collectively form an integrated set of feature parameters that effectively characterize the state of the windings. Comparative analysis under both sinusoidal and harmonic-rich non-sinusoidal excitation conditions demonstrates that the extracted feature parameters exhibit strong immunity to grid harmonic interference. Simulation results based on an SVM fault diagnosis model confirm that the proposed method can reliably identify winding deformations exceeding 2% in severity. Furthermore, the model maintains a diagnostic accuracy of ≥96.30% under various practical disturbances—including harmonic distortion with total harmonic distortion (THD) up to 10%, load variations ranging from 40% to 150%, fluctuations in sampling frequency, and the introduction of ±0.5% amplitude errors in both voltage and current signals—indicating high robustness for engineering applications. Compared with conventional diagnostic approaches based on 2D Lissajous curves, the proposed method achieves an approximately 7.41% improvement in diagnostic accuracy. Real-time performance analysis further verifies that the method satisfies the requirements for online monitoring, demonstrating its practical feasibility in engineering contexts.

It is worth noting that the current investigation is primarily based on a simulation platform, and all analyses and conclusions are derived from simulated data. During the modeling process, the structural parameters of the windings and the fault types were simplified, and the influence of electromagnetic interference typically present in actual operational environments was not fully taken into account. These idealized assumptions limit the applicability of the proposed method in practical engineering contexts. Moreover, although this study superimposed third-order harmonics with fixed amplitude and phase to assess the method’s robustness under non-ideal input conditions, the harmonic components in real power grids are considerably more complex, often comprising multiple orders and randomly phased harmonics.

In light of the above limitations, future work will systematically focus on the following three aspects: (1) Further validation of the feasibility and accuracy of the proposed method through laboratory-scale physical simulations or field measurement data; (2) Exploration and integration of noise suppression techniques tailored for strong electromagnetic interference environments to improve the engineering applicability of the method under complex operating conditions; (3) In-depth investigation of the method’s adaptability to transformers of various types and capacities under diverse operating scenarios, with particular emphasis on diagnostic performance under complex harmonic conditions involving multiple orders and random phase characteristics, thereby enhancing its generalizability in engineering practice.

Acknowledgement: Not applicable.

Funding Statement: This research was funded by National Natural Science Foundation of China, grant number 52277012.

Author Contributions: Conceptualization, Zhenhua Li and Xinyu Yue; methodology, Zhenxing Li and Xinyu Yue; software, Tao Zhang and Xinyu Yue; validation, Yanchun Xu, Xiaozhen Zhao and Xinyu Yue; formal analysis, Tao Zhang, Zhenhua Li and Xinyu Yue; investigation, Tao Zhang and Xinyu Yue; resources, Xiaozhen Zhao and Zhenhua Li; data curation, Yanchun Xu and Xinyu Yue; writing—original draft preparation, Zhenhua Li and Xinyu Yue; writing—review and editing, Zhenxing Li and Xinyu Yue; supervision, Xiaozhen Zhao; funding acquisition, Zhenhua Li. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data available on request from the authors.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Abbasi AR. Fault detection and diagnosis in power transformers: a comprehensive review and classification of publications and methods. Electr Power Syst Res. 2022;209(11):107990. doi:10.1016/j.epsr.2022.107990. [Google Scholar] [CrossRef]

2. Liu J, Yao C, Yu L, Dong S, Liu Y. Using MLP to locate transformer winding fault based on digital twin. Front Energy Res. 2023;11:1175808. doi:10.3389/fenrg.2023.1175808. [Google Scholar] [CrossRef]

3. Chen Y, Zhao Z, Liu J, Tan S, Liu C. Application of generative AI-based data augmentation technique in transformer winding deformation fault diagnosis. Eng Fail Anal. 2024;159(4):108115. doi:10.1016/j.engfailanal.2024.108115. [Google Scholar] [CrossRef]

4. Liu J, Chen S, Wang F, Zhu D, Xu K, Chen P, et al. Theoretical evaluation of transformer short-circuit withstand capability and identification method of winding cumulative deformation. IEEE Trans Power Deliv. 2025;40(5):2454–65. doi:10.1109/tpwrd.2025.3575087. [Google Scholar] [CrossRef]

5. Chen X, Zhao Z, Guo F, Tan S, Wang J. Diagnosis method of transformer winding mechanical deformation fault based on sliding correlation of FRA and series transfer learning. Electr Power Syst Res. 2024;229:110173. doi:10.1016/j.epsr.2024.110173. [Google Scholar] [CrossRef]

6. Parkash C, Abbasi AR. Transformer’s frequency response analysis results interpretation using a novel cross entropy based methodology. Sci Rep. 2023;13(1):6604. doi:10.1038/s41598-023-33606-0. [Google Scholar] [PubMed] [CrossRef]

7. Wang X, Ding J. A novel method for winding faults diagnostic of power transformers based on parameter estimation. IEEJ Trans Electr Electron Eng. 2021;16(12):1563–72. doi:10.1002/tee.23463. [Google Scholar] [CrossRef]

8. Shamlou A, Reza Feyzi M, Behjat V. Winding deformation classification in a power transformer based on the time-frequency image of frequency response analysis using Hilbert-Huang transform and evidence theory. Int J Electr Power Energy Syst. 2021;129(5):106854. doi:10.1016/j.ijepes.2021.106854. [Google Scholar] [CrossRef]

9. Lavrinovich VA, Mytnikov AV. Development of pulsed method for diagnostics of transformer windings based on short probe impulse. IEEE Trans Dielect Electr Insul. 2015;22(4):2041–5. doi:10.1109/tdei.2015.004745. [Google Scholar] [CrossRef]

10. Sant’Ana WC, Lambert-Torres G, Bonaldi EL, Gama BR, Zacarias TG, dos Santos Areias IA, et al. Online frequency response analysis of electric machinery through an active coupling system based on power electronics. Sensors. 2021;21(23):8057. doi:10.3390/s21238057. [Google Scholar] [PubMed] [CrossRef]

11. Seifaddini N, Sekongo B, Fofana I, Akre S, Chehri A, Ouhrouche M. Towards winding deformation assessment from vibration signals using an optical sensor. IET Sci Meas Technol. 2025;19(1):e12224. doi:10.1049/smt2.12224. [Google Scholar] [CrossRef]

12. Zhou Y, Liu Y, Wang N, Han X, Li J. Partial discharge ultrasonic signals pattern recognition in transformer using BSO-SVM based on microfiber coupler sensor. Measurement. 2022;201(7):111737. doi:10.1016/j.measurement.2022.111737. [Google Scholar] [CrossRef]

13. Wu S, Ji S, Zhang Y, Wang S, Liu H. A novel vibration frequency response analysis method for mechanical condition detection of converter transformer windings. IEEE Trans Ind Electron. 2024;71(7):8176–80. doi:10.1109/tie.2023.3301514. [Google Scholar] [CrossRef]

14. Cheng Y, Bi J, Chang W, Xu Y, Pan X, Ma X, et al. Proposed methodology for online frequency response analysis based on magnetic coupling to detect winding deformations in transformers. High Volt. 2020;5(3):343–9. doi:10.1049/hve.2019.0134. [Google Scholar] [CrossRef]

15. Liu S, Zhang L, Yang L, Gu C, Wang Z. The identification method of the winding vibration faults of dry-type transformers. Electronics. 2023;12(1):3. doi:10.3390/electronics12010003. [Google Scholar] [CrossRef]

16. Nezhad A, Samimi M. Investigation of transformer vibration characteristics using the finite element method. Sci Iran. 2024;31(5):441–57. doi:10.24200/sci.2022.59006.6012. [Google Scholar] [CrossRef]

17. Wang J, Zhao Z, Zhu J, Li X, Dong F, Wan S. Improved support vector machine for voiceprint diagnosis of typical faults in power transformers. Machines. 2023;11(5):539. doi:10.3390/machines11050539. [Google Scholar] [CrossRef]

18. Ilkhechi HD, Samimi MH. Applications of the acoustic method in partial discharge measurement: a review. IEEE Trans Dielect Electr Insul. 2021;28(1):42–51. doi:10.1109/tdei.2020.008985. [Google Scholar] [CrossRef]

19. Abu-Siada A, Islam S. A novel online technique to detect power transformer winding faults. IEEE Trans Power Deliv. 2012;27(2):849–57. doi:10.1109/tpwrd.2011.2180932. [Google Scholar] [CrossRef]

20. Ftatsi Mbetmi GDP, Tsegaing Tchatchueng F, Andre B, Macine Ngong S, Ngasop Stephane N. Mathematical models of power transformers winding faults diagnosis based on voltage-current characteristics. Int J Eng Technol. 2024;13(1):156–66. doi:10.14419/9qpy4k87. [Google Scholar] [CrossRef]

21. Zhou Y, Ma Y, Dai M, Zhu T, Du J, Wu H, et al. Remote online monitoring system for transformer winding deformation based on current-voltage characterized graph. Energy Rep. 2023;9(9):688–95. doi:10.1016/j.egyr.2022.11.032. [Google Scholar] [CrossRef]

22. Ngong SM, Guy-de-patience FM, Gonza M, Ngasop N. A machine learning based model for classifying incipient faults and degree in power transformer windings using voltage current diagram technique. Meas Energy. 2025;7(1):100056. doi:10.1016/j.meaene.2025.100056. [Google Scholar] [CrossRef]

23. Youssef M, Abdelaziz ES, Saad H, Attia M. New approach for accurate discrimination and location of power transformers with different internal winding faults. PLoS One. 2024;19(10):e0309926. doi:10.1371/journal.pone.0309926. [Google Scholar] [PubMed] [CrossRef]

24. Hong T, Deswal D, de Leon F. An online data-driven technique for the detection of transformer winding deformations. IEEE Trans Power Deliv. 2018;33(2):600–9. doi:10.1109/tpwrd.2017.2707922. [Google Scholar] [CrossRef]

25. Zhang Z. Derivation of transformer winding equivalent circuit by employing the transfer function obtained from frequency response analysis data. IET Electr Power Appl. 2024;18(7):826–40. doi:10.1049/elp2.12436. [Google Scholar] [CrossRef]

26. Eftekhari M, Moallem M, Sadri S, Hsieh MF. Online detection of induction motor’s Stator Winding short-circuit faults. IEEE Syst J. 2014;8(4):1272–82. doi:10.1109/jsyst.2013.2288172. [Google Scholar] [CrossRef]

27. Jain M, Saihjpal V, Singh N, Singh SB. An overview of variants and advancements of PSO algorithm. Appl Sci. 2022;12(17):8392. doi:10.3390/app12178392. [Google Scholar] [CrossRef]

28. Gong Z, Rao T, Wang G, Li Z, Luo Z, Zhu J, et al. Fault diagnosis method of transformer based on improved particle swarm optimization XGBoost. High Voltage Appar. 2023;59(8):61–9. (In Chinese). doi:10.13296/j.1001-1609.hva.2023.08.007. [Google Scholar] [CrossRef]

29. Shan W, He X, Liu H, Heidari AA, Wang M, Cai Z, et al. Cauchy mutation boosted Harris hawk algorithm: optimal performance design and engineering applications. J Comput Des Eng. 2023;10(2):503–26. doi:10.1093/jcde/qwad002. [Google Scholar] [CrossRef]

30. Zahra ST, Imdad SK, Khan S, Khalid S, Baig NA. Power transformer health index and life span assessment: a comprehensive review of conventional and machine learning based approaches. Eng Appl Artif Intell. 2025;139(3):109474. doi:10.1016/j.engappai.2024.109474. [Google Scholar] [CrossRef]

31. Valero-Carreras D, Alcaraz J, Landete M. Comparing two SVM models through different metrics based on the confusion matrix. Comput Oper Res. 2023;152(1):106131. doi:10.1016/j.cor.2022.106131. [Google Scholar] [CrossRef]