Time synchronization network test method, apparatus, device, and medium

Abstract

The application relates to the technical field of vehicle-mounted control, in particular to a time synchronization network test method and device, equipment and a medium, the method comprising the following steps: acquiring a state equation and system parameters of a preset nonlinear dynamic disturbance; injecting a phase perturbation sequence into a measured device; capturing a synchronization error response sequence output by the measured device under the perturbation excitation; performing variational modal decomposition on the synchronization error response sequence to obtain K intrinsic modal functions with limited bandwidth; performing Hilbert transformation on each intrinsic modal function to obtain a Hilbert spectrum of the measured time-space synchronization loop; constructing a parameterized reduced-order model of the system based on the intrinsic modal functions, predicting a stable domain boundary and an unstable bifurcation point of the system in a parameter space, and generating a synchronization stability boundary atlas to verify bidirectional synchronization capability of the measured device according to the synchronization stability boundary atlas. The dynamic performance of the device can be deeply evaluated by increasing active excitation, so that automatic time testing is performed, and the testing efficiency is improved.

Description

Technical Field

[0001] This application relates to the technical field of vehicle control, and in particular to time synchronization network testing methods, apparatus, equipment and media. Background Technology

[0002] As the electronic and electrical architecture of intelligent vehicles evolves towards central computing and regional control, Time-Sensitive Networking (TSN) technology has become a key infrastructure supporting advanced autonomous driving, high-precision sensor synchronization, and collaborative in-vehicle infotainment systems. The gPTP protocol (IEEE 802.1AS), through a master-slave clock hierarchy distribution mechanism, can theoretically achieve microsecond-level synchronization accuracy across the entire network. However, in practical engineering applications, synchronization accuracy is affected by multiple factors, including crystal oscillator temperature drift, switch forwarding delay uncertainty, wiring harness electromagnetic interference, and differences in protocol stack implementation, leading to significant deviations between measured accuracy and theoretical values.

[0003] In related technologies, manual measurement using an oscilloscope is typically relied upon, which is inefficient and introduces human reading errors (>10ns). When the master and slave devices use heterogeneous clocks (such as a mixture of GPS synchronization and free-running crystal oscillators), the time base needs to be manually calibrated, which is complex and lacks traceability. Furthermore, traditional methods only support single-machine testing and cannot flexibly verify the master-slave bidirectional synchronization capability of TSN switches. Summary of the Invention

[0004] This application provides a time synchronization network testing method, apparatus, device, and medium, which can perform automatic time testing by adding active stimulation to deeply evaluate the dynamic performance of the device and improve testing efficiency.

[0005] On one hand, embodiments of this application provide a time synchronization network testing method, the method comprising: In response to the access command of the device under test, the state equation and system parameters of the preset nonlinear dynamic disturbance are obtained; The modally matched phase perturbation sequence is injected into the device under test, and the synchronization error response sequence output by the device under test under perturbation excitation is continuously captured. The phase perturbation sequence is constructed based on the state equation and system parameters. Variational mode decomposition is performed on the synchronization error response sequence to obtain K eigenmode functions with finite bandwidth, each eigenmode function corresponding to a dynamic mode component of the test system; A Hilbert transform is performed on each intrinsic mode function to obtain the Hilbert spectrum of the spatiotemporal synchronization loop under test. The Hilbert spectrum reflects the variation of the system frequency components with time. Based on the intrinsic mode functions, a parameterized reduced-order model of the system is constructed. The stability domain boundary and instability bifurcation point of the system in the parameter space are predicted by the numerical extension algorithm. A synchronization stability boundary map of the device under test is generated, and the bidirectional synchronization capability of the device under test is verified according to the synchronization stability boundary map.

[0006] Optionally, obtaining the state equation and system parameters of the preset nonlinear dynamic disturbance includes: In response to the access command of the device under test, the type and protocol version of the device under test are identified, and the corresponding nonlinear dynamic system template is loaded; The initial state equation form is read from the system template, wherein the state variables include at least the phase difference, frequency difference and digital loop filter state, and the system parameters include at least the loop gain, crystal oscillator aging coefficient and filter coefficient. The system continuously captures the synchronization error sequence output by the device under test in a non-perturbation state, and estimates the initial values ​​of the system parameters online based on the subspace identification method, so that the state equation matches the actual dynamic characteristics of the device under test. The estimated system parameters are written to the parameter storage area in real time for subsequent perturbation sequence generation modules to use.

[0007] Optionally, before injecting the mode-matched phase perturbation sequence into the device under test and continuously capturing the synchronization error response sequence output by the device under test under perturbation excitation, the method further includes: Calculate the system's natural oscillation frequency and damping ratio based on the loop gain and crystal oscillator aging coefficient in the current system parameters; Based on the natural oscillation frequency and damping ratio, a sinusoidal perturbation sequence is generated as an intrinsic mode perturbation to excite the main oscillation mode of the system; A step perturbation sequence based on the natural oscillation frequency and damping ratio generation amplitude increasing exponentially from the initial value to a preset upper limit is used as a bifurcation detection perturbation. The bifurcation detection perturbation is used to track the stability changes of the system equilibrium point and locate the bifurcation point. Based on the natural oscillation frequency and damping ratio, pulse train perturbations with adjustable pulse width and interval are generated as chaotic edge perturbations to detect the system's sensitivity to sudden disturbances and the conditions for the appearance of chaotic attractors. The intrinsic mode perturbation, bifurcation detection perturbation, and chaotic edge perturbation are combined in a preset order to form a multi-mode phase perturbation sequence, and the parameters of subsequent perturbations are dynamically adjusted according to the response analysis results of the previous round to form an adaptive closed-loop detection.

[0008] Optionally, the variational mode decomposition of the synchronization error response sequence to obtain K finite-bandwidth eigenmode functions includes: Construct a constrained variational problem that minimizes the sum of the estimated bandwidths of each mode and equals the sum of all modes to the original signal; By introducing a quadratic penalty factor and Lagrange multipliers, the constrained variational problem is transformed into an unconstrained variational problem, which is then solved iteratively in the frequency domain using the alternating direction multiplier method. During the iteration process, the number of modes K is dynamically adjusted according to the degree of overlap between the center frequency of the current mode and the frequency band of the adjacent modes, until the center frequencies of each mode are reasonably separated and there is no mode aliasing. Output K intrinsic mode functions, each representing a dynamic mode component of the system, and sorted by center frequency from low to high.

[0009] Optionally, the step of performing a Hilbert transform on each intrinsic mode function to obtain the Hilbert spectrum of the spatiotemporal synchronization loop under test includes: Perform a Hilbert transform on each intrinsic mode function to construct an analytic signal; Calculate the instantaneous amplitude and instantaneous phase from the analyzed signal; The instantaneous frequency is obtained by differentiating the instantaneous phase. The instantaneous amplitudes and frequencies of all modes are mapped onto the time-frequency plane, and the Hilbert spectrum is obtained by superposition.

[0010] Optionally, constructing a parameterized reduced-order model of the system based on the intrinsic mode functions includes: A damping term model is established by fitting nonlinear damping coefficients from the instantaneous amplitude decay patterns of each intrinsic mode function; A Duffing-type frequency response model is established by fitting nonlinear stiffness coefficients from the instantaneous frequency-instantaneous amplitude relationship of each intrinsic mode function. Calculate the cross-correlation function between the perturbation input signal and the response output signal, and use the time delay corresponding to the peak value of the cross-correlation as the estimate of the system time delay; A third-order nonlinear state-space model is constructed using nonlinear damping coefficient, nonlinear stiffness coefficient, and system time delay as state variables. The third-order nonlinear state-space model is determined as the parameterized reduced-order model.

[0011] Optionally, the prediction of the stability domain boundary and instability bifurcation point of the system in the parameter space using the numerical extension algorithm includes: A parameter plane is constructed by using loop gain and crystal oscillator aging rate as continuously varying parameters and fixing other parameters as the current identification values. The pseudo-arc length extension algorithm is used to track the trajectory of the system equilibrium point as the parameters change, and the eigenvalues ​​of the Jacobian matrix of the equilibrium point are calculated at each extension step. When the real part of the eigenvalue changes from negative to positive, record this point as a bifurcation point and mark it on the parametric plane; Connecting the bifurcation points under different initial conditions forms the boundary curve of the stable region, with the inside of the curve being the stable region and the outside being the unstable region; Based on the stable and unstable regions, a two-dimensional stability domain map is generated, and the stability margin of the current operating point from the boundary is calculated and presented as a percentage. Based on the stability margin and the two-dimensional stability domain map, a synchronization stability boundary map is generated.

[0012] On the other hand, embodiments of this application provide a time synchronization network testing apparatus, the apparatus comprising: The acquisition module is used to acquire the state equation and system parameters of the preset nonlinear dynamic disturbance in response to the access command of the device under test. The capture module is used to inject the mode-matched phase perturbation sequence into the device under test and continuously capture the synchronization error response sequence output by the device under test under perturbation excitation. The phase perturbation sequence is constructed based on the state equation and system parameters. The decomposition module is used to perform variational mode decomposition on the synchronization error response sequence to obtain K eigenmode functions with finite bandwidth, each eigenmode function corresponding to a dynamic mode component of the test system; The transformation module is used to perform Hilbert transformation on each intrinsic mode function to obtain the Hilbert spectrum of the tested spatiotemporal synchronization loop. The Hilbert spectrum reflects the variation of the system frequency components with time. The generation module is used to construct a parameterized reduced-order model of the system based on the intrinsic mode functions, predict the stable domain boundary and instability bifurcation point of the system in the parameter space through a numerical extension algorithm, and generate a synchronization stability boundary map of the device under test, so as to verify the bidirectional synchronization capability of the device under test based on the synchronization stability boundary map.

[0013] In another aspect, embodiments of this application provide an electronic device, the device including: a processor and a memory storing computer program instructions; When the processor executes the computer program instructions, it implements the time synchronization network testing method as described in the first aspect.

[0014] In another aspect, embodiments of this application provide a computer storage medium storing computer program instructions, which, when executed by a processor, implement the time synchronization network testing method as described in the first aspect.

[0015] The time synchronization network testing method, apparatus, device, and computer storage medium of this application can deeply analyze the intrinsic dynamic modes of the spatiotemporal synchronization loop of the device under test by acquiring the nonlinear dynamic model of the device under test, injecting a modally matched phase perturbation sequence, and performing variational mode decomposition and Hilbert transform on the synchronization error response sequence. Furthermore, by constructing a parameterized reduced-order model and combining it with a numerical extension algorithm, the stability domain boundary and instability bifurcation point of the system in the parameter space can be predicted, generating a synchronization stability boundary map. Therefore, the bidirectional synchronization capability of the device under test can be verified intuitively and quantitatively, overcoming the limitations of traditional manual measurement, which is inefficient, has large errors, and cannot comprehensively evaluate the complex synchronization performance of time synchronization network (TSN) devices. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating a time synchronization network testing method provided in an embodiment of this application; Figure 2 This is a structural block diagram of a time synchronization network testing device provided in an embodiment of this application; Figure 3 This is a schematic diagram of the structure of an electronic device provided in an embodiment of this application; Detailed Implementation The features and exemplary embodiments of various aspects of this application will be described in detail below. To make the objectives, technical solutions, and advantages of this application clearer, the application will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are only intended to explain this application and not to limit it. For those skilled in the art, this application can be implemented without some of these specific details. The following description of the embodiments is merely to provide a better understanding of this application by illustrating examples.

[0017] It should be noted that, in this document, relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising..." does not exclude the presence of additional identical elements in the process, method, article, or apparatus that includes the element.

[0018] For ease of understanding, the following explains some key terms in this embodiment: State equations and system parameters of nonlinear dynamic disturbances: State equations are mathematical models describing the dynamic behavior of a system, usually represented by a set of differential or difference equations, whose state variables evolve over time. System parameters are constants or coefficients in the state equations, which determine the specific characteristics and response of the system. Nonlinear dynamic disturbances refer to factors existing inside or outside the system that have a nonlinear impact on the system's synchronization performance, such as the nonlinear characteristics of a crystal oscillator or the nonlinear response of a loop filter.

[0019] Phase perturbation sequence for modal matching: A phase perturbation sequence refers to a small, time-varying disturbance signal applied to the clock phase of the device under test (DUT). Modal matching means that the characteristics (such as frequency, amplitude, and waveform) of the perturbation sequence correspond to the inherent dynamic modes within the DUT, aiming to effectively elicit a specific response from the system.

[0020] Synchronization error response sequence: Synchronization error refers to the phase or frequency deviation between the clock signal output by the device under test (DUT) and the ideal reference clock signal. When the DUT is excited by a phase perturbation sequence, its internal synchronization loop will make corresponding adjustments. The record of the resulting synchronization error changing over time is the synchronization error response sequence.

[0021] Variational Mode Decomposition (VMD): Variational mode decomposition is an adaptive, non-recursive signal decomposition method. It decomposes a complex non-stationary signal into several eigenmode functions with different center frequencies and finite bandwidths, each eigenmode function representing an inherent oscillation mode of the signal.

[0022] Intrinsic Mode Function (IMF): The intrinsic mode function is the output of variational mode decomposition. Each IMF is a narrowband signal with good locality and time-frequency characteristics, which can reflect a specific dynamic mode component in the original signal.

[0023] Hilbert Transform: The Hilbert transform is a linear transform that converts a real signal into an analytic signal. Through the Hilbert transform, information such as instantaneous amplitude, instantaneous phase, and instantaneous frequency can be extracted from a real signal, thereby analyzing the signal's local time-frequency characteristics.

[0024] Hilbert spectrum: A Hilbert spectrum is a two-dimensional spectrum obtained by mapping the instantaneous amplitude and instantaneous frequency of a signal onto the time-frequency plane. It can intuitively show the variation of the system's frequency components over time, revealing the non-stationary characteristics and energy distribution of the signal.

[0025] Parametric reduced-order models: Parametric reduced-order models refer to system models constructed by simplifying mathematical expressions or reducing the number of state variables while retaining the main dynamic characteristics of the system. These models typically include adjustable parameters and can reflect the system's behavior under different operating conditions.

[0026] Numerical continuation algorithms are numerical methods used to track the trajectory of equilibrium points or periodic solutions in nonlinear systems as parameters change. They work by progressively changing system parameters and utilizing a predictive-correction mechanism to calculate new equilibrium points or periodic solutions, thereby plotting the system's bifurcation diagram.

[0027] Stability domain boundary and instability bifurcation point: The stability domain boundary refers to the critical line in the system's parameter space where the system transitions from a stable state to an unstable state. An instability bifurcation point is the point where the stability of the equilibrium point or periodic solution changes during parameter variations; for example, it may change from a stable equilibrium point to an unstable equilibrium point, or a new periodic solution may appear.

[0028] Synchronization stability boundary map: A synchronization stability boundary map is a graphical representation in two-dimensional or multi-dimensional parameter space that identifies the stable and unstable regions of a system. It visually demonstrates the range of stability within which the device under test maintains synchronization under different combinations of system parameters.

[0029] Two-way synchronization capability: Two-way synchronization capability refers to the ability of a device in a time-sensitive network to not only receive and synchronize with the master clock, but also act as the master clock to provide synchronization signals to other devices, while maintaining its own stability and accuracy in the process.

[0030] To address the problems of existing technologies, this application provides a method, apparatus, device, and medium for testing time synchronization networks. In this application, by acquiring the nonlinear dynamic model of the device under test (DUT), injecting a mode-matched phase perturbation sequence, and performing variational mode decomposition and Hilbert transform on the synchronization error response sequence, the intrinsic dynamic modes of the DUT's spatiotemporal synchronization loop can be analyzed in depth. Furthermore, by constructing a parameterized reduced-order model and combining it with a numerical extension algorithm, the stability domain boundary and instability bifurcation point of the system in the parameter space can be predicted, generating a synchronization stability boundary map. Thus, the bidirectional synchronization capability of the DUT can be verified intuitively and quantitatively, overcoming the limitations of traditional manual measurements, which are inefficient, have large errors, and cannot comprehensively evaluate the complex synchronization performance of time synchronization network (TSN) devices.

[0031] The time synchronization network testing method provided in the embodiments of this application will be introduced first below.

[0032] Figure 1 A flowchart illustrating a time synchronization network testing method according to an embodiment of this application is shown. Figure 1As shown, the time synchronization network testing method may include S101-S105: S101, in response to the access command of the device under test, obtains the state equation and system parameters of the preset nonlinear dynamic disturbance.

[0033] In this embodiment, the state equation and system parameters can be obtained by manual input or by selecting a general model similar to the type of device under test from a preset database. For example, the operator can manually configure parameters such as the gain and filter type of the clock loop according to the specification manual of the device under test, and select a general nonlinear oscillator model as the state equation.

[0034] S102, inject the modally matched phase perturbation sequence into the device under test, and continuously capture the synchronization error response sequence output by the device under test under perturbation excitation.

[0035] In this embodiment, the phase perturbation sequence is constructed based on the state equation and system parameters. It can be injected into the device under test (DUT) by pre-setting a single-form phase perturbation sequence, such as a sine wave with a fixed frequency and amplitude. The parameters of this perturbation sequence can be calculated based on the initially obtained state equation and system parameters to ensure that it can produce an observable response to the DUT. During the perturbation injection, the phase difference between the clock signal output by the DUT and the reference clock signal is continuously recorded using a high-precision time interval meter or oscilloscope, forming a synchronization error response sequence.

[0036] S103, variational mode decomposition is performed on the synchronization error response sequence to obtain K eigenmode functions with finite bandwidth.

[0037] In this embodiment, each intrinsic mode function (EMF) corresponds to a dynamic mode component of the test system, and the captured synchronization error response sequence can be processed using a variational mode decomposition algorithm. In implementation, a fixed number K of EMFs can be preset, for example, K can be set to 3 or 4 based on experience or preliminary analysis. The decomposition process iteratively optimizes the original signal into these EMFs with finite bandwidth, each representing a major oscillation mode or dynamic component in the synchronization loop of the device under test.

[0038] S104. Perform Hilbert transform on each intrinsic mode function to obtain the Hilbert spectrum of the spatiotemporal synchronization loop under test.

[0039] In this embodiment, the Hilbert spectrum reflects the variation of the system's frequency components over time. A Hilbert transform can be performed on each decomposed intrinsic mode function to obtain its instantaneous frequency and amplitude information. By visualizing the relationship between these instantaneous frequencies and amplitudes over time, a time-frequency spectrum can be generated. This spectrum can preliminarily reflect how the frequency components of different dynamic modes in the tested spatiotemporal synchronization loop evolve over time.

[0040] S105. Based on the intrinsic mode functions, a parameterized reduced-order model of the system is constructed. The stability domain boundary and instability bifurcation point of the system in the parameter space are predicted by the numerical extension algorithm. The synchronization stability boundary map of the device under test is generated to verify the bidirectional synchronization capability of the device under test based on the synchronization stability boundary map.

[0041] In the embodiments of this application, the dynamic characteristics of the intrinsic mode functions can be approximated by curve fitting the time-domain waveforms, such as using polynomial fitting or exponential decay fitting. Based on these fitting results, a simplified second- or third-order linear model can be constructed, which contains a few key parameters to characterize the main oscillations and decay behavior of the system, thus forming a parameterized reduced-order model.

[0042] The stability domain boundary and instability bifurcation point of a system in parameter space are predicted using a numerical extension algorithm. Specifically, one or two key parameters of the system, such as loop gain, can be selected as variable parameters. Discrete sampling is performed within the range of these parameters, and stability analysis is conducted on the reduced-order model at each sampling point, for example, by calculating eigenvalues, to determine whether the system is stable. When the system stability changes, the corresponding parameter values ​​are recorded to roughly determine the stability domain boundary and instability bifurcation point.

[0043] A synchronization stability boundary map of the device under test (DUT) is generated to verify its bidirectional synchronization capability. Specifically, based on the predicted stability domain boundaries and instability bifurcation points, a two-dimensional map can be plotted, with one axis representing the selected key parameters and the other axis representing the system's stability state. This map visually displays the stable and unstable regions of the system under different parameter conditions. By mapping the device's actual operating parameter points onto this map, its synchronization stability can be preliminarily determined, and its ability to maintain bidirectional synchronization in different operating modes (e.g., as a master clock or slave clock) can be evaluated.

[0044] In this embodiment, by acquiring the nonlinear dynamic model of the device under test (DUT), injecting a modally matched phase perturbation sequence, and performing variational mode decomposition and Hilbert transform on the synchronization error response sequence, the intrinsic dynamic modes of the DUT's spatiotemporal synchronization loop can be analyzed in depth. Furthermore, by constructing a parameterized reduced-order model and combining it with a numerical extension algorithm, the stability domain boundary and instability bifurcation point of the system in the parameter space can be predicted, generating a synchronization stability boundary map. Thus, the bidirectional synchronization capability of the DUT can be verified intuitively and quantitatively, overcoming the limitations of traditional manual measurements, which are inefficient, have large errors, and cannot comprehensively evaluate the complex synchronization performance of TSN devices.

[0045] In some other embodiments, S101 may specifically include: In response to the access command of the device under test, the device under test is identified in terms of type and protocol version, and the corresponding nonlinear dynamic system template is loaded. Read the initial state equation form from the system template, where the state variables include at least the phase difference, frequency difference, and digital loop filter state, and the system parameters include at least the loop gain, crystal oscillator aging factor, and filter coefficient. The system continuously captures the synchronization error sequence output by the device under test in a non-perturbation state, and estimates the initial values ​​of the system parameters online based on the subspace identification method, so that the state equation matches the actual dynamic characteristics of the device under test. The estimated system parameters are written to the parameter storage area in real time for subsequent perturbation sequence generation modules to use.

[0046] In this embodiment, firstly, in response to the access command from the device under test (DUT), the system identifies the type and protocol version of the DUT and loads the corresponding nonlinear dynamic system template. When the test system receives the access command from the DUT, it can trigger a device information query or handshake protocol. By parsing the identifiers sent by the device (such as MAC address, device model, firmware version, etc.) or through manual input by the user, the system can identify the specific type of the DUT (e.g., PTP master clock, slave clock, boundary clock, transparent clock, etc.) and the version of the time synchronization protocol it follows (e.g., IEEE 1588v2, NTP, etc.). Based on these identification results, the system retrieves and loads the most matching nonlinear dynamic system template from a pre-established template library. This template contains a representative state equation structure and parameter range preset for a specific device type and protocol version, ensuring that the initial model matches the basic properties of the DUT.

[0047] Next, the initial state equation form is read from the system template. The state variables include at least the phase difference, frequency difference, and digital loop filter state, while the system parameters include at least the loop gain, crystal aging factor, and filter coefficients. The loaded system template not only provides a general state equation framework but also specifically defines the core state variables and system parameters within that framework. For example, for a typical clock synchronization loop, its state variables typically include the phase difference between the device under test (DUT) and the reference clock (reflecting synchronization accuracy), the frequency difference (reflecting clock drift), and the internal state variables of the digital loop filter (such as a PID controller) (reflecting the controller's cumulative error or integral term). The system parameters cover key factors affecting the loop's dynamic response, such as the loop gain (determining loop response speed and stability), the crystal aging factor (reflecting the long-term drift characteristics of the clock source), and the filter coefficients (determining the filter's dynamic characteristics). The initial forms and ranges of these variables and parameters are preset from the template, laying the foundation for subsequent parameter estimation.

[0048] Subsequently, the synchronization error sequence output by the device under test (DUT) is continuously captured under no-perturbation conditions. Initial values ​​of system parameters are estimated online based on the subspace identification method, ensuring the state equation matches the actual dynamic characteristics of the DUT. Without injecting any external perturbation signals into the DUT, the test system continuously monitors and captures the synchronization error sequence output by the DUT under normal operating conditions. These error sequences reflect the dynamic response of the device under natural disturbances (such as internal noise, environmental changes, etc.). To extract the intrinsic dynamic characteristics of the device from these observations, the subspace identification method can be used. Subspace identification is an advanced system identification technique that directly estimates the state-space model parameters of the system in the subspace by analyzing the covariance matrix of the input and output data, without pre-setting the model structure. Specifically, algorithms such as N4SID, MOESP, or CVA can be used to take the captured synchronization error sequence as output data and, combined with the assumed no-perturbation input (or processed through an autoregressive model), estimate the initial values ​​of system parameters such as loop gain, crystal oscillator aging coefficient, and filter coefficients online. This online estimation solves the problem of mismatch between preset parameters and actual equipment, ensuring that the constructed state equation can match the actual dynamic characteristics of the equipment under test to the greatest extent, laying the foundation for subsequent accurate testing.

[0049] Finally, the estimated system parameters are written to the parameter storage area in real time for use by the subsequent perturbation sequence generation module. The system parameters (such as loop gain, crystal oscillator aging coefficient, and filter coefficients) estimated online using the subspace identification method accurately reflect the actual dynamic characteristics of the device under test. To ensure that subsequent test steps can utilize these accurate parameters, the system stores these estimates in real time in a dedicated parameter storage area (e.g., a specific region in memory, a database, or a configuration file). This storage area serves as a shared resource and can be accessed by other modules in the test system. In particular, the module responsible for generating the mode-matched phase perturbation sequence reads the latest identified system parameters from this parameter storage area when constructing the perturbation signal, thereby ensuring that the injected perturbation sequence accurately excites the specific dynamic modes of the device under test, achieving more precise testing.

[0050] In this embodiment, a matching system template can be loaded based on the specific type and protocol version of the device under test (DUT), and system parameters that highly match the actual dynamic characteristics of the DUT can be estimated online using a subspace identification method. This significantly improves the accuracy and specificity of the initial state equations and system parameters, avoiding model mismatch problems that may occur when using general preset parameters. Therefore, the subsequent construction of the mode-matching phase perturbation sequence will be more accurate and the excitation of the dynamic modes of the DUT will be more effective, thus enabling the capture of the synchronization error response sequence, variational mode decomposition, Hilbert transform, and the construction of the parameterized reduced-order model to be based on more reliable and realistic data. Ultimately, the generated synchronization stability boundary map can more accurately reflect the bidirectional synchronization capability of the DUT, providing more solid technical support for the performance evaluation and optimization of the device.

[0051] In some other embodiments, prior to S102, the method may further include: Calculate the system's natural oscillation frequency and damping ratio based on the loop gain and crystal oscillator aging coefficient in the current system parameters; A sinusoidal perturbation sequence is generated based on the natural oscillation frequency and damping ratio, which is used as an intrinsic mode perturbation to excite the main oscillation mode of the system. A step perturbation sequence based on the natural oscillation frequency and damping ratio generation amplitude increasing exponentially from the initial value to a preset upper limit is used as a bifurcation detection perturbation. The bifurcation detection perturbation is used to track the stability changes of the system equilibrium point and locate the bifurcation point. Based on the natural oscillation frequency and damping ratio, pulse train perturbations with adjustable pulse width and interval are generated as chaotic edge perturbations to detect the system's sensitivity to sudden disturbances and the conditions for the appearance of chaotic attractors. The intrinsic mode perturbation, bifurcation detection perturbation, and chaotic edge perturbation are combined in a preset order to form a multi-mode phase perturbation sequence. The parameters of subsequent perturbations are dynamically adjusted based on the response analysis results of the previous round to form an adaptive closed-loop detection.

[0052] In this embodiment, loop gain and crystal oscillator aging coefficient are key parameters affecting the dynamic characteristics of the synchronization loop. The natural oscillation frequency and damping ratio are fundamental indicators describing the dynamic response of second-order or higher-order linear systems, determining the system's response speed and stability to disturbances. For nonlinear systems, these parameters remain important in local linearization analysis. Calculations of these parameters are typically based on known system state equations or their linearized approximations, determined through the roots of the characteristic equation. For example, for a typical second-order phase-locked loop (PLL) model, expressions for the natural oscillation frequency and damping ratio related to loop gain and filter parameters can be directly derived from its transfer function or state-space representation. These calculations provide the foundation for subsequently generating targeted perturbation sequences.

[0053] Intrinsic mode perturbation is an excitation method targeting the inherent oscillation modes of a system. The frequency and amplitude of the sinusoidal perturbation sequence can be set based on the calculated natural oscillation frequency and damping ratio. For example, the frequency of the sinusoidal perturbation sequence can be set to be close to the system's natural oscillation frequency to maximize the excitation of the system's main oscillation mode, thereby observing the system's response characteristics under resonance conditions. Its amplitude can start from a small value and gradually increase to avoid causing excessive impact on the system in the initial stage. This perturbation helps identify the system's gain and phase response at specific frequencies, providing crucial data for subsequent mode decomposition and model construction.

[0054] Bifurcation detection perturbations aim to systematically explore the stability boundaries of a system's equilibrium point. The amplitude of a step perturbation sequence increases exponentially, gradually increasing the intensity of the disturbance with small step sizes, thus precisely tracking the transition from stability to instability. When system parameters (such as loop gain or crystal oscillator aging factor) reach a critical value, the system's equilibrium point may bifurcate, causing the system behavior to change from stable synchronization to periodic oscillation or chaos. By observing the system's response to this increasing step perturbation, the critical point of instability, i.e., the bifurcation point, can be identified, which is crucial for evaluating the system's robustness and designing a stable operating range.

[0055] Chaotic edge perturbation is used to probe the complex dynamic behavior of systems in strongly nonlinear regions. Pulse train perturbation, characterized by its high transience and concentrated energy, with adjustable pulse width and interval, can simulate various sudden disturbances. By adjusting the pulse parameters, the system's response to different types of transient shocks can be systematically explored. When a system is on the edge of chaos, even small perturbations can lead to significant changes in system behavior. This perturbation helps detect the system's sensitivity to sudden disturbances and identify the critical conditions that cause the system to enter a chaotic state. This is crucial for assessing the system's robustness to interference in complex network environments and predicting its long-term stability.

[0056] Multimodal phase perturbation sequences, by combining different types of perturbations, achieve comprehensive detection of the system's dynamic characteristics. The preset sequence can be optimized according to the test objectives; for example, intrinsic mode perturbation can be performed first to obtain a linear response, followed by bifurcation detection perturbation to find stability boundaries, and finally chaotic edge perturbation to explore extreme nonlinear behavior. More importantly, this combination is adaptive. This means that the test system can adjust the parameters (such as amplitude, frequency, pulse width, and interval) of subsequent perturbation sequences in real time based on the analysis results of the synchronization error response sequence output by the device under test under the previous round of perturbation excitation (e.g., whether bifurcation points were detected, whether chaotic signs appeared, etc.). This adaptive closed-loop detection mechanism makes the testing process more efficient and intelligent, focusing on key areas of system behavior, avoiding invalid detection, and ensuring the acquisition of the most comprehensive information on the system's dynamic characteristics within a limited test time.

[0057] First, the natural oscillation frequency and damping ratio of the device under test (DUT) are calculated based on the current system parameters, providing a precise physical basis for the generation of subsequent perturbation sequences. Building upon this, this application introduces three different types of perturbation sequences: intrinsic mode perturbation, bifurcation detection perturbation, and chaotic edge perturbation. Intrinsic mode perturbation effectively excites the system's main oscillation mode, revealing its linear or weakly nonlinear response characteristics. Bifurcation detection perturbation systematically tracks the stability changes of the system's equilibrium point through an exponentially increasing step sequence, thereby accurately locating the instability bifurcation point. Chaotic edge perturbation detects the system's sensitivity to sudden disturbances and the conditions for the appearance of chaotic attractors through pulse trains with adjustable pulse width and interval. These multimodal perturbation sequences are combined in a preset order, and combined with an adaptive closed-loop detection mechanism, the parameters of subsequent perturbations can be dynamically adjusted based on the results of the previous response analysis. This design allows the testing process to comprehensively excite various dynamic modes of the DUT, from linear response to nonlinear bifurcation and even chaotic edges, avoiding the risk of missing key dynamic behaviors. Meanwhile, the adaptive closed-loop detection mechanism significantly improves testing efficiency and intelligence, concentrating testing resources on key areas of system behavior to ensure the acquisition of the most comprehensive information on system dynamic characteristics within a limited testing time. By detecting bifurcation points and chaotic edges, it can more accurately assess the stability boundaries of the device under test under different disturbance intensities and its sensitivity to sudden disturbances. This provides high-quality input data for subsequent variational mode decomposition, Hilbert transform, and parameterized reduced-order model construction, thereby improving the accuracy and reliability of the synchronization stability boundary map and ultimately achieving comprehensive, accurate, and efficient verification of the bidirectional synchronization capability of the device under test.

[0058] In some other embodiments, S103 may include: Construct a constrained variational problem that minimizes the sum of the estimated bandwidths of each mode and equals the sum of all modes to the original signal; By introducing a quadratic penalty factor and Lagrange multipliers, the constrained variational problem is transformed into an unconstrained variational problem, which is then solved iteratively in the frequency domain using the alternating direction multiplier method. During the iteration process, the number of modes K is dynamically adjusted according to the degree of overlap between the center frequency of the current mode and the frequency band of the adjacent modes, until the center frequencies of each mode are reasonably separated and there is no mode aliasing. Output K intrinsic mode functions, each representing a dynamic mode component of the system, and sorted by center frequency from low to high.

[0059] In this embodiment, a constrained variational problem is first constructed to minimize the sum of the estimated bandwidths of each mode, and the sum of all modes equals the original signal. The core of this constrained variational problem lies in its aim to decompose the complex synchronization error response sequence into a series of eigenmode functions with specific physical meanings. Each mode should be as compact as possible, meaning its frequency components are concentrated within a narrow bandwidth, while ensuring that the sum of all decomposed modes can accurately reconstruct the original signal, thus guaranteeing the integrity and accuracy of the decomposition. Specifically, this is typically achieved by defining an energy functional, which consists of two parts: the H1 norm of the demodulated signal for each mode (used to measure bandwidth), and the L2 norm of the difference between the sum of all modes and the original signal (used to measure reconstruction error). By minimizing this functional, a mathematical balance between bandwidth minimization and signal reconstruction can be achieved.

[0060] Secondly, to effectively solve the aforementioned constrained variational problem, a quadratic penalty factor and Lagrange multipliers are introduced to transform the constrained variational problem into an unconstrained variational problem, which is then solved iteratively in the frequency domain using the alternating direction multiplier method. By introducing the quadratic penalty factor and Lagrange multipliers, the original complex constrained problem can be transformed into a series of more manageable subproblems. The alternating direction multiplier method (ADMM) provides an efficient iterative optimization framework that gradually approximates the optimal solution by alternately updating the modes, center frequencies, and Lagrange multipliers in the frequency domain. Solving in the frequency domain leverages the advantages of the Fourier transform, converting convolution operations in the time domain into multiplications in the frequency domain, thereby significantly improving computational efficiency and stability. In practice, a Fourier transform is first performed on the original signal and the modes to be decomposed. Then, in the frequency domain, the center frequency, the mode itself, and the Lagrange multipliers of each mode are iteratively updated. Each iteration involves a Wiener-filter-like update of the current mode to minimize its bandwidth and adjusting the Lagrange multipliers to satisfy the reconstruction constraints.

[0061] Furthermore, during the iteration process, the number of modes K is dynamically adjusted based on the overlap between the center frequency of the current mode and the frequency bands of adjacent modes, until the center frequencies of each mode are reasonably separated and there is no mode aliasing. This adaptive adjustment mechanism is one of the key innovations of this implementation method. Traditional variational mode decomposition methods often require a pre-set number of modes K, which is difficult to guarantee the decomposition effect when facing unknown or complex signals. By monitoring the frequency relationship and overlap between modes in real time, the system can intelligently increase or decrease the number of modes K, ensuring that each intrinsic mode function corresponds to an independent dynamic mode component, effectively avoiding under-decomposition or over-decomposition caused by improper K value setting, as well as the problem of mode aliasing. For example, the distance between the center frequencies of adjacent modes can be calculated, or the degree of overlap between their frequency bands can be evaluated. If the center frequencies of two modes are found to be very close, or their frequency bands are highly overlapping, this may indicate that the K value is too large, with redundant modes or mode aliasing, and in this case, K can be reduced. Conversely, if there is still a large amount of energy in the signal that has not been effectively decomposed, and the interval between the existing modes is large, it may be necessary to increase K.

[0062] Finally, K eigenmode functions are output, each representing a dynamic mode component of the system, and ordered from low to high center frequency. The eigenmode functions obtained after the precise decomposition and adaptive adjustment are pure representations of different dynamic modes in the original synchronization error response sequence. Sort by center frequency not only makes subsequent analysis and visualization more intuitive but also provides a structured data foundation for further Hilbert transform and parameterized reduced-order model construction. After iterative convergence, the K modes in the frequency domain are transformed back to the time domain using an inverse Fourier transform to obtain the final K eigenmode functions. Then, the instantaneous frequency of each mode or its center frequency in the frequency domain is calculated, and the modes are arranged in ascending or descending order based on these frequency values.

[0063] In this embodiment, accurate and adaptive variational mode decomposition (VMD) of the synchronization error response sequence can be achieved. By constructing a constrained variational problem and combining a quadratic penalty factor and Lagrange multipliers, the decomposition process maximizes bandwidth compression for each mode while ensuring signal reconstruction integrity. The iterative solution in the frequency domain using the alternating direction multiplier method significantly improves computational efficiency and robustness. Crucially, the mechanism of dynamically adjusting the number of modes K during iteration adaptively optimizes the decomposition results based on the frequency relationships and overlap between modes, effectively avoiding mode aliasing and ensuring reasonable separation of intrinsic mode functions, thereby accurately capturing the inherent dynamic modes in the synchronization error response sequence of the device under test. This accurate and adaptive VMD provides high-quality input for subsequent Hilbert transform analysis and the construction of parameterized reduced-order models, greatly improving the accuracy and reliability of the dynamic behavior analysis of the spatiotemporal synchronization loop under test, making the generation of synchronization stability boundary maps more accurate, and thus more effectively verifying the bidirectional synchronization capability of the device under test.

[0064] In some other embodiments, S104 may include: Perform a Hilbert transform on each intrinsic mode function to construct an analytic signal; Calculate the instantaneous amplitude and instantaneous phase from the analytic signal; The instantaneous frequency is obtained by differentiating the instantaneous phase. The instantaneous amplitudes and frequencies of all modes are mapped onto the time-frequency plane, and the Hilbert spectrum is obtained by superposition.

[0065] In this embodiment, the Hilbert transform is a mathematical operation that converts a real signal into a complex signal, resulting in an analytic signal. The real part of this analytic signal is the original signal, and the imaginary part is the Hilbert transform of the original signal. By constructing an analytic signal, the amplitude and phase information of the signal can be unified into a complex form, facilitating the subsequent extraction of instantaneous features. In practical implementation, the eigenmode functions can be transformed to the frequency domain using a Fast Fourier Transform (FFT), then the negative frequency components can be set to zero, and the Hilbert transform result can be obtained through an Inverse Fourier Transform (IFFT), thereby constructing the analytic signal.

[0066] Secondly, the instantaneous amplitude and instantaneous phase are calculated from the analytic signal. For the constructed analytic signal, its magnitude is the instantaneous amplitude, reflecting the energy or intensity of the signal at a certain moment; its argument is the instantaneous phase, describing the change of the signal's phase over time. These instantaneous characteristics are fundamental to understanding the dynamic behavior of non-stationary signals, revealing the energy distribution and phase state of the eigenmode functions at different time points.

[0067] Next, the instantaneous frequency is obtained by differentiating the instantaneous phase with respect to time. The instantaneous frequency is the derivative of the instantaneous phase with respect to time; it accurately describes the frequency components of a signal at each moment, especially important for non-stationary signals whose frequencies may dynamically change over time. In discrete data processing, the instantaneous frequency can be approximated by performing a difference operation on the instantaneous phase sequence, such as using the central difference method. Accurate calculation of the instantaneous frequency is crucial for revealing the frequency evolution of the system's dynamic modes and is a core element in constructing a time-frequency spectrum.

[0068] Finally, the instantaneous amplitudes and frequencies of all modes are mapped to the time-frequency plane, and the Hilbert spectrum is obtained by superposition. The Hilbert spectrum is a time-frequency representation that distributes the energy or amplitude of a signal on the time-frequency plane. For each intrinsic mode function (EMF), its calculated instantaneous amplitude and instantaneous frequency constitute a point on the time-frequency plane, and its corresponding "intensity" is the instantaneous amplitude. By collecting the instantaneous amplitude and instantaneous frequency information of all K EMFs and superimposing them on the time-frequency plane, an energy density map can be formed. This map is the Hilbert spectrum, which visually displays the intensity distribution of different frequency components in the measured spatiotemporal synchronization loop over time, thus reflecting the dynamic characteristics of the system and the distribution of energy in different frequency modes.

[0069] In this embodiment, firstly, a Hilbert transform is performed on each intrinsic mode function to construct an analytic signal, laying the mathematical foundation for subsequent extraction of instantaneous features. Next, the instantaneous amplitude and instantaneous phase are directly calculated from the analytic signal, ensuring the accuracy of these features. By differentiating the instantaneous phase, the instantaneous frequency can be obtained, accurately capturing the change in signal frequency over time, which is crucial for analyzing non-stationary synchronization error responses. Finally, the instantaneous amplitudes and frequencies of all modes are mapped and superimposed onto the time-frequency plane, generating a Hilbert spectrum that intuitively and meticulously reflects how each frequency component in the tested spatiotemporal synchronization loop evolves over time. This detailed spectrum not only clearly reveals the energy distribution and dominant frequency modes of the system at different times but also provides more refined and reliable input data for subsequent construction of parameterized reduced-order models based on the intrinsic mode functions. This improves the accuracy of predicting the system's stability domain boundaries and instability bifurcation points, making the verification of the bidirectional synchronization capability of the tested equipment more comprehensive and reliable.

[0070] In some other embodiments, S105 may include: A damping term model is established by fitting nonlinear damping coefficients from the instantaneous amplitude decay patterns of each intrinsic mode function; A Duffing-type frequency response model is established by fitting nonlinear stiffness coefficients from the instantaneous frequency-instantaneous amplitude relationship of each intrinsic mode function. Calculate the cross-correlation function between the perturbation input signal and the response output signal, and use the time delay corresponding to the peak value of the cross-correlation as the estimate of the system time delay; A third-order nonlinear state-space model is constructed using nonlinear damping coefficient, nonlinear stiffness coefficient, and system time delay as state variables. The third-order nonlinear state-space model is determined to be a parameterized reduced-order model.

[0071] Specifically, the nonlinear damping coefficient refers to a parameter that describes the nonlinear relationship between energy dissipation and factors such as vibration amplitude and velocity during system vibration. Establishing a damping term model aims to accurately describe the energy decay of the system over time. Damping characteristics related to amplitude can be extracted from the instantaneous amplitude time series of each eigenmode function obtained from variational mode decomposition using nonlinear regression analysis methods, such as exponential decay fitting, polynomial fitting, or neural network-based fitting. For example, for a mode, its instantaneous amplitude may decay nonlinearly over time. Fitting can yield a functional relationship that describes how the damping force changes with vibration amplitude or velocity, thus forming a damping term model that reflects the energy dissipation mechanism of the actual system.

[0072] Meanwhile, the nonlinear stiffness coefficient refers to a parameter that exhibits a nonlinear relationship between the system's restoring force and displacement. Establishing a Duffing-type frequency response model aims to capture the characteristics of the system's vibration frequency changing with vibration amplitude. It allows analysis of the relationship between the instantaneous frequency and instantaneous amplitude of each eigenmode function. In many nonlinear systems, the vibration frequency is not constant but shifts with the increase or decrease of vibration amplitude; this phenomenon is called frequency-amplitude coupling. By curve fitting the instantaneous frequency-instantaneous amplitude data, such as using polynomial fitting or a specific form of fitting based on the Duffing equation, the nonlinear stiffness coefficient can be obtained. The Duffing-type frequency response model is a typical nonlinear oscillator model that can effectively describe systems with nonlinear restoring forces. Its characteristic feature is that the frequency response curve may exhibit multivaluedness and jump phenomena, thus more accurately reflecting the nonlinear dynamic behavior of the measured spatiotemporal synchronization loop.

[0073] Furthermore, system time delay refers to the time delay required for the system's output signal to respond after an input signal is applied. To accurately estimate this key parameter, the cross-correlation function between the perturbation input signal injected into the device under test (DUT) and the synchronization error response signal output by the DUT can be calculated. The cross-correlation function quantifies the similarity between two signals at different time delays. When the cross-correlation function reaches its peak, the corresponding time delay is considered the best estimate of the system time delay. This method can effectively filter out noise interference and accurately identify the causal time relationship between input and output, which is crucial for understanding and modeling the dynamic response of synchronization loops.

[0074] Building upon this, the state-space model is a mathematical framework that represents the dynamic behavior of a system as a set of first-order differential equations, where state variables describe the internal state of the system. Using the fitted nonlinear damping coefficient, nonlinear stiffness coefficient, and estimated system time delay as core parameters, a third-order nonlinear state-space model is constructed. The third-order model is chosen to fully capture key dynamic variables such as phase difference, frequency difference, and loop filter states in the synchronization loop, as well as their interactions, while incorporating the identified nonlinear characteristics and time delay effects. This model is typically represented as a set of nonlinear differential equations, capable of describing the system's evolution trajectory under different input and initial conditions, providing a precise mathematical foundation for subsequent stability analysis. Finally, the constructed third-order nonlinear state-space model is determined to be a parameterized reduced-order model. Here, "reduced-order" means that the model simplifies the mathematical description of the original complex system while preserving its key nonlinear dynamic characteristics, making it easier to analyze and compute. "Parameterization" emphasizes that the coefficients in the model (such as nonlinear damping coefficients, nonlinear stiffness coefficients, system time delays, etc.) are obtained based on actual measurement data and the fitting process. These parameters can be adjusted and updated according to the specific characteristics of the device under test, thus giving the model a high degree of adaptability and accuracy. This model provides a solid foundation for subsequent prediction of the system's stability domain boundary and instability bifurcation points using numerical extension algorithms.

[0075] By accurately fitting the nonlinear damping coefficient from the instantaneous amplitude decay law of the intrinsic mode functions and establishing a Duffing-type frequency response model from the instantaneous frequency-instantaneous amplitude relationship to capture nonlinear stiffness characteristics, the constructed model can more realistically reflect the energy dissipation and frequency shift phenomena of the device under test under different vibration amplitudes. Simultaneously, by calculating the cross-correlation function between the perturbation input signal and the response output signal to accurately estimate the system time delay, the model ensures accurate characterization of this crucial factor. Finally, using these precisely identified nonlinear damping coefficients, nonlinear stiffness coefficients, and system time delays as state variables, a third-order nonlinear state-space model is constructed and defined as a parameterized reduced-order model, thus providing a mathematical description that effectively simplifies system complexity while fully capturing its nonlinear dynamic characteristics. This significantly improves the accuracy and reliability of predicting the synchronization stability boundary and instability bifurcation point of the device under test, providing more solid and refined technical support for verifying the bidirectional synchronization capability of the device under test.

[0076] In other embodiments, the stability domain boundary and instability bifurcation point of the system in the parameter space are predicted using a numerical extension algorithm, including: A parameter plane is constructed by using loop gain and crystal oscillator aging rate as continuously varying parameters and fixing other parameters as the current identification values. The pseudo-arc length extension algorithm is used to track the trajectory of the system equilibrium point as the parameters change, and the eigenvalues ​​of the Jacobian matrix of the equilibrium point are calculated at each extension step. When the real part of the eigenvalue changes from negative to positive, record this point as a bifurcation point and mark it on the parametric plane; Connecting the bifurcation points under different initial conditions forms the boundary curve of the stable region, with the inside of the curve being the stable region and the outside being the unstable region; Based on the stable and unstable regions, a two-dimensional stability domain map is generated, and the stability margin of the current operating point from the boundary is calculated and presented as a percentage. Based on the stability margin and the two-dimensional stability domain map, a synchronous stability boundary map is generated.

[0077] In this embodiment, loop gain and crystal oscillator aging rate are two key parameters affecting system stability in the time synchronization network. Loop gain determines the response strength of the synchronization loop to error signals, while crystal oscillator aging rate reflects the degree of frequency drift of the local oscillator over time. Selecting these two parameters as continuously varying parameters means that all possible combinations of their values ​​within a certain range will be examined when analyzing system stability. By fixing other parameters as the current identified values, the complexity of the analysis can be simplified, projecting the multidimensional parameter space onto a two-dimensional plane, thereby constructing a parameter plane that can intuitively display system stability. This parameter plane provides a clear analytical foundation for subsequent numerical extension algorithms, allowing changes in system stability to be visually tracked.

[0078] The pseudo-arc-length extension algorithm is a numerical method used to track the trajectory of equilibrium points in a nonlinear system when parameters change continuously. In time-synchronous systems, equilibrium points typically correspond to the phase and frequency differences at which the system reaches a stable synchronized state. This algorithm introduces an arc-length parameter, enabling stable tracking even when the equilibrium trajectory changes direction or bifurcates, avoiding the numerical instability near critical points inherent in traditional parametric extension methods. At each extension step, calculating the eigenvalues ​​of the Jacobian matrix at the current equilibrium point is crucial for determining the system's local stability. The Jacobian matrix describes the linearized dynamic behavior of the system near the equilibrium point, and the sign of the real parts of its eigenvalues ​​directly indicates the stability of the equilibrium point: if all eigenvalues ​​have negative real parts, the equilibrium point is stable; if any eigenvalue has a positive real part, the equilibrium point is unstable.

[0079] During pseudo-arc-length extension, when the real part of the eigenvalues ​​of the Jacobian matrix at the system equilibrium point changes from negative to positive, it signifies the loss of local stability, i.e., a bifurcation has occurred. Bifurcation points are critical points where the system's dynamic behavior undergoes a qualitative change, such as transitioning from a stable synchronous state to a periodic oscillation or chaotic state. By accurately recording these bifurcation points and labeling them on the previously constructed parameter plane, the stable and unstable regions of the system under different parameter combinations can be clearly defined. These bifurcation points constitute the boundaries of the system's stability domain, which is crucial for understanding and predicting system behavior.

[0080] By repeating the pseudo-arc length extension and bifurcation point detection process described above, and starting from different initial parameter values ​​or initial equilibrium points, multiple bifurcation points on the parameter plane can be identified. Connecting these bifurcation points forms the boundary curve of the stability region. This curve divides the parameter plane into two main regions: the inner side of the curve represents the parameter combination in which the system can maintain stable synchronization, i.e., the stable region; the outer side of the curve represents the parameter combination in which the system will lose synchronization or enter other unstable states, i.e., the unstable region. This division intuitively reveals the stability range of the system in the parameter space.

[0081] After identifying the stable and unstable regions, a two-dimensional stability domain map can be generated. This map clearly displays the stability boundaries of the system using loop gain and crystal oscillator aging rate as coordinate axes. Based on this, for the current operating point of the device under test (i.e., the current loop gain and crystal oscillator aging rate) in actual operation, its distance to the nearest stability domain boundary can be calculated and converted into a stability margin. The stability margin is presented as a percentage, intuitively quantifying the "safe distance" from instability from the current operating point, providing a quantitative indicator for evaluating the robustness of the system and predicting potential instability risks.

[0082] Finally, the calculated stability margin information is combined with the two-dimensional stability domain map to generate the final synchronous stability boundary map. This map not only shows the stable and unstable regions of the system in the parameter space, but also further refines the degree of stability within the stable regions through the stability margin.

[0083] The generated comprehensive and intuitive synchronization stability boundary map not only clearly defines the stable and unstable regions of the system under different combinations of loop gain and crystal oscillator aging rate, but also quantifies the "safe distance" of the current operating point from the instability boundary. This allows the verification of the bidirectional synchronization capability of the device under test to go beyond discrete point testing, and to comprehensively evaluate its robustness and stability under continuous parameter changes, thereby significantly improving the accuracy and depth of testing and providing solid data support for system design optimization and fault prediction.

[0084] Based on the time synchronization network testing method provided in the above embodiments, this application also provides a specific implementation of the time synchronization network testing apparatus 200. Please refer to the following embodiments.

[0085] First see Figure 2 The time synchronization network testing device 200 provided in this application embodiment may include: The acquisition module 201 is used to acquire the state equation and system parameters of the preset nonlinear dynamic disturbance in response to the access command of the device under test. The capture module 202 is used to inject the modally matched phase perturbation sequence into the device under test and continuously capture the synchronization error response sequence output by the device under test under perturbation excitation. The phase perturbation sequence is constructed based on the state equation and system parameters. The decomposition module 203 is used to perform variational mode decomposition on the synchronization error response sequence to obtain K eigenmode functions with finite bandwidth, each eigenmode function corresponding to a dynamic mode component of the test system; The transformation module 204 is used to perform Hilbert transformation on each intrinsic mode function to obtain the Hilbert spectrum of the spatiotemporal synchronization loop under test. The Hilbert spectrum reflects the variation of the system frequency components with time. The generation module 205 is used to construct a parameterized reduced-order model of the system based on the intrinsic mode function, predict the stable domain boundary and instability bifurcation point of the system in the parameter space through a numerical extension algorithm, and generate a synchronization stability boundary map of the device under test, so as to verify the bidirectional synchronization capability of the device under test based on the synchronization stability boundary map.

[0086] As an alternative implementation, the acquisition module 201 can also be used for: In response to the access command of the device under test, the device under test is identified in terms of type and protocol version, and the corresponding nonlinear dynamic system template is loaded. Read the initial state equation form from the system template, where the state variables include at least the phase difference, frequency difference, and digital loop filter state, and the system parameters include at least the loop gain, crystal oscillator aging factor, and filter coefficient. The system continuously captures the synchronization error sequence output by the device under test in a non-perturbation state, and estimates the initial values ​​of the system parameters online based on the subspace identification method, so that the state equation matches the actual dynamic characteristics of the device under test. The estimated system parameters are written to the parameter storage area in real time for subsequent perturbation sequence generation modules to use.

[0087] As an alternative implementation, the capture module 202 can also be used for: Calculate the system's natural oscillation frequency and damping ratio based on the loop gain and crystal oscillator aging coefficient in the current system parameters; A sinusoidal perturbation sequence is generated based on the natural oscillation frequency and damping ratio, which is used as an intrinsic mode perturbation to excite the main oscillation mode of the system. A step perturbation sequence based on the natural oscillation frequency and damping ratio generation amplitude increasing exponentially from the initial value to a preset upper limit is used as a bifurcation detection perturbation. The bifurcation detection perturbation is used to track the stability changes of the system equilibrium point and locate the bifurcation point. Based on the natural oscillation frequency and damping ratio, pulse train perturbations with adjustable pulse width and interval are generated as chaotic edge perturbations to detect the system's sensitivity to sudden disturbances and the conditions for the appearance of chaotic attractors. The intrinsic mode perturbation, bifurcation detection perturbation, and chaotic edge perturbation are combined in a preset order to form a multi-mode phase perturbation sequence. The parameters of subsequent perturbations are dynamically adjusted based on the response analysis results of the previous round to form an adaptive closed-loop detection.

[0088] As an alternative implementation, the decomposition module 203 can also be used for: Construct a constrained variational problem that minimizes the sum of the estimated bandwidths of each mode and equals the sum of all modes to the original signal; By introducing a quadratic penalty factor and Lagrange multipliers, the constrained variational problem is transformed into an unconstrained variational problem, which is then solved iteratively in the frequency domain using the alternating direction multiplier method. During the iteration process, the number of modes K is dynamically adjusted according to the degree of overlap between the center frequency of the current mode and the frequency band of the adjacent modes, until the center frequencies of each mode are reasonably separated and there is no mode aliasing. Output K intrinsic mode functions, each representing a dynamic mode component of the system, and sorted by center frequency from low to high.

[0089] As an alternative implementation, the transformation module 204 can also be used for: Perform a Hilbert transform on each intrinsic mode function to construct an analytic signal; Calculate the instantaneous amplitude and instantaneous phase from the analytic signal; The instantaneous frequency is obtained by differentiating the instantaneous phase. The instantaneous amplitudes and frequencies of all modes are mapped onto the time-frequency plane, and the Hilbert spectrum is obtained by superposition.

[0090] As an alternative implementation, the generation module 205 can also be used for: A damping term model is established by fitting nonlinear damping coefficients from the instantaneous amplitude decay patterns of each intrinsic mode function; A Duffing-type frequency response model is established by fitting nonlinear stiffness coefficients from the instantaneous frequency-instantaneous amplitude relationship of each intrinsic mode function. Calculate the cross-correlation function between the perturbation input signal and the response output signal, and use the time delay corresponding to the peak value of the cross-correlation as the estimate of the system time delay; A third-order nonlinear state-space model is constructed using nonlinear damping coefficient, nonlinear stiffness coefficient, and system time delay as state variables. The third-order nonlinear state-space model is determined to be a parameterized reduced-order model.

[0091] As an alternative implementation, the generation module 205 can also be used for: A parameter plane is constructed by using loop gain and crystal oscillator aging rate as continuously varying parameters and fixing other parameters as the current identification values. The pseudo-arc length extension algorithm is used to track the trajectory of the system equilibrium point as the parameters change, and the eigenvalues ​​of the Jacobian matrix of the equilibrium point are calculated at each extension step. When the real part of the eigenvalue changes from negative to positive, record this point as a bifurcation point and mark it on the parametric plane; Connecting the bifurcation points under different initial conditions forms the boundary curve of the stable region, with the inside of the curve being the stable region and the outside being the unstable region; Based on the stable and unstable regions, a two-dimensional stability domain map is generated, and the stability margin of the current operating point from the boundary is calculated and presented as a percentage. Based on the stability margin and the two-dimensional stability domain map, a synchronous stability boundary map is generated.

[0092] Figure 3 A schematic diagram of the hardware structure of the electronic device provided in an embodiment of this application is shown.

[0093] An electronic device may include a processor 301 and a memory 302 storing computer program instructions.

[0094] Specifically, the processor 301 may include a central processing unit (CPU), an application specific integrated circuit (ASIC), or one or more integrated circuits that can be configured to implement the embodiments of this application.

[0095] Memory 302 may include mass storage for data or instructions. For example, and not limitingly, memory 302 may include a hard disk drive (HDD), floppy disk drive, flash memory, optical disk, magneto-optical disk, magnetic tape, or Universal Serial Bus (USB) drive, or a combination of two or more of these. In one instance, memory 302 may include removable or non-removable (or fixed) media, or memory 302 may be non-volatile solid-state memory. Memory 302 may be internal or external to the integrated gateway disaster recovery device.

[0096] In one instance, memory 302 may be read-only memory (ROM). In one instance, the ROM may be a mask-programmed ROM, a programmable ROM (PROM), an erasable PROM (EPROM), an electrically erasable PROM (EEPROM), an electrically rewritable ROM (EAROM), or flash memory, or a combination of two or more of these.

[0097] Memory 302 may include read-only memory (ROM), random access memory (RAM), disk storage media device, optical storage media device, flash memory device, electrical, optical, or other physical / tangible memory storage device. Therefore, generally, memory includes one or more tangible (non-transitory) computer-readable storage media (e.g., memory devices) encoded with software including computer-executable instructions, and when the software is executed (e.g., by one or more processors), it is operable to perform the operations described with reference to the time synchronization network testing method according to the first aspect of this disclosure.

[0098] The processor 301 reads and executes computer program instructions stored in the memory 302 to achieve... Figure 1 A time synchronization network testing method is shown in the embodiment.

[0099] In one example, the electronic device may also include a communication interface 303 and a bus 304. For example, Figure 3 As shown, the processor 301, memory 302, and communication interface 303 are connected through bus 304 and complete communication with each other.

[0100] The communication interface 303 is mainly used to realize communication between various modules, devices, units and / or equipment in the embodiments of this application.

[0101] Bus 304 includes hardware, software, or both, that couples components of an electronic device together. For example, and not as a limitation, the bus may include an Accelerated Graphics Port (AGP) or other graphics bus, an Extended Industry Standard Architecture (EISA) bus, a Front Side Bus (FSB), a Hyper Transport (HT) interconnect, an Industry Standard Architecture (ISA) bus, an Infinite Bandwidth Interconnect, a Low Pin Count (LPC) bus, a memory bus, a Microchannel Architecture (MCA) bus, a Peripheral Component Interconnect (PCI) bus, a PCI-Express (PCI-X) bus, a Serial Advanced Technology Attachment (SATA) bus, a Video Electronics Standards Association Local (VLB) bus, or other suitable buses, or a combination of two or more of these. Where appropriate, bus 304 may include one or more buses. Although specific buses are described and illustrated in embodiments of this application, this application contemplates any suitable bus or interconnect.

[0102] The electronic device can execute the time synchronization network testing method described in the embodiments of this application, thereby achieving a combination Figures 1-2 The method and apparatus for testing time synchronization networks are described.

[0103] Furthermore, in conjunction with the time synchronization network testing methods in the above embodiments, this application embodiment can provide a computer storage medium for implementation. This computer storage medium stores computer program instructions; when these computer program instructions are executed by a processor, they implement any of the time synchronization network testing methods in the above embodiments.

[0104] In an optional embodiment, in conjunction with the time synchronization network testing methods in the above embodiments, this application embodiment can provide a computer program product to implement the method. The instructions in the computer program product are executed by the processor of the electronic device, enabling the electronic device to implement any of the time synchronization network testing methods in the above embodiments.

[0105] It should be clarified that this application is not limited to the specific configurations and processes described above and shown in the figures. For the sake of brevity, detailed descriptions of known methods are omitted here. In the above embodiments, several specific steps are described and shown as examples. However, the method process of this application is not limited to the specific steps described and shown. Those skilled in the art can make various changes, modifications, and additions, or change the order of steps, after understanding the spirit of this application.

[0106] The functional blocks shown in the above block diagram can be implemented as hardware, software, firmware, or a combination thereof. When implemented in hardware, they can be, for example, electronic circuits, application-specific integrated circuits (ASICs), appropriate firmware, plug-ins, function cards, etc. When implemented in software, the elements of this application are programs or code segments used to perform the required tasks. Programs or code segments can be stored on a machine-readable medium or transmitted over a transmission medium or communication link via data signals carried on a carrier wave. "Machine-readable medium" can include any medium capable of storing or transmitting information. Examples of machine-readable media include electronic circuits, semiconductor memory devices, ROM, flash memory, erasable ROM (EROM), floppy disks, CD-ROMs, optical disks, hard disks, fiber optic media, radio frequency (RF) links, etc. Code segments can be downloaded via computer networks such as the Internet, intranets, etc.

[0107] It should also be noted that the exemplary embodiments mentioned in this application describe methods or systems based on a series of steps or apparatus. However, this application is not limited to the order of the above steps; that is, the steps can be performed in the order mentioned in the embodiments, or in a different order, or several steps can be performed simultaneously.

[0108] The aspects of this disclosure have been described above with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this disclosure. It should be understood that each block in the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, or other programmable data processing apparatus to produce a machine such that these instructions, executable via the processor of the computer or other programmable data processing apparatus, enable the implementation of the functions / actions specified in one or more blocks of the flowchart illustrations and / or block diagrams. Such a processor can be, but is not limited to, a general-purpose processor, a special-purpose processor, a special application processor, or a field-programmable logic circuit. It is also understood that each block in the block diagrams and / or flowcharts, and combinations of blocks in the block diagrams and / or flowcharts, can also be implemented by special-purpose hardware performing the specified functions or actions, or can be implemented by a combination of special-purpose hardware and computer instructions.

[0109] The above description is merely a specific implementation of this application. Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the systems, modules, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here. It should be understood that the protection scope of this application is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in this application, and these modifications or substitutions should all be covered within the protection scope of this application.

Claims

1. A method for testing time synchronization networks, characterized in that, include: In response to the access command of the device under test, the state equation and system parameters of the preset nonlinear dynamic disturbance are obtained; The modally matched phase perturbation sequence is injected into the device under test, and the synchronization error response sequence output by the device under test under perturbation excitation is continuously captured. The phase perturbation sequence is constructed based on the state equation and system parameters. Variational mode decomposition is performed on the synchronization error response sequence to obtain K eigenmode functions with finite bandwidth, each eigenmode function corresponding to a dynamic mode component of the test system; A Hilbert transform is performed on each intrinsic mode function to obtain the Hilbert spectrum of the spatiotemporal synchronization loop under test. The Hilbert spectrum reflects the variation of the system frequency components with time. Based on the intrinsic mode functions, a parameterized reduced-order model of the system is constructed. The stability domain boundary and instability bifurcation point of the system in the parameter space are predicted by the numerical extension algorithm. A synchronization stability boundary map of the device under test is generated, and the bidirectional synchronization capability of the device under test is verified according to the synchronization stability boundary map.

2. The method according to claim 1, characterized in that, The process of obtaining the state equation and system parameters of the preset nonlinear dynamic disturbance includes: In response to the access command of the device under test, the type and protocol version of the device under test are identified, and the corresponding nonlinear dynamic system template is loaded; The initial state equation form is read from the system template, wherein the state variables include at least the phase difference, frequency difference and digital loop filter state, and the system parameters include at least the loop gain, crystal oscillator aging coefficient and filter coefficient. The system continuously captures the synchronization error sequence output by the device under test in a non-perturbation state, and estimates the initial values ​​of the system parameters online based on the subspace identification method, so that the state equation matches the actual dynamic characteristics of the device under test. The estimated system parameters are written to the parameter storage area in real time for subsequent perturbation sequence generation modules to use.

3. The method according to claim 1, characterized in that, Before injecting the mode-matched phase perturbation sequence into the device under test and continuously capturing the synchronization error response sequence output by the device under test under perturbation excitation, the method further includes: Calculate the system's natural oscillation frequency and damping ratio based on the loop gain and crystal oscillator aging coefficient in the current system parameters; Based on the natural oscillation frequency and damping ratio, a sinusoidal perturbation sequence is generated as an intrinsic mode perturbation to excite the main oscillation mode of the system; A step perturbation sequence based on the natural oscillation frequency and damping ratio generation amplitude increasing exponentially from the initial value to a preset upper limit is used as a bifurcation detection perturbation. The bifurcation detection perturbation is used to track the stability changes of the system equilibrium point and locate the bifurcation point. Based on the natural oscillation frequency and damping ratio, pulse train perturbations with adjustable pulse width and interval are generated as chaotic edge perturbations to detect the system's sensitivity to sudden disturbances and the conditions for the appearance of chaotic attractors. The intrinsic mode perturbation, bifurcation detection perturbation, and chaotic edge perturbation are combined in a preset order to form a multi-mode phase perturbation sequence, and the parameters of subsequent perturbations are dynamically adjusted according to the response analysis results of the previous round to form an adaptive closed-loop detection.

4. The method according to claim 1, characterized in that, The variational mode decomposition of the synchronization error response sequence yields K finite-bandwidth eigenmode functions, including: Construct a constrained variational problem that minimizes the sum of the estimated bandwidths of each mode and equals the sum of all modes to the original signal; By introducing a quadratic penalty factor and Lagrange multipliers, the constrained variational problem is transformed into an unconstrained variational problem, which is then solved iteratively in the frequency domain using the alternating direction multiplier method. During the iteration process, the number of modes K is dynamically adjusted according to the degree of overlap between the center frequency of the current mode and the frequency band of the adjacent modes, until the center frequencies of each mode are reasonably separated and there is no mode aliasing. Output K intrinsic mode functions, each representing a dynamic mode component of the system, and sorted by center frequency from low to high.

5. The method according to claim 1, characterized in that, The process of performing a Hilbert transform on each intrinsic mode function to obtain the Hilbert spectrum of the spatiotemporal synchronization loop under test includes: Perform a Hilbert transform on each intrinsic mode function to construct an analytic signal; Calculate the instantaneous amplitude and instantaneous phase from the analyzed signal; The instantaneous frequency is obtained by differentiating the instantaneous phase. The instantaneous amplitudes and frequencies of all modes are mapped onto the time-frequency plane, and the Hilbert spectrum is obtained by superposition.

6. The method according to claim 5, characterized in that, The construction of a parameterized reduced-order model of the system based on the intrinsic mode functions includes: A damping term model is established by fitting nonlinear damping coefficients from the instantaneous amplitude decay patterns of each intrinsic mode function; A Duffing-type frequency response model is established by fitting nonlinear stiffness coefficients from the instantaneous frequency-instantaneous amplitude relationship of each intrinsic mode function. Calculate the cross-correlation function between the perturbation input signal and the response output signal, and use the time delay corresponding to the peak value of the cross-correlation as the estimate of the system time delay; A third-order nonlinear state-space model is constructed using nonlinear damping coefficient, nonlinear stiffness coefficient, and system time delay as state variables. The third-order nonlinear state-space model is determined as the parameterized reduced-order model.

7. The method according to claim 1, characterized in that, The method of predicting the stability domain boundary and instability bifurcation point of the system in the parameter space using numerical extension algorithms includes: A parameter plane is constructed by using loop gain and crystal oscillator aging rate as continuously varying parameters and fixing other parameters as the current identification values. The pseudo-arc length extension algorithm is used to track the trajectory of the system equilibrium point as the parameters change, and the eigenvalues ​​of the Jacobian matrix of the equilibrium point are calculated at each extension step. When the real part of the eigenvalue changes from negative to positive, record this point as a bifurcation point and mark it on the parametric plane; Connecting the bifurcation points under different initial conditions forms the boundary curve of the stable region, with the inside of the curve being the stable region and the outside being the unstable region; Based on the stable and unstable regions, a two-dimensional stability domain map is generated, and the stability margin of the current operating point from the boundary is calculated and presented as a percentage. Based on the stability margin and the two-dimensional stability domain map, a synchronization stability boundary map is generated.

8. A time synchronization network testing device, characterized in that, The device includes: The acquisition module is used to acquire the state equation and system parameters of the preset nonlinear dynamic disturbance in response to the access command of the device under test. The capture module is used to inject the mode-matched phase perturbation sequence into the device under test and continuously capture the synchronization error response sequence output by the device under test under perturbation excitation. The phase perturbation sequence is constructed based on the state equation and system parameters. The decomposition module is used to perform variational mode decomposition on the synchronization error response sequence to obtain K eigenmode functions with finite bandwidth, each eigenmode function corresponding to a dynamic mode component of the test system; The transformation module is used to perform Hilbert transformation on each intrinsic mode function to obtain the Hilbert spectrum of the tested spatiotemporal synchronization loop. The Hilbert spectrum reflects the variation of the system frequency components with time. The generation module is used to construct a parameterized reduced-order model of the system based on the intrinsic mode functions, predict the stable domain boundary and instability bifurcation point of the system in the parameter space through a numerical extension algorithm, and generate a synchronization stability boundary map of the device under test, so as to verify the bidirectional synchronization capability of the device under test based on the synchronization stability boundary map.

9. An electronic device, characterized in that, The device includes: a processor and a memory storing computer program instructions; When the processor executes the computer program instructions, it implements the time synchronization network testing method as described in any one of claims 1-7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores computer program instructions, which, when executed by a processor, implement the time synchronization network testing method as described in any one of claims 1-7.

Images